/[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.3 by dimitri, Fri May 30 22:24:25 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 450  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bi Line 531  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bi
531  #else  #else
532           Kgm_tmp = GM_isopycK           Kgm_tmp = GM_isopycK
533  #endif  #endif
534  #ifdef ALLOW_KAPGM_CONTROL  #if (defined (ALLOW_AUTODIFF) && defined (ALLOW_KAPGM_CONTROL))
535       &           + GM_skewflx*kapgm(i,j,k,bi,bj)       &           + GM_skewflx*kapgm(i,j,k,bi,bj)
536  #else  #else
537       &           + GM_skewflx*GM_background_K       &           + GM_skewflx*GM_background_K
# 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 609  CADJ STORE taperFct(:,:)     = comlev1_b Line 715  CADJ STORE taperFct(:,:)     = comlev1_b
715  #else  #else
716       &     ( GM_isopycK       &     ( GM_isopycK
717  #endif  #endif
718  #ifdef ALLOW_KAPGM_CONTROL  #if (defined (ALLOW_AUTODIFF) && defined (ALLOW_KAPGM_CONTROL))
719       &     - GM_skewflx*kapgm(i,j,k,bi,bj)       &     - GM_skewflx*kapgm(i,j,k,bi,bj)
720  #else  #else
721       &     - GM_skewflx*GM_background_K       &     - GM_skewflx*GM_background_K
# Line 722  C     Gradient of Sigma at V points Line 828  C     Gradient of Sigma at V points
828           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 )
829       &                      +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1       &                      +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1
830       &                      )*_maskS(i,j,k,bi,bj)       &                      )*_maskS(i,j,k,bi,bj)
831    
832    #ifdef GM_SUBMESO
833    C--     Depth average of SigmaY at V points
834    C       compute depth average from surface down to the MixLayer depth
835             IF (k.GT.1) THEN
836              IF (-rC(k-1).LT.locMixLayer(i,j) ) THEN
837               IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
838                integrDepth = -rC( k )
839    C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments
840                integrDepth = MIN( integrDepth, locMixLayer(i,j) )
841    C       Distance between level center above and the integration depth
842                deltaH = integrDepth + rC(k-1)
843    C       If negative then we are below the integration level
844    C       (cannot be the case with 2 conditions on maskC & -rC(k-1))
845    C       If positive we limit this to the distance from center above
846                deltaH = MIN( deltaH, drC(k) )
847    C       Now we convert deltaH to a non-dimensional fraction
848                deltaH = deltaH/( integrDepth+rC(1) )
849                dBdyAV(i,j) = dBdyAV(i,j)
850         &            +dSigmaDy(i,j)*deltaH*recip_rhoConst*gravity
851               ENDIF
852              ENDIF
853             ENDIF
854    #endif
855    
856          ENDDO          ENDDO
857         ENDDO         ENDDO
858    
# Line 797  CADJ STORE taperFct(:,:)     = comlev1_b Line 928  CADJ STORE taperFct(:,:)     = comlev1_b
928  #else  #else
929       &     ( GM_isopycK       &     ( GM_isopycK
930  #endif  #endif
931  #ifdef ALLOW_KAPGM_CONTROL  #if (defined (ALLOW_AUTODIFF) && defined (ALLOW_KAPGM_CONTROL))
932       &     - GM_skewflx*kapgm(i,j,k,bi,bj)       &     - GM_skewflx*kapgm(i,j,k,bi,bj)
933  #else  #else
934       &     - GM_skewflx*GM_background_K       &     - GM_skewflx*GM_background_K
# Line 830  C         store in tmp1k Kvz_Redi Line 961  C         store in tmp1k Kvz_Redi
961          ENDDO          ENDDO
962          DO j=1,sNy+1          DO j=1,sNy+1
963           DO i=1,sNx           DO i=1,sNx
964  C-        Vertical gradients interpolated to U points  C-        Vertical gradients interpolated to V points
965            dTdz = (            dTdz =   op5*(
966       &     +recip_drC(k)*       &   +op5*recip_drC(k)*
967       &       ( maskC(i,j-1,k,bi,bj)*       &       ( maskC(i,j-1,k,bi,bj)*
968       &           (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))
969       &        +maskC(i, j ,k,bi,bj)*       &        +maskC(i, j ,k,bi,bj)*
970       &           (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))       &           (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))
971       &       )       &       )
972       &     +recip_drC(kp1)*       &   +op5*recip_drC(kp1)*
973       &       ( maskC(i,j-1,kp1,bi,bj)*       &       ( maskC(i,j-1,kp1,bi,bj)*
974       &           (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))
975       &        +maskC(i, j ,kp1,bi,bj)*       &        +maskC(i, j ,kp1,bi,bj)*
976       &           (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))       &           (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))
977       &       )      ) * 0.25 _d 0       &       )      )
978             tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)             tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)
979       &                * _hFacS(i,j,k,bi,bj)       &                * _hFacS(i,j,k,bi,bj)
980       &                * tmp1k(i,j) * dTdz       &                * tmp1k(i,j) * dTdz
# Line 869  C-- end 3rd  loop on vertical level inde Line 1000  C-- end 3rd  loop on vertical level inde
1000        ENDIF        ENDIF
1001  #endif  #endif
1002    
1003    
1004    #ifdef GM_SUBMESO
1005    CBFK Add the submesoscale contribution, in a 4th k loop
1006          DO k=1,Nr
1007           km1=max(1,k-1)
1008           IF ((k.gt.1).and.(-rF(k-1) .lt. mlmax)) THEN
1009            kp1 = MIN(k+1,Nr)
1010    CBFK Add in the mu vertical structure factor
1011            DO j=1-Oly+1,sNy+Oly-1
1012             DO i=1-Olx+1,sNx+Olx-1
1013    #ifdef ALLOW_KPP
1014              hml=KPPhbl(i,j,bi,bj)
1015    #else
1016              hml=hMixLayer(i,j,bi,bj)
1017    #endif
1018              IF (hml.gt.0 _d 0) THEN
1019               recip_hml=1 _d 0/hml
1020              ELSE
1021               recip_hml=0 _d 0
1022              ENDIF
1023    CBFK We calculate the h^2 mu(z) factor only on w points.  
1024    CBFK It is possible that we might need to calculate it
1025    CBFK on Psi or u,v points independently to prevent spurious
1026    CBFK entrainment.  Unlikely that this will be major
1027    CBFK (it wasnt in offline testing).
1028              qfac=(2*rf(k)*recip_hml+1 _d 0)**2
1029              hsqmu=(1 _d 0-qfac)*(1 _d 0+(5 _d 0)*qfac/21 _d 0)
1030              hsqmu=max(0 _d 0, hsqmu)*hml**2
1031              SM_Lf(i,j)=SM_Lf(i,j)*hsqmu
1032             ENDDO
1033            ENDDO
1034    CBFK Now interpolate to match locations
1035            DO j=1-Oly+1,sNy+Oly-1
1036             DO i=1-Olx+1,sNx+Olx-1
1037    C SM_Lf coefficients are on rVel points
1038    C Psix are on faces above U
1039              SM_PsiX(i,j)=op5*(SM_Lf(i+1,j)+SM_Lf(i,j))*dBdxAV(i,j)
1040         &                        *_maskW(i,j,k,bi,bj)
1041    C Psiy are on faces above V
1042              SM_PsiY(i,j)=op5*(SM_Lf(i,j+1)+SM_Lf(i,j))*dBdyAV(i,j)
1043         &                        *_maskS(i,j,k,bi,bj)
1044    
1045    #ifndef GM_BOLUS_ADVEC
1046    C Kwx,Kwy are on rVel Points
1047              Kwx(i,j,k,bi,bj)  = Kwx(i,j,k,bi,bj)
1048         &                      +op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))
1049              Kwy(i,j,k,bi,bj)  = Kwy(i,j,k,bi,bj)
1050         &                      +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))
1051    #ifdef GM_EXTRA_DIAGONAL
1052              IF (GM_ExtraDiag) THEN
1053    C Kuz,Kvz are on u,v Points
1054               Kuz(i,j,k,bi,bj)  = Kuz(i,j,k,bi,bj)
1055         &                    -op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1056               Kvz(i,j,k,bi,bj)  = Kvz(i,j,k,bi,bj)
1057         &                    -op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1058              ENDIF
1059    #endif
1060    #else
1061              IF (GM_AdvForm) THEN
1062               GM_PsiX(i,j,k,bi,bj)=GM_PsiX(i,j,k,bi,bj)+SM_PsiX(i,j)
1063               GM_PsiY(i,j,k,bi,bj)=GM_PsiY(i,j,k,bi,bj)+SM_PsiY(i,j)
1064              ENDIF
1065    #endif
1066             ENDDO
1067            ENDDO
1068           ELSE
1069            DO j=1-Oly+1,sNy+Oly-1
1070             DO i=1-Olx+1,sNx+Olx-1
1071              SM_PsiX(i,j)=0. _d 0
1072              SM_PsiY(i,j)=0. _d 0
1073             ENDDO
1074            ENDDO
1075           ENDIF  
1076    
1077    #ifdef ALLOW_DIAGNOSTICS
1078           IF ( useDiagnostics ) THEN
1079            IF ( DIAGNOSTICS_IS_ON('SM_PsiX ',myThid) ) THEN
1080             CALL DIAGNOSTICS_FILL(SM_PsiX,'SM_PsiX ',k,1,2,bi,bj,myThid)
1081            ENDIF
1082            IF ( DIAGNOSTICS_IS_ON('SM_PsiY ',myThid) ) THEN
1083             CALL DIAGNOSTICS_FILL(SM_PsiY,'SM_PsiY ',k,1,2,bi,bj,myThid)
1084            ENDIF
1085    
1086    CBFK Note:  for comparision, you can diagnose the bolus form
1087    CBFK or the Kappa form in the same simulation, regardless of other
1088    CBFK settings
1089            IF ( DIAGNOSTICS_IS_ON('SM_ubT  ',myThid) ) THEN
1090             DO j=jMin,jMax
1091              DO i=iMin,iMax
1092               tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiX(i,j)
1093         &                                -SM_PsiXm1(i,j) )
1094         &                               *maskW(i,j,km1,bi,bj)
1095         &            *op5*(Theta(i,j,km1,bi,bj)+Theta(i-1,j,km1,bi,bj))
1096              ENDDO
1097             ENDDO
1098             CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT  ', km1,1,2,bi,bj,myThid)
1099            ENDIF
1100    
1101            IF ( DIAGNOSTICS_IS_ON('SM_vbT  ',myThid) ) THEN
1102             DO j=jMin,jMax
1103              DO i=iMin,iMax
1104               tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiY(i,j)
1105         &                                -SM_PsiYm1(i,j) )
1106         &                               *maskS(i,j,km1,bi,bj)
1107         &            *op5*(Theta(i,j,km1,bi,bj)+Theta(i,j-1,km1,bi,bj))
1108              ENDDO
1109             ENDDO
1110             CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT  ', km1,1,2,bi,bj,myThid)
1111            ENDIF
1112    
1113            IF ( DIAGNOSTICS_IS_ON('SM_wbT  ',myThid) ) THEN
1114             DO j=jMin,jMax
1115              DO i=iMin,iMax
1116               tmp1k(i,j) =
1117         &        (dyG(i+1,j,bi,bj)*SM_PsiX(i+1,j)
1118         &        -dyG( i ,j,bi,bj)*SM_PsiX( i ,j)
1119         &        +dxG(i,j+1,bi,bj)*SM_PsiY(i,j+1)
1120         &        -dxG(i, j ,bi,bj)*SM_PsiY(i, j ))
1121         &             *op5*(Theta(i,j,k,bi,bj)+Theta(i,j,km1,bi,bj))
1122              ENDDO
1123             ENDDO
1124             CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT  ', k,1,2,bi,bj,myThid)
1125    C         print *,'SM_wbT',k,tmp1k
1126            ENDIF
1127    
1128            IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1129             DO j=1,sNy
1130              DO i=1,sNx+1
1131    C-        Vertical gradients interpolated to U points
1132               dTdz = (
1133         &      +recip_drC(k)*
1134         &       ( maskC(i-1,j,k,bi,bj)*
1135         &           (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
1136         &        +maskC( i ,j,k,bi,bj)*
1137         &           (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
1138         &       )
1139         &      +recip_drC(kp1)*
1140         &       ( maskC(i-1,j,kp1,bi,bj)*
1141         &           (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
1142         &        +maskC( i ,j,kp1,bi,bj)*
1143         &           (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
1144         &       )      ) * 0.25 _d 0
1145                tmp1k(i,j) = - dyG(i,j,bi,bj)*drF(k)
1146         &                 * _hFacW(i,j,k,bi,bj)
1147         &                 *op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1148         &                 * dTdz
1149              ENDDO
1150             ENDDO
1151             CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KuzTz', k,1,2,bi,bj,myThid)
1152    C         print *,'SM_KuzTz',k,tmp1k
1153            ENDIF
1154    
1155            IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1156             DO j=1,sNy+1
1157              DO i=1,sNx
1158    C-        Vertical gradients interpolated to V points
1159               dTdz =   op5*(
1160         &      +op5*recip_drC(k)*
1161         &       ( maskC(i,j-1,k,bi,bj)*
1162         &           (Theta(i,j-1,km1,bi,bj)-Theta(i,j-1,k,bi,bj))
1163         &        +maskC(i, j ,k,bi,bj)*
1164         &           (Theta(i, j ,km1,bi,bj)-Theta(i, j ,k,bi,bj))
1165         &       )
1166         &      +op5*recip_drC(kp1)*
1167         &       ( maskC(i,j-1,kp1,bi,bj)*
1168         &           (Theta(i,j-1,k,bi,bj)-Theta(i,j-1,kp1,bi,bj))
1169         &        +maskC(i, j ,kp1,bi,bj)*
1170         &           (Theta(i, j ,k,bi,bj)-Theta(i, j ,kp1,bi,bj))
1171         &       )      )
1172               tmp1k(i,j) = - dxG(i,j,bi,bj)*drF(k)
1173         &                * _hFacS(i,j,k,bi,bj)
1174         &                *op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1175         &                * dTdz
1176              ENDDO
1177             ENDDO
1178             CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KvzTz', k,1,2,bi,bj,myThid)
1179    C         print *,'SM_KvzTz',k,tmp1k
1180            ENDIF
1181    
1182            IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1183             DO j=jMin,jMax
1184              DO i=iMin,iMax
1185    C-      Horizontal gradients interpolated to W points
1186               dTdx = op5*(
1187         &           +op5*(_maskW(i+1,j,k,bi,bj)
1188         &            *_recip_dxC(i+1,j,bi,bj)*
1189         &             (Theta(i+1,j,k,bi,bj)-Theta(i,j,k,bi,bj))
1190         &              +_maskW(i,j,k,bi,bj)
1191         &            *_recip_dxC(i,j,bi,bj)*
1192         &             (Theta(i,j,k,bi,bj)-Theta(i-1,j,k,bi,bj)))
1193         &           +op5*(_maskW(i+1,j,k-1,bi,bj)
1194         &            *_recip_dxC(i+1,j,bi,bj)*
1195         &             (Theta(i+1,j,k-1,bi,bj)-Theta(i,j,k-1,bi,bj))
1196         &              +_maskW(i,j,k-1,bi,bj)
1197         &            *_recip_dxC(i,j,bi,bj)*
1198         &             (Theta(i,j,k-1,bi,bj)-Theta(i-1,j,k-1,bi,bj)))
1199         &       )
1200    
1201               dTdy = op5*(
1202         &           +op5*(_maskS(i,j,k,bi,bj)
1203         &            *_recip_dyC(i,j,bi,bj)*
1204         &             (Theta(i,j,k,bi,bj)-Theta(i,j-1,k,bi,bj))
1205         &              +_maskS(i,j+1,k,bi,bj)
1206         &            *_recip_dyC(i,j+1,bi,bj)*
1207         &             (Theta(i,j+1,k,bi,bj)-Theta(i,j,k,bi,bj)))
1208         &           +op5*(_maskS(i,j,k-1,bi,bj)
1209         &            *_recip_dyC(i,j,bi,bj)*
1210         &             (Theta(i,j,k-1,bi,bj)-Theta(i,j-1,k-1,bi,bj))
1211         &              +_maskS(i,j+1,k-1,bi,bj)
1212         &            *_recip_dyC(i,j+1,bi,bj)*
1213         &             (Theta(i,j+1,k-1,bi,bj)-Theta(i,j,k-1,bi,bj)))
1214         &       )
1215    
1216               tmp1k(i,j) = - _rA(i,j,bi,bj)
1217         &                     *(op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))*dTdx
1218         &                       +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))*dTdy)
1219              ENDDO
1220             ENDDO
1221             CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', k,1,2,bi,bj,myThid)
1222    C         print *,'SM_KrddT',k,tmp1k
1223            ENDIF
1224           ENDIF
1225    #endif
1226           DO j=1-Oly+1,sNy+Oly-1
1227            DO i=1-Olx+1,sNx+Olx-1
1228             SM_PsiXm1(i,j)=SM_PsiX(i,j)
1229             SM_PsiYm1(i,j)=SM_PsiY(i,j)
1230             tmp1k(i,j)=0 _d 0
1231            ENDDO
1232           ENDDO
1233          ENDDO
1234    
1235    CBFK Always Zero at the bottom.
1236          IF ( DIAGNOSTICS_IS_ON('SM_ubT  ',myThid) ) THEN
1237           CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT  ', Nr,1,2,bi,bj,myThid)
1238          ENDIF
1239          IF ( DIAGNOSTICS_IS_ON('SM_vbT  ',myThid) ) THEN
1240           CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT  ', Nr,1,2,bi,bj,myThid)
1241          ENDIF
1242          IF ( DIAGNOSTICS_IS_ON('SM_wbT  ',myThid) ) THEN
1243           CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT  ', Nr,1,2,bi,bj,myThid)
1244          ENDIF
1245          IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1246           CALL DIAGNOSTICS_FILL(tmp1k,'SM_KuzTz', Nr,1,2,bi,bj,myThid)
1247          ENDIF
1248          IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1249           CALL DIAGNOSTICS_FILL(tmp1k,'SM_KvzTz', Nr,1,2,bi,bj,myThid)
1250          ENDIF
1251          IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1252           CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', Nr,1,2,bi,bj,myThid)
1253          ENDIF
1254    #endif
1255    
1256  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
1257  C--   Time-average  C--   Time-average
1258        IF ( taveFreq.GT.0. ) THEN        IF ( taveFreq.GT.0. ) THEN
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22