[MITgcm-devel] more seaice

Martin Losch mlosch at awi-bremerhaven.de
Fri Mar 10 03:52:32 EST 2006 

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

More information about the MITgcm-devel mailing list