C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_calc_tensor.F,v 1.45 2015/10/15 23:06:36 gforget Exp $
C $Name: $
#include "GMREDI_OPTIONS.h"
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif
#ifdef ALLOW_CTRL
# include "CTRL_OPTIONS.h"
#endif
#ifdef ALLOW_KPP
# include "KPP_OPTIONS.h"
#endif
CBOP
C !ROUTINE: GMREDI_CALC_TENSOR
C !INTERFACE:
SUBROUTINE GMREDI_CALC_TENSOR(
I iMin, iMax, jMin, jMax,
I sigmaX, sigmaY, sigmaR,
I bi, bj, myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE GMREDI_CALC_TENSOR
C | o Calculate tensor elements for GM/Redi tensor.
C *==========================================================*
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GMREDI.h"
#include "GMREDI_TAVE.h"
#ifdef ALLOW_CTRL
# include "CTRL_FIELDS.h"
#endif
#ifdef ALLOW_KPP
# include "KPP.h"
#endif
#ifdef ALLOW_AUTODIFF_TAMC
#include "tamc.h"
#include "tamc_keys.h"
#endif /* ALLOW_AUTODIFF_TAMC */
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C bi, bj :: tile indices
C myTime :: Current time in simulation
C myIter :: Current iteration number in simulation
C myThid :: My Thread Id. number
C
INTEGER iMin,iMax,jMin,jMax
_RL sigmaX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL sigmaY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
INTEGER bi, bj
_RL myTime
INTEGER myIter
INTEGER myThid
CEOP
#ifdef ALLOW_GMREDI
C !LOCAL VARIABLES:
C == Local variables ==
INTEGER i,j,k
_RL SlopeX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL SlopeY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dSigmaDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dSigmaDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dSigmaDr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL SlopeSqr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL taperFct(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL ldd97_LrhoC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL ldd97_LrhoW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL ldd97_LrhoS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL Cspd, LrhoInf, LrhoSup, fCoriLoc
_RL Kgm_tmp, isopycK, bolus_K
INTEGER kLow_W (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
INTEGER kLow_S (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL locMixLayer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL baseSlope (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL hTransLay (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL recipLambda(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
INTEGER km1
#if ( defined (GM_NON_UNITY_DIAGONAL) defined (GM_EXTRA_DIAGONAL) )
INTEGER kp1
_RL maskp1
#endif
#ifdef GM_VISBECK_VARIABLE_K
#ifdef OLD_VISBECK_CALC
_RL Ssq(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#else
_RL dSigmaH, dSigmaR
_RL Sloc, M2loc
#endif
_RL recipMaxSlope
_RL deltaH, integrDepth
_RL N2loc, SNloc
#endif /* GM_VISBECK_VARIABLE_K */
#ifdef ALLOW_DIAGNOSTICS
LOGICAL doDiagRediFlx
LOGICAL DIAGNOSTICS_IS_ON
EXTERNAL
#if ( defined (GM_NON_UNITY_DIAGONAL) defined (GM_EXTRA_DIAGONAL) )
_RL dTdz
_RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#endif
#endif /* ALLOW_DIAGNOSTICS */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
act1 = bi - myBxLo(myThid)
max1 = myBxHi(myThid) - myBxLo(myThid) + 1
act2 = bj - myByLo(myThid)
max2 = myByHi(myThid) - myByLo(myThid) + 1
act3 = myThid - 1
max3 = nTx*nTy
act4 = ikey_dynamics - 1
igmkey = (act1 + 1) + act2*max1
& + act3*max1*max2
& + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef ALLOW_DIAGNOSTICS
doDiagRediFlx = .FALSE.
IF ( useDiagnostics ) THEN
doDiagRediFlx = DIAGNOSTICS_IS_ON('GM_KuzTz', myThid )
doDiagRediFlx = doDiagRediFlx .OR.
& DIAGNOSTICS_IS_ON('GM_KvzTz', myThid )
ENDIF
#endif
#ifdef GM_VISBECK_VARIABLE_K
recipMaxSlope = 0. _d 0
IF ( GM_Visbeck_maxSlope.GT.0. _d 0 ) THEN
recipMaxSlope = 1. _d 0 / GM_Visbeck_maxSlope
ENDIF
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
VisbeckK(i,j,bi,bj) = 0. _d 0
ENDDO
ENDDO
#endif
C-- set ldd97_Lrho (for tapering scheme ldd97):
IF ( GM_taper_scheme.EQ.'ldd97' .OR.
& GM_taper_scheme.EQ.'fm07' ) THEN
Cspd = 2. _d 0
LrhoInf = 15. _d 3
LrhoSup = 100. _d 3
C- Tracer point location (center):
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
IF (fCori(i,j,bi,bj).NE.0.) THEN
ldd97_LrhoC(i,j) = Cspd/ABS(fCori(i,j,bi,bj))
ELSE
ldd97_LrhoC(i,j) = LrhoSup
ENDIF
ldd97_LrhoC(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoC(i,j),LrhoSup))
ENDDO
ENDDO
C- U point location (West):
DO j=1-OLy,sNy+OLy
kLow_W(1-OLx,j) = 0
ldd97_LrhoW(1-OLx,j) = LrhoSup
DO i=1-OLx+1,sNx+OLx
kLow_W(i,j) = MIN(kLowC(i-1,j,bi,bj),kLowC(i,j,bi,bj))
fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))
IF (fCoriLoc.NE.0.) THEN
ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)
ELSE
ldd97_LrhoW(i,j) = LrhoSup
ENDIF
ldd97_LrhoW(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoW(i,j),LrhoSup))
ENDDO
ENDDO
C- V point location (South):
DO i=1-OLx+1,sNx+OLx
kLow_S(i,1-OLy) = 0
ldd97_LrhoS(i,1-OLy) = LrhoSup
ENDDO
DO j=1-OLy+1,sNy+OLy
DO i=1-OLx,sNx+OLx
kLow_S(i,j) = MIN(kLowC(i,j-1,bi,bj),kLowC(i,j,bi,bj))
fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))
IF (fCoriLoc.NE.0.) THEN
ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)
ELSE
ldd97_LrhoS(i,j) = LrhoSup
ENDIF
ldd97_LrhoS(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoS(i,j),LrhoSup))
ENDDO
ENDDO
ELSE
C- Just initialize to zero (not use anyway)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
ldd97_LrhoC(i,j) = 0. _d 0
ldd97_LrhoW(i,j) = 0. _d 0
ldd97_LrhoS(i,j) = 0. _d 0
ENDDO
ENDDO
ENDIF
#ifdef GM_BOLUS_ADVEC
DO k=1,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
GM_PsiX(i,j,k,bi,bj) = 0. _d 0
GM_PsiY(i,j,k,bi,bj) = 0. _d 0
ENDDO
ENDDO
ENDDO
#endif /* GM_BOLUS_ADVEC */
#ifdef ALLOW_AUTODIFF
DO k=1,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
Kwx(i,j,k,bi,bj) = 0. _d 0
Kwy(i,j,k,bi,bj) = 0. _d 0
Kwz(i,j,k,bi,bj) = 0. _d 0
# ifdef GM_NON_UNITY_DIAGONAL
Kux(i,j,k,bi,bj) = 0. _d 0
Kvy(i,j,k,bi,bj) = 0. _d 0
# endif
# ifdef GM_EXTRA_DIAGONAL
Kuz(i,j,k,bi,bj) = 0. _d 0
Kvz(i,j,k,bi,bj) = 0. _d 0
# endif
ENDDO
ENDDO
ENDDO
#endif /* ALLOW_AUTODIFF */
C-- Initialise Mixed Layer related array:
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
hTransLay(i,j) = R_low(i,j,bi,bj)
baseSlope(i,j) = 0. _d 0
recipLambda(i,j) = 0. _d 0
locMixLayer(i,j) = 0. _d 0
ENDDO
ENDDO
#ifdef ALLOW_KPP
IF ( useKPP ) THEN
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
locMixLayer(i,j) = KPPhbl(i,j,bi,bj)
ENDDO
ENDDO
ELSE
#else
IF ( .TRUE. ) THEN
#endif
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
locMixLayer(i,j) = hMixLayer(i,j,bi,bj)
ENDDO
ENDDO
ENDIF
#ifdef GM_K3D
IF (GM_useK3D) THEN
C Calculate the 3D diffusivity as per Bates et al. (2013)
CALL TIMER_START('GMREDI_K3D [GMREDI_CALC_TENSOR]',
& myThid)
CALL GMREDI_K3D(
I iMin, iMax, jMin, jMax,
I sigmaX, sigmaY, sigmaR,
I bi, bj, myTime, myThid)
CALL TIMER_STOP('GMREDI_K3D [GMREDI_CALC_TENSOR]',
& myThid)
ENDIF
#endif
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- 1rst loop on k : compute Tensor Coeff. at W points.
DO k=Nr,2,-1
#ifdef ALLOW_AUTODIFF
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
SlopeX(i,j) = 0. _d 0
SlopeY(i,j) = 0. _d 0
dSigmaDx(i,j) = 0. _d 0
dSigmaDy(i,j) = 0. _d 0
dSigmaDr(i,j) = 0. _d 0
SlopeSqr(i,j) = 0. _d 0
taperFct(i,j) = 0. _d 0
ENDDO
ENDDO
#endif /* ALLOW_AUTODIFF */
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
C Gradient of Sigma at rVel points
dSigmaDx(i,j)=op25*( sigmaX(i+1,j,k-1)+sigmaX(i,j,k-1)
& +sigmaX(i+1,j, k )+sigmaX(i,j, k )
& )*maskC(i,j,k,bi,bj)
dSigmaDy(i,j)=op25*( sigmaY(i,j+1,k-1)+sigmaY(i,j,k-1)
& +sigmaY(i,j+1, k )+sigmaY(i,j, k )
& )*maskC(i,j,k,bi,bj)
c dSigmaDr(i,j)=sigmaR(i,j,k)
ENDDO
ENDDO
#ifdef GM_VISBECK_VARIABLE_K
#ifndef OLD_VISBECK_CALC
IF ( GM_Visbeck_alpha.GT.0. .AND.
& -rC(k-1).LT.GM_Visbeck_depth ) THEN
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
dSigmaDr(i,j) = MIN( sigmaR(i,j,k), 0. _d 0 )
ENDDO
ENDDO
C-- Depth average of f/sqrt(Ri) = M^2/N^2 * N
C M^2 and N^2 are horizontal & vertical gradient of buoyancy.
C Calculate terms for mean Richardson number which is used
C in the "variable K" parameterisaton:
C compute depth average from surface down to the bottom or
C GM_Visbeck_depth, whatever is the shallower.
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
integrDepth = -rC( kLowC(i,j,bi,bj) )
C- in 2 steps to avoid mix of RS & RL type in min fct. arguments
integrDepth = MIN( integrDepth, GM_Visbeck_depth )
C- to recover "old-visbeck" form with Visbeck_minDepth = Visbeck_depth
integrDepth = MAX( integrDepth, GM_Visbeck_minDepth )
C Distance between level center above and the integration depth
deltaH = integrDepth + rC(k-1)
C If negative then we are below the integration level
C (cannot be the case with 2 conditions on maskC & -rC(k-1))
C If positive we limit this to the distance from center above
deltaH = MIN( deltaH, drC(k) )
C Now we convert deltaH to a non-dimensional fraction
deltaH = deltaH/( integrDepth+rC(1) )
C-- compute: ( M^2 * S )^1/2 (= S*N since S=M^2/N^2 )
C a 5 points average gives a more "homogeneous" formulation
C (same stencil and same weights as for dSigmaH calculation)
dSigmaR = ( dSigmaDr(i,j)*4. _d 0
& + dSigmaDr(i-1,j)
& + dSigmaDr(i+1,j)
& + dSigmaDr(i,j-1)
& + dSigmaDr(i,j+1)
& )/( 4. _d 0
& + maskC(i-1,j,k,bi,bj)
& + maskC(i+1,j,k,bi,bj)
& + maskC(i,j-1,k,bi,bj)
& + maskC(i,j+1,k,bi,bj)
& )
dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)
& + dSigmaDy(i,j)*dSigmaDy(i,j)
IF ( dSigmaH .GT. 0. _d 0 ) THEN
dSigmaH = SQRT( dSigmaH )
C- compute slope, limited by GM_Visbeck_maxSlope:
IF ( -dSigmaR.GT.dSigmaH*recipMaxSlope ) THEN
Sloc = dSigmaH / ( -dSigmaR )
ELSE
Sloc = GM_Visbeck_maxSlope
ENDIF
M2loc = gravity*recip_rhoConst*dSigmaH
c SNloc = SQRT( Sloc*M2loc )
N2loc = -gravity*recip_rhoConst*dSigmaR
c N2loc = -gravity*recip_rhoConst*dSigmaDr(i,j)
IF ( N2loc.GT.0. _d 0 ) THEN
SNloc = Sloc*SQRT(N2loc)
ELSE
SNloc = 0. _d 0
ENDIF
ELSE
SNloc = 0. _d 0
ENDIF
VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
& +deltaH*GM_Visbeck_alpha
& *GM_Visbeck_length*GM_Visbeck_length*SNloc
ENDIF
ENDDO
ENDDO
ENDIF
#endif /* ndef OLD_VISBECK_CALC */
#endif /* GM_VISBECK_VARIABLE_K */
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
dSigmaDr(i,j)=sigmaR(i,j,k)
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
kkey = (igmkey-1)*Nr + k
CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE baseSlope(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE hTransLay(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE recipLambda(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C Calculate slopes for use in tensor, taper and/or clip
CALL GMREDI_SLOPE_LIMIT(
O SlopeX, SlopeY,
O SlopeSqr, taperFct,
U hTransLay, baseSlope, recipLambda,
U dSigmaDr,
I dSigmaDx, dSigmaDy,
I ldd97_LrhoC, locMixLayer, rF,
I kLowC(1-OLx,1-OLy,bi,bj),
I k, bi, bj, myTime, myIter, myThid )
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
C Mask Iso-neutral slopes
SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)
SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)
SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE SlopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE SlopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C Components of Redi/GM tensor
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
ENDDO
ENDDO
#ifdef GM_VISBECK_VARIABLE_K
#ifdef OLD_VISBECK_CALC
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K
C but do not know if *taperFct (or **2 ?) is necessary
Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)
C-- Depth average of M^2/N^2 * N
C Calculate terms for mean Richardson number
C which is used in the "variable K" parameterisaton.
C Distance between interface above layer and the integration depth
deltaH=abs(GM_Visbeck_depth)-abs(rF(k))
C If positive we limit this to the layer thickness
integrDepth = drF(k)
deltaH=min(deltaH,integrDepth)
C If negative then we are below the integration level
deltaH=max(deltaH, 0. _d 0)
C Now we convert deltaH to a non-dimensional fraction
deltaH=deltaH/GM_Visbeck_depth
IF ( Ssq(i,j).NE.0. .AND. dSigmaDr(i,j).NE.0. ) THEN
N2loc = -gravity*recip_rhoConst*dSigmaDr(i,j)
SNloc = SQRT(Ssq(i,j)*N2loc )
VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
& +deltaH*GM_Visbeck_alpha
& *GM_Visbeck_length*GM_Visbeck_length*SNloc
ENDIF
ENDDO
ENDDO
#endif /* OLD_VISBECK_CALC */
#endif /* GM_VISBECK_VARIABLE_K */
C-- end 1rst loop on vertical level index k
ENDDO
#ifdef GM_VISBECK_VARIABLE_K
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
#endif
IF ( GM_Visbeck_alpha.GT.0. ) THEN
C- Limit range that KapGM can take
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
VisbeckK(i,j,bi,bj)=
& MIN( MAX( VisbeckK(i,j,bi,bj), GM_Visbeck_minVal_K ),
& GM_Visbeck_maxVal_K )
ENDDO
ENDDO
ENDIF
cph( NEW
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
#endif
cph)
#endif /* GM_VISBECK_VARIABLE_K */
C- express the Tensor in term of Diffusivity (= m**2 / s )
DO k=1,Nr
#ifdef ALLOW_AUTODIFF_TAMC
kkey = (igmkey-1)*Nr + k
# if (defined (GM_NON_UNITY_DIAGONAL)
defined (GM_VISBECK_VARIABLE_K))
CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
# endif
#endif
km1 = MAX(k-1,1)
isopycK = GM_isopycK
& *(GM_isoFac1d(km1)+GM_isoFac1d(k))*op5
bolus_K = GM_background_K
& *(GM_bolFac1d(km1)+GM_bolFac1d(k))*op5
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
#ifdef ALLOW_KAPREDI_CONTROL
# ifdef ALLOW_KAPREDI_CONTROL_OLD
Kgm_tmp = kapRedi(i,j,k,bi,bj)
else
Kgm_tmp = op5*(kapRedi(i,j,k,bi,bj)+kapRedi(i,j,km1,bi,bj))
# endif
#else
Kgm_tmp = isopycK*GM_isoFac2d(i,j,bi,bj)
#endif
#ifdef ALLOW_KAPGM_CONTROL
# ifdef ALLOW_KAPGM_CONTROL_OLD
& + GM_skewflx*kapGM(i,j,k,bi,bj)
else
& + GM_skewflx*op5*(kapGM(i,j,k,bi,bj)+kapGM(i,j,km1,bi,bj))
# endif
#else
& + GM_skewflx*bolus_K*GM_bolFac2d(i,j,bi,bj)
#endif
#ifdef GM_VISBECK_VARIABLE_K
& + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)
#endif
#if ((defined GM_K3D) ! (defined GM_K3D_PASSIVE))
& + op5*(K3D(i,j,k,bi,bj)+K3D(i,j,km1,bi,bj))
& *(1. _d 0 + GM_skewflx)
#endif
Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)
Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)
#ifdef ALLOW_KAPREDI_CONTROL
# ifdef ALLOW_KAPREDI_CONTROL_OLD
Kwz(i,j,k,bi,bj)= ( kapRedi(i,j,k,bi,bj)
else
Kwz(i,j,k,bi,bj)= ( op5*(kapRedi(i,j,k,bi,bj)
& +kapRedi(i,j,km1,bi,bj))
# endif
#else
Kwz(i,j,k,bi,bj)= ( isopycK*GM_isoFac2d(i,j,bi,bj)
#endif
#ifdef GM_VISBECK_VARIABLE_K
& + VisbeckK(i,j,bi,bj)
#endif
#if ((defined GM_K3D) ! (defined GM_K3D_PASSIVE))
& + op5*(K3D(i,j,k,bi,bj)+K3D(i,j,km1,bi,bj))
#endif
& )*Kwz(i,j,k,bi,bj)
ENDDO
ENDDO
ENDDO
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics .AND. GM_taper_scheme.EQ.'fm07' ) THEN
CALL DIAGNOSTICS_FILL( hTransLay, 'GM_hTrsL', 0,1,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL( baseSlope, 'GM_baseS', 0,1,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(recipLambda,'GM_rLamb', 0,1,2,bi,bj,myThid)
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- Calculate Stream-Functions used in Advective Form:
#ifdef GM_BOLUS_ADVEC
IF (GM_AdvForm) THEN
#ifdef GM_BOLUS_BVP
IF (GM_UseBVP) THEN
CALL GMREDI_CALC_PSI_BVP(
I bi, bj, iMin, iMax, jMin, jMax,
I sigmaX, sigmaY, sigmaR,
I myThid )
ELSE
#endif
#ifndef GM_K3D_PASSIVE
IF (.NOT. GM_useK3D) THEN
#endif
C If using GM_K3D PsiX and PsiY are calculated in gmredi_k3d
CALL GMREDI_CALC_PSI_B(
I bi, bj, iMin, iMax, jMin, jMax,
I sigmaX, sigmaY, sigmaR,
I ldd97_LrhoW, ldd97_LrhoS,
I myThid )
#ifndef GM_K3D_PASSIVE
ENDIF
#endif
#ifdef GM_BOLUS_BVP
ENDIF
#endif
ENDIF
#endif
#ifndef GM_EXCLUDE_SUBMESO
IF ( GM_useSubMeso .AND. GM_AdvForm ) THEN
CALL SUBMESO_CALC_PSI(
I bi, bj, iMin, iMax, jMin, jMax,
I sigmaX, sigmaY, sigmaR,
I locMixLayer,
I myIter, myThid )
ENDIF
#endif /* ndef GM_EXCLUDE_SUBMESO */
#if ( defined (GM_NON_UNITY_DIAGONAL) defined (GM_EXTRA_DIAGONAL) )
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- 2nd k loop : compute Tensor Coeff. at U point
#ifdef ALLOW_KPP
IF ( useKPP ) THEN
DO j=1-OLy,sNy+OLy
DO i=2-OLx,sNx+OLx
locMixLayer(i,j) = ( KPPhbl(i-1,j,bi,bj)
& + KPPhbl( i ,j,bi,bj) )*op5
ENDDO
ENDDO
ELSE
#else
IF ( .TRUE. ) THEN
#endif
DO j=1-OLy,sNy+OLy
DO i=2-OLx,sNx+OLx
locMixLayer(i,j) = ( hMixLayer(i-1,j,bi,bj)
& + hMixLayer( i ,j,bi,bj) )*op5
ENDDO
ENDDO
ENDIF
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
hTransLay(i,j) = 0.
baseSlope(i,j) = 0.
recipLambda(i,j)= 0.
ENDDO
DO i=2-OLx,sNx+OLx
hTransLay(i,j) = MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
ENDDO
ENDDO
DO k=Nr,1,-1
kp1 = MIN(Nr,k+1)
maskp1 = 1. _d 0
IF (k.GE.Nr) maskp1 = 0. _d 0
#ifdef ALLOW_AUTODIFF_TAMC
kkey = (igmkey-1)*Nr + k
#endif
C Gradient of Sigma at U points
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
dSigmaDx(i,j)=sigmaX(i,j,k)
& *_maskW(i,j,k,bi,bj)
dSigmaDy(i,j)=op25*( sigmaY(i-1,j+1,k)+sigmaY(i,j+1,k)
& +sigmaY(i-1, j ,k)+sigmaY(i, j ,k)
& )*_maskW(i,j,k,bi,bj)
dSigmaDr(i,j)=op25*( sigmaR(i-1,j, k )+sigmaR(i,j, k )
& +(sigmaR(i-1,j,kp1)+sigmaR(i,j,kp1))*maskp1
& )*_maskW(i,j,k,bi,bj)
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE locMixLayer(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE baseSlope(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE hTransLay(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE recipLambda(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C Calculate slopes for use in tensor, taper and/or clip
CALL GMREDI_SLOPE_LIMIT(
O SlopeX, SlopeY,
O SlopeSqr, taperFct,
U hTransLay, baseSlope, recipLambda,
U dSigmaDr,
I dSigmaDx, dSigmaDy,
I ldd97_LrhoW, locMixLayer, rC,
I kLow_W,
I k, bi, bj, myTime, myIter, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef GM_NON_UNITY_DIAGONAL
c IF ( GM_nonUnitDiag ) THEN
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
Kux(i,j,k,bi,bj) =
#ifdef ALLOW_KAPREDI_CONTROL
# ifdef ALLOW_KAPREDI_CONTROL_OLD
& ( kapRedi(i,j,k,bi,bj)
else
& ( op5*(kapRedi(i,j,k,bi,bj)+kapRedi(i-1,j,k,bi,bj))
# endif
#else
& ( GM_isopycK*GM_isoFac1d(k)
& *op5*(GM_isoFac2d(i-1,j,bi,bj)+GM_isoFac2d(i,j,bi,bj))
#endif
#ifdef GM_VISBECK_VARIABLE_K
& +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
#endif
#if ((defined GM_K3D) ! (defined GM_K3D_PASSIVE))
& +op5*(K3D(i,j,k,bi,bj)+K3D(i-1,j,k,bi,bj))
#endif
& )*taperFct(i,j)
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
# ifdef GM_EXCLUDE_CLIPPING
CADJ STORE Kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
# endif
#endif
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
Kux(i,j,k,bi,bj) = MAX( Kux(i,j,k,bi,bj), GM_Kmin_horiz )
ENDDO
ENDDO
c ENDIF
#endif /* GM_NON_UNITY_DIAGONAL */
#ifdef GM_EXTRA_DIAGONAL
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE SlopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
IF ( GM_ExtraDiag ) THEN
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
Kuz(i,j,k,bi,bj) =
#ifdef ALLOW_KAPREDI_CONTROL
# ifdef ALLOW_KAPREDI_CONTROL_OLD
& ( kapRedi(i,j,k,bi,bj)
else
& ( op5*(kapRedi(i,j,k,bi,bj)+kapRedi(i-1,j,k,bi,bj))
# endif
#else
& ( GM_isopycK*GM_isoFac1d(k)
& *op5*(GM_isoFac2d(i-1,j,bi,bj)+GM_isoFac2d(i,j,bi,bj))
#endif
#ifdef ALLOW_KAPGM_CONTROL
# ifdef ALLOW_KAPGM_CONTROL_OLD
& - GM_skewflx*kapGM(i,j,k,bi,bj)
else
& - GM_skewflx*op5*(kapGM(i,j,k,bi,bj)+kapGM(i-1,j,k,bi,bj))
# endif
#else
& - GM_skewflx*GM_background_K*GM_bolFac1d(k)
& *op5*(GM_bolFac2d(i-1,j,bi,bj)+GM_bolFac2d(i,j,bi,bj))
#endif
#ifdef GM_VISBECK_VARIABLE_K
& +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*GM_advect
#endif
#if ((defined GM_K3D) ! (defined GM_K3D_PASSIVE))
& +op5*(K3D(i,j,k,bi,bj)+K3D(i-1,j,k,bi,bj))*GM_advect
#endif
& )*SlopeX(i,j)*taperFct(i,j)
ENDDO
ENDDO
ENDIF
#endif /* GM_EXTRA_DIAGONAL */
#ifdef ALLOW_DIAGNOSTICS
IF (doDiagRediFlx) THEN
km1 = MAX(k-1,1)
DO j=1,sNy
DO i=1,sNx+1
C store in tmp1k Kuz_Redi
#ifdef ALLOW_KAPREDI_CONTROL
# ifdef ALLOW_KAPREDI_CONTROL_OLD
tmp1k(i,j) = ( kapRedi(i,j,k,bi,bj)
else
tmp1k(i,j) = ( op5*(kapRedi(i-1,j,k,bi,bj)
& +kapRedi(i,j,k,bi,bj))
# endif
#else
tmp1k(i,j) = ( GM_isopycK*GM_isoFac1d(k)
& *op5*(GM_isoFac2d(i-1,j,bi,bj)+GM_isoFac2d(i,j,bi,bj))
#endif
#ifdef GM_VISBECK_VARIABLE_K
& +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0
#endif
#if ((defined GM_K3D) ! (defined GM_K3D_PASSIVE))
& +op5*(K3D(i,j,k,bi,bj)+K3D(i-1,j,k,bi,bj))
#endif
& )*SlopeX(i,j)*taperFct(i,j)
ENDDO
ENDDO
DO j=1,sNy
DO i=1,sNx+1
C- Vertical gradients interpolated to U points
dTdz = (
& +recip_drC(k)*
& ( maskC(i-1,j,k,bi,bj)*
& (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
& +maskC( i ,j,k,bi,bj)*
& (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
& )
& +recip_drC(kp1)*
& ( maskC(i-1,j,kp1,bi,bj)*
& (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
& +maskC( i ,j,kp1,bi,bj)*
& (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
& ) ) * 0.25 _d 0
tmp1k(i,j) = dyG(i,j,bi,bj)*drF(k)
& * _hFacW(i,j,k,bi,bj)
& * tmp1k(i,j) * dTdz
ENDDO
ENDDO
CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KuzTz', k,1,2,bi,bj,myThid)
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
C-- end 2nd loop on vertical level index k
ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- 3rd k loop : compute Tensor Coeff. at V point
#ifdef ALLOW_KPP
IF ( useKPP ) THEN
DO j=2-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
locMixLayer(i,j) = ( KPPhbl(i,j-1,bi,bj)
& + KPPhbl(i, j ,bi,bj) )*op5
ENDDO
ENDDO
ELSE
#else
IF ( .TRUE. ) THEN
#endif
DO j=2-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
locMixLayer(i,j) = ( hMixLayer(i,j-1,bi,bj)
& + hMixLayer(i, j ,bi,bj) )*op5
ENDDO
ENDDO
ENDIF
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
hTransLay(i,j) = 0.
baseSlope(i,j) = 0.
recipLambda(i,j)= 0.
ENDDO
ENDDO
DO j=2-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
hTransLay(i,j) = MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
ENDDO
ENDDO
C Gradient of Sigma at V points
DO k=Nr,1,-1
kp1 = MIN(Nr,k+1)
maskp1 = 1. _d 0
IF (k.GE.Nr) maskp1 = 0. _d 0
#ifdef ALLOW_AUTODIFF_TAMC
kkey = (igmkey-1)*Nr + k
#endif
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
& +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k)
& )*_maskS(i,j,k,bi,bj)
dSigmaDy(i,j)=sigmaY(i,j,k)
& *_maskS(i,j,k,bi,bj)
dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k )+sigmaR(i,j, k )
& +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1
& )*_maskS(i,j,k,bi,bj)
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE baseSlope(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE hTransLay(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE recipLambda(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C Calculate slopes for use in tensor, taper and/or clip
CALL GMREDI_SLOPE_LIMIT(
O SlopeX, SlopeY,
O SlopeSqr, taperFct,
U hTransLay, baseSlope, recipLambda,
U dSigmaDr,
I dSigmaDx, dSigmaDy,
I ldd97_LrhoS, locMixLayer, rC,
I kLow_S,
I k, bi, bj, myTime, myIter, myThid )
cph(
#ifdef ALLOW_AUTODIFF_TAMC
cph(
CADJ STORE taperfct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
cph)
#endif /* ALLOW_AUTODIFF_TAMC */
cph)
#ifdef GM_NON_UNITY_DIAGONAL
c IF ( GM_nonUnitDiag ) THEN
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
Kvy(i,j,k,bi,bj) =
#ifdef ALLOW_KAPREDI_CONTROL
# ifdef ALLOW_KAPREDI_CONTROL_OLD
& ( kapRedi(i,j,k,bi,bj)
else
& ( op5*(kapRedi(i,j,k,bi,bj)+kapRedi(i,j-1,k,bi,bj))
# endif
#else
& ( GM_isopycK*GM_isoFac1d(k)
& *op5*(GM_isoFac2d(i,j-1,bi,bj)+GM_isoFac2d(i,j,bi,bj))
#endif
#ifdef GM_VISBECK_VARIABLE_K
& +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
#endif
#if ((defined GM_K3D) ! (defined GM_K3D_PASSIVE))
& +op5*(K3D(i,j,k,bi,bj)+K3D(i,j-1,k,bi,bj))
#endif
& )*taperFct(i,j)
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
# ifdef GM_EXCLUDE_CLIPPING
CADJ STORE Kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
# endif
#endif
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
Kvy(i,j,k,bi,bj) = MAX( Kvy(i,j,k,bi,bj), GM_Kmin_horiz )
ENDDO
ENDDO
c ENDIF
#endif /* GM_NON_UNITY_DIAGONAL */
#ifdef GM_EXTRA_DIAGONAL
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE SlopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
IF ( GM_ExtraDiag ) THEN
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
Kvz(i,j,k,bi,bj) =
#ifdef ALLOW_KAPREDI_CONTROL
# ifdef ALLOW_KAPREDI_CONTROL_OLD
& ( kapRedi(i,j,k,bi,bj)
else
& ( op5*(kapRedi(i,j,k,bi,bj)+kapRedi(i,j-1,k,bi,bj))
# endif
#else
& ( GM_isopycK*GM_isoFac1d(k)
& *op5*(GM_isoFac2d(i,j-1,bi,bj)+GM_isoFac2d(i,j,bi,bj))
#endif
#ifdef ALLOW_KAPGM_CONTROL
# ifdef ALLOW_KAPGM_CONTROL_OLD
& - GM_skewflx*kapGM(i,j,k,bi,bj)
else
& - GM_skewflx*op5*(kapGM(i,j,k,bi,bj)+kapGM(i,j-1,k,bi,bj))
# endif
#else
& - GM_skewflx*GM_background_K*GM_bolFac1d(k)
& *op5*(GM_bolFac2d(i,j-1,bi,bj)+GM_bolFac2d(i,j,bi,bj))
#endif
#ifdef GM_VISBECK_VARIABLE_K
& +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*GM_advect
#endif
#if ((defined GM_K3D) ! (defined GM_K3D_PASSIVE))
& +op5*(K3D(i,j,k,bi,bj)+K3D(i,j-1,k,bi,bj))*GM_advect
#endif
& )*SlopeY(i,j)*taperFct(i,j)
ENDDO
ENDDO
ENDIF
#endif /* GM_EXTRA_DIAGONAL */
#ifdef ALLOW_DIAGNOSTICS
IF (doDiagRediFlx) THEN
km1 = MAX(k-1,1)
DO j=1,sNy+1
DO i=1,sNx
C store in tmp1k Kvz_Redi
#ifdef ALLOW_KAPREDI_CONTROL
# ifdef ALLOW_KAPREDI_CONTROL_OLD
tmp1k(i,j) = ( kapRedi(i,j,k,bi,bj)
else
tmp1k(i,j) = ( op5*(kapRedi(i,j-1,k,bi,bj)
& +kapRedi(i,j,k,bi,bj))
# endif
#else
tmp1k(i,j) = ( GM_isopycK*GM_isoFac1d(k)
& *op5*(GM_isoFac2d(i,j-1,bi,bj)+GM_isoFac2d(i,j,bi,bj))
#endif
#ifdef GM_VISBECK_VARIABLE_K
& +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0
#endif
#if ((defined GM_K3D) ! (defined GM_K3D_PASSIVE))
& +op5*(K3D(i,j,k,bi,bj)+K3D(i,j-1,k,bi,bj))
#endif
& )*SlopeY(i,j)*taperFct(i,j)
ENDDO
ENDDO
DO j=1,sNy+1
DO i=1,sNx
C- Vertical gradients interpolated to U points
dTdz = (
& +recip_drC(k)*
& ( maskC(i,j-1,k,bi,bj)*
& (theta(i,j-1,km1,bi,bj)-theta(i,j-1,k,bi,bj))
& +maskC(i, j ,k,bi,bj)*
& (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))
& )
& +recip_drC(kp1)*
& ( maskC(i,j-1,kp1,bi,bj)*
& (theta(i,j-1,k,bi,bj)-theta(i,j-1,kp1,bi,bj))
& +maskC(i, j ,kp1,bi,bj)*
& (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))
& ) ) * 0.25 _d 0
tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)
& * _hFacS(i,j,k,bi,bj)
& * tmp1k(i,j) * dTdz
ENDDO
ENDDO
CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KvzTz', k,1,2,bi,bj,myThid)
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
C-- end 3rd loop on vertical level index k
ENDDO
#endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
#ifdef ALLOW_TIMEAVE
C-- Time-average
IF ( taveFreq.GT.0. ) THEN
CALL TIMEAVE_CUMULATE( GM_Kwx_T, Kwx, Nr,
& deltaTclock, bi, bj, myThid )
CALL TIMEAVE_CUMULATE( GM_Kwy_T, Kwy, Nr,
& deltaTclock, bi, bj, myThid )
CALL TIMEAVE_CUMULATE( GM_Kwz_T, Kwz, Nr,
& deltaTclock, bi, bj, myThid )
#ifdef GM_VISBECK_VARIABLE_K
IF ( GM_Visbeck_alpha.NE.0. ) THEN
CALL TIMEAVE_CUMULATE( Visbeck_K_T, VisbeckK, 1,
& deltaTclock, bi, bj, myThid )
ENDIF
#endif
#ifdef GM_BOLUS_ADVEC
IF ( GM_AdvForm ) THEN
CALL TIMEAVE_CUMULATE( GM_PsiXtave, GM_PsiX, Nr,
& deltaTclock, bi, bj, myThid )
CALL TIMEAVE_CUMULATE( GM_PsiYtave, GM_PsiY, Nr,
& deltaTclock, bi, bj, myThid )
ENDIF
#endif
GM_timeAve(bi,bj) = GM_timeAve(bi,bj)+deltaTclock
ENDIF
#endif /* ALLOW_TIMEAVE */
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
CALL GMREDI_DIAGNOSTICS_FILL(bi,bj,myThid)
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
#endif /* ALLOW_GMREDI */
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: GMREDI_CALC_TENSOR_DUMMY
C !INTERFACE:
SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
I iMin, iMax, jMin, jMax,
I sigmaX, sigmaY, sigmaR,
I bi, bj, myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE GMREDI_CALC_TENSOR_DUMMY
C | o Calculate tensor elements for GM/Redi tensor.
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "GMREDI.h"
C !INPUT/OUTPUT PARAMETERS:
_RL sigmaX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL sigmaY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
INTEGER iMin,iMax,jMin,jMax
INTEGER bi, bj
_RL myTime
INTEGER myIter
INTEGER myThid
CEOP
#ifdef ALLOW_GMREDI
C !LOCAL VARIABLES:
INTEGER i, j, k
DO k=1,Nr
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
Kwx(i,j,k,bi,bj) = 0.0
Kwy(i,j,k,bi,bj) = 0.0
Kwz(i,j,k,bi,bj) = 0.0
ENDDO
ENDDO
ENDDO
#endif /* ALLOW_GMREDI */
RETURN
END