/[MITgcm]/MITgcm_contrib/submesoscale/code/gmredi_calc_tensor.F
ViewVC logotype

Diff of /MITgcm_contrib/submesoscale/code/gmredi_calc_tensor.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by dimitri, Fri May 30 21:51:23 2008 UTC revision 1.2 by dimitri, Fri May 30 22:13:42 2008 UTC
# Line 73  C     == Local variables == Line 73  C     == Local variables ==
73        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
74        _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
75        _RL maskp1, Kgm_tmp        _RL maskp1, Kgm_tmp
76          _RL deltaH, integrDepth
77        _RL ldd97_LrhoC(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL ldd97_LrhoC(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
78        _RL ldd97_LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL ldd97_LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
79        _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
# Line 85  C     == Local variables == Line 86  C     == Local variables ==
86        _RL hTransLay  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL hTransLay  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
87        _RL recipLambda(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL recipLambda(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
88    
89    #ifdef GM_SUBMESO
90          _RL dBdxAV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
91          _RL dBdyAV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
92          _RL SM_Lf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
93          _RL SM_PsiX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
94          _RL SM_PsiY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
95          _RL SM_PsiXm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
96          _RL SM_PsiYm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
97          _RL hsqmu, hml, recip_hml, qfac, dS, mlmax
98    #endif
99    
100  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
101  #ifdef OLD_VISBECK_CALC  #ifdef OLD_VISBECK_CALC
102        _RL deltaH,zero_rs        _RL zero_rs
103        PARAMETER(zero_rs=0.D0)        PARAMETER(zero_rs=0.D0)
104        _RL N2,SN        _RL N2,SN
105        _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
106  #else  #else
107        _RL dSigmaH        _RL dSigmaH
       _RL deltaH, integrDepth  
108        _RL Sloc, M2loc, SNloc        _RL Sloc, M2loc, SNloc
109  #endif  #endif
110  #endif  #endif
# Line 103  C     == Local variables == Line 114  C     == Local variables ==
114        LOGICAL  DIAGNOSTICS_IS_ON        LOGICAL  DIAGNOSTICS_IS_ON
115        EXTERNAL DIAGNOSTICS_IS_ON        EXTERNAL DIAGNOSTICS_IS_ON
116        INTEGER  km1        INTEGER  km1
117        _RL dTdz        _RL dTdz, dTdx, dTdy
118        _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119  #endif  #endif
120    
# Line 210  C-- 1rst loop on k : compute Tensor Coef Line 221  C-- 1rst loop on k : compute Tensor Coef
221           locMixLayer(i,j) = 0. _d 0           locMixLayer(i,j) = 0. _d 0
222         ENDDO         ENDDO
223        ENDDO        ENDDO
224          mlmax=0. _d 0
225  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
226        IF ( useKPP ) THEN        IF ( useKPP ) THEN
227         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
228          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
229           locMixLayer(i,j) = KPPhbl(i,j,bi,bj)           locMixLayer(i,j) = KPPhbl(i,j,bi,bj)
230             mlmax=max(mlmax,locMixLayer(i,j))
231          ENDDO          ENDDO
232         ENDDO         ENDDO
233        ELSE        ELSE
# Line 224  C-- 1rst loop on k : compute Tensor Coef Line 237  C-- 1rst loop on k : compute Tensor Coef
237         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
238          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
239           locMixLayer(i,j) = hMixLayer(i,j,bi,bj)           locMixLayer(i,j) = hMixLayer(i,j,bi,bj)
240             mlmax=max(mlmax,locMixLayer(i,j))
241          ENDDO          ENDDO
242         ENDDO         ENDDO
243        ENDIF        ENDIF
244    
245    #ifdef GM_SUBMESO
246          DO j=1-Oly,sNy+Oly
247           DO i=1-Olx,sNx+Olx
248            dBdxAV(i,j)     = 0. _d 0
249            dBdyAV(i,j)     = 0. _d 0
250            SM_Lf(i,j)        = 0. _d 0
251            SM_PsiX(i,j)      = 0. _d 0
252            SM_PsiY(i,j)      = 0. _d 0
253            SM_PsiXm1(i,j)    = 0. _d 0
254            SM_PsiXm1(i,j)    = 0. _d 0
255           ENDDO
256          ENDDO
257    #endif
258    
259        DO k=Nr,2,-1        DO k=Nr,2,-1
260    
261  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 270  C      Gradient of Sigma at rVel points Line 298  C      Gradient of Sigma at rVel points
298       &                       +sigmaY(i,j+1, k )+sigmaY(i,j, k )       &                       +sigmaY(i,j+1, k )+sigmaY(i,j, k )
299       &                      )*maskC(i,j,k,bi,bj)       &                      )*maskC(i,j,k,bi,bj)
300           dSigmaDr(i,j)=sigmaR(i,j,k)           dSigmaDr(i,j)=sigmaR(i,j,k)
301    
302    #ifdef GM_SUBMESO
303    #ifdef GM_SUBMESO_VARYLf
304    C--     Depth average of SigmaR at W points
305    C       compute depth average from surface down to the MixLayer depth
306              IF (-rC(k-1).LT.locMixLayer(i,j) ) THEN
307               IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
308                integrDepth = -rC( k )
309    C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments
310                integrDepth = MIN( integrDepth, locMixLayer(i,j) )
311    C       Distance between level center above and the integration depth
312                deltaH = integrDepth + rC(k-1)
313    C       If negative then we are below the integration level
314    C       (cannot be the case with 2 conditions on maskC & -rC(k-1))
315    C       If positive we limit this to the distance from center above
316                deltaH = MIN( deltaH, drC(k) )
317    C       Now we convert deltaH to a non-dimensional fraction
318                deltaH = deltaH/( integrDepth+rC(1) )
319    C--     Store db/dr in SM_Lf for now.
320                SM_Lf(i,j) = SM_Lf(i,j)
321         &                   -gravity*recip_rhoConst*dSigmaDr(i,j)*deltaH
322               ENDIF
323              ENDIF
324    #endif
325    #endif
326          ENDDO          ENDDO
327         ENDDO         ENDDO
328    
# Line 411  C       Now we convert deltaH to a non-d Line 464  C       Now we convert deltaH to a non-d
464  C-- end 1rst loop on vertical level index k  C-- end 1rst loop on vertical level index k
465        ENDDO        ENDDO
466    
467    #ifdef GM_SUBMESO
468    CBFK-- Use the dsigmadr average to construct the coefficients of the SM param
469          DO j=1-Oly+1,sNy+Oly-1
470           DO i=1-Olx+1,sNx+Olx-1
471    #ifdef GM_SUBMESO_VARYLf
472    
473            IF (SM_Lf(i,j).gt.0) THEN
474    CBFK ML def. rad. as Lf if available and not too small
475              SM_Lf(i,j)=max(sqrt(SM_Lf(i,j))*locMixLayer(i,j)
476         &                        /abs(fCori(i,j,bi,bj))
477         &                   ,GM_SM_Lf)
478            ELSE
479    #else
480            IF (.TRUE.) THEN
481    #endif
482    CBFK Otherwise, store just the fixed number
483              SM_Lf(i,j)=GM_SM_Lf
484            ENDIF
485    CBFK Now do the rest of the coefficient
486            dS=2*dxC(i,j,bi,bj)*dyC(i,j,bi,bj)/
487         &              (dxC(i,j,bi,bj)+dyC(i,j,bi,bj))
488    CBFK Scaling only works up to 1 degree.
489            dS=min(dS,GM_SM_Lmax)
490            deltaH=sqrt(fCori(i,j,bi,bj)**2+1 _d 0/(GM_SM_tau**2))
491            SM_Lf(i,j)=GM_SM_Ce*dS/(deltaH*SM_Lf(i,j))
492           ENDDO
493          ENDDO
494    #endif
495    
496  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
497  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 481  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bi Line 562  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bi
562        ENDIF        ENDIF
563  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
564    
   
565  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
566  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
567  C-- 2nd  k loop : compute Tensor Coeff. at U point  C-- 2nd  k loop : compute Tensor Coeff. at U point
# Line 535  C     Gradient of Sigma at U points Line 615  C     Gradient of Sigma at U points
615           dSigmaDr(i,j)=op25*( sigmaR(i-1,j, k )+sigmaR(i,j, k )           dSigmaDr(i,j)=op25*( sigmaR(i-1,j, k )+sigmaR(i,j, k )
616       &                      +(sigmaR(i-1,j,kp1)+sigmaR(i,j,kp1))*maskp1       &                      +(sigmaR(i-1,j,kp1)+sigmaR(i,j,kp1))*maskp1
617       &                      )*_maskW(i,j,k,bi,bj)       &                      )*_maskW(i,j,k,bi,bj)
618    
619    #ifdef GM_SUBMESO
620    C--     Depth average of SigmaX at U points
621    C       compute depth average from surface down to the MixLayer depth
622             IF (k.GT.1) THEN
623              IF (-rC(k-1).LT.locMixLayer(i,j) ) THEN
624               IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
625                integrDepth = -rC( k )
626    C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments
627                integrDepth = MIN( integrDepth, locMixLayer(i,j) )
628    C       Distance between level center above and the integration depth
629                deltaH = integrDepth + rC(k-1)
630    C       If negative then we are below the integration level
631    C       (cannot be the case with 2 conditions on maskC & -rC(k-1))
632    C       If positive we limit this to the distance from center above
633                deltaH = MIN( deltaH, drC(k) )
634    C       Now we convert deltaH to a non-dimensional fraction
635                deltaH = deltaH/( integrDepth+rC(1) )
636    C--      compute: ( M^2 * S )^1/2   (= M^2 / N since S=M^2/N^2 )
637                dBdxAV(i,j) = dBdxAV(i,j)
638         &            +dSigmaDx(i,j)*deltaH*recip_rhoConst*gravity
639               ENDIF
640              ENDIF
641             ENDIF
642    #endif
643    
644          ENDDO          ENDDO
645         ENDDO         ENDDO
646    
# Line 670  C-- end 2nd  loop on vertical level inde Line 776  C-- end 2nd  loop on vertical level inde
776    
777  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
778  C-- 3rd  k loop : compute Tensor Coeff. at V point  C-- 3rd  k loop : compute Tensor Coeff. at V point
   
779  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
780        IF ( useKPP ) THEN        IF ( useKPP ) THEN
781         DO j=2-Oly,sNy+Oly         DO j=2-Oly,sNy+Oly
# Line 722  C     Gradient of Sigma at V points Line 827  C     Gradient of Sigma at V points
827           dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k )+sigmaR(i,j, k )           dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k )+sigmaR(i,j, k )
828       &                      +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1       &                      +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1
829       &                      )*_maskS(i,j,k,bi,bj)       &                      )*_maskS(i,j,k,bi,bj)
830    
831    #ifdef GM_SUBMESO
832    C--     Depth average of SigmaY at V points
833    C       compute depth average from surface down to the MixLayer depth
834             IF (k.GT.1) THEN
835              IF (-rC(k-1).LT.locMixLayer(i,j) ) THEN
836               IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
837                integrDepth = -rC( k )
838    C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments
839                integrDepth = MIN( integrDepth, locMixLayer(i,j) )
840    C       Distance between level center above and the integration depth
841                deltaH = integrDepth + rC(k-1)
842    C       If negative then we are below the integration level
843    C       (cannot be the case with 2 conditions on maskC & -rC(k-1))
844    C       If positive we limit this to the distance from center above
845                deltaH = MIN( deltaH, drC(k) )
846    C       Now we convert deltaH to a non-dimensional fraction
847                deltaH = deltaH/( integrDepth+rC(1) )
848                dBdyAV(i,j) = dBdyAV(i,j)
849         &            +dSigmaDy(i,j)*deltaH*recip_rhoConst*gravity
850               ENDIF
851              ENDIF
852             ENDIF
853    #endif
854    
855          ENDDO          ENDDO
856         ENDDO         ENDDO
857    
# Line 830  C         store in tmp1k Kvz_Redi Line 960  C         store in tmp1k Kvz_Redi
960          ENDDO          ENDDO
961          DO j=1,sNy+1          DO j=1,sNy+1
962           DO i=1,sNx           DO i=1,sNx
963  C-        Vertical gradients interpolated to U points  C-        Vertical gradients interpolated to V points
964            dTdz = (            dTdz =   op5*(
965       &     +recip_drC(k)*       &   +op5*recip_drC(k)*
966       &       ( maskC(i,j-1,k,bi,bj)*       &       ( maskC(i,j-1,k,bi,bj)*
967       &           (theta(i,j-1,km1,bi,bj)-theta(i,j-1,k,bi,bj))       &           (theta(i,j-1,km1,bi,bj)-theta(i,j-1,k,bi,bj))
968       &        +maskC(i, j ,k,bi,bj)*       &        +maskC(i, j ,k,bi,bj)*
969       &           (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))       &           (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))
970       &       )       &       )
971       &     +recip_drC(kp1)*       &   +op5*recip_drC(kp1)*
972       &       ( maskC(i,j-1,kp1,bi,bj)*       &       ( maskC(i,j-1,kp1,bi,bj)*
973       &           (theta(i,j-1,k,bi,bj)-theta(i,j-1,kp1,bi,bj))       &           (theta(i,j-1,k,bi,bj)-theta(i,j-1,kp1,bi,bj))
974       &        +maskC(i, j ,kp1,bi,bj)*       &        +maskC(i, j ,kp1,bi,bj)*
975       &           (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))       &           (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))
976       &       )      ) * 0.25 _d 0       &       )      )
977             tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)             tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)
978       &                * _hFacS(i,j,k,bi,bj)       &                * _hFacS(i,j,k,bi,bj)
979       &                * tmp1k(i,j) * dTdz       &                * tmp1k(i,j) * dTdz
# Line 869  C-- end 3rd  loop on vertical level inde Line 999  C-- end 3rd  loop on vertical level inde
999        ENDIF        ENDIF
1000  #endif  #endif
1001    
1002    
1003    #ifdef GM_SUBMESO
1004    CBFK Add the submesoscale contribution, in a 4th k loop
1005          DO k=1,Nr
1006           km1=max(1,k-1)
1007           IF ((k.gt.1).and.(-rF(k-1) .lt. mlmax)) THEN
1008            kp1 = MIN(k+1,Nr)
1009    CBFK Add in the mu vertical structure factor
1010            DO j=1-Oly+1,sNy+Oly-1
1011             DO i=1-Olx+1,sNx+Olx-1
1012    #ifdef ALLOW_KPP
1013              hml=KPPhbl(i,j,bi,bj)
1014    #else
1015              hml=hMixLayer(i,j,bi,bj)
1016    #endif
1017              IF (hml.gt.0 _d 0) THEN
1018               recip_hml=1 _d 0/hml
1019              ELSE
1020               recip_hml=0 _d 0
1021              ENDIF
1022    CBFK We calculate the h^2 mu(z) factor only on w points.  
1023    CBFK It is possible that we might need to calculate it
1024    CBFK on Psi or u,v points independently to prevent spurious
1025    CBFK entrainment.  Unlikely that this will be major
1026    CBFK (it wasnt in offline testing).
1027              qfac=(2*rf(k)*recip_hml+1 _d 0)**2
1028              hsqmu=(1 _d 0-qfac)*(1 _d 0+(5 _d 0)*qfac/21 _d 0)
1029              hsqmu=max(0 _d 0, hsqmu)*hml**2
1030              SM_Lf(i,j)=SM_Lf(i,j)*hsqmu
1031             ENDDO
1032            ENDDO
1033    CBFK Now interpolate to match locations
1034            DO j=1-Oly+1,sNy+Oly-1
1035             DO i=1-Olx+1,sNx+Olx-1
1036    C SM_Lf coefficients are on rVel points
1037    C Psix are on faces above U
1038              SM_PsiX(i,j)=op5*(SM_Lf(i+1,j)+SM_Lf(i,j))*dBdxAV(i,j)
1039         &                        *_maskW(i,j,k,bi,bj)
1040    C Psiy are on faces above V
1041              SM_PsiY(i,j)=op5*(SM_Lf(i,j+1)+SM_Lf(i,j))*dBdyAV(i,j)
1042         &                        *_maskS(i,j,k,bi,bj)
1043    
1044    #ifndef GM_BOLUS_ADVEC
1045    C Kwx,Kwy are on rVel Points
1046              Kwx(i,j,k,bi,bj)  = Kwx(i,j,k,bi,bj)
1047         &                      +op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))
1048              Kwy(i,j,k,bi,bj)  = Kwy(i,j,k,bi,bj)
1049         &                      +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))
1050    #ifdef GM_EXTRA_DIAGONAL
1051              IF (GM_ExtraDiag) THEN
1052    C Kuz,Kvz are on u,v Points
1053               Kuz(i,j,k,bi,bj)  = Kuz(i,j,k,bi,bj)
1054         &                    -op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1055               Kvz(i,j,k,bi,bj)  = Kvz(i,j,k,bi,bj)
1056         &                    -op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1057              ENDIF
1058    #endif
1059    #else
1060              IF (GM_AdvForm) THEN
1061               GM_PsiX(i,j,k,bi,bj)=GM_PsiX(i,j,k,bi,bj)+SM_PsiX(i,j)
1062               GM_PsiY(i,j,k,bi,bj)=GM_PsiY(i,j,k,bi,bj)+SM_PsiY(i,j)
1063              ENDIF
1064    #endif
1065             ENDDO
1066            ENDDO
1067           ELSE
1068            DO j=1-Oly+1,sNy+Oly-1
1069             DO i=1-Olx+1,sNx+Olx-1
1070              SM_PsiX(i,j)=0. _d 0
1071              SM_PsiY(i,j)=0. _d 0
1072             ENDDO
1073            ENDDO
1074           ENDIF  
1075    
1076    #ifdef ALLOW_DIAGNOSTICS
1077           IF ( useDiagnostics ) THEN
1078            IF ( DIAGNOSTICS_IS_ON('SM_PsiX ',myThid) ) THEN
1079             CALL DIAGNOSTICS_FILL(SM_PsiX,'SM_PsiX ',k,1,2,bi,bj,myThid)
1080            ENDIF
1081            IF ( DIAGNOSTICS_IS_ON('SM_PsiY ',myThid) ) THEN
1082             CALL DIAGNOSTICS_FILL(SM_PsiY,'SM_PsiY ',k,1,2,bi,bj,myThid)
1083            ENDIF
1084    
1085    CBFK Note:  for comparision, you can diagnose the bolus form
1086    CBFK or the Kappa form in the same simulation, regardless of other
1087    CBFK settings
1088            IF ( DIAGNOSTICS_IS_ON('SM_ubT  ',myThid) ) THEN
1089             DO j=jMin,jMax
1090              DO i=iMin,iMax
1091               tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiX(i,j)
1092         &                                -SM_PsiXm1(i,j) )
1093         &                               *maskW(i,j,km1,bi,bj)
1094         &            *op5*(Theta(i,j,km1,bi,bj)+Theta(i-1,j,km1,bi,bj))
1095              ENDDO
1096             ENDDO
1097             CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT  ', km1,1,2,bi,bj,myThid)
1098            ENDIF
1099    
1100            IF ( DIAGNOSTICS_IS_ON('SM_vbT  ',myThid) ) THEN
1101             DO j=jMin,jMax
1102              DO i=iMin,iMax
1103               tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiY(i,j)
1104         &                                -SM_PsiYm1(i,j) )
1105         &                               *maskS(i,j,km1,bi,bj)
1106         &            *op5*(Theta(i,j,km1,bi,bj)+Theta(i,j-1,km1,bi,bj))
1107              ENDDO
1108             ENDDO
1109             CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT  ', km1,1,2,bi,bj,myThid)
1110            ENDIF
1111    
1112            IF ( DIAGNOSTICS_IS_ON('SM_wbT  ',myThid) ) THEN
1113             DO j=jMin,jMax
1114              DO i=iMin,iMax
1115               tmp1k(i,j) =
1116         &        (dyG(i+1,j,bi,bj)*SM_PsiX(i+1,j)
1117         &        -dyG( i ,j,bi,bj)*SM_PsiX( i ,j)
1118         &        +dxG(i,j+1,bi,bj)*SM_PsiY(i,j+1)
1119         &        -dxG(i, j ,bi,bj)*SM_PsiY(i, j ))
1120         &             *op5*(Theta(i,j,k,bi,bj)+Theta(i,j,km1,bi,bj))
1121              ENDDO
1122             ENDDO
1123             CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT  ', k,1,2,bi,bj,myThid)
1124    C         print *,'SM_wbT',k,tmp1k
1125            ENDIF
1126    
1127            IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1128             DO j=1,sNy
1129              DO i=1,sNx+1
1130    C-        Vertical gradients interpolated to U points
1131               dTdz = (
1132         &      +recip_drC(k)*
1133         &       ( maskC(i-1,j,k,bi,bj)*
1134         &           (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
1135         &        +maskC( i ,j,k,bi,bj)*
1136         &           (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
1137         &       )
1138         &      +recip_drC(kp1)*
1139         &       ( maskC(i-1,j,kp1,bi,bj)*
1140         &           (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
1141         &        +maskC( i ,j,kp1,bi,bj)*
1142         &           (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
1143         &       )      ) * 0.25 _d 0
1144                tmp1k(i,j) = - dyG(i,j,bi,bj)*drF(k)
1145         &                 * _hFacW(i,j,k,bi,bj)
1146         &                 *op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1147         &                 * dTdz
1148              ENDDO
1149             ENDDO
1150             CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KuzTz', k,1,2,bi,bj,myThid)
1151    C         print *,'SM_KuzTz',k,tmp1k
1152            ENDIF
1153    
1154            IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1155             DO j=1,sNy+1
1156              DO i=1,sNx
1157    C-        Vertical gradients interpolated to V points
1158               dTdz =   op5*(
1159         &      +op5*recip_drC(k)*
1160         &       ( maskC(i,j-1,k,bi,bj)*
1161         &           (Theta(i,j-1,km1,bi,bj)-Theta(i,j-1,k,bi,bj))
1162         &        +maskC(i, j ,k,bi,bj)*
1163         &           (Theta(i, j ,km1,bi,bj)-Theta(i, j ,k,bi,bj))
1164         &       )
1165         &      +op5*recip_drC(kp1)*
1166         &       ( maskC(i,j-1,kp1,bi,bj)*
1167         &           (Theta(i,j-1,k,bi,bj)-Theta(i,j-1,kp1,bi,bj))
1168         &        +maskC(i, j ,kp1,bi,bj)*
1169         &           (Theta(i, j ,k,bi,bj)-Theta(i, j ,kp1,bi,bj))
1170         &       )      )
1171               tmp1k(i,j) = - dxG(i,j,bi,bj)*drF(k)
1172         &                * _hFacS(i,j,k,bi,bj)
1173         &                *op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1174         &                * dTdz
1175              ENDDO
1176             ENDDO
1177             CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KvzTz', k,1,2,bi,bj,myThid)
1178    C         print *,'SM_KvzTz',k,tmp1k
1179            ENDIF
1180    
1181            IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1182             DO j=jMin,jMax
1183              DO i=iMin,iMax
1184    C-      Horizontal gradients interpolated to W points
1185               dTdx = op5*(
1186         &           +op5*(_maskW(i+1,j,k,bi,bj)
1187         &            *_recip_dxC(i+1,j,bi,bj)*
1188         &             (Theta(i+1,j,k,bi,bj)-Theta(i,j,k,bi,bj))
1189         &              +_maskW(i,j,k,bi,bj)
1190         &            *_recip_dxC(i,j,bi,bj)*
1191         &             (Theta(i,j,k,bi,bj)-Theta(i-1,j,k,bi,bj)))
1192         &           +op5*(_maskW(i+1,j,k-1,bi,bj)
1193         &            *_recip_dxC(i+1,j,bi,bj)*
1194         &             (Theta(i+1,j,k-1,bi,bj)-Theta(i,j,k-1,bi,bj))
1195         &              +_maskW(i,j,k-1,bi,bj)
1196         &            *_recip_dxC(i,j,bi,bj)*
1197         &             (Theta(i,j,k-1,bi,bj)-Theta(i-1,j,k-1,bi,bj)))
1198         &       )
1199    
1200               dTdy = op5*(
1201         &           +op5*(_maskS(i,j,k,bi,bj)
1202         &            *_recip_dyC(i,j,bi,bj)*
1203         &             (Theta(i,j,k,bi,bj)-Theta(i,j-1,k,bi,bj))
1204         &              +_maskS(i,j+1,k,bi,bj)
1205         &            *_recip_dyC(i,j+1,bi,bj)*
1206         &             (Theta(i,j+1,k,bi,bj)-Theta(i,j,k,bi,bj)))
1207         &           +op5*(_maskS(i,j,k-1,bi,bj)
1208         &            *_recip_dyC(i,j,bi,bj)*
1209         &             (Theta(i,j,k-1,bi,bj)-Theta(i,j-1,k-1,bi,bj))
1210         &              +_maskS(i,j+1,k-1,bi,bj)
1211         &            *_recip_dyC(i,j+1,bi,bj)*
1212         &             (Theta(i,j+1,k-1,bi,bj)-Theta(i,j,k-1,bi,bj)))
1213         &       )
1214    
1215               tmp1k(i,j) = - _rA(i,j,bi,bj)
1216         &                     *(op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))*dTdx
1217         &                       +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))*dTdy)
1218              ENDDO
1219             ENDDO
1220             CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', k,1,2,bi,bj,myThid)
1221    C         print *,'SM_KrddT',k,tmp1k
1222            ENDIF
1223           ENDIF
1224    #endif
1225           DO j=1-Oly+1,sNy+Oly-1
1226            DO i=1-Olx+1,sNx+Olx-1
1227             SM_PsiXm1(i,j)=SM_PsiX(i,j)
1228             SM_PsiYm1(i,j)=SM_PsiY(i,j)
1229             tmp1k(i,j)=0 _d 0
1230            ENDDO
1231           ENDDO
1232          ENDDO
1233    
1234    CBFK Always Zero at the bottom.
1235          IF ( DIAGNOSTICS_IS_ON('SM_ubT  ',myThid) ) THEN
1236           CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT  ', Nr,1,2,bi,bj,myThid)
1237          ENDIF
1238          IF ( DIAGNOSTICS_IS_ON('SM_vbT  ',myThid) ) THEN
1239           CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT  ', Nr,1,2,bi,bj,myThid)
1240          ENDIF
1241          IF ( DIAGNOSTICS_IS_ON('SM_wbT  ',myThid) ) THEN
1242           CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT  ', Nr,1,2,bi,bj,myThid)
1243          ENDIF
1244          IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1245           CALL DIAGNOSTICS_FILL(tmp1k,'SM_KuzTz', Nr,1,2,bi,bj,myThid)
1246          ENDIF
1247          IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1248           CALL DIAGNOSTICS_FILL(tmp1k,'SM_KvzTz', Nr,1,2,bi,bj,myThid)
1249          ENDIF
1250          IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1251           CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', Nr,1,2,bi,bj,myThid)
1252          ENDIF
1253    #endif
1254    
1255  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
1256  C--   Time-average  C--   Time-average
1257        IF ( taveFreq.GT.0. ) THEN        IF ( taveFreq.GT.0. ) THEN
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22