[MITgcm-devel] lab_sea.hb87 restart problem
Jean-Michel Campin
jmc at ocean.mit.edu
Tue Nov 13 13:04:31 EST 2007
Hi Martin,
On Tue, Nov 13, 2007 at 11:36:53AM +0100, Martin Losch wrote:
> Hi Jean-Michel,
>
> nice work, now it makes sense:
> for free slip e12, seaice_shear, and seaice_sigma12 are zero on grid
> (corner) points that are part of a dry/land cell. In lab_sea this is
> the case for j=17 and j=1, so e12(:,1)=e12(:,17)=0.
> The same is (should be) true for i=1 and i=21. For no slip e12 is non-
> zero on cornerpoints belonging to land and the problem you describe
> (inappropriate exchanges of seaice_sigma12 from j=1 to j=sNx+1, etc.)
> occurs.
>
You are right.
> What we really need is a pickup that stores seaice_sigma12(1:sNx
> +1,1:sNy+1,:,:), right?
> Actually the exchange of seaice_sigma12 in seaice_pickup is wrong! It
> should be a Z-exchange. Could you tell me what the correct macro is.
> I can then fix this.
We don't have the proper IO S/R to do
the correct pickup (up to sNy+1). Note that we will
have also a problem on the cube (the 2 missing corners !)
unless sea-ice is not floating around the corners.
And when I want to find plenty of exotic EXCH call,
I look to ini_curvilinear_grid.F. No macro, just:
CALL EXCH_Z_3D_RL( seaice_sigma12, 1, myThid )
should do it.
> The quick solution is your suggestion: only use domains that have
> their solid boundaries in the North and East within the domain. I
> didn't like the exchange solution, either (and I knew that e12 needs
> a Z-exchange, was too lazy to look it up), and in fact, it's wrong.
>
> Go ahead and add extra parameters to seaice_calc_strainrates
I think I will do the same opration for seaice_calc_viscosities,
since it's the same set of S/R to modify.
> Yes, in ocean_stress, uIce/vIce should be copied to 4D fields which
> should then be passed to calc_strain rates. Alternatively one could
> move the bi/bj-loops out of calc strain rates and then pass only 2D
> arrays. Sorry about that. What do you prefer?
>
> Martin
>
It's up to you. There is also a 3rd possibility, no copy, but
require to add the number of "levels" of uIce,vIce as argument
to S/R seaice_calc_strainrates.F (e.g.: kSize)
(and inside seaice_calc_strainrates.F, declare those 2 array with
kSize but only refer to the level 1). This way you don't have
to change the bi,bj loops, just need to call seaice_calc_strainrates
with the right kSize.
And regarding the usefullness of a "good" script (like tst1+1)
to test the restart, I agree, and I think it should be possible
to avoid the addition of a "data.tst" file (my bad suggestion of
the other day). If we agree to drop few options (globalFiles,
useSingleCpuIO, startTime, endTime, mdsioLocalDir, pickupSuff,
mnc pickups) for this restart test, should be able to generate
automatically a data.tst and do all the steps in 1.
Cheers,
Jean-Michel
>
> On 13 Nov 2007, at 02:35, Jean-Michel Campin wrote:
>
> >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
> >_______________________________________________
> >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