/[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.4 by dimitri, Fri May 30 23:18:08 2008 UTC revision 1.7 by zhc, Fri Mar 19 19:29:47 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 531  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bi Line 580  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bi
580  #else  #else
581           Kgm_tmp = GM_isopycK           Kgm_tmp = GM_isopycK
582  #endif  #endif
583  #if (defined (ALLOW_AUTODIFF) && defined (ALLOW_KAPGM_CONTROL))  #ifdef ALLOW_KAPGM_CONTROL
584       &           + GM_skewflx*kapgm(i,j,k,bi,bj)       &           + GM_skewflx*kapgm(i,j,k,bi,bj)
585  #else  #else
586       &           + GM_skewflx*GM_background_K       &           + GM_skewflx*GM_background_K
# 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 715  CADJ STORE taperFct(:,:)     = comlev1_b Line 764  CADJ STORE taperFct(:,:)     = comlev1_b
764  #else  #else
765       &     ( GM_isopycK       &     ( GM_isopycK
766  #endif  #endif
767  #if (defined (ALLOW_AUTODIFF) && defined (ALLOW_KAPGM_CONTROL))  #ifdef ALLOW_KAPGM_CONTROL
768       &     - GM_skewflx*kapgm(i,j,k,bi,bj)       &     - GM_skewflx*kapgm(i,j,k,bi,bj)
769  #else  #else
770       &     - GM_skewflx*GM_background_K       &     - GM_skewflx*GM_background_K
# Line 852  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 926  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  #if (defined (ALLOW_AUTODIFF) && defined (ALLOW_KAPGM_CONTROL))  #ifdef ALLOW_KAPGM_CONTROL
980       &     - GM_skewflx*kapgm(i,j,k,bi,bj)       &     - GM_skewflx*kapgm(i,j,k,bi,bj)
981  #else  #else
982       &     - GM_skewflx*GM_background_K       &     - GM_skewflx*GM_background_K
# Line 961  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 1000  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 1038  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 1273  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 1291  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.4  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22