/[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.5 by zhc, Fri Mar 12 18:31:00 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    #ifndef GM_BOLUS_ADVEC
1089    C Kwx,Kwy are on rVel Points
1090              Kwx(i,j,k,bi,bj)  = Kwx(i,j,k,bi,bj)
1091         &                      +op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))
1092              Kwy(i,j,k,bi,bj)  = Kwy(i,j,k,bi,bj)
1093         &                      +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))
1094    #ifdef GM_EXTRA_DIAGONAL
1095              IF (GM_ExtraDiag) THEN
1096    C Kuz,Kvz are on u,v Points
1097               Kuz(i,j,k,bi,bj)  = Kuz(i,j,k,bi,bj)
1098         &                    -op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1099               Kvz(i,j,k,bi,bj)  = Kvz(i,j,k,bi,bj)
1100         &                    -op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1101              ENDIF
1102    #endif
1103    #else
1104              IF (GM_AdvForm) THEN
1105               GM_PsiX(i,j,k,bi,bj)=GM_PsiX(i,j,k,bi,bj)+SM_PsiX(i,j)
1106               GM_PsiY(i,j,k,bi,bj)=GM_PsiY(i,j,k,bi,bj)+SM_PsiY(i,j)
1107              ENDIF
1108    #endif
1109             ENDDO
1110            ENDDO
1111           ELSE
1112            DO j=1-Oly+1,sNy+Oly-1
1113             DO i=1-Olx+1,sNx+Olx-1
1114              SM_PsiX(i,j)=0. _d 0
1115              SM_PsiY(i,j)=0. _d 0
1116             ENDDO
1117            ENDDO
1118           ENDIF  
1119    
1120    #ifdef ALLOW_DIAGNOSTICS
1121           IF ( useDiagnostics ) THEN
1122            IF ( DIAGNOSTICS_IS_ON('SM_PsiX ',myThid) ) THEN
1123             CALL DIAGNOSTICS_FILL(SM_PsiX,'SM_PsiX ',k,1,2,bi,bj,myThid)
1124            ENDIF
1125            IF ( DIAGNOSTICS_IS_ON('SM_PsiY ',myThid) ) THEN
1126             CALL DIAGNOSTICS_FILL(SM_PsiY,'SM_PsiY ',k,1,2,bi,bj,myThid)
1127            ENDIF
1128    
1129    CBFK Note:  for comparision, you can diagnose the bolus form
1130    CBFK or the Kappa form in the same simulation, regardless of other
1131    CBFK settings
1132            IF ( DIAGNOSTICS_IS_ON('SM_ubT  ',myThid) ) THEN
1133             DO j=jMin,jMax
1134              DO i=iMin,iMax
1135               tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiX(i,j)
1136         &                                -SM_PsiXm1(i,j) )
1137         &                               *maskW(i,j,km1,bi,bj)
1138         &            *op5*(Theta(i,j,km1,bi,bj)+Theta(i-1,j,km1,bi,bj))
1139              ENDDO
1140             ENDDO
1141             CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT  ', km1,1,2,bi,bj,myThid)
1142            ENDIF
1143    
1144            IF ( DIAGNOSTICS_IS_ON('SM_vbT  ',myThid) ) THEN
1145             DO j=jMin,jMax
1146              DO i=iMin,iMax
1147               tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiY(i,j)
1148         &                                -SM_PsiYm1(i,j) )
1149         &                               *maskS(i,j,km1,bi,bj)
1150         &            *op5*(Theta(i,j,km1,bi,bj)+Theta(i,j-1,km1,bi,bj))
1151              ENDDO
1152             ENDDO
1153             CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT  ', km1,1,2,bi,bj,myThid)
1154            ENDIF
1155    
1156            IF ( DIAGNOSTICS_IS_ON('SM_wbT  ',myThid) ) THEN
1157             DO j=jMin,jMax
1158              DO i=iMin,iMax
1159               tmp1k(i,j) =
1160         &        (dyG(i+1,j,bi,bj)*SM_PsiX(i+1,j)
1161         &        -dyG( i ,j,bi,bj)*SM_PsiX( i ,j)
1162         &        +dxG(i,j+1,bi,bj)*SM_PsiY(i,j+1)
1163         &        -dxG(i, j ,bi,bj)*SM_PsiY(i, j ))
1164         &             *op5*(Theta(i,j,k,bi,bj)+Theta(i,j,km1,bi,bj))
1165              ENDDO
1166             ENDDO
1167             CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT  ', k,1,2,bi,bj,myThid)
1168    C         print *,'SM_wbT',k,tmp1k
1169            ENDIF
1170    
1171            IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1172             DO j=1,sNy
1173              DO i=1,sNx+1
1174    C-        Vertical gradients interpolated to U points
1175               dTdz = (
1176         &      +recip_drC(k)*
1177         &       ( maskC(i-1,j,k,bi,bj)*
1178         &           (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
1179         &        +maskC( i ,j,k,bi,bj)*
1180         &           (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
1181         &       )
1182         &      +recip_drC(kp1)*
1183         &       ( maskC(i-1,j,kp1,bi,bj)*
1184         &           (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
1185         &        +maskC( i ,j,kp1,bi,bj)*
1186         &           (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
1187         &       )      ) * 0.25 _d 0
1188                tmp1k(i,j) = - dyG(i,j,bi,bj)*drF(k)
1189         &                 * _hFacW(i,j,k,bi,bj)
1190         &                 *op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1191         &                 * dTdz
1192              ENDDO
1193             ENDDO
1194             CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KuzTz', k,1,2,bi,bj,myThid)
1195    C         print *,'SM_KuzTz',k,tmp1k
1196            ENDIF
1197    
1198            IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1199             DO j=1,sNy+1
1200              DO i=1,sNx
1201    C-        Vertical gradients interpolated to V points
1202               dTdz =   op5*(
1203         &      +op5*recip_drC(k)*
1204         &       ( maskC(i,j-1,k,bi,bj)*
1205         &           (Theta(i,j-1,km1,bi,bj)-Theta(i,j-1,k,bi,bj))
1206         &        +maskC(i, j ,k,bi,bj)*
1207         &           (Theta(i, j ,km1,bi,bj)-Theta(i, j ,k,bi,bj))
1208         &       )
1209         &      +op5*recip_drC(kp1)*
1210         &       ( maskC(i,j-1,kp1,bi,bj)*
1211         &           (Theta(i,j-1,k,bi,bj)-Theta(i,j-1,kp1,bi,bj))
1212         &        +maskC(i, j ,kp1,bi,bj)*
1213         &           (Theta(i, j ,k,bi,bj)-Theta(i, j ,kp1,bi,bj))
1214         &       )      )
1215               tmp1k(i,j) = - dxG(i,j,bi,bj)*drF(k)
1216         &                * _hFacS(i,j,k,bi,bj)
1217         &                *op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1218         &                * dTdz
1219              ENDDO
1220             ENDDO
1221             CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KvzTz', k,1,2,bi,bj,myThid)
1222    C         print *,'SM_KvzTz',k,tmp1k
1223            ENDIF
1224    
1225            IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1226             DO j=jMin,jMax
1227              DO i=iMin,iMax
1228    C-      Horizontal gradients interpolated to W points
1229               dTdx = op5*(
1230         &           +op5*(_maskW(i+1,j,k,bi,bj)
1231         &            *_recip_dxC(i+1,j,bi,bj)*
1232         &             (Theta(i+1,j,k,bi,bj)-Theta(i,j,k,bi,bj))
1233         &              +_maskW(i,j,k,bi,bj)
1234         &            *_recip_dxC(i,j,bi,bj)*
1235         &             (Theta(i,j,k,bi,bj)-Theta(i-1,j,k,bi,bj)))
1236         &           +op5*(_maskW(i+1,j,k-1,bi,bj)
1237         &            *_recip_dxC(i+1,j,bi,bj)*
1238         &             (Theta(i+1,j,k-1,bi,bj)-Theta(i,j,k-1,bi,bj))
1239         &              +_maskW(i,j,k-1,bi,bj)
1240         &            *_recip_dxC(i,j,bi,bj)*
1241         &             (Theta(i,j,k-1,bi,bj)-Theta(i-1,j,k-1,bi,bj)))
1242         &       )
1243    
1244               dTdy = op5*(
1245         &           +op5*(_maskS(i,j,k,bi,bj)
1246         &            *_recip_dyC(i,j,bi,bj)*
1247         &             (Theta(i,j,k,bi,bj)-Theta(i,j-1,k,bi,bj))
1248         &              +_maskS(i,j+1,k,bi,bj)
1249         &            *_recip_dyC(i,j+1,bi,bj)*
1250         &             (Theta(i,j+1,k,bi,bj)-Theta(i,j,k,bi,bj)))
1251         &           +op5*(_maskS(i,j,k-1,bi,bj)
1252         &            *_recip_dyC(i,j,bi,bj)*
1253         &             (Theta(i,j,k-1,bi,bj)-Theta(i,j-1,k-1,bi,bj))
1254         &              +_maskS(i,j+1,k-1,bi,bj)
1255         &            *_recip_dyC(i,j+1,bi,bj)*
1256         &             (Theta(i,j+1,k-1,bi,bj)-Theta(i,j,k-1,bi,bj)))
1257         &       )
1258    
1259               tmp1k(i,j) = - _rA(i,j,bi,bj)
1260         &                     *(op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))*dTdx
1261         &                       +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))*dTdy)
1262              ENDDO
1263             ENDDO
1264             CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', k,1,2,bi,bj,myThid)
1265    C         print *,'SM_KrddT',k,tmp1k
1266            ENDIF
1267           ENDIF
1268    #endif
1269           DO j=1-Oly+1,sNy+Oly-1
1270            DO i=1-Olx+1,sNx+Olx-1
1271             SM_PsiXm1(i,j)=SM_PsiX(i,j)
1272             SM_PsiYm1(i,j)=SM_PsiY(i,j)
1273             tmp1k(i,j)=0 _d 0
1274            ENDDO
1275           ENDDO
1276          ENDDO
1277    
1278    CBFK Always Zero at the bottom.
1279          IF ( DIAGNOSTICS_IS_ON('SM_ubT  ',myThid) ) THEN
1280           CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT  ', Nr,1,2,bi,bj,myThid)
1281          ENDIF
1282          IF ( DIAGNOSTICS_IS_ON('SM_vbT  ',myThid) ) THEN
1283           CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT  ', Nr,1,2,bi,bj,myThid)
1284          ENDIF
1285          IF ( DIAGNOSTICS_IS_ON('SM_wbT  ',myThid) ) THEN
1286           CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT  ', Nr,1,2,bi,bj,myThid)
1287          ENDIF
1288          IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1289           CALL DIAGNOSTICS_FILL(tmp1k,'SM_KuzTz', Nr,1,2,bi,bj,myThid)
1290          ENDIF
1291          IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1292           CALL DIAGNOSTICS_FILL(tmp1k,'SM_KvzTz', Nr,1,2,bi,bj,myThid)
1293          ENDIF
1294          IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1295           CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', Nr,1,2,bi,bj,myThid)
1296          ENDIF
1297    #endif
1298    
1299  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
1300  C--   Time-average  C--   Time-average
1301        IF ( taveFreq.GT.0. ) THEN        IF ( taveFreq.GT.0. ) THEN
# Line 893  C--   Time-average Line 1320  C--   Time-average
1320       &                          deltaTclock, bi, bj, myThid )       &                          deltaTclock, bi, bj, myThid )
1321         ENDIF         ENDIF
1322  #endif  #endif
1323         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  
1324    
1325        ENDIF        ENDIF
1326  #endif /* ALLOW_TIMEAVE */  #endif /* ALLOW_TIMEAVE */
# Line 911  C--   Time-average Line 1336  C--   Time-average
1336        RETURN        RETURN
1337        END        END
1338    
1339    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1340    
1341        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
1342       I             iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22