/[MITgcm]/MITgcm/pkg/gmredi/gmredi_calc_tensor.F
ViewVC logotype

Diff of /MITgcm/pkg/gmredi/gmredi_calc_tensor.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.26 by jmc, Thu May 24 22:34:38 2007 UTC revision 1.37 by jmc, Tue Jan 11 00:54:45 2011 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "GMREDI_OPTIONS.h"  #include "GMREDI_OPTIONS.h"
5  #undef OLD_VISBECK_CALC  #ifdef ALLOW_KPP
6    # include "KPP_OPTIONS.h"
7    #endif
8    
9  CStartOfInterface  CBOP
10    C     !ROUTINE: GMREDI_CALC_TENSOR
11    C     !INTERFACE:
12        SUBROUTINE GMREDI_CALC_TENSOR(        SUBROUTINE GMREDI_CALC_TENSOR(
13       I             bi, bj, iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
14       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
15       I             myThid )       I             bi, bj, myTime, myIter, myThid )
16  C     /==========================================================\  
17  C     | SUBROUTINE GMREDI_CALC_TENSOR                            |  C     !DESCRIPTION: \bv
18  C     | o Calculate tensor elements for GM/Redi tensor.          |  C     *==========================================================*
19  C     |==========================================================|  C     | SUBROUTINE GMREDI_CALC_TENSOR
20  C     \==========================================================/  C     | o Calculate tensor elements for GM/Redi tensor.
21    C     *==========================================================*
22    C     *==========================================================*
23    C     \ev
24    
25    C     !USES:
26        IMPLICIT NONE        IMPLICIT NONE
27    
28  C     == Global variables ==  C     == Global variables ==
# Line 24  C     == Global variables == Line 33  C     == Global variables ==
33  #include "PARAMS.h"  #include "PARAMS.h"
34  #include "GMREDI.h"  #include "GMREDI.h"
35  #include "GMREDI_TAVE.h"  #include "GMREDI_TAVE.h"
36    #ifdef ALLOW_KPP
37    # include "KPP.h"
38    #endif
39    
40  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
41  #include "tamc.h"  #include "tamc.h"
42  #include "tamc_keys.h"  #include "tamc_keys.h"
43  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
44    
45    C     !INPUT/OUTPUT PARAMETERS:
46  C     == Routine arguments ==  C     == Routine arguments ==
47    C     bi, bj    :: tile indices
48    C     myTime    :: Current time in simulation
49    C     myIter    :: Current iteration number in simulation
50    C     myThid    :: My Thread Id. number
51  C  C
52          INTEGER iMin,iMax,jMin,jMax
53        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
54        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
55        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
56        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi, bj
57          _RL     myTime
58          INTEGER myIter
59        INTEGER myThid        INTEGER myThid
60  CEndOfInterface  CEOP
61    
62  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
63    
64    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 50  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)
       _RL maskp1, Kgm_tmp  
