C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/submesoscale/code/gmredi_calc_tensor.F,v 1.7 2010/03/19 19:29:47 zhc Exp $ C $Name: $ #include "GMREDI_OPTIONS.h" #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_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 Kgm_tmp _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 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) c#if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) ) INTEGER kp1 _RL maskp1 c#endif #ifdef GM_SUBMESO _RL dBdxAV(1-Olx:sNx+Olx,1-Oly:sNy+Oly) _RL dBdyAV(1-Olx:sNx+Olx,1-Oly:sNy+Oly) _RL SM_Lf(1-Olx:sNx+Olx,1-Oly:sNy+Oly) _RL SM_PsiX(1-Olx:sNx+Olx,1-Oly:sNy+Oly) _RL SM_PsiY(1-Olx:sNx+Olx,1-Oly:sNy+Oly) _RL SM_PsiXm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly) _RL SM_PsiYm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly) _RL hsqmu, hml, recip_hml, qfac, dS, mlmax #endif c#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 c#endif /* GM_VISBECK_VARIABLE_K */ #ifdef ALLOW_DIAGNOSTICS LOGICAL doDiagRediFlx LOGICAL DIAGNOSTICS_IS_ON EXTERNAL DIAGNOSTICS_IS_ON c#if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) ) INTEGER km1 _RL dTdz, dTdx, dTdy _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy) c#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 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- 1rst loop on k : compute Tensor Coeff. at W points. 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 C SM(1) mlmax=0. _d 0 #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) C SM(1) mlmax=max(mlmax,locMixLayer(i,j)) 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) C SM(1) mlmax=max(mlmax,locMixLayer(i,j)) ENDDO ENDDO ENDIF #ifdef GM_SUBMESO DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx dBdxAV(i,j) = 0. _d 0 dBdyAV(i,j) = 0. _d 0 SM_Lf(i,j) = 0. _d 0 SM_PsiX(i,j) = 0. _d 0 SM_PsiY(i,j) = 0. _d 0 SM_PsiXm1(i,j) = 0. _d 0 SM_PsiYm1(i,j) = 0. _d 0 ENDDO ENDDO #endif DO k=Nr,2,-1 #ifdef ALLOW_AUTODIFF_TAMC kkey = (igmkey-1)*Nr + k 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 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 # ifdef GM_BOLUS_ADVEC GM_PsiX(i,j,k,bi,bj) = 0. _d 0 GM_PsiY(i,j,k,bi,bj) = 0. _d 0 # endif ENDDO ENDDO #endif /* ALLOW_AUTODIFF_TAMC */ 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) dSigmaDr(i,j)=sigmaR(i,j,k) #ifdef GM_SUBMESO #ifdef GM_SUBMESO_VARYLf C-- Depth average of SigmaR at W points C compute depth average from surface down to the MixLayer depth IF (-rC(k-1).LT.locMixLayer(i,j) ) THEN IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN integrDepth = -rC( k ) C- in 2 steps to avoid mix of RS & RL type in min fct. arguments integrDepth = MIN( integrDepth, locMixLayer(i,j) ) 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-- Store db/dr in SM_Lf for now. SM_Lf(i,j) = SM_Lf(i,j) & -gravity*recip_rhoConst*dSigmaDr(i,j)*deltaH ENDIF ENDIF #endif #endif 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 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_SUBMESO CBFK-- Use the dsigmadr average to construct the coefficients of the SM param DO j=1-Oly+1,sNy+Oly-1 DO i=1-Olx+1,sNx+Olx-1 #ifdef GM_SUBMESO_VARYLf IF (SM_Lf(i,j).gt.0) THEN CBFK ML def. rad. as Lf if available and not too small SM_Lf(i,j)=max(sqrt(SM_Lf(i,j))*locMixLayer(i,j) & /abs(fCori(i,j,bi,bj)) & ,GM_SM_Lf) ELSE #else IF (.TRUE.) THEN #endif CBFK Otherwise, store just the fixed number SM_Lf(i,j)=GM_SM_Lf ENDIF CBFK Now do the rest of the coefficient dS=2*dxC(i,j,bi,bj)*dyC(i,j,bi,bj)/ & (dxC(i,j,bi,bj)+dyC(i,j,bi,bj)) CBFK Scaling only works up to 1 degree. dS=min(dS,GM_SM_Lmax) deltaH=sqrt(fCori(i,j,bi,bj)**2+1 _d 0/(GM_SM_tau**2)) SM_Lf(i,j)=GM_SM_Ce*dS/(deltaH*SM_Lf(i,j)) ENDDO ENDDO #endif #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 DO j=1-Oly+1,sNy+Oly-1 DO i=1-Olx+1,sNx+Olx-1 #ifdef ALLOW_KAPREDI_CONTROL Kgm_tmp = kapredi(i,j,k,bi,bj) #else Kgm_tmp = GM_isopycK #endif #ifdef ALLOW_KAPGM_CONTROL & + GM_skewflx*kapgm(i,j,k,bi,bj) #else & + GM_skewflx*GM_background_K #endif #ifdef GM_VISBECK_VARIABLE_K & + VisbeckK(i,j,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 Kwz(i,j,k,bi,bj)= ( kapredi(i,j,k,bi,bj) #else Kwz(i,j,k,bi,bj)= ( GM_isopycK #endif #ifdef GM_VISBECK_VARIABLE_K & + VisbeckK(i,j,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 */ #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) #ifdef GM_SUBMESO C-- Depth average of SigmaX at U points C compute depth average from surface down to the MixLayer depth IF (k.GT.1) THEN IF (-rC(k-1).LT.locMixLayer(i,j) ) THEN IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN integrDepth = -rC( k ) C- in 2 steps to avoid mix of RS & RL type in min fct. arguments integrDepth = MIN( integrDepth, locMixLayer(i,j) ) 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 (= M^2 / N since S=M^2/N^2 ) dBdxAV(i,j) = dBdxAV(i,j) & +dSigmaDx(i,j)*deltaH*recip_rhoConst*gravity ENDIF ENDIF ENDIF #endif 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 & ( kapredi(i,j,k,bi,bj) #else & ( GM_isopycK #endif #ifdef GM_VISBECK_VARIABLE_K & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,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 & ( kapredi(i,j,k,bi,bj) #else & ( GM_isopycK #endif #ifdef ALLOW_KAPGM_CONTROL & - GM_skewflx*kapgm(i,j,k,bi,bj) #else & - GM_skewflx*GM_background_K #endif #ifdef GM_VISBECK_VARIABLE_K & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,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 tmp1k(i,j) = ( kapredi(i,j,k,bi,bj) #else tmp1k(i,j) = ( GM_isopycK #endif #ifdef GM_VISBECK_VARIABLE_K & +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0 #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) #ifdef GM_SUBMESO C-- Depth average of SigmaY at V points C compute depth average from surface down to the MixLayer depth IF (k.GT.1) THEN IF (-rC(k-1).LT.locMixLayer(i,j) ) THEN IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN integrDepth = -rC( k ) C- in 2 steps to avoid mix of RS & RL type in min fct. arguments integrDepth = MIN( integrDepth, locMixLayer(i,j) ) 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) ) dBdyAV(i,j) = dBdyAV(i,j) & +dSigmaDy(i,j)*deltaH*recip_rhoConst*gravity ENDIF ENDIF ENDIF #endif 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 & ( kapredi(i,j,k,bi,bj) #else & ( GM_isopycK #endif #ifdef GM_VISBECK_VARIABLE_K & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,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 & ( kapredi(i,j,k,bi,bj) #else & ( GM_isopycK #endif #ifdef ALLOW_KAPGM_CONTROL & - GM_skewflx*kapgm(i,j,k,bi,bj) #else & - GM_skewflx*GM_background_K #endif #ifdef GM_VISBECK_VARIABLE_K & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,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 tmp1k(i,j) = ( kapredi(i,j,k,bi,bj) #else tmp1k(i,j) = ( GM_isopycK #endif #ifdef GM_VISBECK_VARIABLE_K & +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0 #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 GM_BOLUS_ADVEC IF (GM_AdvForm) THEN CALL GMREDI_CALC_PSI_B( I bi, bj, iMin, iMax, jMin, jMax, I sigmaX, sigmaY, sigmaR, I ldd97_LrhoW, ldd97_LrhoS, I myThid ) ENDIF #endif #ifdef GM_SUBMESO CBFK Add the submesoscale contribution, in a 4th k loop DO k=1,Nr km1=max(1,k-1) IF ((k.gt.1).and.(-rF(k-1) .lt. mlmax)) THEN kp1 = MIN(k+1,Nr) CBFK Add in the mu vertical structure factor DO j=1-Oly+1,sNy+Oly-1 DO i=1-Olx+1,sNx+Olx-1 hml=hMixLayer(i,j,bi,bj) IF (hml.gt.0 _d 0) THEN recip_hml=1 _d 0/hml ELSE recip_hml=0 _d 0 ENDIF CBFK We calculate the h^2 mu(z) factor only on w points. CBFK It is possible that we might need to calculate it CBFK on Psi or u,v points independently to prevent spurious CBFK entrainment. Unlikely that this will be major CBFK (it wasnt in offline testing). qfac=(2*rf(k)*recip_hml+1 _d 0)**2 hsqmu=(1 _d 0-qfac)*(1 _d 0+(5 _d 0)*qfac/21 _d 0) hsqmu=max(0 _d 0, hsqmu)*hml**2 SM_Lf(i,j)=SM_Lf(i,j)*hsqmu ENDDO ENDDO CBFK Now interpolate to match locations DO j=1-Oly+1,sNy+Oly-1 DO i=1-Olx+1,sNx+Olx-1 C SM_Lf coefficients are on rVel points C Psix are on faces above U SM_PsiX(i,j)=op5*(SM_Lf(i+1,j)+SM_Lf(i,j))*dBdxAV(i,j) & *_maskW(i,j,k,bi,bj) C Psiy are on faces above V SM_PsiY(i,j)=op5*(SM_Lf(i,j+1)+SM_Lf(i,j))*dBdyAV(i,j) & *_maskS(i,j,k,bi,bj) c hzhang: clipping here, ref to Baylor paper 3 appendix A MOM hml=hMixLayer(i,j,bi,bj) IF (hml .lt. (delR(1)+delR(2)+delR(3)+delR(4))) THEN SM_PsiX(i,j)=0 _d 0 SM_PsiY(i,j)=0 _d 0 ENDIF hml=.5 * delR(k) IF ( abs(SM_PsiX(i,j)) .gt. hml ) THEN SM_PsiX(i,j)=SIGN( hml, SM_PsiX(i,j) ) ENDIF IF ( abs(SM_PsiY(i,j)) .gt. hml ) THEN SM_PsiY(i,j)=SIGN( hml, SM_PsiY(i,j) ) ENDIF c hzhang: clipping done #ifndef GM_BOLUS_ADVEC C Kwx,Kwy are on rVel Points Kwx(i,j,k,bi,bj) = Kwx(i,j,k,bi,bj) & +op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j)) Kwy(i,j,k,bi,bj) = Kwy(i,j,k,bi,bj) & +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j)) #ifdef GM_EXTRA_DIAGONAL IF (GM_ExtraDiag) THEN C Kuz,Kvz are on u,v Points Kuz(i,j,k,bi,bj) = Kuz(i,j,k,bi,bj) & -op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j)) Kvz(i,j,k,bi,bj) = Kvz(i,j,k,bi,bj) & -op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j)) ENDIF #endif #else IF (GM_AdvForm) THEN GM_PsiX(i,j,k,bi,bj)=GM_PsiX(i,j,k,bi,bj)+SM_PsiX(i,j) GM_PsiY(i,j,k,bi,bj)=GM_PsiY(i,j,k,bi,bj)+SM_PsiY(i,j) ENDIF #endif ENDDO ENDDO ELSE DO j=1-Oly+1,sNy+Oly-1 DO i=1-Olx+1,sNx+Olx-1 SM_PsiX(i,j)=0. _d 0 SM_PsiY(i,j)=0. _d 0 ENDDO ENDDO ENDIF #ifdef ALLOW_DIAGNOSTICS IF ( useDiagnostics ) THEN IF ( DIAGNOSTICS_IS_ON('SM_PsiX ',myThid) ) THEN CALL DIAGNOSTICS_FILL(SM_PsiX,'SM_PsiX ',k,1,2,bi,bj,myThid) ENDIF IF ( DIAGNOSTICS_IS_ON('SM_PsiY ',myThid) ) THEN CALL DIAGNOSTICS_FILL(SM_PsiY,'SM_PsiY ',k,1,2,bi,bj,myThid) ENDIF CBFK Note: for comparision, you can diagnose the bolus form CBFK or the Kappa form in the same simulation, regardless of other CBFK settings IF ( DIAGNOSTICS_IS_ON('SM_ubT ',myThid) ) THEN DO j=jMin,jMax DO i=iMin,iMax tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiX(i,j) & -SM_PsiXm1(i,j) ) & *maskW(i,j,km1,bi,bj) & *op5*(Theta(i,j,km1,bi,bj)+Theta(i-1,j,km1,bi,bj)) ENDDO ENDDO CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT ', km1,1,2,bi,bj,myThid) ENDIF IF ( DIAGNOSTICS_IS_ON('SM_vbT ',myThid) ) THEN DO j=jMin,jMax DO i=iMin,iMax tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiY(i,j) & -SM_PsiYm1(i,j) ) & *maskS(i,j,km1,bi,bj) & *op5*(Theta(i,j,km1,bi,bj)+Theta(i,j-1,km1,bi,bj)) ENDDO ENDDO CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT ', km1,1,2,bi,bj,myThid) ENDIF IF ( DIAGNOSTICS_IS_ON('SM_wbT ',myThid) ) THEN DO j=jMin,jMax DO i=iMin,iMax tmp1k(i,j) = & (dyG(i+1,j,bi,bj)*SM_PsiX(i+1,j) & -dyG( i ,j,bi,bj)*SM_PsiX( i ,j) & +dxG(i,j+1,bi,bj)*SM_PsiY(i,j+1) & -dxG(i, j ,bi,bj)*SM_PsiY(i, j )) & *op5*(Theta(i,j,k,bi,bj)+Theta(i,j,km1,bi,bj)) ENDDO ENDDO CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT ', k,1,2,bi,bj,myThid) C print *,'SM_wbT',k,tmp1k ENDIF IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN 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) & *op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j)) & * dTdz ENDDO ENDDO CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KuzTz', k,1,2,bi,bj,myThid) C print *,'SM_KuzTz',k,tmp1k ENDIF IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN DO j=1,sNy+1 DO i=1,sNx C- Vertical gradients interpolated to V points dTdz = op5*( & +op5*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)) & ) & +op5*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)) & ) ) tmp1k(i,j) = - dxG(i,j,bi,bj)*drF(k) & * _hFacS(i,j,k,bi,bj) & *op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j)) & * dTdz ENDDO ENDDO CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KvzTz', k,1,2,bi,bj,myThid) C print *,'SM_KvzTz',k,tmp1k ENDIF IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN DO j=jMin,jMax DO i=iMin,iMax C- Horizontal gradients interpolated to W points dTdx = op5*( & +op5*(_maskW(i+1,j,k,bi,bj) & *_recip_dxC(i+1,j,bi,bj)* & (Theta(i+1,j,k,bi,bj)-Theta(i,j,k,bi,bj)) & +_maskW(i,j,k,bi,bj) & *_recip_dxC(i,j,bi,bj)* & (Theta(i,j,k,bi,bj)-Theta(i-1,j,k,bi,bj))) & +op5*(_maskW(i+1,j,k-1,bi,bj) & *_recip_dxC(i+1,j,bi,bj)* & (Theta(i+1,j,k-1,bi,bj)-Theta(i,j,k-1,bi,bj)) & +_maskW(i,j,k-1,bi,bj) & *_recip_dxC(i,j,bi,bj)* & (Theta(i,j,k-1,bi,bj)-Theta(i-1,j,k-1,bi,bj))) & ) dTdy = op5*( & +op5*(_maskS(i,j,k,bi,bj) & *_recip_dyC(i,j,bi,bj)* & (Theta(i,j,k,bi,bj)-Theta(i,j-1,k,bi,bj)) & +_maskS(i,j+1,k,bi,bj) & *_recip_dyC(i,j+1,bi,bj)* & (Theta(i,j+1,k,bi,bj)-Theta(i,j,k,bi,bj))) & +op5*(_maskS(i,j,k-1,bi,bj) & *_recip_dyC(i,j,bi,bj)* & (Theta(i,j,k-1,bi,bj)-Theta(i,j-1,k-1,bi,bj)) & +_maskS(i,j+1,k-1,bi,bj) & *_recip_dyC(i,j+1,bi,bj)* & (Theta(i,j+1,k-1,bi,bj)-Theta(i,j,k-1,bi,bj))) & ) tmp1k(i,j) = - _rA(i,j,bi,bj) & *(op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))*dTdx & +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))*dTdy) ENDDO ENDDO CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', k,1,2,bi,bj,myThid) C print *,'SM_KrddT',k,tmp1k ENDIF ENDIF #endif DO j=1-Oly+1,sNy+Oly-1 DO i=1-Olx+1,sNx+Olx-1 SM_PsiXm1(i,j)=SM_PsiX(i,j) SM_PsiYm1(i,j)=SM_PsiY(i,j) tmp1k(i,j)=0 _d 0 ENDDO ENDDO ENDDO CBFK Always Zero at the bottom. IF ( DIAGNOSTICS_IS_ON('SM_ubT ',myThid) ) THEN CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT ', Nr,1,2,bi,bj,myThid) ENDIF IF ( DIAGNOSTICS_IS_ON('SM_vbT ',myThid) ) THEN CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT ', Nr,1,2,bi,bj,myThid) ENDIF IF ( DIAGNOSTICS_IS_ON('SM_wbT ',myThid) ) THEN CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT ', Nr,1,2,bi,bj,myThid) ENDIF IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN CALL DIAGNOSTICS_FILL(tmp1k,'SM_KuzTz', Nr,1,2,bi,bj,myThid) ENDIF IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN CALL DIAGNOSTICS_FILL(tmp1k,'SM_KvzTz', Nr,1,2,bi,bj,myThid) ENDIF IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', Nr,1,2,bi,bj,myThid) ENDIF #endif #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-|--+----| SUBROUTINE GMREDI_CALC_TENSOR_DUMMY( I iMin, iMax, jMin, jMax, I sigmaX, sigmaY, sigmaR, I bi, bj, myTime, myIter, myThid ) C /==========================================================\ C | SUBROUTINE GMREDI_CALC_TENSOR | C | o Calculate tensor elements for GM/Redi tensor. | C |==========================================================| C \==========================================================/ IMPLICIT NONE C == Global variables == #include "SIZE.h" #include "EEPARAMS.h" #include "GMREDI.h" C == Routine arguments == C _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 CEndOfInterface #ifdef ALLOW_GMREDI 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