C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_diagnostics_state.F,v 1.34 2017/05/23 16:23:46 mlosch Exp $
C $Name: $
#include "SEAICE_OPTIONS.h"
CBOP
C !ROUTINE: SEAICE_DIAGNOSTICS_STATE
C !INTERFACE:
SUBROUTINE SEAICE_DIAGNOSTICS_STATE(
I myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | S/R SEAICE_DIAGNOSTICS_STATE
C | o fill-in diagnostics array for SEAICE state variables
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "SEAICE_SIZE.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE.h"
#include "SEAICE_TRACER.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine Arguments ==
C myTime :: time counter for this thread
C myIter :: iteration counter for this thread
C bi,bj :: tile indices
C myThid :: thread number for this instance of the routine.
_RL myTime
INTEGER myIter
INTEGER myThid
CEOP
#ifdef ALLOW_DIAGNOSTICS
C == Local variables ==
INTEGER i,j,bi,bj
_RL tIce (1-oLx:sNx+oLx,1-oLy:sNy+oLy,nSx,nSy)
#ifdef SEAICE_CGRID
_RL sig1(1-oLx:sNx+oLx,1-oLy:sNy+oLy)
_RL sig2(1-oLx:sNx+oLx,1-oLy:sNy+oLy)
_RL sig11(1-oLx:sNx+oLx,1-oLy:sNy+oLy)
_RL sig22(1-oLx:sNx+oLx,1-oLy:sNy+oLy)
_RL sig12(1-oLx:sNx+oLx,1-oLy:sNy+oLy)
_RL sigp, sigm, sigTmp, recip_prs
#endif
INTEGER k
_RL recip_multDim
_RL tmp
#ifdef ALLOW_SITRACER
INTEGER iTracer
CHARACTER*8 diagName
#endif
LOGICAL DIAGNOSTICS_IS_ON
EXTERNAL
IF ( useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL( AREA, 'SIarea ', 0, 1, 0, 1,1, myThid )
CALL DIAGNOSTICS_FILL( HEFF, 'SIheff ', 0, 1, 0, 1,1, myThid )
CALL DIAGNOSTICS_FILL( UICE, 'SIuice ', 0, 1, 0, 1,1, myThid )
CALL DIAGNOSTICS_FILL( VICE, 'SIvice ', 0, 1, 0, 1,1, myThid )
IF ( DIAGNOSTICS_IS_ON('SItices ',myThid) ) THEN
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
C use tIce as a temporary field, as it is done in seaice_growth
DO j=1,sNy
DO i=1,sNx
tIce(I,J,bi,bj) = 0. _d 0
ENDDO
ENDDO
C division by zero is not possible
recip_multDim = 1. _d 0/SEAICE_multDim
DO k=1,SEAICE_multDim
DO j=1,sNy
DO i=1,sNx
tmp = 1. _d 0
#ifdef SEAICE_ITD
IF (AREA(I,J,bi,bj) .GT. ZERO)
& tmp=AREAITD(I,J,K,bi,bj)/AREA(I,J,bi,bj)
#endif /* SEAICE_ITD */
tIce(I,J,bi,bj) = tIce(I,J,bi,bj)
& + TICES(I,J,K,bi,bj)*tmp*recip_multDim
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
CALL DIAGNOSTICS_FRACT_FILL(
I tIce, AREA, 1. _d 0, 1, 'SItices ',
I 0, 1, 0, 1, 1, myThid )
ENDIF
CALL DIAGNOSTICS_FILL(HSNOW ,'SIhsnow ',0,1 ,0,1,1,myThid)
#ifdef SEAICE_ITD
CALL DIAGNOSTICS_FILL(HEFFITD ,'SIheffN ',0,nITD ,0,1,1,myThid)
CALL DIAGNOSTICS_FILL(AREAITD ,'SIareaN ',0,nITD ,0,1,1,myThid)
CALL DIAGNOSTICS_FILL(HSNOWITD,'SIhsnowN',0,nITD ,0,1,1,myThid)
#endif
#ifdef ALLOW_SITRACER
DO iTracer = 1, SItrNumInUse
WRITE(diagName,'(A4,I2.2,A2)') 'SItr',iTracer,' '
if (SItrMate(iTracer).EQ.'HEFF') then
CALL DIAGNOSTICS_FRACT_FILL(
I SItracer(1-OLx,1-OLy,1,1,iTracer), HEFF,
I ONE, 1, diagName, 0, 1, 0, 1, 1, myThid )
else
CALL DIAGNOSTICS_FRACT_FILL(
I SItracer(1-OLx,1-OLy,1,1,iTracer), AREA,
I ONE, 1, diagName, 0, 1, 0, 1, 1, myThid )
endif
ENDDO
#endif
#ifdef SEAICE_VARIABLE_SALINITY
CALL DIAGNOSTICS_FILL(HSALT ,'SIhsalt ',0,1 ,0,1,1,myThid)
#endif
CALL DIAGNOSTICS_FILL(zeta ,'SIzeta ',0,1 ,0,1,1,myThid)
CALL DIAGNOSTICS_FILL(eta ,'SIeta ',0,1 ,0,1,1,myThid)
CALL DIAGNOSTICS_FILL(press ,'SIpress ',0,1 ,0,1,1,myThid)
CALL DIAGNOSTICS_FILL(deltaC ,'SIdelta ',0,1 ,0,1,1,myThid)
#ifdef SEAICE_CGRID
IF ( DIAGNOSTICS_IS_ON('SItensil',myThid) ) THEN
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
C use tIce as a temporary field
DO j=1,sNy
DO i=1,sNx
IF ( tensileStrFac(I,J,bi,bj) .EQ. 1. _d 0) THEN
C This special case of tensile strength equal to compressive strength
C is not very physical and should actually not happen but you never know;
C in this case, press = P-T = P*(1-k) = 0. and we have to use press0 to
C get something
tIce(I,J,bi,bj) = press0(I,J,bi,bj)
ELSE
C This is more complicated than you think because press = P-T = P*(1-k),
C but we are looking for T = k*P = k*press/(1-k)
tIce(I,J,bi,bj) = tensileStrFac(I,J,bi,bj)
& *press(I,J,bi,bj)/(1. _d 0 - tensileStrFac(I,J,bi,bj))
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
CALL DIAGNOSTICS_FILL(tIce,'SItensil',0,1 ,0,1,1,myThid)
ENDIF
IF ( DIAGNOSTICS_IS_ON('SIsig1 ',myThid) .OR.
& DIAGNOSTICS_IS_ON('SIsig2 ',myThid) .OR.
& DIAGNOSTICS_IS_ON('SIshear ',myThid) ) THEN
#ifdef SEAICE_ALLOW_EVP
IF ( SEAICEuseEVP ) THEN
C for EVP compute principle stress components from recent
C stress state and normalize with latest
C PRESS = PRESS(n-1), n = number of sub-cycling steps
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO j=1,sNy
DO i=1,sNx
sigp = seaice_sigma1(I,J,bi,bj)
sigm = seaice_sigma2(I,J,bi,bj)
sig12(I,J) = 0.25 _d 0 *
& ( seaice_sigma12(I, J, bi,bj)
& + seaice_sigma12(I+1,J, bi,bj)
& + seaice_sigma12(I+1,J+1,bi,bj)
& + seaice_sigma12(I ,J+1,bi,bj) )
sigTmp = SQRT( sigm*sigm + 4. _d 0*sig12(I,J)*sig12(I,J) )
recip_prs = 0. _d 0
IF ( press0(I,J,bi,bj) .GT. 1. _d -13 )
& recip_prs = 1./press0(I,J,bi,bj)
sig1(I,J) = 0.5*(sigp + sigTmp)*recip_prs
sig2(I,J) = 0.5*(sigp - sigTmp)*recip_prs
ENDDO
ENDDO
CALL DIAGNOSTICS_FILL(sig1,'SIsig1 ',0,1,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(sig2,'SIsig2 ',0,1,2,bi,bj,myThid)
ENDDO
ENDDO
ELSE
#else
IF ( .TRUE. ) THEN
#endif /* SEAICE_ALLOW_EVP */
C recompute strainrates from up-to-date velocities
CALL SEAICE_CALC_STRAINRATES(
I uIce, vIce,
O e11, e22, e12,
I 0, myTime, myIter, myThid )
CML CALL SEAICE_CALC_VISCOSITIES(
CML I e11, e22, e12, zMin, zMax, hEffM, press0, tensileStrFac,
CML O eta, etaZ, zeta, zetaZ, press, deltaC,
CML 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)
CALL SEAICE_CALC_STRESS(
I e11, e22, e12, press, zeta, eta, etaZ,
O sig11, sig22, sig12,
I bi, bj, myTime, myIter, myThid )
DO j=1,sNy
DO i=1,sNx
sigp = sig11(I,J) + sig22(I,J)
sigm = sig11(I,J) - sig22(I,J)
C This should be the way of computing sig12 at C-points,
C sigTmp = 0.25 _d 0 *
C & ( sig12(I ,J ) + sig12(I+1,J )
C & + sig12(I ,J+1) + sig12(I+1,J+1) )
C but sig12 = 2*etaZ*e12, and because of strong gradients in eta,
C etaZ can be very large for a cell with small eta and the straightforward
C way of averaging mixes large etaZ with small press0, so we have to do it
C in different way to get meaningfull sig12C (=sigTmp):
sigTmp = 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( sigm*sigm + 4.*sigTmp*sigTmp )
recip_prs = 0. _d 0
IF ( press0(I,J,bi,bj) .GT. 1. _d -13 )
& recip_prs = 1./press0(I,J,bi,bj)
sig1(I,J) = 0.5*(sigp + sigTmp)*recip_prs
sig2(I,J) = 0.5*(sigp - sigTmp)*recip_prs
ENDDO
ENDDO
CALL DIAGNOSTICS_FILL(sig1,'SIsig1 ',0,1,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(sig2,'SIsig2 ',0,1,2,bi,bj,myThid)
ENDDO
ENDDO
C endif SEAICEuseEVP
ENDIF
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO j=1,sNy
DO i=1,sNx
sigm = e11(I,J,bi,bj) - e22(I,J,bi,bj)
sigTmp =
& ( e12(I, J, bi,bj)**2 + e12(I+1,J, bi,bj)**2
& + e12(I+1,J+1,bi,bj)**2 + e12(I ,J+1,bi,bj)**2 )
C shear deformation as sqrt((e11-e22)**2 + 4*e12**2); the 4 pulled into
C the average
sig1(I,J) = sqrt(sigm*sigm + sigTmp)
ENDDO
ENDDO
CALL DIAGNOSTICS_FILL(sig1,'SIshear ',0,1,2,bi,bj,myThid)
ENDDO
ENDDO
C endif DIAGNOSTICS_IS_ON(SIsig1/2)
ENDIF
#endif /* SEAICE_CGRID */
C abuse press as a temporary field
IF ( DIAGNOSTICS_IS_ON('SIuheff ',myThid) ) THEN
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO j = 1,sNy
DO i = 1,sNx+1
press(i,j,bi,bj) =
#ifdef SEAICE_CGRID
& UICE(i,j,bi,bj)
#else
C average B-grid velocities to C-grid points
& 0.5 _d 0*(UICE(i,j,bi,bj)+UICE(i,j+1,bi,bj))
#endif /* SEAICE_CGRID */
& *0.5 _d 0*(HEFF(i,j,bi,bj)+HEFF(i-1,j,bi,bj))
ENDDO
ENDDO
ENDDO
ENDDO
CALL DIAGNOSTICS_FILL(press,'SIuheff ',0,1,0,1,1,myThid)
ENDIF
IF ( DIAGNOSTICS_IS_ON('SIvheff ',myThid) ) THEN
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO j = 1,sNy+1
DO i = 1,sNx
press(i,j,bi,bj) =
#ifdef SEAICE_CGRID
& VICE(i,j,bi,bj)
#else
C average B-grid velocities to C-grid points
& 0.5 _d 0*(VICE(i,j,bi,bj)+VICE(i+1,j,bi,bj))
#endif /* SEAICE_CGRID */
& *0.5 _d 0*(HEFF(i,j,bi,bj)+HEFF(i,j-1,bi,bj))
ENDDO
ENDDO
ENDDO
ENDDO
CALL DIAGNOSTICS_FILL(press,'SIvheff ',0,1,0,1,1,myThid)
ENDIF
C endif useDiagnostics
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
RETURN
END