[MITgcm-support] sea ice diagnostics
Georgy Manucharayan
gmanuch at caltech.edu
Fri Sep 9 12:39:45 EDT 2016
Hi Martin,
I am using the default (LSR) solver and code version c65e. I looked over
the code and found that seaice_diagnostics_state.F outputs the available
diagnostic SIsigI and SIsigII for principle components (code variables
are sigI and sigII). It does so by first calculating the stress
components sig1,sig2,sig12 which I managed to add to the sea ice
diagnostics output. The file seaice_calc_lhs.F has a calculation of
stress divergence (variables stressDivergenceX,stressDivergenceY) but I
wasn't able to
simply do call diagnostics_fill to output it from the main loop of
seaice_diagnostics_state.F. Instead, I recalculated stress divergence
in seaice_diagnostics_state.F using the same code as in
seaice_calc_lhs.F. I defined new variables sigma11,sigma22,sigma12 that
store the stress tensor elements so that I can call diagnostics_fill on
them. Here is the part of the code that I modified, can you please check
if this is correct?
In file seaice_diagnostics_state.F :
...
#ifdef SEAICE_CGRID
C adding new variables to store data
_RL sigma11 (1-oLx:sNx+oLx,1-oLy:sNy+oLy)
_RL sigma12 (1-oLx:sNx+oLx,1-oLy:sNy+oLy)
_RL sigma22 (1-oLx:sNx+oLx,1-oLy:sNy+oLy)
_RL dsigmaX (1-oLx:sNx+oLx,1-oLy:sNy+oLy)
_RL dsigmaY (1-oLx:sNx+oLx,1-oLy:sNy+oLy)
....
#endif /* SEAICE_ALLOW_EVP */
C entering the case of default seaice dynamics
C recompute strainrates from up-to-date velocities
CALL SEAICE_CALC_STRAINRATES(
I uIce, vIce,
O e11, e22, e12,
I 0, myTime, myIter, myThid )
C but use old viscosities and pressure for the
C principle stress components
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO j=1,sNy
DO i=1,sNx
sig1 = 2.*zeta(I,J,bi,bj)
& * (e11(I,J,bi,bj) + e22(I,J,bi,bj))
& - press(I,J,bi,bj)
sig2 = 2.* eta(I,J,bi,bj)
& * (e11(I,J,bi,bj) - e22(I,J,bi,bj))
sig12 = 2.*eta(I,J,bi,bj) * 0.25 _d 0 *
& ( e12(I ,J ,bi,bj) + e12(I+1,J ,bi,bj)
& + e12(I ,J+1,bi,bj) + e12(I+1,J+1,bi,bj) )
sigTmp = SQRT( sig2*sig2 + 4. _d 0*sig12*sig12 )
C Saving stress components for output
sigma11(I,J)=sig1
sigma12(I,J)=sig12
sigma22(I,J)=sig2
recip_prs = 0. _d 0
IF ( press(I,J,bi,bj) .GT. 1. _d -13 )
& recip_prs = 1./press(I,J,bi,bj)
sigI (I,J) = 0.5*(sig1 + sigTmp)*recip_prs
sigII(I,J) = 0.5*(sig1 - sigTmp)*recip_prs
ENDDO
ENDDO
CALL DIAGNOSTICS_FILL(sigI ,'SIsigI ',0,1,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(sigII,'SIsigII ',0,1,2,bi,bj,myThid)
C output of stress tensor
CALL DIAGNOSTICS_FILL(sigma11,'SIsigma11',0,1,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(sigma22,'SIsigma22',0,1,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(sigma12,'SIsigma12',0,1,2,bi,bj,myThid)
C compute divergence of stress tensor
C
DO J=1,sNy
DO I=1,sNx
dsigmaX(I,J) =
& ( sigma11(I ,J ) * _dyF(I ,J,bi,bj)
& - sigma11(I-1,J ) * _dyF(I-1,J,bi,bj)
& + sigma12(I ,J+1) * _dxV(I,J+1,bi,bj)
& - sigma12(I ,J ) * _dxV(I,J ,bi,bj)
& ) * recip_rAw(I,J,bi,bj)
dsigmaY(I,J) =
& ( sigma22(I ,J ) * _dxF(I,J ,bi,bj)
& - sigma22(I ,J-1) * _dxF(I,J-1,bi,bj)
& + sigma12(I+1,J ) * _dyU(I+1,J,bi,bj)
& - sigma12(I ,J ) * _dyU(I ,J,bi,bj)
& ) * recip_rAs(I,J,bi,bj)
ENDDO
C output of stress divergence components
CALL DIAGNOSTICS_FILL(dsigmaX,'SIdsigmaX',0,1,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(dsigmaY,'SIdsigmaY',0,1,2,bi,bj,myThid)
Also in file seaice_diagnostics_init.F I added corresponding definitions
for the new fields SIdsigmaX/Y, sigma11/12/22.
These modifications do compile and output the data. One problem with the
stress divergence output is that it shows discontinuities at the
boundaries of the core sub-domains. Am I doing something wrong here?
Also, I was able to output FORCEX and FORCEY but I am not sure what is
its physical meaning.
Thanks again,
Georgy
On 09/09/2016 07:10 AM, mitgcm-support-request at mitgcm.org wrote:
> ------------------------------ Message: 2 Date: Thu, 8 Sep 2016
> 17:08:00 -0600 From: Martin Losch <martin.losch at awi.de> To: MITgcm
> Support <mitgcm-support at mitgcm.org> Subject: Re: [MITgcm-support] sea
> ice diagnostics Message-ID:
> <2CF4D43A-A7C0-464D-AC38-885427F3BA35 at awi.de> Content-Type:
> text/plain; charset="utf-8" Hi George, all this depends a little on
> which solution technique you?ll use. ?F? including stress divergence
> is never explicitly computed for the default ?LSR? solver, but for the
> ?EVP? solver and even for the JFNK solver (which you will probably not
> want to use, because it?s the most expensive one). If you want to
> avoid the dependence on solvers, you?ll have to recompute F before
> saving it (reusing code to do so). Which solver are you using? Martin
>> >On 08 Sep 2016, at 10:41, Georgy Manucharayan<gmanuch at caltech.edu> wrote:
>> >
>> >Hi all,
>> >
>> >I'm interested in outputting the dynamic/rheology diagnostics for the sea ice package; in particular, the force F in the ice momentum equations. As I understand one can calculate it as a divergence of stress tensor, in which case I'd like to output the stress tensor components (sigma11,sigma12,sigma22). The default outputs only its principle components (sigI, sigII) and I'm not sure if this is sufficient to calculate the force vector. Eqns. (5,6) in Hibler 1979 paper provide a recipe for calculating F based on strain rates, viscosities, and pressure; however, I would like to have an output of F exactly as it is used in the model (i.e. consistent with the used numerical scheme).
>> >
>> >To create the user defined diagnostic I need to know which variable names in the code denote F and sigma (and also a term -mg gradH if there is one). I guess that those might be FORCEX, FORCEY, sig1,sig2,sig12 but can someone confirm this?
>> >
>> >Also is SIpress the same as P (pressure) in ice momentum equations?
>> >
>> >Thanks,
>> >Georgy
>> >
>> >--
>> >Georgy Manucharyan,
>> >Stanback Postdoctoral Scholar,
>> >Environmental Science and Engineering,
>> >California Institute of Technology,
>> >MS 131-24, 1200 California Blvd,
>> >Pasadena, CA, 91125
>> >
>> >Office: L+R, Room 226; Tel: (626) 395-8715
>> >Web:http://www.its.caltech.edu/~gmanuch
>> >
>> >
>> >
>> >_______________________________________________
>> >MITgcm-support mailing list
>> >MITgcm-support at mitgcm.org
>> >http://mitgcm.org/mailman/listinfo/mitgcm-support
--
Georgy Manucharyan,
Stanback Postdoctoral Scholar,
Environmental Science and Engineering,
California Institute of Technology,
MS 131-24, 1200 California Blvd,
Pasadena, CA, 91125
Office: L+R, Room 226; Tel: (626) 395-8715
Web: http://www.its.caltech.edu/~gmanuch
More information about the MITgcm-support
mailing list