[MITgcm-devel] seaice adjoint and EVP

Martin Losch Martin.Losch at awi.de
Thu May 24 09:06:21 EDT 2007 

I have a new feeble theory for the noise in the evp solver over open  
ocean:

heff = 0 over open ocean, therefore seaiceMassU/V = 0.
momentum equation in seaice_evp is discretized (in time) as
m*duice/dt = -m dphi/dx + tau_air + cd*(uice-uocean) + m*f*vice +  
\nabla\sigma
m*(uice(n+1)-uice(n))/dt = -m dphi(n)/dx + tau_air(n) - cd*(uice(n+1)- 
uocean(n)) + m*f*vice(n) + \nabla\sigma(n),
so coriolis is explicit, ice-ocean stress is implicit. if the mass m  
is zero (and zetaMin=0, so that zeta=eta=press = 0 over open ocen)  
this reduces to
cd*uice(n+1) = tau_air(n) + cd*uocean(N)
so that uice ist a purely diagnostic quantity and not time stepped.  
cd is a function of uice-uocean at the nth time step, averaged to  
center points and the averaged back to velocity points.

Dimitris, could that be the problem, somehow I don't think so, but  
you can try by putting a minimum seaiceMassU/V in seaice_dynsolver.F,  
say seaiceMassU = max(seaiceMassU,SEAICE_rhoIce*0.05)

Martin

On 22 May 2007, at 18:59, Jinlun Zhang wrote:

