[MITgcm-devel] seaice adjoint and EVP

[Jinlun Zhang]

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
>>>
>>> _______________________________________________
>>