[MITgcm-devel] seaice_calc_viscosities

[Martin Losch]

Hi Jean-Michel,

in seaice_calc_viscosities, there’s a piece of code (the computation of deltaC) that causes problems.

this is how it should look, delta=sqrt( (eps11+eps22)**2 + (eps11-eps22)**2/e**2 + 4 eps12**2 /e**2 )
for some reason I expanded the squares and then I get this formulation:
          deltaC(I,J,bi,bj) = 
     &         (e11(i,j,bi,bj)**2+e22(i,j,bi,bj)**2)*(ONE+ecm2)
     &         + 4. _d 0*ecm2*e12Csqr(i,j)
     &         + 2. _d 0*e11(i,j,bi,bj)*e22(i,j,bi,bj)*(ONE-ecm2)

Unfortunately, this can be less that zero in some rare cases (we have one at very high resolution) and then the sqrt(deltaC) produces NaN. I would like to change the code, but the many tests suffer from the truncation error. by up to 7digits (see below, I didn’t try the adjont yet). Should I include the modification within ifdefs for now, so that we can reproduce old results, or should I just update the tests in this case? I have it in the back of my head that I preferred the current version because it avoids evaluating e11-e22. Was that a good reason?

Martin

Thu Apr 30 18:21:36 CEST 2015
run: ./testreport -t 'global_ocean.cs32x15 lab_sea offline_exf_seaice seaice_itd seaice_obcs' -j 8 -q
on : Linux bkli04l006 3.13.0-49-generic #83-Ubuntu SMP Fri Apr 10 20:11:33 UTC 2015 x86_64 x86_64 x86_64 GNU/Linux

No "OPTFILE" was specified ; genmake2 found and uses:
                         OPTFILE=../../../tools/build_options/linux_amd64_gfortran

default 10  ----T-----  ----S-----  ----U-----  ----V-----  --PTR 01--  --PTR 02--  --PTR 03--  --PTR 04--  --PTR 05--
G D M    c        m  s        m  s        m  s        m  s        m  s        m  s        m  s        m  s        m  s
e p a R  g  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .
n n k u  2  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d
2 d e n  d  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .

Y Y Y Y>16<16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 pass  global_ocean.cs32x15
Y Y Y Y>11<16 16 16 16 16 16 16 14 16 16 14 16 16 16 14 16 16 16 16 16 16 16 pass  global_ocean.cs32x15.icedyn
Y Y Y Y>11<16 16 16 16 16 16 16 16 16 16 14 13 16 16 14 16 22 16 16 16 22 16 16 16 pass  global_ocean.cs32x15.seaice
Y Y Y Y>16<16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 13 16 13 16 16 16 pass  global_ocean.cs32x15.thsice
Y Y Y Y>16<16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 pass  global_ocean.cs32x15.viscA4
Y Y Y Y>12<16 16 16 16 16 16 16 16 16 16 14 16 16 16 13 14 22 16 16 16 22 16 16 16 pass  lab_sea
Y Y Y Y>16<16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 22 16 16 16 22 16 16 16 pass  lab_sea.fd
Y Y Y Y>16<16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 22 16 16 16 22 16 16 16 pass  lab_sea.hb87
Y Y Y Y>12<16 16 16 16 16 16 16 16 16 16 13 16 16 14 13 16 22 16 16 16 22 16 16 16 pass  lab_sea.salt_plume
Y Y Y Y -- 16 16 16 13 16 16 16 16 12 13 13>12<13 13 12 13 pass  offline_exf_seaice
Y Y Y Y -- 13 14 14 12 13 12 13 16 12 12 12>12<12 11 11 11 pass  offline_exf_seaice.dyn_jfnk
Y Y Y Y -- 13 16 16 13 13 14 16 13 12 13 13>13<13 13 11 13 pass  offline_exf_seaice.dyn_lsr
Y Y Y Y -- 16 16 16 16 16 16 16 22 22 16 16 16 22 16 16>16<pass  offline_exf_seaice.thermo
Y Y Y Y -- 16 16 16 16 16 16 16 22 16 16 16 16 16 16 16>16<pass  offline_exf_seaice.thsice
Y Y Y Y -- 16 16 16 16 16 16 16 16 12 14 12>12<13 16 12 13 pass  seaice_itd
Y Y Y Y -- 16 16 16 13 16 16 16 16 12 14 13>12<13 16 13 13 pass  seaice_itd.lipscomb07
Y Y Y Y -- 16 11  8  8 16 11 12  9  7  7  8> 7< 6  9  7  7 FAIL  seaice_itd.thermo
Y Y Y Y>14<16 16 16 16 16 16 16 16 16 16 14 16 16 16 14 16 22 16 16 16 22 16 16 16 pass  seaice_obcs
Y Y Y Y>13<16 16 16 16 16 16 16 16 16 16 13 14 16 16 16 16 22 16 16 16 22 16 16 16 pass  seaice_obcs.seaiceSponge
Y Y Y Y>14<16 16 16 16 16 16 16 16 13 16 16 16 16 16 16 16 22 16 16 16 22 16 16 16 pass  seaice_obcs.tides
Start time:  Thu Apr 30 18:21:36 CEST 2015
End time:    Thu Apr 30 18:27:26 CEST 2015