/[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.4 by dimitri, Fri May 30 23:18:08 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              hml=hMixLayer(i,j,bi,bj)
1014              IF (hml.gt.0 _d 0) THEN
1015               recip_hml=1 _d 0/hml
1016              ELSE
1017               recip_hml=0 _d 0
1018              ENDIF
1019    CBFK We calculate the h^2 mu(z) factor only on w points.  
1020    CBFK It is possible that we might need to calculate it
1021    CBFK on Psi or u,v points independently to prevent spurious
1022    CBFK entrainment.  Unlikely that this will be major
1023    CBFK (it wasnt in offline testing).
1024              qfac=(2*rf(k)*recip_hml+1 _d 0)**2
1025              hsqmu=(1 _d 0-qfac)*(1 _d 0+(5 _d 0)*qfac/21 _d 0)
1026              hsqmu=max(0 _d 0, hsqmu)*hml**2
1027              SM_Lf(i,j)=SM_Lf(i,j)*hsqmu
1028             ENDDO
1029            ENDDO
1030    CBFK Now interpolate to match locations
1031            DO j=1-Oly+1,sNy+Oly-1
1032             DO i=1-Olx+1,sNx+Olx-1
1033    C SM_Lf coefficients are on rVel points
1034    C Psix are on faces above U
1035              SM_PsiX(i,j)=op5*(SM_Lf(i+1,j)+SM_Lf(i,j))*dBdxAV(i,j)
1036         &                        *_maskW(i,j,k,bi,bj)
1037    C Psiy are on faces above V
1038              SM_PsiY(i,j)=op5*(SM_Lf(i,j+1)+SM_Lf(i,j))*dBdyAV(i,j)
1039         &                        *_maskS(i,j,k,bi,bj)
1040    
1041    #ifndef GM_BOLUS_ADVEC
1042    C Kwx,Kwy are on rVel Points
1043              Kwx(i,j,k,bi,bj)  = Kwx(i,j,k,bi,bj)
1044         &                      +op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))
1045              Kwy(i,j,k,bi,bj)  = Kwy(i,j,k,bi,bj)
1046         &                      +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))
1047    #ifdef GM_EXTRA_DIAGONAL
1048              IF (GM_ExtraDiag) THEN
1049    C Kuz,Kvz are on u,v Points
1050               Kuz(i,j,k,bi,bj)  = Kuz(i,j,k,bi,bj)
1051         &                    -op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1052               Kvz(i,j,k,bi,bj)  = Kvz(i,j,k,bi,bj)
1053         &                    -op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1054              ENDIF
1055    #endif
1056    #else
1057              IF (GM_AdvForm) THEN
1058               GM_PsiX(i,j,k,bi,bj)=GM_PsiX(i,j,k,bi,bj)+SM_PsiX(i,j)
1059               GM_PsiY(i,j,k,bi,bj)=GM_PsiY(i,j,k,bi,bj)+SM_PsiY(i,j)
1060              ENDIF
1061    #endif
1062             ENDDO
1063            ENDDO
1064           ELSE
1065            DO j=1-Oly+1,sNy+Oly-1
1066             DO i=1-Olx+1,sNx+Olx-1
1067              SM_PsiX(i,j)=0. _d 0
1068              SM_PsiY(i,j)=0. _d 0
1069             ENDDO
1070            ENDDO
1071           ENDIF  
1072    
1073    #ifdef ALLOW_DIAGNOSTICS
1074           IF ( useDiagnostics ) THEN
1075            IF ( DIAGNOSTICS_IS_ON('SM_PsiX ',myThid) ) THEN
1076             CALL DIAGNOSTICS_FILL(SM_PsiX,'SM_PsiX ',k,1,2,bi,bj,myThid)
1077            ENDIF
1078            IF ( DIAGNOSTICS_IS_ON('SM_PsiY ',myThid) ) THEN
1079             CALL DIAGNOSTICS_FILL(SM_PsiY,'SM_PsiY ',k,1,2,bi,bj,myThid)
1080            ENDIF
1081    
1082    CBFK Note:  for comparision, you can diagnose the bolus form
1083    CBFK or the Kappa form in the same simulation, regardless of other
1084    CBFK settings
1085            IF ( DIAGNOSTICS_IS_ON('SM_ubT  ',myThid) ) THEN
1086             DO j=jMin,jMax
1087              DO i=iMin,iMax
1088               tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiX(i,j)
1089         &                                -SM_PsiXm1(i,j) )
1090         &                               *maskW(i,j,km1,bi,bj)
1091         &            *op5*(Theta(i,j,km1,bi,bj)+Theta(i-1,j,km1,bi,bj))
1092              ENDDO
1093             ENDDO
1094             CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT  ', km1,1,2,bi,bj,myThid)
1095            ENDIF
1096    
1097            IF ( DIAGNOSTICS_IS_ON('SM_vbT  ',myThid) ) THEN
1098             DO j=jMin,jMax
1099              DO i=iMin,iMax
1100               tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiY(i,j)
1101         &                                -SM_PsiYm1(i,j) )
1102         &                               *maskS(i,j,km1,bi,bj)
1103         &            *op5*(Theta(i,j,km1,bi,bj)+Theta(i,j-1,km1,bi,bj))
1104              ENDDO
1105             ENDDO
1106             CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT  ', km1,1,2,bi,bj,myThid)
1107            ENDIF
1108    
1109            IF ( DIAGNOSTICS_IS_ON('SM_wbT  ',myThid) ) THEN
1110             DO j=jMin,jMax
1111              DO i=iMin,iMax
1112               tmp1k(i,j) =
1113         &        (dyG(i+1,j,bi,bj)*SM_PsiX(i+1,j)
1114         &        -dyG( i ,j,bi,bj)*SM_PsiX( i ,j)
1115         &        +dxG(i,j+1,bi,bj)*SM_PsiY(i,j+1)
1116         &        -dxG(i, j ,bi,bj)*SM_PsiY(i, j ))
1117         &             *op5*(Theta(i,j,k,bi,bj)+Theta(i,j,km1,bi,bj))
1118              ENDDO
1119             ENDDO
1120             CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT  ', k,1,2,bi,bj,myThid)
1121    C         print *,'SM_wbT',k,tmp1k
1122            ENDIF
1123    
1124            IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1125             DO j=1,sNy
1126              DO i=1,sNx+1
1127    C-        Vertical gradients interpolated to U points
1128               dTdz = (
1129         &      +recip_drC(k)*
1130         &       ( maskC(i-1,j,k,bi,bj)*
1131         &           (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
1132         &        +maskC( i ,j,k,bi,bj)*
1133         &           (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
1134         &       )
1135         &      +recip_drC(kp1)*
1136         &       ( maskC(i-1,j,kp1,bi,bj)*
1137         &           (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
1138         &        +maskC( i ,j,kp1,bi,bj)*
1139         &           (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
1140         &       )      ) * 0.25 _d 0
1141                tmp1k(i,j) = - dyG(i,j,bi,bj)*drF(k)
1142         &                 * _hFacW(i,j,k,bi,bj)
1143         &                 *op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1144         &                 * dTdz
1145              ENDDO
1146             ENDDO
1147             CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KuzTz', k,1,2,bi,bj,myThid)
1148    C         print *,'SM_KuzTz',k,tmp1k
1149            ENDIF
1150    
1151            IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1152             DO j=1,sNy+1
1153              DO i=1,sNx
1154    C-        Vertical gradients interpolated to V points
1155               dTdz =   op5*(
1156         &      +op5*recip_drC(k)*
1157         &       ( maskC(i,j-1,k,bi,bj)*
1158         &           (Theta(i,j-1,km1,bi,bj)-Theta(i,j-1,k,bi,bj))
1159         &        +maskC(i, j ,k,bi,bj)*
1160         &           (Theta(i, j ,km1,bi,bj)-Theta(i, j ,k,bi,bj))
1161         &       )
1162         &      +op5*recip_drC(kp1)*
1163         &       ( maskC(i,j-1,kp1,bi,bj)*
1164         &           (Theta(i,j-1,k,bi,bj)-Theta(i,j-1,kp1,bi,bj))
1165         &        +maskC(i, j ,kp1,bi,bj)*
1166         &           (Theta(i, j ,k,bi,bj)-Theta(i, j ,kp1,bi,bj))
1167         &       )      )
1168               tmp1k(i,j) = - dxG(i,j,bi,bj)*drF(k)
1169         &                * _hFacS(i,j,k,bi,bj)
1170         &                *op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1171         &                * dTdz
1172              ENDDO
1173             ENDDO
1174             CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KvzTz', k,1,2,bi,bj,myThid)
1175    C         print *,'SM_KvzTz',k,tmp1k
1176            ENDIF
1177    
1178            IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1179             DO j=jMin,jMax
1180              DO i=iMin,iMax
1181    C-      Horizontal gradients interpolated to W points
1182               dTdx = op5*(
1183         &           +op5*(_maskW(i+1,j,k,bi,bj)
1184         &            *_recip_dxC(i+1,j,bi,bj)*
1185         &             (Theta(i+1,j,k,bi,bj)-Theta(i,j,k,bi,bj))
1186         &              +_maskW(i,j,k,bi,bj)
1187         &            *_recip_dxC(i,j,bi,bj)*
1188         &             (Theta(i,j,k,bi,bj)-Theta(i-1,j,k,bi,bj)))
1189         &           +op5*(_maskW(i+1,j,k-1,bi,bj)
1190         &            *_recip_dxC(i+1,j,bi,bj)*
1191         &             (Theta(i+1,j,k-1,bi,bj)-Theta(i,j,k-1,bi,bj))
1192         &              +_maskW(i,j,k-1,bi,bj)
1193         &            *_recip_dxC(i,j,bi,bj)*
1194         &             (Theta(i,j,k-1,bi,bj)-Theta(i-1,j,k-1,bi,bj)))
1195         &       )
1196    
1197               dTdy = op5*(
1198         &           +op5*(_maskS(i,j,k,bi,bj)
1199         &            *_recip_dyC(i,j,bi,bj)*
1200         &             (Theta(i,j,k,bi,bj)-Theta(i,j-1,k,bi,bj))
1201         &              +_maskS(i,j+1,k,bi,bj)
1202         &            *_recip_dyC(i,j+1,bi,bj)*
1203         &             (Theta(i,j+1,k,bi,bj)-Theta(i,j,k,bi,bj)))
1204         &           +op5*(_maskS(i,j,k-1,bi,bj)
1205         &            *_recip_dyC(i,j,bi,bj)*
1206         &             (Theta(i,j,k-1,bi,bj)-Theta(i,j-1,k-1,bi,bj))
1207         &              +_maskS(i,j+1,k-1,bi,bj)
1208         &            *_recip_dyC(i,j+1,bi,bj)*
1209         &             (Theta(i,j+1,k-1,bi,bj)-Theta(i,j,k-1,bi,bj)))
1210         &       )
1211    
1212               tmp1k(i,j) = - _rA(i,j,bi,bj)
1213         &                     *(op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))*dTdx
1214         &                       +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))*dTdy)
1215              ENDDO
1216             ENDDO
1217             CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', k,1,2,bi,bj,myThid)
1218    C         print *,'SM_KrddT',k,tmp1k
1219            ENDIF
1220           ENDIF
1221    #endif
1222           DO j=1-Oly+1,sNy+Oly-1
1223            DO i=1-Olx+1,sNx+Olx-1
1224             SM_PsiXm1(i,j)=SM_PsiX(i,j)
1225             SM_PsiYm1(i,j)=SM_PsiY(i,j)
1226             tmp1k(i,j)=0 _d 0
1227            ENDDO
1228           ENDDO
1229          ENDDO
1230    
1231    CBFK Always Zero at the bottom.
1232          IF ( DIAGNOSTICS_IS_ON('SM_ubT  ',myThid) ) THEN
1233           CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT  ', Nr,1,2,bi,bj,myThid)
1234          ENDIF
1235          IF ( DIAGNOSTICS_IS_ON('SM_vbT  ',myThid) ) THEN
1236           CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT  ', Nr,1,2,bi,bj,myThid)
1237          ENDIF
1238          IF ( DIAGNOSTICS_IS_ON('SM_wbT  ',myThid) ) THEN
1239           CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT  ', Nr,1,2,bi,bj,myThid)
1240          ENDIF
1241          IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1242           CALL DIAGNOSTICS_FILL(tmp1k,'SM_KuzTz', Nr,1,2,bi,bj,myThid)
1243          ENDIF
1244          IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1245           CALL DIAGNOSTICS_FILL(tmp1k,'SM_KvzTz', Nr,1,2,bi,bj,myThid)
1246          ENDIF
1247          IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1248           CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', Nr,1,2,bi,bj,myThid)
1249          ENDIF
1250    #endif
1251    
1252  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
1253  C--   Time-average  C--   Time-average
1254        IF ( taveFreq.GT.0. ) THEN        IF ( taveFreq.GT.0. ) THEN
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22