74        _RL ldd97_LrhoC(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL ldd97_LrhoC(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
75        _RL ldd97_LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL ldd97_LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
76        _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
77        _RL Cspd, LrhoInf, LrhoSup, fCoriLoc        _RL Cspd, LrhoInf, LrhoSup, fCoriLoc
78          _RL Kgm_tmp, isopycK, bolus_K
79    
80          INTEGER kLow_W (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
81          INTEGER kLow_S (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
82          _RL locMixLayer(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
83          _RL baseSlope  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
84          _RL hTransLay  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
85          _RL recipLambda(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
86          INTEGER  km1
87    #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
88          INTEGER kp1
89          _RL maskp1
90    #endif
91    
92  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
93  #ifdef OLD_VISBECK_CALC  #ifdef OLD_VISBECK_CALC
       _RL deltaH,zero_rs  
       PARAMETER(zero_rs=0.D0)  
       _RL N2,SN  
94        _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
95  #else  #else
96        _RL dSigmaH        _RL dSigmaH, dSigmaR
97        _RL deltaH, integrDepth        _RL Sloc, M2loc
       _RL Sloc, M2loc, SNloc  
 #endif  
98  #endif  #endif
99          _RL recipMaxSlope
100          _RL deltaH, integrDepth
101          _RL N2loc, SNloc
102    #endif /* GM_VISBECK_VARIABLE_K */
103    
104  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
105        LOGICAL  doDiagRediFlx        LOGICAL  doDiagRediFlx
106        LOGICAL  DIAGNOSTICS_IS_ON        LOGICAL  DIAGNOSTICS_IS_ON
107        EXTERNAL DIAGNOSTICS_IS_ON        EXTERNAL DIAGNOSTICS_IS_ON
108        INTEGER  km1  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
109        _RL dTdz        _RL dTdz
110        _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111  #endif  #endif
112    #endif /* ALLOW_DIAGNOSTICS */
113    
114  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
115    
# Line 103  C---+----1----+----2----+----3----+----4 Line 136  C---+----1----+----2----+----3----+----4
136  #endif  #endif
137    
138  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
139          recipMaxSlope = 0. _d 0
140          IF ( GM_Visbeck_maxSlope.GT.0. _d 0 ) THEN
141            recipMaxSlope = 1. _d 0 / GM_Visbeck_maxSlope
142          ENDIF
143        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
144         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
145          VisbeckK(i,j,bi,bj) = 0. _d 0          VisbeckK(i,j,bi,bj) = 0. _d 0
# Line 111  C---+----1----+----2----+----3----+----4 Line 148  C---+----1----+----2----+----3----+----4
148  #endif  #endif
149    
150  C--   set ldd97_Lrho (for tapering scheme ldd97):  C--   set ldd97_Lrho (for tapering scheme ldd97):
151        IF (GM_taper_scheme.EQ.'ldd97') THEN        IF ( GM_taper_scheme.EQ.'ldd97' .OR.
152         &     GM_taper_scheme.EQ.'fm07' ) THEN
153         Cspd = 2. _d 0         Cspd = 2. _d 0
154         LrhoInf = 15. _d 3         LrhoInf = 15. _d 3
155         LrhoSup = 100. _d 3         LrhoSup = 100. _d 3
# Line 128  C-     Tracer point location (center): Line 166  C-     Tracer point location (center):
166         ENDDO         ENDDO
167  C-     U point location (West):  C-     U point location (West):
168         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
169            kLow_W(1-Olx,j) = 0
170          ldd97_LrhoW(1-Olx,j) = LrhoSup          ldd97_LrhoW(1-Olx,j) = LrhoSup
171          DO i=1-Olx+1,sNx+Olx          DO i=1-Olx+1,sNx+Olx
172             kLow_W(i,j) = MIN(kLowC(i-1,j,bi,bj),kLowC(i,j,bi,bj))
173           fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))           fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))
174           IF (fCoriLoc.NE.0.) THEN           IF (fCoriLoc.NE.0.) THEN
175             ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)             ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)
# Line 141  C-     U point location (West): Line 181  C-     U point location (West):
181         ENDDO         ENDDO
182  C-     V point location (South):  C-     V point location (South):
183         DO i=1-Olx+1,sNx+Olx         DO i=1-Olx+1,sNx+Olx
184             kLow_S(i,1-Oly) = 0
185           ldd97_LrhoS(i,1-Oly) = LrhoSup           ldd97_LrhoS(i,1-Oly) = LrhoSup
186         ENDDO         ENDDO
187         DO j=1-Oly+1,sNy+Oly         DO j=1-Oly+1,sNy+Oly
188          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
189             kLow_S(i,j) = MIN(kLowC(i,j-1,bi,bj),kLowC(i,j,bi,bj))
190           fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))           fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))
191           IF (fCoriLoc.NE.0.) THEN           IF (fCoriLoc.NE.0.) THEN
192             ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)             ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)
# Line 164  C-    Just initialize to zero (not use a Line 206  C-    Just initialize to zero (not use a
206          ENDDO          ENDDO
207         ENDDO         ENDDO
208        ENDIF        ENDIF
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
209    
210        DO k=2,Nr  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
211  C-- 1rst loop on k : compute Tensor Coeff. at W points.  C-- 1rst loop on k : compute Tensor Coeff. at W points.
212    
213          DO j=1-Oly,sNy+Oly
214           DO i=1-Olx,sNx+Olx
215             hTransLay(i,j) = R_low(i,j,bi,bj)
216             baseSlope(i,j) =  0. _d 0
217             recipLambda(i,j) = 0. _d 0
218             locMixLayer(i,j) = 0. _d 0
219           ENDDO
220          ENDDO
221    #ifdef ALLOW_KPP
222          IF ( useKPP ) THEN
223           DO j=1-Oly,sNy+Oly
224            DO i=1-Olx,sNx+Olx
225             locMixLayer(i,j) = KPPhbl(i,j,bi,bj)
226            ENDDO
227           ENDDO
228          ELSE
229    #else
230          IF ( .TRUE. ) THEN
231    #endif
232           DO j=1-Oly,sNy+Oly
233            DO i=1-Olx,sNx+Olx
234             locMixLayer(i,j) = hMixLayer(i,j,bi,bj)
235            ENDDO
236           ENDDO
237          ENDIF
238    
239          DO k=Nr,2,-1
240    
241  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
242         kkey = (igmkey-1)*Nr + k         kkey = (igmkey-1)*Nr + k
243         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
# Line 197  C-- 1rst loop on k : compute Tensor Coef Line 266  C-- 1rst loop on k : compute Tensor Coef
266  # endif  # endif
267          ENDDO          ENDDO
268         ENDDO         ENDDO
269  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
270    
271         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
272          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
# Line 208  C      Gradient of Sigma at rVel points Line 277  C      Gradient of Sigma at rVel points
277           dSigmaDy(i,j)=op25*( sigmaY(i,j+1,k-1)+sigmaY(i,j,k-1)           dSigmaDy(i,j)=op25*( sigmaY(i,j+1,k-1)+sigmaY(i,j,k-1)
278       &                       +sigmaY(i,j+1, k )+sigmaY(i,j, k )       &                       +sigmaY(i,j+1, k )+sigmaY(i,j, k )
279       &                      )*maskC(i,j,k,bi,bj)       &                      )*maskC(i,j,k,bi,bj)
280           dSigmaDr(i,j)=sigmaR(i,j,k)  c        dSigmaDr(i,j)=sigmaR(i,j,k)
281          ENDDO          ENDDO
282         ENDDO         ENDDO
283    
 #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  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
