[MITgcm-support] sea ice diagnostics
Martin Losch
martin.losch at awi.de
Sun Sep 11 01:34:11 EDT 2016
Hi George,
all of this is a little tricky and my comments in the code may not be very good (or even wrong). None of the quantities that you are interested in are explicitly computed in seaice_lsr, so you’ll have to do this in the diagnostics routine, as you have started to do.
The proper nomenclature is probably this: the “cartesian” stress tensor components are sigma11, sigma22, sigma12. In the code I often use sig1=sig11+sig22 and sig2=sig11-sig22. These are just abbreviations. You can use the code in seaice_calc_lhs.F to compute the tensor components and the stress divergence, lines 135-170. Your code computes sig1 and sig2, which is not quite what you want. (sigma11 = 0.5*(sig1+sig2), sigma22=0.5*(sig1-sig2))
The lines in the diagnostics fields may actually not have anything to do with your code, but with the linear solver not converging very well. To make them go away try increasing the tolerance from default 1e-4 to 1e-5 or 1e-6 (LSR_ERROR). This will make your solutions more expensive. An alternative is to use the Krylov solver (but you need more up to date code for that) or the JFNK solver (very expensive if you are not interested in converged solutions), both of which are not sentitive to parallelisation.
FORCEX/Y contain all explicit terms of the forcing, but not the implicit part, so it does not contain grad(sigma) nor the part of the ice-ocean stress that is treated implicitly. You’ll have to add these contributions for a meaningful field, I guess, like this:
FORCEX - 0.5 _d 0 * ( DWATN(I,J,bi,bj)+DWATN(I-1,J,bi,bj) ) * COSWAT * uIce(I,J,,bi,bj) (+ dsigmaX)
Does that make sense?
Martin
> On 09 Sep 2016, at 10:39, Georgy Manucharayan <gmanuch at caltech.edu> wrote:
>
> 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
>
>
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support
More information about the MITgcm-support
mailing list