[MITgcm-devel] lab_sea.hb87 restart problem
Jean-Michel Campin
jmc at ocean.mit.edu
Mon Nov 12 20:35:48 EST 2007
One more thing:
seaice_sigma12 is located at a corner point,
and some grid distance , area or tanPhiAtU ...
are really function of latitude, are not exchanged
and have different value @ j=1 & j=Ny+1.
This might be enough to explain why, with noSlip BC,
seaice_sigma12 for j=1 & j=Ny+1 diverge.
If I close the lab_sea domanin at the Northern boundary
(adding a wall for j=16), the 1+1 =2 test pass with no diff,
even when I split the domain in 4 tiles (nSx=nSy=2).
May be this is the solution (if noSlip=TRUE) which, from
my point of view, is better than adding EXCH (tends to copy
a wrong calculation at the Northern edge but don't really fix it).
Jean-Michel
On Mon, Nov 12, 2007 at 07:10:59PM -0500, Jean-Michel Campin wrote:
> Hi Martin,
>
> I think I know (well, part of it) why the restart does not work
> in this lab_sea test:
> I printed seaice_sigma12(i,j,bi,bj) for i=3,
> and clearly the value for j=1 and j=sNy+1 tend to diverge
> (but both are stepped forward in seaice_evp.F, but not exchanged).
> And the fact they diverge might be particular to the no-slip BC
> (why ? not sure yet).
> When we restart, we only write/read [1:sNy] and do an exch (1 -> j+1),
> and this explains why we get differences.
>
> 2 little things:
> 1) I was thinking of adding iEVPstep,myTimer & myIter as argument
> to S/R seaice_calc_strainrates.F . It really helps to debug.
> I've done it in my own version, but could be worth to check-in.
> 2) The call to seaice_calc_strainrates in seaice_ocean_stress.F
> is wrong. Did you already wrote about this in a previous e-mail ?
>
> Cheers,
> Jean-Michel
>
> On Mon, Nov 12, 2007 at 11:46:10AM -0500, Jean-Michel Campin wrote:
> > Hi Martin,
> >
> > I started to investigate into this direction (overlap & loop range)
> > but did not go very far.
> > Trying to do some debub print & plots, but it's a pain that
> > myTime,myIter are not passed to all those S/R.
> >
> > If we add some EXCH, it will still be good to know why we
> > need them. And also to call the right EXCH (e.g., e12 is
> > at a corner point, so an EXCH_Z would be better).
> >
> > Cheers,
> > Jean-Michel
> >
> > On Mon, Nov 12, 2007 at 10:20:05AM +0100, Martin Losch wrote:
> > > Jean-Michel,
> > > here's the solution:
> > > I need to exchange the strain rates e11/e22/e12 (all three of them)
> > > after they have been comuted in seaice_calc_strainrates.F
> > > _EXCH_XY_R8( e11, myThid )
> > > _EXCH_XY_R8( e22, myThid )
> > > _EXCH_XY_R8( e12, myThid )
> > > then the pickup is OK.
> > > This changes the testreport results of lab_sea.hb87 (but only that
> > > one and only the exchanges of e11 and e12). I do not see any reason,
> > > why these extra exchanges should change lab_sea.hb87 but not any of
> > > the other experiments. The computational stencil in the IF
> > > (SEAICE_no_slip) THEN/ENDIF block is absolutely the same as in the
> > > previous lines of code.
> > >
> > > Absolutely clueless,
> > > Martin
> > >
> > > On 9 Nov 2007, at 20:44, Martin Losch wrote:
> > >
> > > >Hi Jean-Michel,
> > > >
> > > >I am stuck. With the following parameters in data.seaice
> > > >>&SEAICE_PARM01
> > > >> SEAICEwriteState = .TRUE.,
> > > >> SEAICE_initialHEFF = 1.0,
> > > >> SEAICE_deltaTtherm = 3600.,
> > > >> SEAICE_deltaTdyn = 3600.,
> > > >># useHB87stressCoupling = .true.,
> > > >> SEAICE_no_slip = .TRUE.,
> > > >> SEAICE_deltaTevp = 60.,
> > > >I can get a perfect pickup, if I set hFacU to zero in the if
> > > >(SEAICE_no_slip) block in seaice_calc_strainrates.F like this:
> > > >> hFacU = _maskW(i,j,k,bi,bj) - _maskW(i,j-1,k,bi,bj)
> > > >> hFacV = _maskS(i,j,k,bi,bj) - _maskS(i-1,j,k,bi,bj)
> > > >> hFacU = 0. _d 0
> > > >>CML hFacV = 0. _d 0
> > > >>
> > > >> e12(I,J,bi,bj) = e12(I,J,bi,bj)
> > > >> & + recip_rAz(i,j,bi,bj) *
> > > >> & ( hFacU * ( _dxC(i,j-1,bi,bj)*uFld(i,j ,bi,bj)
> > > >> & + _dxC(i,j, bi,bj)*uFld(i,j-1,bi,bj) )
> > > >> & + hFacV * ( _dyC(i-1,j,bi,bj)*vFld(i ,j,bi,bj)
> > > >> & + _dyC(i, j,bi,bj)*vFld(i-1,j,bi,bj) ) )
> > > >> & - hFacU
> > > >> & * 0.25 _d 0 * (uFld(I,J,bi,bj)+uFld(I ,J-1,bi,bj))
> > > >> & * ( _tanPhiAtU(I,J,bi,bj) + _tanPhiAtU(I,J-1,bi,bj) )
> > > >> & *recip_rSphere
> > > >This means, that at least the no-slip problem is related to uIce.
> > > >But it cannot be directly related to uice, because the pickup (and
> > > >this piece of code) is ok for many other cases (e.g. LSR).
> > > >Somewhere the true problem is masked. I don't have a clue. I don't
> > > >even know what to test next ...
> > > >
> > > >Martin
> > > >
> > > >
> > > >On 9 Nov 2007, at 17:40, Jean-Michel Campin wrote:
> > > >
> > > >>Hi Martin,
> > > >>
> > > >>I don't know much about the cause of the problem,
> > > >>but you narrow down quiet well where it sould be.
> > > >>
> > > >>I agree that it should not be too difficult to
> > > >>make this tst1+1 script more tolerent regarding blanks
> > > >>upper-case/lower case.
> > > >>
> > > >>I have a general question:
> > > >>What about adding few "data.tst" in some verification/*/input
> > > >>experiments
> > > >>to facilitate those kind of test ? I am not sure it's a good
> > > >>idea, but would interesting to know what other people think.
> > > >>
> > > >>Cheers,
> > > >>Jean-Michel
> > > >>
> > > >>On Fri, Nov 09, 2007 at 05:21:37PM +0100, Martin Losch wrote:
> > > >>>OK, there error was somewhere between keyboard and chair: your
> > > >>>script
> > > >>>is case sensitive and I was sloppy in editing data.tst.
> > > >>>
> > > >>>Now I have repeated the tests with your script and this is what I
> > > >>>get:
> > > >>>original lab_sea.hb87
> > > >>>>csysm3::tr_run.hb87> ./tst1+1 2
> > > >>>>10 12 14
> > > >>>>0000000010 0000000012 0000000014
> > > >>>>gcmExc=mitgcmuv
> > > >>>>-- compare cg2d_init_res :
> > > >>>>run 1iA:
> > > >>>> 6.08551932893086E-02
> > > >>>> 6.55586521348751E-02
> > > >>>>run 1iB:
> > > >>>> 7.02066252452855E-02
> > > >>>> 7.29377130933022E-02
> > > >>>>run 2it:
> > > >>>> 6.08551932893086E-02
> > > >>>> 6.55586521348751E-02
> > > >>>> 7.01960854205512E-02
> > > >>>> 7.29481919712837E-02
> > > >>>lab_sea.hb87 with useHB87stressCoupling = .false.,
> > > >>>>csysm3::tr_run.hb87> ./tst1+1 2
> > > >>>>10 12 14
> > > >>>>0000000010 0000000012 0000000014
> > > >>>>gcmExc=mitgcmuv
> > > >>>>-- compare cg2d_init_res :
> > > >>>>run 1iA:
> > > >>>> 6.08586878524103E-02
> > > >>>> 6.55715983076715E-02
> > > >>>>run 1iB:
> > > >>>> 7.02079681169673E-02
> > > >>>> 7.29487338830253E-02
> > > >>>>run 2it:
> > > >>>> 6.08586878524103E-02
> > > >>>> 6.55715983076715E-02
> > > >>>> 7.02078133431693E-02
> > > >>>> 7.29489594851236E-02
> > > >>>with useHB87stressCoupling = .false., and evp turned off (#
> > > >>>SEAICE_deltaTevp = 60.,)
> > > >>>>10 12 14
> > > >>>>0000000010 0000000012 0000000014
> > > >>>>gcmExc=mitgcmuv
> > > >>>>-- compare cg2d_init_res :
> > > >>>>run 1iA:
> > > >>>> 6.08355827055190E-02
> > > >>>> 1.16426574689837E-01
> > > >>>>run 1iB:
> > > >>>> 7.80912267704927E-02
> > > >>>> 7.36159469063700E-02
> > > >>>>run 2it:
> > > >>>> 6.08355827055190E-02
> > > >>>> 1.16426574689837E-01
> > > >>>> 7.80912267704927E-02
> > > >>>> 7.36159469063700E-02
> > > >>>with useHB87stressCoupling = .true., and evp turned off (#
> > > >>>SEAICE_deltaTevp = 60.,)
> > > >>>and with useSeaice=.false.,
> > > >>>>10 12 14
> > > >>>>0000000010 0000000012 0000000014
> > > >>>>gcmExc=mitgcmuv
> > > >>>>-- compare cg2d_init_res :
> > > >>>>run 1iA:
> > > >>>> 6.13113158130375E-02
> > > >>>> 6.57463167891239E-02
> > > >>>>run 1iB:
> > > >>>> 7.04816593980862E-02
> > > >>>> 7.33615824658603E-02
> > > >>>>run 2it:
> > > >>>> 6.13113158130375E-02
> > > >>>> 6.57463167891239E-02
> > > >>>> 7.04816593980862E-02
> > > >>>> 7.33615824658603E-02
> > > >>>
> > > >>>My conclusions from these numbers:
> > > >>>1. Your script is great!!!
> > > >>>2. The restart problem is related to EVP.
> > > >>>However, I can reproduce your result in lab_sea/run, that the
> > > >>>restart
> > > >>>does work with EVP. Puzzling. Therefore I reran the test in lab_sea/
> > > >>>run with SEAICE_no_slip = .true. (which is also used in
> > > >>>tr_run.hb87),
> > > >>>and then:
> > > >>>>-- compare cg2d_init_res :
> > > >>>>run 1iA:
> > > >>>> 6.08981197561294E-02
> > > >>>> 6.55972956406261E-02
> > > >>>>run 1iB:
> > > >>>> 7.02007081962388E-02
> > > >>>> 7.29011820001111E-02
> > > >>>>run 2it:
> > > >>>> 6.08981197561294E-02
> > > >>>> 6.55972956406261E-02
> > > >>>> 7.02007451116285E-02
> > > >>>> 7.29010160470697E-02
> > > >>>For tr_run.lsr everything is OK even with SEAICE_no_slip = .true.
> > > >>>
> > > >>>So the problem is related to a combination of EVP and no-slip
> > > >>>boundary conditions (and not the hb87 stress formulation). Now, the
> > > >>>flag SEAICE_no_slip is not used anywhere in the EVP code, except in
> > > >>>seaice_calc_strainrates.F, which is also used by the LSR code. I
> > > >>>don't have no clue where this problem comes from. Why should break
> > > >>>this flag the pickup for the evp code only, if it is not used in the
> > > >>>code that is exclusive to EVP?
> > > >>>
> > > >>>Any idea?
> > > >>>
> > > >>>Martin
> > > >>>
> > > >>>
> > > >>>
> > > >>>
> > > >>>
> > > >>>_______________________________________________
> > > >>>MITgcm-devel mailing list
> > > >>>MITgcm-devel at mitgcm.org
> > > >>>http://mitgcm.org/mailman/listinfo/mitgcm-devel
> > > >>_______________________________________________
> > > >>MITgcm-devel mailing list
> > > >>MITgcm-devel at mitgcm.org
> > > >>http://mitgcm.org/mailman/listinfo/mitgcm-devel
> > > >
> > > >_______________________________________________
> > > >MITgcm-devel mailing list
> > > >MITgcm-devel at mitgcm.org
> > > >http://mitgcm.org/mailman/listinfo/mitgcm-devel
> > >
> > > _______________________________________________
> > > MITgcm-devel mailing list
> > > MITgcm-devel at mitgcm.org
> > > http://mitgcm.org/mailman/listinfo/mitgcm-devel
> > _______________________________________________
> > MITgcm-devel mailing list
> > MITgcm-devel at mitgcm.org
> > http://mitgcm.org/mailman/listinfo/mitgcm-devel
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel
More information about the MITgcm-devel
mailing list