284  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
285  #ifndef OLD_VISBECK_CALC  #ifndef OLD_VISBECK_CALC
286         IF ( GM_Visbeck_alpha.GT.0. .AND.         IF ( GM_Visbeck_alpha.GT.0. .AND.
287       &      -rC(k-1).LT.GM_Visbeck_depth ) THEN       &      -rC(k-1).LT.GM_Visbeck_depth ) THEN
288    
289            DO j=1-Oly,sNy+Oly
290             DO i=1-Olx,sNx+Olx
291              dSigmaDr(i,j) = MIN( sigmaR(i,j,k), 0. _d 0 )
292             ENDDO
293            ENDDO
294    
295  C--     Depth average of f/sqrt(Ri) = M^2/N^2 * N  C--     Depth average of f/sqrt(Ri) = M^2/N^2 * N
296  C       M^2 and N^2 are horizontal & vertical gradient of buoyancy.  C       M^2 and N^2 are horizontal & vertical gradient of buoyancy.
297    
# Line 237  C       GM_Visbeck_depth, whatever is th Line 306  C       GM_Visbeck_depth, whatever is th
306             integrDepth = -rC( kLowC(i,j,bi,bj) )             integrDepth = -rC( kLowC(i,j,bi,bj) )
307  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
308             integrDepth = MIN( integrDepth, GM_Visbeck_depth )             integrDepth = MIN( integrDepth, GM_Visbeck_depth )
309    C-      to recover "old-visbeck" form with Visbeck_minDepth = Visbeck_depth
310               integrDepth = MAX( integrDepth, GM_Visbeck_minDepth )
311  C       Distance between level center above and the integration depth  C       Distance between level center above and the integration depth
312             deltaH = integrDepth + rC(k-1)             deltaH = integrDepth + rC(k-1)
313  C       If negative then we are below the integration level  C       If negative then we are below the integration level
# Line 246  C       If positive we limit this to the Line 317  C       If positive we limit this to the
317  C       Now we convert deltaH to a non-dimensional fraction  C       Now we convert deltaH to a non-dimensional fraction
318             deltaH = deltaH/( integrDepth+rC(1) )             deltaH = deltaH/( integrDepth+rC(1) )
319    
320  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 )
321    C        a 5 points average gives a more "homogeneous" formulation
322    C        (same stencil and same weights as for dSigmaH calculation)
323               dSigmaR = ( dSigmaDr(i,j)*4. _d 0
324         &               + dSigmaDr(i-1,j)
325         &               + dSigmaDr(i+1,j)
326         &               + dSigmaDr(i,j-1)
327         &               + dSigmaDr(i,j+1)
328         &               )/( 4. _d 0
329         &                 + maskC(i-1,j,k,bi,bj)
330         &                 + maskC(i+1,j,k,bi,bj)
331         &                 + maskC(i,j-1,k,bi,bj)
332         &                 + maskC(i,j+1,k,bi,bj)
333         &                 )
334             dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)             dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)
335       &             + dSigmaDy(i,j)*dSigmaDy(i,j)       &             + dSigmaDy(i,j)*dSigmaDy(i,j)
336             IF ( dSigmaH .GT. 0. _d 0 ) THEN             IF ( dSigmaH .GT. 0. _d 0 ) THEN
337               dSigmaH = SQRT( dSigmaH )               dSigmaH = SQRT( dSigmaH )
338  C-       compute slope, limited by GM_maxSlope:  C-       compute slope, limited by GM_Visbeck_maxSlope:
339               IF ( -dSigmaDr(i,j).GT.dSigmaH*GM_rMaxSlope ) THEN               IF ( -dSigmaR.GT.dSigmaH*recipMaxSlope ) THEN
340                Sloc = dSigmaH / ( -dSigmaDr(i,j) )                Sloc = dSigmaH / ( -dSigmaR )
341               ELSE               ELSE
342                Sloc = GM_maxSlope                Sloc = GM_Visbeck_maxSlope
343                 ENDIF
344                 M2loc = gravity*recip_rhoConst*dSigmaH
345    c            SNloc = SQRT( Sloc*M2loc )
346                 N2loc = -gravity*recip_rhoConst*dSigmaR
347    c            N2loc = -gravity*recip_rhoConst*dSigmaDr(i,j)
348                 IF ( N2loc.GT.0. _d 0 ) THEN
349                   SNloc = Sloc*SQRT(N2loc)
350                 ELSE
351                   SNloc = 0. _d 0
352               ENDIF               ENDIF
              M2loc = Gravity*recip_RhoConst*dSigmaH  
              SNloc = SQRT( Sloc*M2loc )  
