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

[Martin Losch]

Hi Mark,

has anyone answered you yet?

Thanks for pointing out this problem. I checked in your fix.

Martin

On 26 Apr 2007, at 06:15, Mark Hadfield wrote:

> 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)