[MITgcm-support] OBS_CALC bug: divide by zero at internal boundaries

[Mark Hadfield]

There's a problem in obs_calc.F that manifests itself as a crash with 
some compilers in serial runs when nSx or nSy is greater than one and 
useOBCSbalance is true.

The code in question is at the end of subroutine OBS_CALC and enclosed 
in an IF block

      IF ( useOBCSbalance) THEN
...
      ENDIF

For the east, west, north and south boundaries in turn the code 
calculates the boundary-integrated area (Ar_T) and transport (Tr_T). It 
then calculates a correction term

      Tr_T = (0. - Tr_T)/Ar_T

and adds this to the velocities at that boundary.

Consider the case where the boundaries are at the edges of the domain, eg:

  OB_Iwest = 160*1
  OB_Ieast = 160*-1
  OB_Jsouth = 198*1
  OB_Jnorth = 198*-1

If, say nSx = 2 and nSy = 1, then tile 0 has no eastern boundary points 
and tile 1 has no western boundary points. In this case, Ar_T = Tr_T = 0 
and the expression above is 0/0. With some compilers this evaluates to 
zero and the code does what it should (i.e., nothing). With other 
compilers (eg. g95) the division by zero leads to a crash.

In the attached file the problem is fixed by checking for Ar_T > 0.

-- 
Mark Hadfield          "Ka puwaha te tai nei, Hoea tahi tatou"
m.hadfield at niwa.co.nz
National Institute for Water and Atmospheric Research (NIWA)


-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: obcs_calc.F
URL: <http://mitgcm.org/pipermail/mitgcm-support/attachments/20070426/d854cca7/attachment.el>