[MITgcm-devel] more seaice
Martin Losch
mlosch at awi-bremerhaven.de
Fri Mar 10 03:14:40 EST 2006
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
More information about the MITgcm-devel
mailing list