/[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.2 by dimitri, Fri May 30 22:13:42 2008 UTC revision 1.6 by zhc, Fri Mar 19 19:17:45 2010 UTC
# Line 5  C $Name$ Line 5  C $Name$
5  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
6  # include "KPP_OPTIONS.h"  # include "KPP_OPTIONS.h"
7  #endif  #endif
 #undef OLD_VISBECK_CALC  
8    
9  CBOP  CBOP
10  C     !ROUTINE: GMREDI_CALC_TENSOR  C     !ROUTINE: GMREDI_CALC_TENSOR
# Line 64  CEOP Line 63  CEOP
63    
64  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
65  C     == Local variables ==  C     == Local variables ==
66        INTEGER i,j,k,kp1        INTEGER i,j,k
67        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
68        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
69        _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
# Line 72  C     == Local variables == Line 71  C     == Local variables ==
71        _RL dSigmaDr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
72        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
73        _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
74        _RL maskp1, Kgm_tmp        _RL Kgm_tmp
       _RL deltaH, integrDepth  
75        _RL ldd97_LrhoC(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL ldd97_LrhoC(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
76        _RL ldd97_LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL ldd97_LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
77        _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 83  C     == Local variables ==
83        _RL baseSlope  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL baseSlope  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
84        _RL hTransLay  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL hTransLay  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
85        _RL recipLambda(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL recipLambda(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
86    c#if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
87          INTEGER kp1
88          _RL maskp1
89    c#endif
90    
91  #ifdef GM_SUBMESO  #ifdef GM_SUBMESO
92        _RL dBdxAV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dBdxAV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
# Line 96  C     == Local variables == Line 98  C     == Local variables ==
98        _RL SM_PsiYm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SM_PsiYm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
99        _RL hsqmu, hml, recip_hml, qfac, dS, mlmax        _RL hsqmu, hml, recip_hml, qfac, dS, mlmax
100  #endif  #endif
101                                                    
102    
103  #ifdef GM_VISBECK_VARIABLE_K  c#ifdef GM_VISBECK_VARIABLE_K
104  #ifdef OLD_VISBECK_CALC  #ifdef OLD_VISBECK_CALC
       _RL zero_rs  
       PARAMETER(zero_rs=0.D0)  
       _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, dSigmaR
108        _RL Sloc, M2loc, SNloc        _RL Sloc, M2loc
 #endif  
109  #endif  #endif
110          _RL recipMaxSlope
111          _RL deltaH, integrDepth
112          _RL N2loc, SNloc
113    c#endif /* GM_VISBECK_VARIABLE_K */
114    
115  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
116        LOGICAL  doDiagRediFlx        LOGICAL  doDiagRediFlx
117        LOGICAL  DIAGNOSTICS_IS_ON        LOGICAL  DIAGNOSTICS_IS_ON
118        EXTERNAL DIAGNOSTICS_IS_ON        EXTERNAL DIAGNOSTICS_IS_ON
119    c#if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
120        INTEGER  km1        INTEGER  km1
121        _RL dTdz, dTdx, dTdy        _RL dTdz, dTdx, dTdy
122        _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123  #endif  c#endif
124    #endif /* ALLOW_DIAGNOSTICS */
125    
126  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
127    
# Line 143  C---+----1----+----2----+----3----+----4 Line 148  C---+----1----+----2----+----3----+----4
148  #endif  #endif
149    
150  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
151          recipMaxSlope = 0. _d 0
152          IF ( GM_Visbeck_maxSlope.GT.0. _d 0 ) THEN
153            recipMaxSlope = 1. _d 0 / GM_Visbeck_maxSlope
154          ENDIF
155        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
156         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
157          VisbeckK(i,j,bi,bj) = 0. _d 0          VisbeckK(i,j,bi,bj) = 0. _d 0
# Line 221  C-- 1rst loop on k : compute Tensor Coef Line 230  C-- 1rst loop on k : compute Tensor Coef
230           locMixLayer(i,j) = 0. _d 0           locMixLayer(i,j) = 0. _d 0
231         ENDDO         ENDDO
232        ENDDO        ENDDO
233        mlmax=0. _d 0  C SM(1)
234            mlmax=0. _d 0
235  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
236        IF ( useKPP ) THEN        IF ( useKPP ) THEN
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) = KPPhbl(i,j,bi,bj)           locMixLayer(i,j) = KPPhbl(i,j,bi,bj)
240           mlmax=max(mlmax,locMixLayer(i,j))  C SM(1)  
241             mlmax=max(mlmax,locMixLayer(i,j))
242          ENDDO          ENDDO
243         ENDDO         ENDDO
244        ELSE        ELSE
# Line 237  C-- 1rst loop on k : compute Tensor Coef Line 248  C-- 1rst loop on k : compute Tensor Coef
248         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
249          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
250           locMixLayer(i,j) = hMixLayer(i,j,bi,bj)           locMixLayer(i,j) = hMixLayer(i,j,bi,bj)
251           mlmax=max(mlmax,locMixLayer(i,j))  C SM(1)
252             mlmax=max(mlmax,locMixLayer(i,j))
253          ENDDO          ENDDO
254         ENDDO         ENDDO
255        ENDIF        ENDIF
256    
257  #ifdef GM_SUBMESO  #ifdef GM_SUBMESO
258        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
259         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
260          dBdxAV(i,j)     = 0. _d 0          dBdxAV(i,j)     = 0. _d 0
261          dBdyAV(i,j)     = 0. _d 0          dBdyAV(i,j)     = 0. _d 0
# Line 251  C-- 1rst loop on k : compute Tensor Coef Line 263  C-- 1rst loop on k : compute Tensor Coef
263          SM_PsiX(i,j)      = 0. _d 0          SM_PsiX(i,j)      = 0. _d 0
264          SM_PsiY(i,j)      = 0. _d 0          SM_PsiY(i,j)      = 0. _d 0
265          SM_PsiXm1(i,j)    = 0. _d 0          SM_PsiXm1(i,j)    = 0. _d 0
266          SM_PsiXm1(i,j)    = 0. _d 0          SM_PsiYm1(i,j)    = 0. _d 0
267         ENDDO         ENDDO
268        ENDDO        ENDDO
269  #endif  #endif
270                                                                                      
271    
272        DO k=Nr,2,-1        DO k=Nr,2,-1
273    
# Line 286  C-- 1rst loop on k : compute Tensor Coef Line 299  C-- 1rst loop on k : compute Tensor Coef
299  # endif  # endif
300          ENDDO          ENDDO
301         ENDDO         ENDDO
302  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
303    
304         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
305          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
# Line 298  C      Gradient of Sigma at rVel points Line 311  C      Gradient of Sigma at rVel points
311       &                       +sigmaY(i,j+1, k )+sigmaY(i,j, k )       &                       +sigmaY(i,j+1, k )+sigmaY(i,j, k )
312       &                      )*maskC(i,j,k,bi,bj)       &                      )*maskC(i,j,k,bi,bj)
313           dSigmaDr(i,j)=sigmaR(i,j,k)           dSigmaDr(i,j)=sigmaR(i,j,k)
   
314  #ifdef GM_SUBMESO  #ifdef GM_SUBMESO
315  #ifdef GM_SUBMESO_VARYLf  #ifdef GM_SUBMESO_VARYLf
316  C--     Depth average of SigmaR at W points  C--     Depth average of SigmaR at W points
317  C       compute depth average from surface down to the MixLayer depth  C       compute depth average from surface down to the MixLayer depth
318            IF (-rC(k-1).LT.locMixLayer(i,j) ) THEN            IF (-rC(k-1).LT.locMixLayer(i,j) ) THEN
319             IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN             IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
320              integrDepth = -rC( k )              integrDepth = -rC( k )
321  C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments  C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments
322              integrDepth = MIN( integrDepth, locMixLayer(i,j) )              integrDepth = MIN( integrDepth, locMixLayer(i,j) )
323  C       Distance between level center above and the integration depth  C       Distance between level center above and the integration depth
# Line 326  C--     Store db/dr in SM_Lf for now. Line 338  C--     Store db/dr in SM_Lf for now.
338          ENDDO          ENDDO
339         ENDDO         ENDDO
340    
 #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 */  
   
341  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
342  #ifndef OLD_VISBECK_CALC  #ifndef OLD_VISBECK_CALC
343         IF ( GM_Visbeck_alpha.GT.0. .AND.         IF ( GM_Visbeck_alpha.GT.0. .AND.
344       &      -rC(k-1).LT.GM_Visbeck_depth ) THEN       &      -rC(k-1).LT.GM_Visbeck_depth ) THEN
345    
346            DO j=1-Oly,sNy+Oly
347             DO i=1-Olx,sNx+Olx
348              dSigmaDr(i,j) = MIN( sigmaR(i,j,k), 0. _d 0 )
349             ENDDO
350            ENDDO
351    
352  C--     Depth average of f/sqrt(Ri) = M^2/N^2 * N  C--     Depth average of f/sqrt(Ri) = M^2/N^2 * N
353  C       M^2 and N^2 are horizontal & vertical gradient of buoyancy.  C       M^2 and N^2 are horizontal & vertical gradient of buoyancy.
354    
# Line 354  C       GM_Visbeck_depth, whatever is th Line 363  C       GM_Visbeck_depth, whatever is th
363             integrDepth = -rC( kLowC(i,j,bi,bj) )             integrDepth = -rC( kLowC(i,j,bi,bj) )
364  C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments  C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments
365             integrDepth = MIN( integrDepth, GM_Visbeck_depth )             integrDepth = MIN( integrDepth, GM_Visbeck_depth )
366    C-      to recover "old-visbeck" form with Visbeck_minDepth = Visbeck_depth
367               integrDepth = MAX( integrDepth, GM_Visbeck_minDepth )
368  C       Distance between level center above and the integration depth  C       Distance between level center above and the integration depth
369             deltaH = integrDepth + rC(k-1)             deltaH = integrDepth + rC(k-1)
370  C       If negative then we are below the integration level  C       If negative then we are below the integration level
# Line 363  C       If positive we limit this to the Line 374  C       If positive we limit this to the
374  C       Now we convert deltaH to a non-dimensional fraction  C       Now we convert deltaH to a non-dimensional fraction
375             deltaH = deltaH/( integrDepth+rC(1) )             deltaH = deltaH/( integrDepth+rC(1) )
376    
377  C--      compute: ( M^2 * S )^1/2   (= M^2 / N since S=M^2/N^2 )  C--      compute: ( M^2 * S )^1/2   (= S*N since S=M^2/N^2 )
378    C        a 5 points average gives a more "homogeneous" formulation
379    C        (same stencil and same weights as for dSigmaH calculation)
380               dSigmaR = ( dSigmaDr(i,j)*4. _d 0
381         &               + dSigmaDr(i-1,j)
382         &               + dSigmaDr(i+1,j)
383         &               + dSigmaDr(i,j-1)
384         &               + dSigmaDr(i,j+1)
385         &               )/( 4. _d 0
386         &                 + maskC(i-1,j,k,bi,bj)
387         &                 + maskC(i+1,j,k,bi,bj)
388         &                 + maskC(i,j-1,k,bi,bj)
389         &                 + maskC(i,j+1,k,bi,bj)
390         &                 )
391             dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)             dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)
392       &             + dSigmaDy(i,j)*dSigmaDy(i,j)       &             + dSigmaDy(i,j)*dSigmaDy(i,j)
393             IF ( dSigmaH .GT. 0. _d 0 ) THEN             IF ( dSigmaH .GT. 0. _d 0 ) THEN
394               dSigmaH = SQRT( dSigmaH )               dSigmaH = SQRT( dSigmaH )
395  C-       compute slope, limited by GM_maxSlope:  C-       compute slope, limited by GM_Visbeck_maxSlope:
396               IF ( -dSigmaDr(i,j).GT.dSigmaH*GM_rMaxSlope ) THEN               IF ( -dSigmaR.GT.dSigmaH*recipMaxSlope ) THEN
397                Sloc = dSigmaH / ( -dSigmaDr(i,j) )                Sloc = dSigmaH / ( -dSigmaR )
398                 ELSE
399                  Sloc = GM_Visbeck_maxSlope
400                 ENDIF
401                 M2loc = gravity*recip_rhoConst*dSigmaH
402    c            SNloc = SQRT( Sloc*M2loc )
403                 N2loc = -gravity*recip_rhoConst*dSigmaR
404    c            N2loc = -gravity*recip_rhoConst*dSigmaDr(i,j)
405                 IF ( N2loc.GT.0. _d 0 ) THEN
406                   SNloc = Sloc*SQRT(N2loc)
407               ELSE               ELSE
408                Sloc = GM_maxSlope                 SNloc = 0. _d 0
409               ENDIF               ENDIF
              M2loc = Gravity*recip_RhoConst*dSigmaH  
              SNloc = SQRT( Sloc*M2loc )  
