[Mitgcm-support] GM c48

[mitgcm-support at dev.mitgcm.org]

Hi guys,

I started to look at c48 version of GM. 

o The fact that the 4 experiments that used GM fail the testscript 
  in both checkpoint47g_post & checkpoint48 is annoying (for us and
  also for others that download the code, run an example and
  cannot reproduce the standard results).

o Things are not simple. gmredi_slope_limit.F is becoming almost
  as hard to read as pieces of code that I usually write myself.
  Let's concentrate on the 2 S/R that are responsible for the output 
  differences: gmredi_slope_limit.F & gmredi_calc_tensor.F

----------------------
gmredi_slope_limit.F 
 a) a 2nd tapering scheme has been added over all the others.
   This is enough to change the results:
Y Y Y Y  8 13 10 10  9 16 12 12  8  7  9  8  7  8  9  7  7 FAIL  global_ocean.90x40x15
Y Y Y Y  8 13 10 10  9 16 11 13  9  7  9  8  8  8  9  7  8 FAIL  global_with_exf  
   This is not a compiler Optimization issue, this is a real change in
   the code. Therefore, it would nice to put the 2nd slope_max (=Large_SlopeSqr)
   in data.gmredi so that one can recover the original code & results
   setting e.g. Large_SlopeSqr = 1.e48 (I checked that it works). 
   And for the adjoint, Large_SlopeSqr= 1.e8.

 b) A 2 steps computation (multiplication by the inverse) has been replaced
   by a single division:
   c48, line 250,251 :
           SlopeX(i,j) = -dSigmaDx(i,j)/dSigmaDrReal(i,j)
           SlopeY(i,j) = -dSigmaDy(i,j)/dSigmaDrReal(i,j) 
   47f_post:
           dRdSigmaLtd(i,j) = 1./(dSigmaDrReal(i,j))
           SlopeX(i,j)=-dSigmaDx(i,j)*dRdSigmaLtd(i,j)
           SlopeY(i,j)=-dSigmaDy(i,j)*dRdSigmaLtd(i,j)
   This is enough to change the results:
Y Y Y Y  9 16 16 16 16 16 16 16 16 13 13 12 16 12 11 10 13 FAIL  front_relax 
Y Y Y Y 12 16 16 16 16 16 16 16 16 13 13 12 13 13 13 12 13 FAIL  lab_sea 
   This is not sensitive to the Optimization level (I get differences also
   if I use no optimization at all).
   We can decide to change the output of those 2 exp. But I would prefer
   to change at the same time both the tapering part and the clipping part
   (right now, only the tapering part has been simplified).
   What I propose is:
    1) return to the old form.
    2) make a tag
    3) replace the *dRdSigmaLtd by /dSigmaDrReal in BOTH the tapering part
       and the clipping part. (if it work with TAF for one, it should work
       also for the other ?). 
    4) change the output that are affected. put some clear comments that
       explain what causes the differences.
    5) make a tag.

 c) GM_Small_Number=1.D-12 (but not used in gmredi_slope_limit.F)
   Small_Number=1.D-20 (before was Small_Number=1.D-12)
   Again, this is enough to change the results: 
Y Y Y Y 12 16 16 13 13 16 16 16 13 12 12 12 12 11 12 12 12 FAIL  global_ocean.90x40x15
Y Y Y Y 12 16 16 13 13 16 16 16 13 12 12 12 12 11 12 12 12 FAIL  global_with_exf 

 d) Cspd=2. _d 3 (before was  Cspd=2.) in the ldd97  tapering scheme.
   This looks like a bug ?
   this affects the results of lab_sea:
Y Y Y Y  1  3  4  5  4  5  5  6  4  3  4  2  3  3  4  2  3 FAIL  lab_sea 
 
----------------------
gmredi_calc_tensor.F
   c48, line 271:
     &          *_maskW(i,j,k,bi,bj)*maskp1
   c48, line 344:
     &          *_maskS(i,j,k,bi,bj)*maskp1
   *maskp1 was not there before and it should not be there.
   I think that this breaks the "adiabatic property' of the Redi tensor.
   Anyway, this is enough to change the results:
Y Y Y Y  6 16  8  8  7 16  9 12  6  6  6  7  5  5  6  5  5 FAIL  global_ocean.90x40x15
Y Y Y Y  6 16  7  8  7 16  9 12  6  6  6  7  5  5  6  5  5 FAIL  global_with_exf
   
   Patrick: could you check that the adjoint is OK with the right 
     form, i.e., with no *maskp1
----------------------

  Concerning the advective form of GM (with GM_AdvSeparate T or F),
  I can reproduce the results of c46. 
 
Question for Patrick:
  In gmredi_calc_tensor.F, the TAF specific initializations are done
  in the 1rst k loop, DO k=2,Nr ; why not in a specific block, before 
  this K=2,Nr loop, and this time from k=1,Nr ? 

what I will do now: try to reach point 2 of note.b above
  
See you,

Jean-Michel