/[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.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
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 84  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
92          _RL dBdxAV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
93          _RL dBdyAV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
94          _RL SM_Lf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
95          _RL SM_PsiX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
96          _RL SM_PsiY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
97          _RL SM_PsiXm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
98          _RL SM_PsiYm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
99          _RL hsqmu, hml, recip_hml, qfac, dS, mlmax
100    #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 deltaH,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 deltaH, integrDepth        _RL Sloc, M2loc
       _RL Sloc, M2loc, SNloc  
 #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        _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 132  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 210  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    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    C SM(1)  
241             mlmax=max(mlmax,locMixLayer(i,j))
242          ENDDO          ENDDO
243         ENDDO         ENDDO
244        ELSE        ELSE
# Line 224  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    C SM(1)
252             mlmax=max(mlmax,locMixLayer(i,j))
253          ENDDO          ENDDO
254         ENDDO         ENDDO
255        ENDIF        ENDIF
256    
257    #ifdef GM_SUBMESO
258          DO j=1-Oly,sNy+Oly
259           DO i=1-Olx,sNx+Olx
260            dBdxAV(i,j)     = 0. _d 0
261            dBdyAV(i,j)     = 0. _d 0
262            SM_Lf(i,j)        = 0. _d 0
263            SM_PsiX(i,j)      = 0. _d 0
264            SM_PsiY(i,j)      = 0. _d 0
265            SM_PsiXm1(i,j)    = 0. _d 0
266            SM_PsiYm1(i,j)    = 0. _d 0
267           ENDDO
268          ENDDO
269    #endif
270                                                                                      
271    
272        DO k=Nr,2,-1        DO k=Nr,2,-1
273    
274  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 258  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 270  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
315    #ifdef GM_SUBMESO_VARYLf
316    C--     Depth average of SigmaR at W points
317    C       compute depth average from surface down to the MixLayer depth
318              IF (-rC(k-1).LT.locMixLayer(i,j) ) THEN
319               IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
320                integrDepth = -rC( k )
321    C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments
322                integrDepth = MIN( integrDepth, locMixLayer(i,j) )
323    C       Distance between level center above and the integration depth
324                deltaH = integrDepth + rC(k-1)
325    C       If negative then we are below the integration level
326    C       (cannot be the case with 2 conditions on maskC & -rC(k-1))
327    C       If positive we limit this to the distance from center above
328                deltaH = MIN( deltaH, drC(k) )
329    C       Now we convert deltaH to a non-dimensional fraction
330                deltaH = deltaH/( integrDepth+rC(1) )
331    C--     Store db/dr in SM_Lf for now.
332                SM_Lf(i,j) = SM_Lf(i,j)
333         &                   -gravity*recip_rhoConst*dSigmaDr(i,j)*deltaH
334               ENDIF
335              ENDIF
336    #endif
337    #endif
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 301  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 310  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               ELSE
399                Sloc = GM_maxSlope                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
408                   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 335  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 389  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 412  C-- end 1rst loop on vertical level inde Line 511  C-- end 1rst loop on vertical level inde
511        ENDDO        ENDDO
512    
513    
514    #ifdef GM_SUBMESO
515    CBFK-- Use the dsigmadr average to construct the coefficients of the SM param
516          DO j=1-Oly+1,sNy+Oly-1
517           DO i=1-Olx+1,sNx+Olx-1
518    #ifdef GM_SUBMESO_VARYLf
519    
520            IF (SM_Lf(i,j).gt.0) THEN
521    CBFK ML def. rad. as Lf if available and not too small
522              SM_Lf(i,j)=max(sqrt(SM_Lf(i,j))*locMixLayer(i,j)
523         &                        /abs(fCori(i,j,bi,bj))
524         &                   ,GM_SM_Lf)
525            ELSE
526    #else
527            IF (.TRUE.) THEN
528    #endif
529    CBFK Otherwise, store just the fixed number
530              SM_Lf(i,j)=GM_SM_Lf
531            ENDIF
532    CBFK Now do the rest of the coefficient
533            dS=2*dxC(i,j,bi,bj)*dyC(i,j,bi,bj)/
534         &              (dxC(i,j,bi,bj)+dyC(i,j,bi,bj))
535    CBFK Scaling only works up to 1 degree.
536            dS=min(dS,GM_SM_Lmax)
537            deltaH=sqrt(fCori(i,j,bi,bj)**2+1 _d 0/(GM_SM_tau**2))
538            SM_Lf(i,j)=GM_SM_Ce*dS/(deltaH*SM_Lf(i,j))
539           ENDDO
540          ENDDO
541    #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 421  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 461  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 535  C     Gradient of Sigma at U points Line 665  C     Gradient of Sigma at U points
665           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 )
666       &                      +(sigmaR(i-1,j,kp1)+sigmaR(i,j,kp1))*maskp1       &                      +(sigmaR(i-1,j,kp1)+sigmaR(i,j,kp1))*maskp1
667       &                      )*_maskW(i,j,k,bi,bj)       &                      )*_maskW(i,j,k,bi,bj)
668    
669    #ifdef GM_SUBMESO
670    C--     Depth average of SigmaX at U points
671    C       compute depth average from surface down to the MixLayer depth
672             IF (k.GT.1) THEN
673              IF (-rC(k-1).LT.locMixLayer(i,j) ) THEN
674               IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
675                integrDepth = -rC( k )
676    C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments
677                integrDepth = MIN( integrDepth, locMixLayer(i,j) )
678    C       Distance between level center above and the integration depth
679                deltaH = integrDepth + rC(k-1)
680    C       If negative then we are below the integration level
681    C       (cannot be the case with 2 conditions on maskC & -rC(k-1))
682    C       If positive we limit this to the distance from center above
683                deltaH = MIN( deltaH, drC(k) )
684    C       Now we convert deltaH to a non-dimensional fraction
685                deltaH = deltaH/( integrDepth+rC(1) )
686    C--      compute: ( M^2 * S )^1/2   (= M^2 / N since S=M^2/N^2 )
687                dBdxAV(i,j) = dBdxAV(i,j)
688         &            +dSigmaDx(i,j)*deltaH*recip_rhoConst*gravity
689               ENDIF
690              ENDIF
691             ENDIF
692    #endif
693          ENDDO          ENDDO
694         ENDDO         ENDDO
695    
# Line 722  C     Gradient of Sigma at V points Line 877  C     Gradient of Sigma at V points
877           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 )
878       &                      +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1       &                      +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1
879       &                      )*_maskS(i,j,k,bi,bj)       &                      )*_maskS(i,j,k,bi,bj)
880    
881    #ifdef GM_SUBMESO
882    C--     Depth average of SigmaY at V points
883    C       compute depth average from surface down to the MixLayer depth
884             IF (k.GT.1) THEN
885              IF (-rC(k-1).LT.locMixLayer(i,j) ) THEN
886               IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
887                integrDepth = -rC( k )
888    C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments
889                integrDepth = MIN( integrDepth, locMixLayer(i,j) )
890    C       Distance between level center above and the integration depth
891                deltaH = integrDepth + rC(k-1)
892    C       If negative then we are below the integration level
893    C       (cannot be the case with 2 conditions on maskC & -rC(k-1))
894    C       If positive we limit this to the distance from center above
895                deltaH = MIN( deltaH, drC(k) )
896    C       Now we convert deltaH to a non-dimensional fraction
897                deltaH = deltaH/( integrDepth+rC(1) )
898                dBdyAV(i,j) = dBdyAV(i,j)
899         &            +dSigmaDy(i,j)*deltaH*recip_rhoConst*gravity
900               ENDIF
901              ENDIF
902             ENDIF
903    #endif
904          ENDDO          ENDDO
905         ENDDO         ENDDO
906    
# Line 795  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 869  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
1052    CBFK Add the submesoscale contribution, in a 4th k loop
1053          DO k=1,Nr
1054           km1=max(1,k-1)
1055           IF ((k.gt.1).and.(-rF(k-1) .lt. mlmax)) THEN
1056            kp1 = MIN(k+1,Nr)
1057    CBFK Add in the mu vertical structure factor
1058            DO j=1-Oly+1,sNy+Oly-1
1059             DO i=1-Olx+1,sNx+Olx-1
1060              hml=hMixLayer(i,j,bi,bj)
1061              IF (hml.gt.0 _d 0) THEN
1062               recip_hml=1 _d 0/hml
1063              ELSE
1064               recip_hml=0 _d 0
1065              ENDIF
1066    CBFK We calculate the h^2 mu(z) factor only on w points.  
1067    CBFK It is possible that we might need to calculate it
1068    CBFK on Psi or u,v points independently to prevent spurious
1069    CBFK entrainment.  Unlikely that this will be major
1070    CBFK (it wasnt in offline testing).
1071              qfac=(2*rf(k)*recip_hml+1 _d 0)**2
1072              hsqmu=(1 _d 0-qfac)*(1 _d 0+(5 _d 0)*qfac/21 _d 0)
1073              hsqmu=max(0 _d 0, hsqmu)*hml**2
1074              SM_Lf(i,j)=SM_Lf(i,j)*hsqmu
1075             ENDDO
1076            ENDDO
1077    CBFK Now interpolate to match locations
1078            DO j=1-Oly+1,sNy+Oly-1
1079             DO i=1-Olx+1,sNx+Olx-1
1080    C SM_Lf coefficients are on rVel points
1081    C Psix are on faces above U
1082              SM_PsiX(i,j)=op5*(SM_Lf(i+1,j)+SM_Lf(i,j))*dBdxAV(i,j)
1083         &                        *_maskW(i,j,k,bi,bj)
1084    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)
1086         &                        *_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
1106    C Kwx,Kwy are on rVel Points
1107              Kwx(i,j,k,bi,bj)  = Kwx(i,j,k,bi,bj)
1108         &                      +op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))
1109              Kwy(i,j,k,bi,bj)  = Kwy(i,j,k,bi,bj)
1110         &                      +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))
1111    #ifdef GM_EXTRA_DIAGONAL
1112              IF (GM_ExtraDiag) THEN
1113    C Kuz,Kvz are on u,v Points
1114               Kuz(i,j,k,bi,bj)  = Kuz(i,j,k,bi,bj)
1115         &                    -op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1116               Kvz(i,j,k,bi,bj)  = Kvz(i,j,k,bi,bj)
1117         &                    -op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1118              ENDIF
1119    #endif
1120    #else
1121              IF (GM_AdvForm) THEN
1122               GM_PsiX(i,j,k,bi,bj)=GM_PsiX(i,j,k,bi,bj)+SM_PsiX(i,j)
1123               GM_PsiY(i,j,k,bi,bj)=GM_PsiY(i,j,k,bi,bj)+SM_PsiY(i,j)
1124              ENDIF
1125    #endif
1126             ENDDO
1127            ENDDO
1128           ELSE
1129            DO j=1-Oly+1,sNy+Oly-1
1130             DO i=1-Olx+1,sNx+Olx-1
1131              SM_PsiX(i,j)=0. _d 0
1132              SM_PsiY(i,j)=0. _d 0
1133             ENDDO
1134            ENDDO
1135           ENDIF  
1136    
1137    #ifdef ALLOW_DIAGNOSTICS
1138           IF ( useDiagnostics ) THEN
1139            IF ( DIAGNOSTICS_IS_ON('SM_PsiX ',myThid) ) THEN
1140             CALL DIAGNOSTICS_FILL(SM_PsiX,'SM_PsiX ',k,1,2,bi,bj,myThid)
1141            ENDIF
1142            IF ( DIAGNOSTICS_IS_ON('SM_PsiY ',myThid) ) THEN
1143             CALL DIAGNOSTICS_FILL(SM_PsiY,'SM_PsiY ',k,1,2,bi,bj,myThid)
1144            ENDIF
1145    
1146    CBFK Note:  for comparision, you can diagnose the bolus form
1147    CBFK or the Kappa form in the same simulation, regardless of other
1148    CBFK settings
1149            IF ( DIAGNOSTICS_IS_ON('SM_ubT  ',myThid) ) THEN
1150             DO j=jMin,jMax
1151              DO i=iMin,iMax
1152               tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiX(i,j)
1153         &                                -SM_PsiXm1(i,j) )
1154         &                               *maskW(i,j,km1,bi,bj)
1155         &            *op5*(Theta(i,j,km1,bi,bj)+Theta(i-1,j,km1,bi,bj))
1156              ENDDO
1157             ENDDO
1158             CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT  ', km1,1,2,bi,bj,myThid)
1159            ENDIF
1160    
1161            IF ( DIAGNOSTICS_IS_ON('SM_vbT  ',myThid) ) THEN
1162             DO j=jMin,jMax
1163              DO i=iMin,iMax
1164               tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiY(i,j)
1165         &                                -SM_PsiYm1(i,j) )
1166         &                               *maskS(i,j,km1,bi,bj)
1167         &            *op5*(Theta(i,j,km1,bi,bj)+Theta(i,j-1,km1,bi,bj))
1168              ENDDO
1169             ENDDO
1170             CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT  ', km1,1,2,bi,bj,myThid)
1171            ENDIF
1172    
1173            IF ( DIAGNOSTICS_IS_ON('SM_wbT  ',myThid) ) THEN
1174             DO j=jMin,jMax
1175              DO i=iMin,iMax
1176               tmp1k(i,j) =
1177         &        (dyG(i+1,j,bi,bj)*SM_PsiX(i+1,j)
1178         &        -dyG( i ,j,bi,bj)*SM_PsiX( i ,j)
1179         &        +dxG(i,j+1,bi,bj)*SM_PsiY(i,j+1)
1180         &        -dxG(i, j ,bi,bj)*SM_PsiY(i, j ))
1181         &             *op5*(Theta(i,j,k,bi,bj)+Theta(i,j,km1,bi,bj))
1182              ENDDO
1183             ENDDO
1184             CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT  ', k,1,2,bi,bj,myThid)
1185    C         print *,'SM_wbT',k,tmp1k
1186            ENDIF
1187    
1188            IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1189             DO j=1,sNy
1190              DO i=1,sNx+1
1191    C-        Vertical gradients interpolated to U points
1192               dTdz = (
1193         &      +recip_drC(k)*
1194         &       ( maskC(i-1,j,k,bi,bj)*
1195         &           (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
1196         &        +maskC( i ,j,k,bi,bj)*
1197         &           (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
1198         &       )
1199         &      +recip_drC(kp1)*
1200         &       ( maskC(i-1,j,kp1,bi,bj)*
1201         &           (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
1202         &        +maskC( i ,j,kp1,bi,bj)*
1203         &           (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
1204         &       )      ) * 0.25 _d 0
1205                tmp1k(i,j) = - dyG(i,j,bi,bj)*drF(k)
1206         &                 * _hFacW(i,j,k,bi,bj)
1207         &                 *op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1208         &                 * dTdz
1209              ENDDO
1210             ENDDO
1211             CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KuzTz', k,1,2,bi,bj,myThid)
1212    C         print *,'SM_KuzTz',k,tmp1k
1213            ENDIF
1214    
1215            IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1216             DO j=1,sNy+1
1217              DO i=1,sNx
1218    C-        Vertical gradients interpolated to V points
1219               dTdz =   op5*(
1220         &      +op5*recip_drC(k)*
1221         &       ( maskC(i,j-1,k,bi,bj)*
1222         &           (Theta(i,j-1,km1,bi,bj)-Theta(i,j-1,k,bi,bj))
1223         &        +maskC(i, j ,k,bi,bj)*
1224         &           (Theta(i, j ,km1,bi,bj)-Theta(i, j ,k,bi,bj))
1225         &       )
1226         &      +op5*recip_drC(kp1)*
1227         &       ( maskC(i,j-1,kp1,bi,bj)*
1228         &           (Theta(i,j-1,k,bi,bj)-Theta(i,j-1,kp1,bi,bj))
1229         &        +maskC(i, j ,kp1,bi,bj)*
1230         &           (Theta(i, j ,k,bi,bj)-Theta(i, j ,kp1,bi,bj))
1231         &       )      )
1232               tmp1k(i,j) = - dxG(i,j,bi,bj)*drF(k)
1233         &                * _hFacS(i,j,k,bi,bj)
1234         &                *op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1235         &                * dTdz
1236              ENDDO
1237             ENDDO
1238             CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KvzTz', k,1,2,bi,bj,myThid)
1239    C         print *,'SM_KvzTz',k,tmp1k
1240            ENDIF
1241    
1242            IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1243             DO j=jMin,jMax
1244              DO i=iMin,iMax
1245    C-      Horizontal gradients interpolated to W points
1246               dTdx = op5*(
1247         &           +op5*(_maskW(i+1,j,k,bi,bj)
1248         &            *_recip_dxC(i+1,j,bi,bj)*
1249         &             (Theta(i+1,j,k,bi,bj)-Theta(i,j,k,bi,bj))
1250         &              +_maskW(i,j,k,bi,bj)
1251         &            *_recip_dxC(i,j,bi,bj)*
1252         &             (Theta(i,j,k,bi,bj)-Theta(i-1,j,k,bi,bj)))
1253         &           +op5*(_maskW(i+1,j,k-1,bi,bj)
1254         &            *_recip_dxC(i+1,j,bi,bj)*
1255         &             (Theta(i+1,j,k-1,bi,bj)-Theta(i,j,k-1,bi,bj))
1256         &              +_maskW(i,j,k-1,bi,bj)
1257         &            *_recip_dxC(i,j,bi,bj)*
1258         &             (Theta(i,j,k-1,bi,bj)-Theta(i-1,j,k-1,bi,bj)))
1259         &       )
1260    
1261               dTdy = op5*(
1262         &           +op5*(_maskS(i,j,k,bi,bj)
1263         &            *_recip_dyC(i,j,bi,bj)*
1264         &             (Theta(i,j,k,bi,bj)-Theta(i,j-1,k,bi,bj))
1265         &              +_maskS(i,j+1,k,bi,bj)
1266         &            *_recip_dyC(i,j+1,bi,bj)*
1267         &             (Theta(i,j+1,k,bi,bj)-Theta(i,j,k,bi,bj)))
1268         &           +op5*(_maskS(i,j,k-1,bi,bj)
1269         &            *_recip_dyC(i,j,bi,bj)*
1270         &             (Theta(i,j,k-1,bi,bj)-Theta(i,j-1,k-1,bi,bj))
1271         &              +_maskS(i,j+1,k-1,bi,bj)
1272         &            *_recip_dyC(i,j+1,bi,bj)*
1273         &             (Theta(i,j+1,k-1,bi,bj)-Theta(i,j,k-1,bi,bj)))
1274         &       )
1275    
1276               tmp1k(i,j) = - _rA(i,j,bi,bj)
1277         &                     *(op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))*dTdx
1278         &                       +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))*dTdy)
1279              ENDDO
1280             ENDDO
1281             CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', k,1,2,bi,bj,myThid)
1282    C         print *,'SM_KrddT',k,tmp1k
1283            ENDIF
1284           ENDIF
1285    #endif
1286           DO j=1-Oly+1,sNy+Oly-1
1287            DO i=1-Olx+1,sNx+Olx-1
1288             SM_PsiXm1(i,j)=SM_PsiX(i,j)
1289             SM_PsiYm1(i,j)=SM_PsiY(i,j)
1290             tmp1k(i,j)=0 _d 0
1291            ENDDO
1292           ENDDO
1293          ENDDO
1294    
1295    CBFK Always Zero at the bottom.
1296          IF ( DIAGNOSTICS_IS_ON('SM_ubT  ',myThid) ) THEN
1297           CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT  ', Nr,1,2,bi,bj,myThid)
1298          ENDIF
1299          IF ( DIAGNOSTICS_IS_ON('SM_vbT  ',myThid) ) THEN
1300           CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT  ', Nr,1,2,bi,bj,myThid)
1301          ENDIF
1302          IF ( DIAGNOSTICS_IS_ON('SM_wbT  ',myThid) ) THEN
1303           CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT  ', Nr,1,2,bi,bj,myThid)
1304          ENDIF
1305          IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1306           CALL DIAGNOSTICS_FILL(tmp1k,'SM_KuzTz', Nr,1,2,bi,bj,myThid)
1307          ENDIF
1308          IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1309           CALL DIAGNOSTICS_FILL(tmp1k,'SM_KvzTz', Nr,1,2,bi,bj,myThid)
1310          ENDIF
1311          IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1312           CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', Nr,1,2,bi,bj,myThid)
1313          ENDIF
1314    #endif
1315    
1316  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
1317  C--   Time-average  C--   Time-average
1318        IF ( taveFreq.GT.0. ) THEN        IF ( taveFreq.GT.0. ) THEN
# Line 893  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 911  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.1  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22