[MITgcm-devel] more seaice

Wed Mar 8 03:03:43 EST 2006

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