[MITgcm-devel] lab_sea.hb87 restart problem
Martin Losch
Martin.Losch at awi.de
Tue Nov 13 05:36:53 EST 2007
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.
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.
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
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
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
More information about the MITgcm-devel
mailing list