[MITgcm-devel] seaice adjoint and EVP

[Martin Losch]

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