> Hi Martin,
> Yeah, we are sort of stuck, but hey it is very interesting and  
> revealing.
> I would vote against masking ice velocities over open water  
> because, as mentioned earlier, the ice velocities would be wrong at  
> ice edge and the ice velocity discontinuity at ice edge will get  
> into ocean. (o:. We don't do the masking with LSR solver, perhaps  
> we can avoid doing that with EVP.
> Jinlun
>
> Martin Losch wrote:
>
>> Hi Jinlun,
>> the evp-solver is only in place for the C-grid. I don't have the  
>> time  to code the solver for the b-grid now. The b-grid code (for  
>> LSR) is  still working, but I have not kept it up to date, so  
>> there may be a  few thing different other than the different grids.
>>
>> In general I though that the c-grid is perfect for evp as all the   
>> discretizations fall in place naturally. Only for this \delta  
>> term  one needs to average from center to corner points and vice  
>> versa  (have a look at seaice_calc_strainrates and seaice_evp).  
>> However,  there may be issues with the coriolis terms (commonly a  
>> problem with  the c-grid).
>>
>> Actually, Elizabeth told us that she masks ice velocities over  
>> open  water in CICE.
>>
>> Now we are a little stuck, aren't we?
>>
>> Martin
>>
>> PS. I need to be able to reproduce these results myself (I  
>> haven't  been able to, yet), maybe I can debug the stuff this way.  
>> Via email  etc. it's quite demanding (o:
>>
>> On 21 May 2007, at 19:15, Jinlun Zhang wrote:
>>
>>> I wouldn't think C-grid is problematic with EVP as we have seen.   
>>> But just to make sure, is it possible to use the original B-grid   
>>> EVP to see if the same things occur? There was a B-grid ice  
>>> model  setup in place that may be used for doing B-grid.
>>> Better not zap out things over open ocean. Otherwise,  
>>> discontinuity  may occur and ocean may be screwed up.
>>> Jinlun
>>>
>>>>> Hi Martin,
>>>>>
>>>>>> 1. Are these figures all with with zMin = 0?
>>>>>
>>>>>
>>>> In this case it may be worth turning of individual terms in the   
>>>> rhs  of the momentum equations
>>>> 1. dphiSurf/dx and dphiSurf/dy (in seaice_dynsolver)
>>>> 2. surface wind stress (taux/y=0 in seaice_get_dynforcing)
>>>> 3. ice-ocean stress (DWATN in seaice_evp)
>>>> 4. Coriolis
>>>> 5. stressDivergence
>>>> 4 and 5 should be zero over open ocean anyway so I do not see  
>>>> how   these terms can lead to the stripes.
>>>> We should get to the bottom of what is causing these stripes.   
>>>> that  way we can probably understand the noise in the ice  
>>>> fields,  too.
>>>>
>>>>>
>>>>> Yes, all the figures and results under
>>>>> http://ecco2.jpl.nasa.gov/data1/arctic/output/tests/
>>>>> (except for the oldtest subdirectory) are with zMin=0.
>>>>>
>>>>>> 2. Do you have an EVP run that does not blow up at all   
>>>>>> (regardless  of noise)?
>>>>>
>>>>>
>>>>> I have not run any of the zMin=0/SEAICEuseFlooding=.true.  
>>>>> tests  out  for very
>>>>> long, but I am almost certain that none of these new   
>>>>> integrations  will crash,
>>>>> including the SEAICE_deltaTevp=60.
>>>>> The crashes had to do with snow accumulation and could happen  
>>>>> to   both LSR or to
>>>>> EVP solutions.
>>>>
>>>>
>>>> That's good news. It mean that we can (in principle) maskRHS  
>>>> flag  and  not worry about the stripes.
>>>>
>>>>>
>>>>>> 3. What's the convergence criterion for LSR, and how many    
>>>>>> interations do you allow/do? In other words how close is the   
>>>>>> LSR  solution to VP?
>>>>>
>>>>>
>>>>>
>>>>>       LSR_ERROR          = 2e-4,
>>>>>       SOLV_MAX_ITERS=1500
>>>>
>>>>
>>>> That's not very much, is it? For an accurate VP solution I  
>>>> would  put  LSR_ERROR = 1e-7 to 1e-13, right?
>>>>
>>>>>
>>>>>> c. the same is true for the wind-ice/ocean-ice stress terms   
>>>>>> which  in involve
>>>>>>  averaging perpendicular to the stripes (unless the turning   
>>>>>> angle  is not
>>>>>> equal to zero, in which case there is also averaging in the   
>>>>>> other  directions,
>>>>>> but you don't do that, do you?).
>>>>>
>>>>>
>>>>>
>>>>> No I use SEAICE_airTurnAngle=SEAICE_waterTurnAngle=0.
>>>>
>>>>
>>>> Good.
>>>>
>>>>>
>>>>>> About question 3 (is it really a VP solution): Could you   
>>>>>> diagnose  SIsigI and SIsigII (snapshots!!!! I guess one is   
>>>>>> enough) for all  (or some) solutions and
>>>>>>  plot them (plot(SIsigII(:),SIsigI(:),'x')? These should be  
>>>>>> the   principle components of sigma normalized by the  
>>>>>> strength/ pressure P.
>>>>>
>>>>>
>>>>>
>>>>> With SEAICE_dumpFreq, SIGMA1, SIGMA2, and SIGMA12 are  
>>>>> diagnosed  by  default for
>>>>> the EVP solutions but not for LSR.  Are these the same as  
>>>>> SIsigI   and SIsigII?
>>>>> Figure for SIGMA1, SIGMA2 for EVP solution is here:
>>>>> http://ecco2.jpl.nasa.gov/data1/arctic/output/tests/figs/  
>>>>> SIGMA2232.ps
>>>>> Does it look as expected?
>>>>
>>>>
>>>> sigma1/2/12 are not the principle stress components. I have  
>>>> added   diagnostics that are called SIsigI and SIsigII, which is  
>>>> what you   want. In principle you could computed them yourself  
>>>> (from snapshots):
>>>> SIsigI = 0.5*(sigma1 + sqrt(sigma2^2 + 4*sigma12^2)/Press
>>>> SIsigII = 0.5*(sigma1 - sqrt(sigma2^2 + 4*sigma12^2)/Press
>>>>
>>>> Press = max(1.e-13,Pstar * HEFF *exp( -20*(1-AREA)));
>>>>
>>>> see seaice_do_diags.F (and seaice_dynsolver.F)
>>>>
>>>>>
>>>>>> I am also a little concerned that the LSR and EVP solutions   
>>>>>> look  so different
>>>>>> in the ice-covered area, can that be attributed to that   
>>>>>> different  boundary
>>>>>> conditons? Can you try a run with no slip for the evp solver?
>>>>>
>>>>>
>>>>>
>>>>> Is LSR no slip by default?  How do you specify no slip for evp   
>>>>> solver?
>>>>
>>>>
>>>> LSR is half slip and that's hardwired. I didn't want to bother   
>>>> this  the boundary conditions if EVP works, because it's so  
>>>> much  simpler to  do that in EVP. But now I may have to  
>>>> reconsider this  decision.
>>>> EVP is free slip by default. SEAICE_no_slip = .true. makes it  
>>>> no  slip.
>>>>
>>>> Martin
>>>>
>>>> _______________________________________________
>>>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel

More information about the MITgcm-devel mailing list