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