353             ELSE             ELSE
354               SNloc = 0. _d 0               SNloc = 0. _d 0
355             ENDIF             ENDIF
# Line 271  C-       compute slope, limited by GM_ma Line 362  C-       compute slope, limited by GM_ma
362         ENDIF         ENDIF
363  #endif /* ndef OLD_VISBECK_CALC */  #endif /* ndef OLD_VISBECK_CALC */
364  #endif /* GM_VISBECK_VARIABLE_K */  #endif /* GM_VISBECK_VARIABLE_K */
365           DO j=1-Oly,sNy+Oly
366            DO i=1-Olx,sNx+Olx
367             dSigmaDr(i,j)=sigmaR(i,j,k)
368            ENDDO
369           ENDDO
370    
371    #ifdef ALLOW_AUTODIFF_TAMC
372    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
373    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
374    CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
375    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
376    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
377    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
378    #endif /* ALLOW_AUTODIFF_TAMC */
379    
380  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
381         CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
382       O             SlopeX, SlopeY,       O             SlopeX, SlopeY,
383       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
384         U             hTransLay, baseSlope, recipLambda,
385       U             dSigmaDr,       U             dSigmaDr,
386       I             dSigmaDx, dSigmaDy,       I             dSigmaDx, dSigmaDy,
387       I             ldd97_LrhoC,rF(k),k,       I             ldd97_LrhoC, locMixLayer, rF,
388       I             bi, bj, myThid )       I             kLowC(1-Olx,1-Oly,bi,bj),
389         I             k, bi, bj, myTime, myIter, myThid )
390    
391         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
392          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
# Line 298  CADJ STORE dSigmaDr(:,:)     = comlev1_b Line 405  CADJ STORE dSigmaDr(:,:)     = comlev1_b
405  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
406  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
407    
408    C      Components of Redi/GM tensor
409         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
410          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
411  C      Components of Redi/GM tensor            Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
412           Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)            Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
413           Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)            Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
414           Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)          ENDDO
415           ENDDO
416    
417  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
418  #ifdef OLD_VISBECK_CALC  #ifdef OLD_VISBECK_CALC
419           DO j=1-Oly+1,sNy+Oly-1
420            DO i=1-Olx+1,sNx+Olx-1
421    
422  C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K  C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K
423  C           but do not know if *taperFct (or **2 ?) is necessary  C           but do not know if *taperFct (or **2 ?) is necessary
424          Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)          Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)
425    
426  C--     Depth average of M^2/N^2 * N  C--     Depth average of M^2/N^2 * N
# Line 319  C       which is used in the "variable K Line 430  C       which is used in the "variable K
430  C       Distance between interface above layer and the integration depth  C       Distance between interface above layer and the integration depth
431          deltaH=abs(GM_Visbeck_depth)-abs(rF(k))          deltaH=abs(GM_Visbeck_depth)-abs(rF(k))
432  C       If positive we limit this to the layer thickness  C       If positive we limit this to the layer thickness
433          deltaH=min(deltaH,drF(k))          integrDepth = drF(k)
434            deltaH=min(deltaH,integrDepth)
435  C       If negative then we are below the integration level  C       If negative then we are below the integration level
436          deltaH=max(deltaH,zero_rs)          deltaH=max(deltaH, 0. _d 0)
437  C       Now we convert deltaH to a non-dimensional fraction  C       Now we convert deltaH to a non-dimensional fraction
438          deltaH=deltaH/GM_Visbeck_depth          deltaH=deltaH/GM_Visbeck_depth
439    
         IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.  