410             ELSE             ELSE
411               SNloc = 0. _d 0               SNloc = 0. _d 0
412             ENDIF             ENDIF
# Line 388  C-       compute slope, limited by GM_ma Line 419  C-       compute slope, limited by GM_ma
419         ENDIF         ENDIF
420  #endif /* ndef OLD_VISBECK_CALC */  #endif /* ndef OLD_VISBECK_CALC */
421  #endif /* GM_VISBECK_VARIABLE_K */  #endif /* GM_VISBECK_VARIABLE_K */
422           DO j=1-Oly,sNy+Oly
423            DO i=1-Olx,sNx+Olx
424             dSigmaDr(i,j)=sigmaR(i,j,k)
425            ENDDO
426           ENDDO
427    
428    #ifdef ALLOW_AUTODIFF_TAMC
429    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
430    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
431    CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
432    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
433    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
434    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
435    #endif /* ALLOW_AUTODIFF_TAMC */
436    
437  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
438         CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
# Line 442  C       which is used in the "variable K Line 487  C       which is used in the "variable K
487  C       Distance between interface above layer and the integration depth  C       Distance between interface above layer and the integration depth
488          deltaH=abs(GM_Visbeck_depth)-abs(rF(k))          deltaH=abs(GM_Visbeck_depth)-abs(rF(k))
489  C       If positive we limit this to the layer thickness  C       If positive we limit this to the layer thickness
490          deltaH=min(deltaH,drF(k))          integrDepth = drF(k)
491            deltaH=min(deltaH,integrDepth)
492  C       If negative then we are below the integration level  C       If negative then we are below the integration level
493          deltaH=max(deltaH,zero_rs)          deltaH=max(deltaH, 0. _d 0)
494  C       Now we convert deltaH to a non-dimensional fraction  C       Now we convert deltaH to a non-dimensional fraction
495          deltaH=deltaH/GM_Visbeck_depth          deltaH=deltaH/GM_Visbeck_depth
496    
         IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.  
