[MITgcm-devel] (not so) funny things happen in seaice_lsr and pickups

[Jinlun Zhang]

Martin Losch wrote:
> The question marathon continues ...
> Jinlun, my observation is, that the metric terms cause all the trouble 
> (no problems so far in cubed sphere integration of any resolution, nor 
> on cartesian grid).  I have the suspicion that the symmetry of the 
> laplacian operator is broken, because the metric terms are distributed 
> unevenly between left and right hand side of the u and v equations. In 
> order to ensure the diagonal dominance of the coefficient matrix, 
> would it be possible to move all metric terms to the rhs? Just 
> wondering ...
Hi Martin, equations governing sea ice physics would be incorrect if 
metric terms are moved to rhs. I doubt metric terms would break  
symmetry under any conditions (I have not seen that). This is because a 
rotation of a grid (or a change of coordinate system) would not change 
the properties of a matrix.
>
> BTW, the system seems to be stable when I have a non-zero ZMIN.
Non-zero ZMIN is not right and should not be used. With non-zero ZMIN, 
the system is heavily damped. I doubt non-symmetry caused trouble would 
be fixed by non-zero ZMIN, it would eventually blow up no matter what 
ZMIN is.

Jinlun
>
> Martin
>
> On Mar 9, 2009, at 5:54 PM, Jinlun Zhang wrote:
>
>> Hi Martin,
>>
>> As you said, for B-grid, the value does not matter because of the 
>> non-slip bc's. I am not sure about C-grid non-slip bc's. If you need 
>> that value at the boundary for C-grid, it seems that your approach 
>> probably makes better sense.
>>
>> Jinlun
>>
>> Martin Losch wrote:
>>> Hi again,
>>>
>>> I have another question for Jinlun. I am now quite confident, that I 
>>> have found all (serious) problems in the metric terms. However, my 
>>> code still blows up. Could the discretization of \eta be a problem? 
>>> On Z-points ("bottom left" corner of C-grid cell, on a B-grid this 
>>> is where the velocities are), I average eta like this:
>>>        DO j=1,sNy+1
>>>         DO i=1,sNx+1
>>>          etaMeanZ (I,J,bi,bj) =
>>>     &         (           ETA (I,J  ,bi,bj)  + ETA (I-1,J  ,bi,bj)
>>>     &         +           ETA (I,J-1,bi,bj)  + ETA (I-1,J-1,bi,bj) )
>>>     &         / MAX(1.D0,maskC(I,J,  k,bi,bj)+maskC(I-1,J,  k,bi,bj)
>>>     &         +          maskC(I,J-1,k,bi,bj)+maskC(I-1,J-1,k,bi,bj) )
>>>         ENDDO
>>>        ENDDO
>>> In the B-grid code I find this:
>>>        DO j=1-Oly+1,sNy+Oly
>>>         DO i=1-Olx+1,sNx+Olx
>>>          ETAMEAN(I,J,bi,bj) =QUART*(
>>>     &          ETA(I,J-1,bi,bj) + ETA(I-1,J-1,bi,bj)
>>>     &         +ETA(I,J  ,bi,bj) + ETA(I-1,J  ,bi,bj))
>>>          ZETAMEAN(I,J,bi,bj)=QUART*(
>>>     &          ZETA(I,J-1,bi,bj) + ZETA(I-1,J-1,bi,bj)
>>>     &         +ZETA(I,J  ,bi,bj) + ZETA(I-1,J  ,bi,bj))
>>>         ENDDO
>>>        ENDDO
>>> This means that \eta/\zeta on the boundaries are only half as large 
>>> as in my C-grid discretization (however, these values may actually 
>>> never be used because of the no-slip BCs on the B-grid). So I am 
>>> wondering, what the "correct" value of \eta on the boundary is. I 
>>> observe, that code with
>>>        DO j=1,sNy+1
>>>         DO i=1,sNx+1
>>>          etaMeanZ (I,J,bi,bj) = 0.25 _d 0 *
>>>     &         (           ETA (I,J  ,bi,bj)  + ETA (I-1,J  ,bi,bj)
>>>     &         +           ETA (I,J-1,bi,bj)  + ETA (I-1,J-1,bi,bj) )
>>>         ENDDO
>>>        ENDDO
>>> tends to be stable (at least in the shorter test that I made). 
>>> Jinlun, what's your opinion on this?
>>>
>>> Martin
> _______________________________________________
> 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