[MITgcm-devel] seaice adjoint and EVP

Wed May 16 20:04:55 EDT 2007

We just need to set zmin=0 for either open water or ice, that is, always.
Jinlun

Martin Losch wrote:

> There you go, I always thought zmin>0 over open water does not make  
> sense. What do you suggest: zMin = 4e8 if AREA>0 and 0 otherwise? Or  
> zMin = 1e8*HEFF? Or just zMin = 0. always?
>
> Dimitris, can you test that? In seaice_dynsolver set zMin to zero if  
> there is no ice (zMax is zero for HEFF=0 already) and test for both  
> LSR and EVP.
>
> (btw, setting e=1 was only meant to be a test to see what's going on,  
> I would never propose to do something like that for good.)
>
> Martin
>
> On 16 May 2007, at 18:14, Jinlun Zhang wrote:
>
>> Martin,
>>
>> I am not sure what is going on, but setting e=1 is no good, too  
>> unphysical. Also, zeta, eta over open water should be zero. I used  
>> the value of 4e8 for testing and forgot to change it back to zero.
>>
>> Jinlun
>>
>> Martin Losch wrote:
>>
>>> Hi Dimitris,
>>>
>>> this thread it getting very confusing, because it talks about many  
>>> different things. Here I will only answer to EVP related things, OK?
>>>
>>> 1. thanks for your plots and data, too bad that the bug fix did  not 
>>> help (updating to latest code will not either, but you can  always 
>>> try).
>>> 2. I plotted the gradients of sigma1/2/12 as they appear in the  
>>> stress divergence and the stripes in the velocity do show up in  
>>> dsigma12/dx (for v) and dsigma12/dy (for u) which is expected as  
>>> sigma12 is really 2*eta*(du/dy + dv/dx)
>>> 3. sigma1 is different from sigma2 because it contains P:
>>> sigma1 = sig11+sig22 = 2*zeta*(du/dx+dv/dy) - P
>>> sigma2 = sig11-sig22 = 2*eta(du/dx-dv/dy)
>>> so the big negative contribution is the ice strength P, eta and  
>>> zeta are the same except for the factor e^2=4. If you diagnose P  
>>> (SIpress) you could have a look at sigma1+P vs. sigma2 and they  
>>> should have a similar structure (may actuallly be worth trying for  
>>> one timestep just to make sure for snapshots it's even trivial to  
>>> compute P = Pstar*Heff*exp(-20*(1-Area), what's your Pstar?).
>>>
>>> I still do not have any clue what's going on. P (and zeta and eta)  
>>> is set to some (large) constant over open what (see seaice_evp,  
>>> seaice_dynsolver, etc) but that should not change the divergence  of 
>>> the stress. Currently I am having a similar problem with  stripes in 
>>> a highly anisotropic grid (dx=2deg, dy=0.5deg), and I  am suspecting 
>>> that the viscosity that Leith computes is too low  for the coarse 
>>> direction. Maybe we have similar problem here:
>>> When you write down the divergence of the stress tensor you'll get  
>>> something like this:
>>> dsig11/dx+dsig12/dy = d/dx[(e^2+1)*eta du/dx] + d/dy[2*eta du/dy]  + 
>>> d/dx[(e^2-1)*eta dv/dy] + d/dy[2*eta dv/dx]
>>> for the u equation and analogously for the v-equation. If you set  
>>> e=SEAICE_ecc = 1 (not elliptic but circular yield curve), the  
>>> second last term will drop out (unfortunately), but zeta=eta and  
>>> the components of the stress divergence will be symmetric (dd/dxx  + 
>>> dd/dyy, except fo the last term). The problem is, that if this  
>>> works, I still don't know how to fix the problem because we do  want 
>>> to have e=2. May we can increase the value of zeta/eta over  open 
>>> ocean (so zMin which is 4e8 right now, in seaice_dynsolver)?  After 
>>> all the stripes appear where velocity is high.
>>>
>>> Do you (especially Jinlun) think that this is worth a  try?
>>>
>>> Martin
>>>
>>> On 15 May 2007, at 20:50, Dimitris Menemenlis wrote:
>>>
>>>> Martin,
>>>>
>>>>> I have corrected the bug in seaice_ocean_stress (I am wondering  
>>>>> if this will
>>>>> fix the problem that Dimitris has with evp. I don't think so,  but 
>>>>> there is
>>>>> always hope o:), and updated lab_sea but not the
>>>>
>>>>
>>>> Unfortunately the stripes are still there.  This is not with very  
>>>> latest code. It is yesterday's pkg/seaice code with
>>>>
>>>>> columbia15: diff seaice_ocean_stress.F ../MITgcm/pkg/seaice/ 
>>>>> seaice_ocean_stress.F 184c184
>>>>> <           fuIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I-1,J,bi,bj) )*
>>>>> ---
>>>>>
>>>>>>           fuIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I,J+1,bi,bj) )*
>>>>>
>>>>> 195c195
>>>>> <           fvIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I,J-1,bi,bj) )*
>>>>> ---
>>>>>
>>>>>>           fvIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I+1,J,bi,bj) )*
>>>>>
>>>>
>>>> So I am assuming EVP will eventually crash as before.
>>>> Figure comparing evp to lsr solution is here
>>>> http://ecco2.jpl.nasa.gov/data1/arctic/output/tests/evp/ 
>>>> evp_vs_lsr.ps.gz
>>>> Directory
>>>> http://ecco2.jpl.nasa.gov/data1/arctic/output/tests/evp/
>>>> contains pkg/seaice output, including sigmas, matlab routine to  
>>>> read/plot and an additional figure
>>>> http://ecco2.jpl.nasa.gov/data1/arctic/output/tests/evp/evp.ps.gz
>>>> of all the evp fields
>>>> sigma* colorscale is chosen to emphasize open water.
>>>> character of sigma1 over open water is much different from sigma2  
>>>> and sigma12, is this what you expect?
>>>>
>>>> D.
>>>> _______________________________________________
>>>> 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

-- 

Jinlun Zhang
Polar Science Center, Applied Physics Laboratory
University of Washington, 1013 NE 40th St, Seattle, WA 98105-6698

Phone: (206)-543-5569;  Fax: (206)-616-3142
zhang at apl.washington.edu
http://psc.apl.washington.edu/pscweb2002/Staff/zhang/zhang.html