[MITgcm-devel] seaice adjoint and EVP

[Martin Losch]

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