[MITgcm-support] sea ice diagnostics
Georgy Manucharayan
gmanuch at caltech.edu
Wed Sep 14 14:25:36 EDT 2016
The diagnostics log file does contain codes for the grids on which
variables are defined, however I don't understand how it is choosing the
grid for the new variables that were manually added. E.g., I want to
output sigma11, 22 which I believe should be defined on a different grid
from sigma12; how do I reflect this in my code? Do I need to modify any
of the DIAGNOSTICS_FILL parameters to reflect the grid of a variable?
At the moment I added this code in seaice_diagnostics_state.F that
outputs sigma11/12/22 :
_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)
DO j=0,sNy
DO i=0,sNx
eplus = e11(I,J,bi,bj) + e22(I,J,bi,bj)
eminus= e11(I,J,bi,bj) - e22(I,J,bi,bj)
sigma11(I,J) = zeta(I,J,bi,bj)*eplus + eta(I,J,bi,bj)*eminus
& - 0.5 _d 0 * PRESS(I,J,bi,bj)
sigma22(I,J) = zeta(I,J,bi,bj)*eplus - eta(I,J,bi,bj)*eminus
& - 0.5 _d 0 * PRESS(I,J,bi,bj)
ENDDO
ENDDO
DO j=1,sNy+1
DO i=1,sNx+1
sigma12(I,J) = 2. _d 0 * e12(I,J,bi,bj) * etaZ(I,J,bi,bj)
ENDDO
ENDDO
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)
Do I need to change a definition for sigma12 to_RL sigma12
(2-oLx:sNx+oLx,2-oLy:sNy+oLy) using one indexes higher by 1?
In general, is there a way to determine the grid of an arbitrary
variable in the code, e.g. stressDivergenceX and stressDivergenceY?
Thanks,
Georgy
On 09/14/2016 09:00 AM, mitgcm-support-request at mitgcm.org wrote:
> Send MITgcm-support mailing list submissions to
> mitgcm-support at mitgcm.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
> http://mitgcm.org/mailman/listinfo/mitgcm-support
> or, via email, send a message with subject or body 'help' to
> mitgcm-support-request at mitgcm.org
>
> You can reach the person managing the list at
> mitgcm-support-owner at mitgcm.org
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of MITgcm-support digest..."
>
>
> Today's Topics:
>
> 1. Re: sea ice diagnostics (Martin Losch)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Wed, 14 Sep 2016 07:27:28 +0200
> 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: <B85A2032-4E95-460E-B579-8996047CAC39 at awi.de>
> Content-Type: text/plain; charset="utf-8"
>
> Hi Gregory,
>
> the diagnostics package prints an ?available_diagnostics.log? into the run directory. That should have some ?secret? code that tells you where individual variables are stored (given that that has been properly specified at the ?init? step); somewhere (online-documentation) you can find what the secret code means. So expect some of the variables on C-points on U and V points. Usually the diagnostics are not averaged between grid points before output.
>
> Martin
>
>> On 13 Sep 2016, at 02:27, Georgy Manucharayan <gmanuch at caltech.edu> wrote:
>>
>> Thanks Martin,
>>
>> Yes, this all makes sense now! I have added appropriate sigmas into the code and hopefully now computing a proper stress divergence term. I am still a little uncertain about how to calculate the remaining forces in the momentum equations but I think those can be computed offline.
>>
>> For an offline calculation of the ice advection terms, tilt, and ocean drag terms from the diagnostic output files do I need to consider that different variables are defined on different grids? or does the diagnostic package output all variables at the center of the grid boxes?
>>
>> Georgy
>>
>>
>>
>> On 09/11/2016 09:00 AM, mitgcm-support-request at mitgcm.org wrote:
>>> Send MITgcm-support mailing list submissions to
>>> mitgcm-support at mitgcm.org
>>>
>>> To subscribe or unsubscribe via the World Wide Web, visit
>>> http://mitgcm.org/mailman/listinfo/mitgcm-support
>>> or, via email, send a message with subject or body 'help' to
>>> mitgcm-support-request at mitgcm.org
>>>
>>> You can reach the person managing the list at
>>> mitgcm-support-owner at mitgcm.org
>>>
>>> When replying, please edit your Subject line so it is more specific
>>> than "Re: Contents of MITgcm-support digest..."
>>>
>>>
>>> Today's Topics:
>>>
>>> 1. Re: sea ice diagnostics (Martin Losch)
>>>
>>>
>>> ----------------------------------------------------------------------
>>>
>>> Message: 1
>>> Date: Sat, 10 Sep 2016 23:34:11 -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: <7F7F691E-1DAC-4E94-B86F-11036AA961CA at awi.de>
>>> Content-Type: text/plain; charset="utf-8"
>>>
>>> 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
>>>
>>>
>>> ------------------------------
>>>
>>> _______________________________________________
>>> MITgcm-support mailing list
>>> MITgcm-support at mitgcm.org
>>> http://mitgcm.org/mailman/listinfo/mitgcm-support
>>>
>>>
>>> End of MITgcm-support Digest, Vol 159, Issue 11
>>> ***********************************************
>> --
>> 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
>
>
>
> ------------------------------
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support
>
>
> End of MITgcm-support Digest, Vol 159, Issue 17
> ***********************************************
--
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