440          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
441           N2= -Gravity*recip_RhoConst*dSigmaDr(i,j)           N2loc = -gravity*recip_rhoConst*dSigmaDr(i,j)
442           SN=sqrt(Ssq(i,j)*N2)           SNloc = SQRT(Ssq(i,j)*N2loc )
443           VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH           VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
444       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN       &       +deltaH*GM_Visbeck_alpha
445         &              *GM_Visbeck_length*GM_Visbeck_length*SNloc
446          ENDIF          ENDIF
447    
 #endif /* OLD_VISBECK_CALC */  
 #endif /* GM_VISBECK_VARIABLE_K */  
448          ENDDO          ENDDO
449         ENDDO         ENDDO
450    #endif /* OLD_VISBECK_CALC */
451    #endif /* GM_VISBECK_VARIABLE_K */
452    
453  C-- end 1rst loop on vertical level index k  C-- end 1rst loop on vertical level index k
454        ENDDO        ENDDO
# Line 346  C-- end 1rst loop on vertical level inde Line 458  C-- end 1rst loop on vertical level inde
458  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
459  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
460  #endif  #endif
461        IF ( GM_Visbeck_alpha.NE.0. ) THEN        IF ( GM_Visbeck_alpha.GT.0. ) THEN
462  C-     Limit range that KapGM can take  C-     Limit range that KapGM can take
463         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
464          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
465           VisbeckK(i,j,bi,bj)=           VisbeckK(i,j,bi,bj)=
466       &       MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)       &       MIN( MAX( VisbeckK(i,j,bi,bj), GM_Visbeck_minVal_K ),
467         &            GM_Visbeck_maxVal_K )
468          ENDDO          ENDDO
469         ENDDO         ENDDO
470        ENDIF        ENDIF
# Line 362  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1 Line 475  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1
475  cph)  cph)
476  #endif /* GM_VISBECK_VARIABLE_K */  #endif /* GM_VISBECK_VARIABLE_K */
477    
478    C-    express the Tensor in term of Diffusivity (= m**2 / s )
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
   
 C-- 2nd loop on k : compute Tensor Coeff. at U,V levels.  
479        DO k=1,Nr        DO k=1,Nr
        kp1 = MIN(Nr,k+1)  
        maskp1 = 1. _d 0  
        IF (k.GE.Nr) maskp1 = 0. _d 0  
   
