[MITgcm-devel] more seaice

[Jinlun Zhang]

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