[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