[MITgcm-devel] more seaice

[Jinlun Zhang]

Hi Martin,

I see. The clipping was left in the code because I forget commenting 
them out. They were there because of some sensitivity tests I made. 
Dimitris, we need to comment them. They are not needed. So if you mark 
URT and VRT for the purpose of replacing clipping, you don't have to. 
The idea of this particular ice dynamics is that it has to be solved 
throughout the whole model domain, with ice or without ice, to avoid 
velocity discontinuity.

I have tried all sorts of schemes for calculating surface stress, all 
would blow up the code, unless spatial averaging is used, which I don't 
want to. I am still working on it......
Jinlun

Martin Losch wrote:

> Forgot to say what I mean by "clipping":
> at the end of dynsolver (and seaice_dynsolver) you have:
> uIce = min(0.40,uIce)
> uIce = max(-0.40,uIce)
> That's what I call clipping.
>
> The next plan is to look into the stress coupling of Hibler and Bryan  
> (1987, JPO): Instead of applying ice-ocean to the ocean model, use  
> the internal stress of the ice (div sigma). It sounds plausible. Do  
> have any experience with that?
>
> Martin
>
> On Mar 10, 2006, at 9:14 AM, Martin Losch wrote:
>
>> Hi Jinlun,
>>
>> first, a quick update on the AWI sea-ice model: the solver is  
>> basically the original Hibler(1979) solver, with some small  
>> modifications (I checked again with Ralph Timmermann, he said that  
>> the Belgian group that he worked in for some time uses or used to  
>> use your code, but they had problems porting it to the ORCA grid of  
>> OPA).
>>
>> Clippling vs. Masking:
>> I understand LSR as follows: it sets up and solvers a linear system  
>> of equations equation:
>> CoefficientMatrix*UICE(k+1) = URT(UICE(k)) and the same for VICE.  So 
>> if I mask the URT with the "area mask" (seaiceMaskU(I,J) = 1, if  
>> AREA(I,J)+AREA(I-1,J) > 0 otherwise = 0), there is no forcing  
>> outside the ice-area and thus the ice velocities go to zero (in  
>> fact, since it is a global solver to very small numbers) in regions  
>> where there is no ice. The ice regions are probably affected by  this 
>> near the ice boundaries, but according to my quick  tests not  very 
>> much, but I need to check that more thoroughly. This way you  will 
>> not advect zero ice thickness around with large velocities and  thus 
>> avoid numerical problems. I can get away without clipping (in  terms 
>> of numerical stability, I still don't like the solution, not  enough 
>> ice )o: )
>>
>> Martin
>>
>> On Mar 10, 2006, at 1:25 AM, Jinlun Zhang wrote:
>>
>>> Hi Martin,
>>>
>>> Thanks for the info about AWI ice modeling.
>>> I haven't looked at  mitgcm ice part for a while, but wow, you  have 
>>> done lots of work. And good you can run it on a cubed sphere  with C 
>>> grid. But I am not sure why you have to mark out URT and  VRT with 
>>> AREA. You said it is because velocity clipping. What is it?
>>>
>>> Cheers, Jinlun
>>>
>>> Martin Losch wrote:
>>>
>>>> Hello Jinlun,
>>>>
>>>> Thanks for your reply.
>>>> Stresses: I talked to Ralph Timmermann who is currently probably  
>>>> the  one who knows most about "the" AWI ice model (there are by  
>>>> now many  flavors around, one coupled to SPEM, one coupled to  
>>>> MOM2, and a stand  alone version which is basically a handy-down  
>>>> over many generations  of PhD students so that the current  
>>>> students do not really now too  much about the details of it. The  
>>>> basis of all of this is a model by  Harder and Lempke:
>>>> Harder, M.(1996). Roughness and drift trajectories of sea ice in   
>>>> large-scale simulations and their use in model verifications,  Ann  
>>>> Glaciol,25.
>>>> Kreyscher, M., Harder, M., Lemke, P.(1997). First results of the  
>>>> Sea  Ice Model Intercomparison Project (SIMIP), Annals of  
>>>> Glaciology , 25,  8-11.
>>>> And as far as I can see it does not use the LSR solver but  
>>>> something  different which has been changed by members of our  
>>>> computer  department and other. Do you know the term "electronic  
>>>> archeology"?  One would need to do a immense amount of that to  
>>>> figure this things  out (o:.
>>>> And NO, it is not parallel, something that hampers the thrive  for  
>>>> efficiency in the coupled ocean-ice models. Plus there is a  
>>>> hibler- type model on finite elements which will become the new  
>>>> "state of the  art"!)
>>>>
>>>> I don't know if you get the CVS-email notification (and if you  
>>>> do,  pay any attention to it), but I have now written and checked  
>>>> in new  versions of dynsolver, lrs and ostres. They are called   
>>>> seaice_dynsolver.F, seaice_lsr.F, and seaice_ocean_stress.F. They  
>>>> can  be called as an alternative to the original code (which is  
>>>> not  original anymore, since I did all these modifications to the  
>>>> grid  information). The main feature of the new code is that it  is 
>>>> on a C- grid for a cartesian and a spherical grid (still NOT  
>>>> fully  curvilinear). I did that because I want to run the seaice  
>>>> model on a  cubed sphere, where the ice is allowed to move across  
>>>> into other  faces of the cube (which was not possible with the  
>>>> "old version").  Along with this I made a few other changes:
>>>> 1. formulated the ice stress in a way that I find consistent:
>>>>     tau_air-ice = cdair*sqrt([wind]^2)*wind (without averaging  
>>>> over  the entire grid cell)
>>>>     tau_ice-ocean = cdiceocean*sqrt([rel vel]^2)*[rel vel]
>>>>     the ocean feels a weighted average (according to ice cover)  
>>>> of  ice-ocean stress and air-ocean stress
>>>> Because uice and u (vice and v) are now on the same grid there is  
>>>> no  averaging involved (between C and B-grid variables) except  for 
>>>> the  Coriols terms and terms involving this turning angle  (which 
>>>> is zero  by default); and of course, the wind field is on  tracer 
>>>> points, so  there some averaging there, but much less  compared to 
>>>> before!
>>>> 2. I have tried masking *NOT* the velocity field, but the rhs of  
>>>> the  solver, that is, URT and VRT with AREA. That way the  
>>>> velocities are  not discontinuous at all. This solution runs  
>>>> stably (at least for  100yrs in a 2x2 degree run and the coarse  
>>>> cubed sphere with 200km  resolution) and I can even get rid of  the 
>>>> velocity clipping in  seaice_dynsolver.F. What do you think  of 
>>>> that? In fact the solutions  are not that different, but it works
>>>> 3. I have substituted the implementation of the tilt by the real  
>>>> tile  (-gravity*gradient(ssh)).
>>>> 4. I have moved parts of the repeated computations in  
>>>> seaice_dysolver  of E11, etc, FORCEX/Y etc into seaice_lsr.F in  
>>>> order to save code  lines and simplify potential future  
>>>> modifications. I realize that  seaice_lsr is now much more than  
>>>> just the lsr solver, but it was more  than that anyway.
>>>> 5. small changes in handling the parameters (moved previously  
>>>> hard  code values into name lists)
>>>>
>>>> As for lsr.F, as long as the new seaice_lsr.F works so well for  
>>>> me, I  am fine with the currently checked-in version of lsr.F.
>>>>
>>>> Martin
>>>>
>>>> On Mar 8, 2006, at 2:20 AM, Jinlun Zhang wrote:
>>>>
>>>>> Hi Martin,
>>>>>
>>>>> Sorry for the delayed reply for your emails. Cool you have been   
>>>>> able to run the version with new indexing and other modifications.
>>>>>
>>>>> As for the calculation of surface stress, as Dimitris  mentioned,  
>>>>> the one implemented right now is a temporary fix. As  a temporary  
>>>>> fix, I prefer the existing treatment. I have been  trying to find 
>>>>> a  permanent solution that would get real stress  at the 
>>>>> ice-ocean  interface into the system, but it is tough to  do.  I 
>>>>> don't know how  AWI ice models calculate stress. I  remember one 
>>>>> AWI paper published  some years back says that some  spatial 
>>>>> average is needed in  calc'ing stress, otherwise, the  code would 
>>>>> blow up. If it is true,  then it is the same problem  we have 
>>>>> here, but I would not want to  use spatial average. By  the way, 
>>>>> would you mind checking AWI ice  models to see if they  use the 
>>>>> LSR solver or use the  Hibler (1979)  solver? I sent out  LSR 
>>>>> model (rectangular version) to some European  groups more  than 10 
>>>>> years ago, I hope it did not lost in  translation (-:.  And I hope 
>>>>> they can use a parallel LSR, if not  already using  some other 
>>>>> parallel code.
>>>>>
>>>>> As for the mask in the LSR solver, better not use AREA=0 as  
>>>>> mask.  This is because it would create velocity discontinuities  
>>>>> between  areas with and without ice.
>>>>>
>>>>> I haven't got deep into the new LSR routine, since you have run  
>>>>> the  code for many years, it is probably ok except that the  
>>>>> metric terms  have not been included, as you mentioned.
>>>>>
>>>>> Cheers, Jinlun
>>>>>
>>>>>
>>>>>
>>>>> Martin Losch wrote:
>>>>>
>>>>>> Hello Dimitris,
>>>>>>
>>>>>> FORCEX/Y at that point are air-ice-stress interactions, the    
>>>>>> coefficient ist call DAIRN. The ice-ocean stresses are put in   
>>>>>> later  after the computation ETA,ZETA,Eij with the coefficent   
>>>>>> DWATN. So I  think that at this point averaging does not make  
>>>>>> much  sense: It would  mean that the ice is driven both by wind  
>>>>>> over ice  and wind over open  water.
>>>>>>
>>>>>> I think in the c-grid version I'll try to do it the way I find  
>>>>>> it   intuitive: separate wind stress over ocean and ice and  
>>>>>> compute  stress  on ocean as the average of ice-ocean stress  and 
>>>>>> wind  stress (on the  ocean). This is how it is done in the  
>>>>>> various  coupled AWI ice models.
>>>>>>
>>>>>> Martin
>>>>>>
>>>>>> On Mar 6, 2006, at 7:54 PM, Dimitris Menemenlis wrote:
>>>>>>
>>>>>>>> The way I understand this is that the stress on the ice is  
>>>>>>>> an   average over ocean stress and ice stress ( the statement  
>>>>>>>> before   FORCEX(I,J,bi,bj) = ...), whereas the stress over  the 
>>>>>>>> ocean is   just the ocean stress not weighted at all. I  find 
>>>>>>>> this a little   inconsistent if not wrong. I would think  that 
>>>>>>>> these
>>>>>>>>  terms should be treated separately, with only the ice  
>>>>>>>> stress   driving the
>>>>>>>> ice, right?
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> Martin, yes it is inconsistent.  As Jinlun mentioned in an   
>>>>>>> earlier  message, the
>>>>>>> original formulation for ocean stress, the one marked by CPP flag
>>>>>>> SEAICE_ORIGINAL_BAD_ICE_STRESS in ostres.F caused model    
>>>>>>> instabilities.  As a
>>>>>>> temporary fix, the presence of ice is ignored in the  
>>>>>>> computation  of  ocean
>>>>>>> surface stress, variables WINDX and WINDY.
>>>>>>>
>>>>>>> Regarding ocean stress at bottom of ice, variables FORCEX  and   
>>>>>>> FORCEY, my understanding is that dynsolver assumes that  thin 
>>>>>>> ice   covers the open ocean everywhere, hence the  weighted sum 
>>>>>>> of ice-  covered and ice-free components in the  computation of 
>>>>>>> FORCEX and   FORCEY.  But I do not know whether  and why this 
>>>>>>> thin-ice  assumption  is required nor what would  be impact of 
>>>>>>> setting the  URT/VRT mask to  zero where AREA=0.
>>>>>>>
>>>>>>> Dimitris
>>>>>>