497          IF ( Ssq(i,j).NE.0. .AND. dSigmaDr(i,j).NE.0. ) THEN          IF ( Ssq(i,j).NE.0. .AND. dSigmaDr(i,j).NE.0. ) THEN
498           N2= -Gravity*recip_RhoConst*dSigmaDr(i,j)           N2loc = -gravity*recip_rhoConst*dSigmaDr(i,j)
499           SN=sqrt(Ssq(i,j)*N2)           SNloc = SQRT(Ssq(i,j)*N2loc )
500           VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH           VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
501       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN       &       +deltaH*GM_Visbeck_alpha
502         &              *GM_Visbeck_length*GM_Visbeck_length*SNloc
503          ENDIF          ENDIF
504    
505          ENDDO          ENDDO
# Line 464  C       Now we convert deltaH to a non-d Line 510  C       Now we convert deltaH to a non-d
510  C-- end 1rst loop on vertical level index k  C-- end 1rst loop on vertical level index k
511        ENDDO        ENDDO
512    
513    
514  #ifdef GM_SUBMESO  #ifdef GM_SUBMESO
515  CBFK-- Use the dsigmadr average to construct the coefficients of the SM param  CBFK-- Use the dsigmadr average to construct the coefficients of the SM param
516        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
# Line 493  CBFK Scaling only works up to 1 degree. Line 540  CBFK Scaling only works up to 1 degree.
540        ENDDO        ENDDO
541  #endif  #endif
542    
543    
544  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
545  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
546  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
# Line 502  C-     Limit range that KapGM can take Line 550  C-     Limit range that KapGM can take
550         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
551          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
552           VisbeckK(i,j,bi,bj)=           VisbeckK(i,j,bi,bj)=
553       &       MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)       &       MIN( MAX( VisbeckK(i,j,bi,bj), GM_Visbeck_minVal_K ),
554         &            GM_Visbeck_maxVal_K )
555          ENDDO          ENDDO
556         ENDDO         ENDDO
557        ENDIF        ENDIF
# Line 542  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bi Line 591  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bi
591           Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)           Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)
592           Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)           Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)
593  #ifdef ALLOW_KAPREDI_CONTROL  #ifdef ALLOW_KAPREDI_CONTROL
594           Kwz(i,j,k,bi,bj)= ( kapredi(i,j,k,bi,bj)           Kwz(i,j,k,bi,bj)= ( kapredi(i,j,k,bi,bj)
595  #else  #else
596           Kwz(i,j,k,bi,bj)= ( GM_isopycK           Kwz(i,j,k,bi,bj)= ( GM_isopycK
597  #endif  #endif
# Line 562  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bi Line 611  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bi
611        ENDIF        ENDIF
612  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
613    
614    
615  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
616  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
617  C-- 2nd  k loop : compute Tensor Coeff. at U point  C-- 2nd  k loop : compute Tensor Coeff. at U point
# Line 640  C--      compute: ( M^2 * S )^1/2   (= M Line 690  C--      compute: ( M^2 * S )^1/2   (= M
690            ENDIF            ENDIF
691           ENDIF           ENDIF
692  #endif  #endif
   
693          ENDDO          ENDDO
694         ENDDO         ENDDO
695    
# Line 776  C-- end 2nd  loop on vertical level inde Line 825  C-- end 2nd  loop on vertical level inde
825    
826  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
827  C-- 3rd  k loop : compute Tensor Coeff. at V point  C-- 3rd  k loop : compute Tensor Coeff. at V point
828    
829  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
830        IF ( useKPP ) THEN        IF ( useKPP ) THEN
831         DO j=2-Oly,sNy+Oly         DO j=2-Oly,sNy+Oly
# Line 851  C       Now we convert deltaH to a non-d Line 901  C       Now we convert deltaH to a non-d
901            ENDIF            ENDIF
902           ENDIF           ENDIF
903  #endif  #endif
   
904          ENDDO          ENDDO
905         ENDDO         ENDDO
906    
# Line 925  CADJ STORE taperFct(:,:)     = comlev1_b Line 974  CADJ STORE taperFct(:,:)     = comlev1_b
974  #ifdef ALLOW_KAPREDI_CONTROL  #ifdef ALLOW_KAPREDI_CONTROL
975       &     ( kapredi(i,j,k,bi,bj)       &     ( kapredi(i,j,k,bi,bj)
976  #else  #else
977       &     ( GM_isopycK       &     ( GM_isopycK
978  #endif  #endif
979  #ifdef ALLOW_KAPGM_CONTROL  #ifdef ALLOW_KAPGM_CONTROL
980       &     - GM_skewflx*kapgm(i,j,k,bi,bj)       &     - GM_skewflx*kapgm(i,j,k,bi,bj)
# Line 960  C         store in tmp1k Kvz_Redi Line 1009  C         store in tmp1k Kvz_Redi
1009          ENDDO          ENDDO
1010          DO j=1,sNy+1          DO j=1,sNy+1
1011           DO i=1,sNx           DO i=1,sNx
1012  C-        Vertical gradients interpolated to V points  C-        Vertical gradients interpolated to U points
1013            dTdz =   op5*(            dTdz = (
1014       &   +op5*recip_drC(k)*       &     +recip_drC(k)*
1015       &       ( maskC(i,j-1,k,bi,bj)*       &       ( maskC(i,j-1,k,bi,bj)*
1016       &           (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))
1017       &        +maskC(i, j ,k,bi,bj)*       &        +maskC(i, j ,k,bi,bj)*
1018       &           (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))       &           (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))
1019       &       )       &       )
1020       &   +op5*recip_drC(kp1)*       &     +recip_drC(kp1)*
1021       &       ( maskC(i,j-1,kp1,bi,bj)*       &       ( maskC(i,j-1,kp1,bi,bj)*
1022       &           (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))
1023       &        +maskC(i, j ,kp1,bi,bj)*       &        +maskC(i, j ,kp1,bi,bj)*
1024       &           (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))       &           (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))
1025       &       )      )       &       )      ) * 0.25 _d 0
1026             tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)             tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)
1027       &                * _hFacS(i,j,k,bi,bj)       &                * _hFacS(i,j,k,bi,bj)
1028       &                * tmp1k(i,j) * dTdz       &                * tmp1k(i,j) * dTdz
# Line 999  C-- end 3rd  loop on vertical level inde Line 1048  C-- end 3rd  loop on vertical level inde
1048        ENDIF        ENDIF
1049  #endif  #endif
1050    
   