480  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
481         kkey = (igmkey-1)*Nr + k         kkey = (igmkey-1)*Nr + k
482  #if (defined (GM_NON_UNITY_DIAGONAL) || \  # if (defined (GM_NON_UNITY_DIAGONAL) || \
483       defined (GM_VISBECK_VARIABLE_K))        defined (GM_VISBECK_VARIABLE_K))
484  CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
485  CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
486  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
487    # endif
488  #endif  #endif
489  #endif         km1 = MAX(k-1,1)
490           isopycK = GM_isopycK
491  C-    express the Tensor in term of Diffusivity (= m**2 / s )       &         *(GM_isoFac1d(km1)+GM_isoFac1d(k))*op5
492           bolus_K = GM_background_K
493         &         *(GM_bolFac1d(km1)+GM_bolFac1d(k))*op5
494         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
495          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
496    #ifdef ALLOW_KAPREDI_CONTROL
497             Kgm_tmp = kapredi(i,j,k,bi,bj)
498    #else
499             Kgm_tmp = isopycK*GM_isoFac2d(i,j,bi,bj)
500    #endif
501  #ifdef ALLOW_KAPGM_CONTROL  #ifdef ALLOW_KAPGM_CONTROL
502           Kgm_tmp = GM_isopycK + GM_skewflx*kapgm(i,j,k,bi,bj)       &           + GM_skewflx*kapgm(i,j,k,bi,bj)
503  #else  #else
504           Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K       &           + GM_skewflx*bolus_K*GM_bolFac2d(i,j,bi,bj)
505  #endif  #endif
506  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
507       &           + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)       &           + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)
508  #endif  #endif
509           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)
510           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)
511           Kwz(i,j,k,bi,bj)= ( GM_isopycK  #ifdef ALLOW_KAPREDI_CONTROL
512             Kwz(i,j,k,bi,bj)= ( kapredi(i,j,k,bi,bj)
513    #else
514             Kwz(i,j,k,bi,bj)= ( GM_isopycK*GM_isoFac2d(i,j,bi,bj)
515    #endif
516  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
517       &                     + VisbeckK(i,j,bi,bj)       &                     + VisbeckK(i,j,bi,bj)
518  #endif  #endif
519       &                     )*Kwz(i,j,k,bi,bj)       &                     )*Kwz(i,j,k,bi,bj)
520          ENDDO          ENDDO
521         ENDDO         ENDDO
522          ENDDO
523    
524    #ifdef ALLOW_DIAGNOSTICS
525          IF ( useDiagnostics .AND. GM_taper_scheme.EQ.'fm07' ) THEN
526           CALL DIAGNOSTICS_FILL( hTransLay, 'GM_hTrsL', 0,1,2,bi,bj,myThid)
527           CALL DIAGNOSTICS_FILL( baseSlope, 'GM_baseS', 0,1,2,bi,bj,myThid)
528           CALL DIAGNOSTICS_FILL(recipLambda,'GM_rLamb', 0,1,2,bi,bj,myThid)
529          ENDIF
530    #endif /* ALLOW_DIAGNOSTICS */
531    
532    
533  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
534    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
535    C-- 2nd  k loop : compute Tensor Coeff. at U point
536    
537    #ifdef ALLOW_KPP
538          IF ( useKPP ) THEN
539           DO j=1-Oly,sNy+Oly
540            DO i=2-Olx,sNx+Olx
541             locMixLayer(i,j) = ( KPPhbl(i-1,j,bi,bj)
542         &                      + KPPhbl( i ,j,bi,bj) )*op5
543            ENDDO
544           ENDDO
545          ELSE
546    #else
547          IF ( .TRUE. ) THEN
548    #endif
549           DO j=1-Oly,sNy+Oly
550            DO i=2-Olx,sNx+Olx
551             locMixLayer(i,j) = ( hMixLayer(i-1,j,bi,bj)
552         &                      + hMixLayer( i ,j,bi,bj) )*op5
553            ENDDO
554           ENDDO
555          ENDIF
556          DO j=1-Oly,sNy+Oly
557           DO i=1-Olx,sNx+Olx
558             hTransLay(i,j) =  0.
559             baseSlope(i,j) =  0.
560             recipLambda(i,j)= 0.
561           ENDDO
562           DO i=2-Olx,sNx+Olx
563             hTransLay(i,j) = MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
564           ENDDO
565          ENDDO
566    
567          DO k=Nr,1,-1
568           kp1 = MIN(Nr,k+1)
569           maskp1 = 1. _d 0
570           IF (k.GE.Nr) maskp1 = 0. _d 0
571    #ifdef ALLOW_AUTODIFF_TAMC
572           kkey = (igmkey-1)*Nr + k
573    #endif
574    
575  C     Gradient of Sigma at U points  C     Gradient of Sigma at U points
576         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
# Line 422  C     Gradient of Sigma at U points Line 590  C     Gradient of Sigma at U points
590  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
591  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
592  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
593  CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
594    CADJ STORE locMixLayer(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
595    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
596    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
597    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
598  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
599    
600  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
601         CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
602       O             SlopeX, SlopeY,       O             SlopeX, SlopeY,
603       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
604         U             hTransLay, baseSlope, recipLambda,
605       U             dSigmaDr,       U             dSigmaDr,
606       I             dSigmaDx, dSigmaDy,       I             dSigmaDx, dSigmaDy,
607       I             ldd97_LrhoW,rC(k),k,       I             ldd97_LrhoW, locMixLayer, rC,
608       I             bi, bj, myThid )       I             kLow_W,
609         I             k, bi, bj, myTime, myIter, myThid )
610    
 cph( NEW  
611  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
 cph(  
612  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
613  CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
 cph)  
614  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
 cph)  
