[MITgcm-devel] bug: divisioin by zero

[Martin Losch]

Hi,

this is probably a bug:

on our new sx8 I get a division by zeros in mom_calc_visc.F for  
curvilinear grids on the cubed sphere/llc and variable viscosity. The  
culprits are recip_dxV and recip_dyU, because they are the only ones  
that are not exchanged in ini_curvilinear_grid.F
Why are they not exchanged? I gather from the comments in  
ini_curvilinear_grid, that for cube sphere exchanges this has not  
been sorted out and some sort of hack has been applied which is now  
catching up with me. What can we do?
1. use the wrong exchanges so that dxV/dyU and their reciprocals have  
values other than zero in the overlaps and keep our fingers crossed  
that no-one will notice (I assume that the viscosity values in the  
overlaps are not that important (o:)
2. modify mom_calc_visc, so that it does not compute anything in the  
overlaps (will that work?)

solution 1 changes the results of testreport for  
global_ocean.cs32x15, if I uncomment these two lines in  
ini_curvilinear_grid.F:
c     CALL EXCH_Z_3D_RS( dxV, 1, myThid )
c     CALL EXCH_Z_3D_RS( dyU, 1, myThid )
I could do this:
#ifdef TARGET_SX_NEC
      CALL EXCH_Z_3D_RS( dxV, 1, myThid )
      CALL EXCH_Z_3D_RS( dyU, 1, myThid )
#endif
which will fix the problem for me but as soon as you run into a  
compiler that is as picky as the sxf90 you'll have the problem again.

(BTW. sxf90 has options to initial all stack/heap with nan. If I do  
that the model stops in mom_vi_hdissip.F, havent' figured out yet, why.)

What do you think?

Martin

PS. I hate going to new machines, all these problems that you  
discover ...