[MITgcm-devel] seaice
Martin Losch
mlosch at awi-bremerhaven.de
Tue Feb 14 03:21:53 EST 2006
Hi Jinlun and others,
I do not know exactly, what the "STANDARD CONSERVATIVE ADVECTION"
scheme is, but to me it looks very much like a 2nd order central
difference scheme of this expression:
(1/(R*cos(phi))*(d (u*H) /dlambda) + (1/(R*cos(phi))*(d (v*H*cos
(phi))/dphi),
except, that first term R*dlamda=dxtice, R*dphi=dytice and cos(phi)
=cstice are not evaluated at the ij points I had expected. I would
have expected this:
dxtice(I+1,J+1)*cstice(I+1,J+1) and dytice(I+1,J+1)*cstice(I+1,J+1)
Instead we have
dxtice(I+1,J)*cstice(I,J+1) and dytice(I,J+1)*cstice(I,J+1)
for the lab_sea experiment in spherical coordinates this does make
any difference because cstice is at the correct latitude and dxtice,
dytice are uniform, but it breaks a cubed_sphere setup.
My formulation (sum of all fluxes through the cell "interfaces"
divided by the cell area) should be completely equivalent to the
first (and I think correct) version:
> CML HEFF(I,J,1,bi,bj)=HEFF(I,J,K3,bi,bj)
> CML & -DELTT * QUART * (
> CML & + _dyG(I+1,J,bi,bj) *
> CML & (HEFF(I ,J ,2,bi,bj)+HEFF(I+1,J ,2,bi,bj))
> CML & *( UI(I+1,J+1, bi,bj)+ UI(I+1,J , bi,bj))
> CML & - _dyG(I,J,bi,bj) *
> CML & (HEFF(I ,J ,2,bi,bj)+HEFF(I-1,J ,2,bi,bj))
> CML & *( UI(I ,J+1, bi,bj)+ UI(I ,J , bi,bj))
> CML & + _dxG(I ,J+1,bi,bj) *
> CML & (HEFF(I ,J ,2,bi,bj)+HEFF(I ,J+1,2,bi,bj))
> CML & *( VI(I ,J+1, bi,bj)+ VI(I+1,J+1, bi,bj))
> CML & - _dxG(I ,J ,bi,bj)*
> CML & (HEFF(I ,J ,2,bi,bj)+HEFF(I ,J-1,2,bi,bj))
> CML & *( VI(I ,J , bi,bj)+ VI(I+1,J , bi,bj))
> CML & )*recip_rA(I,J,bi,bj)
(note that I have reorder the idices here: I+1 -> I, J+1 -> 1 and
changed the loop boundaries to 1, sNx etc), but that doesn't make any
difference difference, I have done that to the current version as well)
You can check by substituting _dxG = dxtice*cstice, dyG = dytice and
recip_rA=1/(_dxG*_dyG), that is is exactly the same as the current
version (for a lat-lon grid dyG=dyF and also dxG=dyF if yC(J) = 0.5*
(yG(J)+yG(J+1)) ).
So the differences I see in the lab_sea experiment are only due to
round off (or bugs introduced by me (o:). In the cubed_sphere
configuration, the two discretiztions are different, I agree
(although I would prefer my discritization, because it is consistent
with other MITgcm code).
But you can see how sensitive the system is by considering the other
two substitution-examples that change the results:
fCoriG is exactly the same as TWO*OMEGA*SINEICE, but gives different
results
rSphere = 6370. _d 3 (in ini_defaults.F) but RADIUS = rSphere changes
the results. (and I have not even used recip_rSphere instead of 1/
RADIUS)
Personally I like "my" (it's really stolen form other parts of the
MITgcm-code) discretization better, but if it changes the results ...
Martin
PS. all of this will not help me with my Weddell Sea problems.
On Feb 14, 2006, at 1:33 AM, Jinlun Zhang wrote:
> The sea ice grid info is derived from ocean grid, so I guess the
> difference must be from the different advection schemes.
> Jinlun
>
> Martin Losch wrote:
>
>> Dimitris and other seaice modelers,
>>
>> I have modified
>> dynsolver.F
>> diffus.F
>> advect.F
>> to replace all grid information (dxtice,dytice,etc) with generic
>> MITgcm grid parameters (dxF,dyF, etc, defined in GRID.h)
>> So far, this has NOT changed the results of lab_sea experiment,
>> except for the advection part.
>> I tried to use a flux-form, similar to the second-order central
>> difference form in pkg/generic_advection/gad_c2_adv_x.F etc. but
>> there are small differences (the model needs a better advection
>> scheme that is positive!!!).
>>
>> Unfortunately on a cubed_sphere grid with seaice, the results do
>> change a lot, but is that surprising?
>>
>> What also changes the results:
>> in dynsolver.F:
>> RADIUS = rSphere (default = 6370. _d 3)
>> TWO*OMEGA*SINEICE(I,J,bi,bj) = _fCoriG(I,J,bi,bj)
>> (=2._d0*omega*sin (_yG(i,j,bi,bj)*deg2rad)
>>
>> I don't see why.
>>
>> how should I proceed?:
>> I could send the routines around for testing or looking at and we
>> could agree on checking them in or not checking them in.
>> Or I could forget about the project altogether, because it is
>> tedious and unrewarding. Is it worth it? Can you deal with the
>> small differences that these changes imply?
>>
>> Martin
>>
>>
>> On Feb 13, 2006, at 3:37 PM, Dimitris Menemenlis wrote:
>>
>>> Martin, I am attaching a draft wish list for MITgcm sea-ice code
>>> for you and others to consider and modify. The list is based on
>>> your recent e-mails and work on ice shelves as well as on
>>> discussions with Jean-Michel and Patrick during ECCO meeting two
>>> weeks ago and with Elizabeth Hunke and Bill Lipscomb at NCAR
>>> last week. I am not sure when these things will happen and in
>>> what order. I am presently trying to hire someone that can help
>>> move the MITgcm sea-ice code forward. Timing also depends on
>>> Jinlun's timetable and availability.
>>>
>>> 1. Clean-up grid parameters as you suggest. Make sure that
>>> there are no duplicate or repeated computations. Use model/inc/
>>> GRID.h parameters where possible, define the rest in pkg/seaice/
>>> SEAICE_GRID.h and seaice_init.F.
>>>
>>> 2. Adapt LSR solver for curvilinear coordinates so that it is
>>> "properly" compatible with cubed-sphere and terraced grids.
>>>
>>> 3. Port LSR solver to C-grid. This is something that Jinlun is
>>> working on for GISS, and he may be able to spend two months next
>>> year to implement on MITgcm.
>>>
>>> 4. Better coupling of pkg/seaice thermodynamics with surface
>>> level, similar to what Jean-Michel has done for pkg/thsice.
>>>
>>> 5. Coupling of pkg/thsice thermodynamics with pkg/seaice dynamics.
>>>
>>> 6. Ice shelves and cavern representation.
>>>
>>> 7. Import more modern, multi-category sea-ice model. The plan
>>> is to import the Los Alamos CICE sea ice model (http://
>>> climate.lanl.gov/Models/CICE/), ideally with interchangeable
>>> dynamics between Jinlun's LSR VP solver and Elizabeth's explicit
>>> EVP solver.
>>>
>>> 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
>
>
> --
>
> Jinlun Zhang
> Polar Science Center, Applied Physics Laboratory
> University of Washington, 1013 NE 40th St, Seattle, WA 98105-6698
>
> Phone: (206)-543-5569; Fax: (206)-616-3142
> zhang at apl.washington.edu
> http://psc.apl.washington.edu/pscweb2002/Staff/zhang/zhang.html
>
>
>
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel
More information about the MITgcm-devel
mailing list