1051  #ifdef GM_SUBMESO  #ifdef GM_SUBMESO
1052  CBFK Add the submesoscale contribution, in a 4th k loop  CBFK Add the submesoscale contribution, in a 4th k loop
1053        DO k=1,Nr        DO k=1,Nr
# Line 1009  CBFK Add the submesoscale contribution, Line 1057  CBFK Add the submesoscale contribution,
1057  CBFK Add in the mu vertical structure factor  CBFK Add in the mu vertical structure factor
1058          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
1059           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
 #ifdef ALLOW_KPP  
           hml=KPPhbl(i,j,bi,bj)  
 #else  
1060            hml=hMixLayer(i,j,bi,bj)            hml=hMixLayer(i,j,bi,bj)
 #endif  
1061            IF (hml.gt.0 _d 0) THEN            IF (hml.gt.0 _d 0) THEN
1062             recip_hml=1 _d 0/hml             recip_hml=1 _d 0/hml
1063            ELSE            ELSE
# Line 1041  C Psiy are on faces above V Line 1085  C Psiy are on faces above V
1085            SM_PsiY(i,j)=op5*(SM_Lf(i,j+1)+SM_Lf(i,j))*dBdyAV(i,j)            SM_PsiY(i,j)=op5*(SM_Lf(i,j+1)+SM_Lf(i,j))*dBdyAV(i,j)
1086       &                        *_maskS(i,j,k,bi,bj)       &                        *_maskS(i,j,k,bi,bj)
1087    
1088    c hzhang: clipping here, ref to Baylor paper 3 appendix A MOM
1089            hml=hMixLayer(i,j,bi,bj)
1090            IF (hml .lt. (delR(1)+delR(2)+delR(3)+delR(4))) THEN
1091                    SM_PsiX(i,j)=0 _d 0
1092                    SM_PsiY(i,j)=0 _d 0
1093            ENDIF
1094    
1095            hml=.5 * delR(k)
1096            IF ( abs(SM_PsiX(i,j) .gt. hml) ) THEN
1097                    SM_PsiX(i,j)=SIGN( hml, SM_PsiX(i,j) )
1098            ENDIF
1099            IF ( abs(SM_PsiY(i,j) .gt. hml) ) THEN
1100                    SM_PsiY(i,j)=SIGN( hml, SM_PsiY(i,j) )
1101            ENDIF
1102    c hzhang: clipping done
1103    
1104    
1105  #ifndef GM_BOLUS_ADVEC  #ifndef GM_BOLUS_ADVEC
1106  C Kwx,Kwy are on rVel Points  C Kwx,Kwy are on rVel Points
1107            Kwx(i,j,k,bi,bj)  = Kwx(i,j,k,bi,bj)            Kwx(i,j,k,bi,bj)  = Kwx(i,j,k,bi,bj)
# Line 1276  C--   Time-average Line 1337  C--   Time-average
1337       &                          deltaTclock, bi, bj, myThid )       &                          deltaTclock, bi, bj, myThid )
1338         ENDIF         ENDIF
1339  #endif  #endif
1340         DO k=1,Nr         GM_timeAve(bi,bj) = GM_timeAve(bi,bj)+deltaTclock
          GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock  
        ENDDO  
1341    
1342        ENDIF        ENDIF
1343  #endif /* ALLOW_TIMEAVE */  #endif /* ALLOW_TIMEAVE */
# Line 1294  C--   Time-average Line 1353  C--   Time-average
1353        RETURN        RETURN
1354        END        END
1355    
1356    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1357    
1358        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
1359       I             iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22