615    
616  #ifdef GM_NON_UNITY_DIAGONAL  #ifdef GM_NON_UNITY_DIAGONAL
617  c      IF ( GM_nonUnitDiag ) THEN  c      IF ( GM_nonUnitDiag ) THEN
618          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
619           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
620            Kux(i,j,k,bi,bj) =            Kux(i,j,k,bi,bj) =
621       &     ( GM_isopycK  #ifdef ALLOW_KAPREDI_CONTROL
622         &     ( kapredi(i,j,k,bi,bj)
623    #else
624         &     ( GM_isopycK*GM_isoFac1d(k)
625         &        *op5*(GM_isoFac2d(i-1,j,bi,bj)+GM_isoFac2d(i,j,bi,bj))
626    #endif
627  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
628       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
629  #endif  #endif
630       &     )       &     )*taperFct(i,j)
      &     *taperFct(i,j)  
631           ENDDO           ENDDO
632          ENDDO          ENDDO
633  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 475  c      ENDIF Line 649  c      ENDIF
649  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
650  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
651  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
652         IF (GM_ExtraDiag) THEN         IF ( GM_ExtraDiag ) THEN
653          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
654           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
655            Kuz(i,j,k,bi,bj) =            Kuz(i,j,k,bi,bj) =
656    #ifdef ALLOW_KAPREDI_CONTROL
657         &     ( kapredi(i,j,k,bi,bj)
658    #else
659         &     ( GM_isopycK*GM_isoFac1d(k)
660         &        *op5*(GM_isoFac2d(i-1,j,bi,bj)+GM_isoFac2d(i,j,bi,bj))
661    #endif
662  #ifdef ALLOW_KAPGM_CONTROL  #ifdef ALLOW_KAPGM_CONTROL
663       &     ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)       &     - GM_skewflx*kapgm(i,j,k,bi,bj)
664  #else  #else
665       &     ( GM_isopycK - GM_skewflx*GM_background_K       &     - GM_skewflx*GM_background_K*GM_bolFac1d(k)
666         &        *op5*(GM_bolFac2d(i-1,j,bi,bj)+GM_bolFac2d(i,j,bi,bj))
667  #endif  #endif
668  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
669       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*GM_advect       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*GM_advect
# Line 499  CADJ STORE taperFct(:,:)     = comlev1_b Line 680  CADJ STORE taperFct(:,:)     = comlev1_b
680          DO j=1,sNy          DO j=1,sNy
681           DO i=1,sNx+1           DO i=1,sNx+1
682  C         store in tmp1k Kuz_Redi  C         store in tmp1k Kuz_Redi
683    #ifdef ALLOW_KAPREDI_CONTROL
684              tmp1k(i,j) = ( kapredi(i,j,k,bi,bj)
685    #else
686            tmp1k(i,j) = ( GM_isopycK            tmp1k(i,j) = ( GM_isopycK
687    #endif
688  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
689       &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0       &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0
690  #endif  #endif
# Line 531  C-        Vertical gradients interpolate Line 716  C-        Vertical gradients interpolate
716         ENDIF         ENDIF
717  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
718    
719    C-- end 2nd  loop on vertical level index k
720          ENDDO
721    
722    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
723    C-- 3rd  k loop : compute Tensor Coeff. at V point
724    
725    #ifdef ALLOW_KPP
726          IF ( useKPP ) THEN
727           DO j=2-Oly,sNy+Oly
728            DO i=1-Olx,sNx+Olx
729             locMixLayer(i,j) = ( KPPhbl(i,j-1,bi,bj)
730         &                      + KPPhbl(i, j ,bi,bj) )*op5
731            ENDDO
732           ENDDO
733          ELSE
734    #else
735          IF ( .TRUE. ) THEN
736    #endif
737           DO j=2-Oly,sNy+Oly
738            DO i=1-Olx,sNx+Olx
739             locMixLayer(i,j) = ( hMixLayer(i,j-1,bi,bj)
740         &                      + hMixLayer(i, j ,bi,bj) )*op5
741            ENDDO
742           ENDDO
743          ENDIF
744          DO j=1-Oly,sNy+Oly
745           DO i=1-Olx,sNx+Olx
746             hTransLay(i,j) =  0.
747             baseSlope(i,j) =  0.
748             recipLambda(i,j)= 0.
749           ENDDO
750          ENDDO
751          DO j=2-Oly,sNy+Oly
752           DO i=1-Olx,sNx+Olx
753             hTransLay(i,j) = MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
754           ENDDO
755          ENDDO
756    
757  C     Gradient of Sigma at V points  C     Gradient of Sigma at V points
758          DO k=Nr,1,-1
759           kp1 = MIN(Nr,k+1)
760           maskp1 = 1. _d 0
761           IF (k.GE.Nr) maskp1 = 0. _d 0
762    #ifdef ALLOW_AUTODIFF_TAMC
763           kkey = (igmkey-1)*Nr + k
764    #endif
765    
766         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
767          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
768           dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)           dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
# Line 548  C     Gradient of Sigma at V points Line 779  C     Gradient of Sigma at V points
779  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
780  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
781  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
782  CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
783    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
784    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
785    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
786  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
787    
788  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
789         CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
790       O             SlopeX, SlopeY,       O             SlopeX, SlopeY,
791       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
792         U             hTransLay, baseSlope, recipLambda,
793       U             dSigmaDr,       U             dSigmaDr,
794       I             dSigmaDx, dSigmaDy,       I             dSigmaDx, dSigmaDy,
795       I             ldd97_LrhoS,rC(k),k,       I             ldd97_LrhoS, locMixLayer, rC,
796       I             bi, bj, myThid )       I             kLow_S,
797         I             k, bi, bj, myTime, myIter, myThid )
798    
799  cph(  cph(
800  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 573  c      IF ( GM_nonUnitDiag ) THEN Line 809  c      IF ( GM_nonUnitDiag ) THEN
809          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
810           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
811            Kvy(i,j,k,bi,bj) =            Kvy(i,j,k,bi,bj) =
812       &     ( GM_isopycK  #ifdef ALLOW_KAPREDI_CONTROL
813         &     ( kapredi(i,j,k,bi,bj)
814    #else
815         &     ( GM_isopycK*GM_isoFac1d(k)
816         &        *op5*(GM_isoFac2d(i,j-1,bi,bj)+GM_isoFac2d(i,j,bi,bj))
817    #endif
818  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
819       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
820  #endif  #endif
821       &     )       &     )*taperFct(i,j)
      &     *taperFct(i,j)  
822           ENDDO           ENDDO
823          ENDDO          ENDDO
824  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 600  c      ENDIF Line 840  c      ENDIF
840  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
841  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
842  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
843         IF (GM_ExtraDiag) THEN         IF ( GM_ExtraDiag ) THEN
844          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
845           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
846            Kvz(i,j,k,bi,bj) =            Kvz(i,j,k,bi,bj) =
847    #ifdef ALLOW_KAPREDI_CONTROL
848         &     ( kapredi(i,j,k,bi,bj)
849    #else
850         &     ( GM_isopycK*GM_isoFac1d(k)
851         &        *op5*(GM_isoFac2d(i,j-1,bi,bj)+GM_isoFac2d(i,j,bi,bj))
852    #endif
853  #ifdef ALLOW_KAPGM_CONTROL  #ifdef ALLOW_KAPGM_CONTROL
854       &     ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)       &     - GM_skewflx*kapgm(i,j,k,bi,bj)
855  #else  #else
856       &     ( GM_isopycK - GM_skewflx*GM_background_K       &     - GM_skewflx*GM_background_K*GM_bolFac1d(k)
857         &        *op5*(GM_bolFac2d(i,j-1,bi,bj)+GM_bolFac2d(i,j,bi,bj))
858  #endif  #endif
859  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
860       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*GM_advect       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*GM_advect
# Line 620  CADJ STORE taperFct(:,:)     = comlev1_b Line 867  CADJ STORE taperFct(:,:)     = comlev1_b
867    
868  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
869         IF (doDiagRediFlx) THEN         IF (doDiagRediFlx) THEN
870  c       km1 = MAX(k-1,1)          km1 = MAX(k-1,1)
871          DO j=1,sNy+1          DO j=1,sNy+1
872           DO i=1,sNx           DO i=1,sNx
873  C         store in tmp1k Kvz_Redi  C         store in tmp1k Kvz_Redi
874    #ifdef ALLOW_KAPREDI_CONTROL
875              tmp1k(i,j) = ( kapredi(i,j,k,bi,bj)
876    #else
877            tmp1k(i,j) = ( GM_isopycK            tmp1k(i,j) = ( GM_isopycK
878    #endif
879  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
880       &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0       &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0
881  #endif  #endif
# Line 656  C-        Vertical gradients interpolate Line 907  C-        Vertical gradients interpolate
907         ENDIF         ENDIF
908  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
909    
910  #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */  C-- end 3rd  loop on vertical level index k
   
 C-- end 2nd loop on vertical level index k  
911        ENDDO        ENDDO
912    
913    #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
914    
915    
916  #ifdef GM_BOLUS_ADVEC  #ifdef GM_BOLUS_ADVEC
917        IF (GM_AdvForm) THEN        IF (GM_AdvForm) THEN
# Line 696  C--   Time-average Line 947  C--   Time-average
947       &                          deltaTclock, bi, bj, myThid )       &                          deltaTclock, bi, bj, myThid )
948         ENDIF         ENDIF
949  #endif  #endif
950         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  
951    
952        ENDIF        ENDIF
953  #endif /* ALLOW_TIMEAVE */  #endif /* ALLOW_TIMEAVE */
# Line 714  C--   Time-average Line 963  C--   Time-average
963        RETURN        RETURN
964        END        END
965    
966    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
967    
968        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
969       I             bi, bj, iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
970       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
971       I             myThid )       I             bi, bj, myTime, myIter, myThid )
972  C     /==========================================================\  C     /==========================================================\
973  C     | SUBROUTINE GMREDI_CALC_TENSOR                            |  C     | SUBROUTINE GMREDI_CALC_TENSOR                            |
974  C     | o Calculate tensor elements for GM/Redi tensor.          |  C     | o Calculate tensor elements for GM/Redi tensor.          |
# Line 736  C Line 986  C
986        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
987        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
988        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
989        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
990          INTEGER bi, bj
991          _RL     myTime
992          INTEGER myIter
993        INTEGER myThid        INTEGER myThid
994  CEndOfInterface  CEndOfInterface
995    
       INTEGER i, j, k  
   
996  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
997    
998          INTEGER i, j, k
999    
1000        DO k=1,Nr        DO k=1,Nr
1001         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
1002          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
Legend:
Removed from v.1.26  
changed lines
  Added in v.1.37
  ViewVC Help
Powered by ViewVC 1.1.22