/[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.14 by heimbach, Fri Jan 10 00:48:39 2003 UTC revision 1.27 by jmc, Thu Jun 21 01:33:01 2007 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    #ifdef ALLOW_KPP
6    # include "KPP_OPTIONS.h"
7    #endif
8    #undef OLD_VISBECK_CALC
9    
10  CStartOfInterface  CBOP
11    C     !ROUTINE: GMREDI_CALC_TENSOR
12    C     !INTERFACE:
13        SUBROUTINE GMREDI_CALC_TENSOR(        SUBROUTINE GMREDI_CALC_TENSOR(
14       I             bi, bj, iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
15       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
16       I             myThid )       I             bi, bj, myTime, myIter, myThid )
17  C     /==========================================================\  
18  C     | SUBROUTINE GMREDI_CALC_TENSOR                            |  C     !DESCRIPTION: \bv
19  C     | o Calculate tensor elements for GM/Redi tensor.          |  C     *==========================================================*
20  C     |==========================================================|  C     | SUBROUTINE GMREDI_CALC_TENSOR
21  C     \==========================================================/  C     | o Calculate tensor elements for GM/Redi tensor.
22    C     *==========================================================*
23    C     *==========================================================*
24    C     \ev
25    
26    C     !USES:
27        IMPLICIT NONE        IMPLICIT NONE
28    
29  C     == Global variables ==  C     == Global variables ==
# Line 22  C     == Global variables == Line 33  C     == Global variables ==
33  #include "EEPARAMS.h"  #include "EEPARAMS.h"
34  #include "PARAMS.h"  #include "PARAMS.h"
35  #include "GMREDI.h"  #include "GMREDI.h"
36  #include "GMREDI_DIAGS.h"  #include "GMREDI_TAVE.h"
37    #ifdef ALLOW_KPP
38    # include "KPP.h"
39    #endif
40    
41  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
42  #include "tamc.h"  #include "tamc.h"
43  #include "tamc_keys.h"  #include "tamc_keys.h"
44  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
45    
46    C     !INPUT/OUTPUT PARAMETERS:
47  C     == Routine arguments ==  C     == Routine arguments ==
48    C     bi, bj    :: tile indices
49    C     myTime    :: Current time in simulation
50    C     myIter    :: Current iteration number in simulation
51    C     myThid    :: My Thread Id. number
52  C  C
53          INTEGER iMin,iMax,jMin,jMax
54        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
55        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
56        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
57        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi, bj
58          _RL     myTime
59          INTEGER myIter
60        INTEGER myThid        INTEGER myThid
61  CEndOfInterface  CEOP
62    
63  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
64    
65    C     !LOCAL VARIABLES:
66  C     == Local variables ==  C     == Local variables ==
67        INTEGER i,j,k,km1,kp1        INTEGER i,j,k,kp1
68        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
69        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
70        _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
71        _RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
72        _RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
73        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
74        _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
75        _RL maskp1, maskm1, Kgm_tmp        _RL maskp1, Kgm_tmp
76          _RL ldd97_LrhoC(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
77          _RL ldd97_LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
78          _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
79          _RL Cspd, LrhoInf, LrhoSup, fCoriLoc
80    
81          INTEGER kLow_W (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
82          INTEGER kLow_S (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
83          _RL locMixLayer(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
84          _RL baseSlope  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
85          _RL hTransLay  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
86          _RL recipLambda(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
87    
88  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
89    #ifdef OLD_VISBECK_CALC
90        _RL deltaH,zero_rs        _RL deltaH,zero_rs
91        PARAMETER(zero_rs=0.D0)        PARAMETER(zero_rs=0.D0)
92        _RL N2,SN        _RL N2,SN
93        _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
94    #else
95          _RL dSigmaH
96          _RL deltaH, integrDepth
97          _RL Sloc, M2loc, SNloc
98    #endif
99  #endif  #endif
100    
101    #ifdef ALLOW_DIAGNOSTICS
102          LOGICAL  doDiagRediFlx
103          LOGICAL  DIAGNOSTICS_IS_ON
104          EXTERNAL DIAGNOSTICS_IS_ON
105          INTEGER  km1
106          _RL dTdz
107          _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108    #endif
109    
110    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
111    
112  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
113            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
114            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
# Line 71  C     == Local variables == Line 122  C     == Local variables ==
122       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
123  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
124    
125    #ifdef ALLOW_DIAGNOSTICS
126          doDiagRediFlx = .FALSE.
127          IF ( useDiagnostics ) THEN
128            doDiagRediFlx = DIAGNOSTICS_IS_ON('GM_KuzTz', myThid )
129            doDiagRediFlx = doDiagRediFlx .OR.
130         &                  DIAGNOSTICS_IS_ON('GM_KvzTz', myThid )
131          ENDIF
132    #endif
133    
134  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
135        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
136         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
# Line 79  C     == Local variables == Line 139  C     == Local variables ==
139        ENDDO        ENDDO
140  #endif  #endif
141    
142        DO k=2,Nr  C--   set ldd97_Lrho (for tapering scheme ldd97):
143          IF ( GM_taper_scheme.EQ.'ldd97' .OR.
144         &     GM_taper_scheme.EQ.'fm07' ) THEN
145           Cspd = 2. _d 0
146           LrhoInf = 15. _d 3
147           LrhoSup = 100. _d 3
148    C-     Tracer point location (center):
149           DO j=1-Oly,sNy+Oly
150            DO i=1-Olx,sNx+Olx
151             IF (fCori(i,j,bi,bj).NE.0.) THEN
152               ldd97_LrhoC(i,j) = Cspd/ABS(fCori(i,j,bi,bj))
153             ELSE
154               ldd97_LrhoC(i,j) = LrhoSup
155             ENDIF
156             ldd97_LrhoC(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoC(i,j),LrhoSup))
157            ENDDO
158           ENDDO
159    C-     U point location (West):
160           DO j=1-Oly,sNy+Oly
161            kLow_W(1-Olx,j) = 0
162            ldd97_LrhoW(1-Olx,j) = LrhoSup
163            DO i=1-Olx+1,sNx+Olx
164             kLow_W(i,j) = MIN(kLowC(i-1,j,bi,bj),kLowC(i,j,bi,bj))
165             fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))
166             IF (fCoriLoc.NE.0.) THEN
167               ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)
168             ELSE
169               ldd97_LrhoW(i,j) = LrhoSup
170             ENDIF
171             ldd97_LrhoW(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoW(i,j),LrhoSup))
172            ENDDO
173           ENDDO
174    C-     V point location (South):
175           DO i=1-Olx+1,sNx+Olx
176             kLow_S(i,1-Oly) = 0
177             ldd97_LrhoS(i,1-Oly) = LrhoSup
178           ENDDO
179           DO j=1-Oly+1,sNy+Oly
180            DO i=1-Olx,sNx+Olx
181             kLow_S(i,j) = MIN(kLowC(i,j-1,bi,bj),kLowC(i,j,bi,bj))
182             fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))
183             IF (fCoriLoc.NE.0.) THEN
184               ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)
185             ELSE
186               ldd97_LrhoS(i,j) = LrhoSup
187             ENDIF
188             ldd97_LrhoS(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoS(i,j),LrhoSup))
189            ENDDO
190           ENDDO
191          ELSE
192    C-    Just initialize to zero (not use anyway)
193           DO j=1-Oly,sNy+Oly
194            DO i=1-Olx,sNx+Olx
195              ldd97_LrhoC(i,j) = 0. _d 0
196              ldd97_LrhoW(i,j) = 0. _d 0
197              ldd97_LrhoS(i,j) = 0. _d 0
198            ENDDO
199           ENDDO
200          ENDIF
201    
202    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
203  C-- 1rst loop on k : compute Tensor Coeff. at W points.  C-- 1rst loop on k : compute Tensor Coeff. at W points.
204         km1 = MAX(1,k-1)  
205         maskm1 = 1. _d 0  #ifdef ALLOW_KPP
206         IF (k.LE.1) maskm1 = 0. _d 0        IF ( useKPP ) THEN
207           DO j=1-Oly,sNy+Oly
208            DO i=1-Olx,sNx+Olx
209             locMixLayer(i,j) = KPPhbl(i,j,bi,bj)
210            ENDDO
211           ENDDO
212          ELSE
213    #else
214          IF ( .TRUE. ) THEN
215    #endif
216           DO j=1-Oly,sNy+Oly
217            DO i=1-Olx,sNx+Olx
218             locMixLayer(i,j) = hMixLayer(i,j,bi,bj)
219            ENDDO
220           ENDDO
221          ENDIF
222          DO j=1-Oly,sNy+Oly
223           DO i=1-Olx,sNx+Olx
224             hTransLay(i,j) = R_low(i,j,bi,bj)
225             baseSlope(i,j) =  0.
226             recipLambda(i,j)= 0.
227           ENDDO
228          ENDDO
229    
230          DO k=Nr,2,-1
231    
232  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
233         kkey = (igmkey-1)*Nr + k         kkey = (igmkey-1)*Nr + k
# Line 93  C-- 1rst loop on k : compute Tensor Coef Line 237  C-- 1rst loop on k : compute Tensor Coef
237           SlopeY(i,j)       = 0. _d 0           SlopeY(i,j)       = 0. _d 0
238           dSigmaDx(i,j)     = 0. _d 0           dSigmaDx(i,j)     = 0. _d 0
239           dSigmaDy(i,j)     = 0. _d 0           dSigmaDy(i,j)     = 0. _d 0
240           dSigmaDrReal(i,j) = 0. _d 0           dSigmaDr(i,j)     = 0. _d 0
241           SlopeSqr(i,j)     = 0. _d 0           SlopeSqr(i,j)     = 0. _d 0
242           taperFct(i,j)     = 0. _d 0           taperFct(i,j)     = 0. _d 0
243           Kwx(i,j,k,bi,bj)  = 0. _d 0           Kwx(i,j,k,bi,bj)  = 0. _d 0
# Line 115  C-- 1rst loop on k : compute Tensor Coef Line 259  C-- 1rst loop on k : compute Tensor Coef
259         ENDDO         ENDDO
260  #endif  #endif
261    
262        DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
263         DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
264  C      Gradient of Sigma at rVel points  C      Gradient of Sigma at rVel points
265          dSigmaDx(i,j)=op25*( sigmaX(i+1, j ,km1) +sigmaX(i,j,km1)           dSigmaDx(i,j)=op25*( sigmaX(i+1,j,k-1)+sigmaX(i,j,k-1)
266       &                    +sigmaX(i+1, j , k ) +sigmaX(i,j, k ) )       &                       +sigmaX(i+1,j, k )+sigmaX(i,j, k )
267       &                  *maskC(i,j,k,bi,bj)*maskm1       &                      )*maskC(i,j,k,bi,bj)
268          dSigmaDy(i,j)=op25*( sigmaY( i ,j+1,km1) +sigmaY(i,j,km1)           dSigmaDy(i,j)=op25*( sigmaY(i,j+1,k-1)+sigmaY(i,j,k-1)
269       &                    +sigmaY( i ,j+1, k ) +sigmaY(i,j, k ) )       &                       +sigmaY(i,j+1, k )+sigmaY(i,j, k )
270       &                  *maskC(i,j,k,bi,bj)*maskm1       &                      )*maskC(i,j,k,bi,bj)
271          dSigmaDrReal(i,j)=sigmaR(i,j,k)*maskm1           dSigmaDr(i,j)=sigmaR(i,j,k)
272            ENDDO
273         ENDDO         ENDDO
       ENDDO  
274    
275  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
276  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
277  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
278  CADJ STORE dsigmadrreal(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
279  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
280    
281  C     Calculate slopes for use in tensor, taper and/or clip  #ifdef GM_VISBECK_VARIABLE_K
282        CALL GMREDI_SLOPE_LIMIT(  #ifndef OLD_VISBECK_CALC
283       U             dSigmadRReal,         IF ( GM_Visbeck_alpha.GT.0. .AND.
284       I             rF(K),K,       &      -rC(k-1).LT.GM_Visbeck_depth ) THEN
285       U             SlopeX, SlopeY,  
286       U             dSigmaDx, dSigmaDy,  C--     Depth average of f/sqrt(Ri) = M^2/N^2 * N
287       O             SlopeSqr, taperFct,  C       M^2 and N^2 are horizontal & vertical gradient of buoyancy.
288       I             bi, bj, myThid )  
289    C       Calculate terms for mean Richardson number which is used
290    C       in the "variable K" parameterisaton:
291    C       compute depth average from surface down to the bottom or
292    C       GM_Visbeck_depth, whatever is the shallower.
293    
294        DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
295         DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
296              IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
297               integrDepth = -rC( kLowC(i,j,bi,bj) )
298    C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments
299               integrDepth = MIN( integrDepth, GM_Visbeck_depth )
300    C       Distance between level center above and the integration depth
301               deltaH = integrDepth + rC(k-1)
302    C       If negative then we are below the integration level
303    C       (cannot be the case with 2 conditions on maskC & -rC(k-1))
304    C       If positive we limit this to the distance from center above
305               deltaH = MIN( deltaH, drC(k) )
306    C       Now we convert deltaH to a non-dimensional fraction
307               deltaH = deltaH/( integrDepth+rC(1) )
308    
309  C       Mask Iso-neutral slopes  C--      compute: ( M^2 * S )^1/2   (= M^2 / N since S=M^2/N^2 )
310          SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)*maskm1             dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)
311          SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)*maskm1       &             + dSigmaDy(i,j)*dSigmaDy(i,j)
312          SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)*maskm1             IF ( dSigmaH .GT. 0. _d 0 ) THEN
313                 dSigmaH = SQRT( dSigmaH )
314    C-       compute slope, limited by GM_maxSlope:
315                 IF ( -dSigmaDr(i,j).GT.dSigmaH*GM_rMaxSlope ) THEN
316                  Sloc = dSigmaH / ( -dSigmaDr(i,j) )
317                 ELSE
318                  Sloc = GM_maxSlope
319                 ENDIF
320                 M2loc = Gravity*recip_RhoConst*dSigmaH
321                 SNloc = SQRT( Sloc*M2loc )
322               ELSE
323                 SNloc = 0. _d 0
324               ENDIF
325               VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
326         &       +deltaH*GM_Visbeck_alpha
327         &              *GM_Visbeck_length*GM_Visbeck_length*SNloc
328              ENDIF
329             ENDDO
330            ENDDO
331           ENDIF
332    #endif /* ndef OLD_VISBECK_CALC */
333    #endif /* GM_VISBECK_VARIABLE_K */
334    
335    C     Calculate slopes for use in tensor, taper and/or clip
336           CALL GMREDI_SLOPE_LIMIT(
337         O             SlopeX, SlopeY,
338         O             SlopeSqr, taperFct,
339         U             hTransLay, baseSlope, recipLambda,
340         U             dSigmaDr,
341         I             dSigmaDx, dSigmaDy,
342         I             ldd97_LrhoC, locMixLayer, rF,
343         I             kLowC(1-Olx,1-Oly,bi,bj),
344         I             k, bi, bj, myTime, myIter, myThid )
345    
346           DO j=1-Oly+1,sNy+Oly-1
347            DO i=1-Olx+1,sNx+Olx-1
348    C      Mask Iso-neutral slopes
349             SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)
350             SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)
351             SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)
352            ENDDO
353         ENDDO         ENDDO
       ENDDO  
354    
355  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
356    CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
357    CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
358  CADJ STORE SlopeSqr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeSqr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
359    CADJ STORE dSigmaDr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
360    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
361  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
362    
363        DO j=1-Oly+1,sNy+Oly-1  C      Components of Redi/GM tensor
364         DO i=1-Olx+1,sNx+Olx-1         DO j=1-Oly+1,sNy+Oly-1
365            DO i=1-Olx+1,sNx+Olx-1
366  C       Components of Redi/GM tensor            Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
367          Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)            Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
368          Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)            Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
369          Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)          ENDDO
370           ENDDO
371    
372  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
373    #ifdef OLD_VISBECK_CALC
374           DO j=1-Oly+1,sNy+Oly-1
375            DO i=1-Olx+1,sNx+Olx-1
376    
377  C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K  C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K
378  C           but don't know if *taperFct (or **2 ?) is necessary  C           but do not know if *taperFct (or **2 ?) is necessary
379          Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)          Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)
380    
381  C--     Depth average of M^2/N^2 * N  C--     Depth average of M^2/N^2 * N
# Line 186  C       Now we convert deltaH to a non-d Line 392  C       Now we convert deltaH to a non-d
392          deltaH=deltaH/GM_Visbeck_depth          deltaH=deltaH/GM_Visbeck_depth
393    
394          IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.          IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.
395          IF ( Ssq(i,j).NE.0. .AND. dSigmaDrReal(i,j).NE.0. ) THEN          IF ( Ssq(i,j).NE.0. .AND. dSigmaDr(i,j).NE.0. ) THEN
396           N2= -Gravity*recip_RhoConst*dSigmaDrReal(i,j)           N2= -Gravity*recip_RhoConst*dSigmaDr(i,j)
397           SN=sqrt(Ssq(i,j)*N2)           SN=sqrt(Ssq(i,j)*N2)
398           VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH           VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH
399       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN
400          ENDIF          ENDIF
401    
402  #endif /* GM_VISBECK_VARIABLE_K */          ENDDO
   
403         ENDDO         ENDDO
404        ENDDO  #endif /* OLD_VISBECK_CALC */
405    #endif /* GM_VISBECK_VARIABLE_K */
406    
407  C-- end 1rst loop on vertical level index k  C-- end 1rst loop on vertical level index k
408        ENDDO        ENDDO
# Line 206  C-- end 1rst loop on vertical level inde Line 412  C-- end 1rst loop on vertical level inde
412  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
413  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
414  #endif  #endif
415        IF ( GM_Visbeck_alpha.NE.0. ) THEN        IF ( GM_Visbeck_alpha.GT.0. ) THEN
416  C-     Limit range that KapGM can take  C-     Limit range that KapGM can take
417         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
418          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
419           VisbeckK(i,j,bi,bj)=           VisbeckK(i,j,bi,bj)=
420       &       MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)       &       MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)
 #ifdef ALLOW_TIMEAVE  
          Visbeck_K_T(i,j,bi,bj)=Visbeck_K_T(i,j,bi,bj)  
      &                         +VisbeckK(i,j,bi,bj)*deltaTclock  
 #endif  
421          ENDDO          ENDDO
422         ENDDO         ENDDO
423        ENDIF        ENDIF
424    cph( NEW
425    #ifdef ALLOW_AUTODIFF_TAMC
426    CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
427    #endif
428    cph)
429  #endif /* GM_VISBECK_VARIABLE_K */  #endif /* GM_VISBECK_VARIABLE_K */
430    
   
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
   
 C-- 2nd loop on k : compute Tensor Coeff. at U,V levels.  
       DO k=1,Nr  
        kp1 = MIN(Nr,k+1)  
        maskp1 = 1. _d 0  
        IF (k.GE.Nr) maskp1 = 0. _d 0  
   
431  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
432         kkey = (igmkey-1)*Nr + k         kkey = (igmkey-1)*Nr + k
433  #ifdef GM_NON_UNITY_DIAGONAL  #if (defined (GM_NON_UNITY_DIAGONAL) || \
434         defined (GM_VISBECK_VARIABLE_K))
435  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
436  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
437  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
# Line 240  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bi Line 439  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bi
439  #endif  #endif
440    
441  C-    express the Tensor in term of Diffusivity (= m**2 / s )  C-    express the Tensor in term of Diffusivity (= m**2 / s )
442        DO j=1-Oly+1,sNy+Oly-1        DO k=1,Nr
443         DO i=1-Olx+1,sNx+Olx-1         DO j=1-Oly+1,sNy+Oly-1
444          Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K          DO i=1-Olx+1,sNx+Olx-1
445    #ifdef ALLOW_KAPGM_CONTROL
446             Kgm_tmp = GM_isopycK + GM_skewflx*kapgm(i,j,k,bi,bj)
447    #else
448             Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K
449    #endif
450  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
451       &          + VisbeckK(i,j,bi,bj)*(1.+GM_skewflx)           &           + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)
452  #endif  #endif
453          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)
454          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)
455          Kwz(i,j,k,bi,bj)= ( GM_isopycK           Kwz(i,j,k,bi,bj)= ( GM_isopycK
456  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
457       &                    + VisbeckK(i,j,bi,bj)       &                     + VisbeckK(i,j,bi,bj)
458  #endif  #endif
459       &                    )*Kwz(i,j,k,bi,bj)       &                     )*Kwz(i,j,k,bi,bj)
460            ENDDO
461         ENDDO         ENDDO
462        ENDDO        ENDDO
463    
464    #ifdef ALLOW_DIAGNOSTICS
465          IF ( useDiagnostics .AND. GM_taper_scheme.EQ.'fm07' ) THEN
466           CALL DIAGNOSTICS_FILL( hTransLay, 'GM_hTrsL', 0,1,2,bi,bj,myThid)
467           CALL DIAGNOSTICS_FILL( baseSlope, 'GM_baseS', 0,1,2,bi,bj,myThid)
468           CALL DIAGNOSTICS_FILL(recipLambda,'GM_rLamb', 0,1,2,bi,bj,myThid)
469          ENDIF
470    #endif /* ALLOW_DIAGNOSTICS */
471    
472    
473  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
474    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
475    C-- 2nd  k loop : compute Tensor Coeff. at U point
476    
477  C     Gradient of Sigma at U points  #ifdef ALLOW_KPP
478        DO j=1-Oly+1,sNy+Oly-1        IF ( useKPP ) THEN
479         DO i=1-Olx+1,sNx+Olx-1         DO j=1-Oly,sNy+Oly
480          dSigmaDx(i,j)=sigmaX(i,j,k)          DO i=2-Olx,sNx+Olx
481       &          *_maskW(i,j,k,bi,bj)           locMixLayer(i,j) = ( KPPhbl(i-1,j,bi,bj)
482          dSigmaDy(i,j)=op25*( sigmaY(i-1,j+1,k) +sigmaY(i,j+1,k)       &                      + KPPhbl( i ,j,bi,bj) )*op5
483       &                      +sigmaY(i-1, j ,k) +sigmaY(i, j ,k) )          ENDDO
484       &          *_maskW(i,j,k,bi,bj)         ENDDO
485          dSigmaDrReal(i,j)=op25*( sigmaR(i-1,j, k ) +sigmaR(i,j, k )        ELSE
486       &                  +maskp1*(sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1)) )  #else
487       &          *_maskW(i,j,k,bi,bj)*maskp1        IF ( .TRUE. ) THEN
488    #endif
489           DO j=1-Oly,sNy+Oly
490            DO i=2-Olx,sNx+Olx
491             locMixLayer(i,j) = ( hMixLayer(i-1,j,bi,bj)
492         &                      + hMixLayer( i ,j,bi,bj) )*op5
493            ENDDO
494           ENDDO
495          ENDIF
496          DO j=1-Oly,sNy+Oly
497           DO i=1-Olx,sNx+Olx
498             hTransLay(i,j) =  0.
499             baseSlope(i,j) =  0.
500             recipLambda(i,j)= 0.
501           ENDDO
502           DO i=2-Olx,sNx+Olx
503             hTransLay(i,j) = MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
504         ENDDO         ENDDO
505        ENDDO        ENDDO
506    
507          DO k=Nr,1,-1
508           kp1 = MIN(Nr,k+1)
509           maskp1 = 1. _d 0
510           IF (k.GE.Nr) maskp1 = 0. _d 0
511    
512    C     Gradient of Sigma at U points
513           DO j=1-Oly+1,sNy+Oly-1
514            DO i=1-Olx+1,sNx+Olx-1
515             dSigmaDx(i,j)=sigmaX(i,j,k)
516         &                       *_maskW(i,j,k,bi,bj)
517             dSigmaDy(i,j)=op25*( sigmaY(i-1,j+1,k)+sigmaY(i,j+1,k)
518         &                       +sigmaY(i-1, j ,k)+sigmaY(i, j ,k)
519         &                      )*_maskW(i,j,k,bi,bj)
520             dSigmaDr(i,j)=op25*( sigmaR(i-1,j, k )+sigmaR(i,j, k )
521         &                      +(sigmaR(i-1,j,kp1)+sigmaR(i,j,kp1))*maskp1
522         &                      )*_maskW(i,j,k,bi,bj)
523            ENDDO
524           ENDDO
525    
526  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
527    CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
528  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
529  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
530  CADJ STORE dsigmadrreal(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte
531  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
532    
533  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
534        CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
535       U             dSigmadRReal,       O             SlopeX, SlopeY,
      I             rF(K),K,  
      U             SlopeX, SlopeY,  
      U             dSigmaDx, dSigmaDy,  
536       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
537       I             bi, bj, myThid )       U             hTransLay, baseSlope, recipLambda,
538         U             dSigmaDr,
539         I             dSigmaDx, dSigmaDy,
540         I             ldd97_LrhoW, locMixLayer, rC,
541         I             kLow_W,
542         I             k, bi, bj, myTime, myIter, myThid )
543    
544    cph( NEW
545    #ifdef ALLOW_AUTODIFF_TAMC
546    cph(
547    CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
548    CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
549    cph)
550    #endif /* ALLOW_AUTODIFF_TAMC */
551    cph)
552    
553  #ifdef GM_NON_UNITY_DIAGONAL  #ifdef GM_NON_UNITY_DIAGONAL
554    c      IF ( GM_nonUnitDiag ) THEN
555          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
556           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
557            Kux(i,j,k,bi,bj) =            Kux(i,j,k,bi,bj) =
# Line 295  C     Calculate slopes for use in tensor Line 559  C     Calculate slopes for use in tensor
559  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
560       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
561  #endif  #endif
562       &     )       &     )*taperFct(i,j)
      &     *taperFct(i,j)  
563           ENDDO           ENDDO
564          ENDDO          ENDDO
565  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
566  # ifndef GM_TAPER_ORIG_CLIPPING  # ifdef GM_EXCLUDE_CLIPPING
567  CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
568  # endif  # endif
569  #endif  #endif
# Line 309  CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_b Line 572  CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_b
572            Kux(i,j,k,bi,bj) = MAX( Kux(i,j,k,bi,bj), GM_Kmin_horiz )            Kux(i,j,k,bi,bj) = MAX( Kux(i,j,k,bi,bj), GM_Kmin_horiz )
573           ENDDO           ENDDO
574          ENDDO          ENDDO
575    c      ENDIF
576  #endif /* GM_NON_UNITY_DIAGONAL */  #endif /* GM_NON_UNITY_DIAGONAL */
577    
578  #ifdef GM_EXTRA_DIAGONAL  #ifdef GM_EXTRA_DIAGONAL
# Line 317  CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_b Line 581  CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_b
581  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
582  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
583  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
584        IF (GM_ExtraDiag) THEN         IF ( GM_ExtraDiag ) THEN
585          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
586           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
587            Kuz(i,j,k,bi,bj) =            Kuz(i,j,k,bi,bj) =
588    #ifdef ALLOW_KAPGM_CONTROL
589         &     ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)
590    #else
591       &     ( GM_isopycK - GM_skewflx*GM_background_K       &     ( GM_isopycK - GM_skewflx*GM_background_K
592    #endif
593  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
594       &     +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
595  #endif  #endif
596       &     )*SlopeX(i,j)*taperFct(i,j)       &     )*SlopeX(i,j)*taperFct(i,j)
597           ENDDO           ENDDO
598          ENDDO          ENDDO
599        ENDIF         ENDIF
600  #endif /* GM_EXTRA_DIAGONAL */  #endif /* GM_EXTRA_DIAGONAL */
601    
602  C     Gradient of Sigma at V points  #ifdef ALLOW_DIAGNOSTICS
603        DO j=1-Oly+1,sNy+Oly-1         IF (doDiagRediFlx) THEN
604         DO i=1-Olx+1,sNx+Olx-1          km1 = MAX(k-1,1)
605          dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)          DO j=1,sNy
606       &                    +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k) )           DO i=1,sNx+1
607       &          *_maskS(i,j,k,bi,bj)  C         store in tmp1k Kuz_Redi
608          dSigmaDy(i,j)=sigmaY(i,j,k)            tmp1k(i,j) = ( GM_isopycK
609       &          *_maskS(i,j,k,bi,bj)  #ifdef GM_VISBECK_VARIABLE_K
610          dSigmaDrReal(i,j)=op25*( sigmaR(i,j-1, k ) +sigmaR(i,j, k )       &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0
611       &                  +maskp1*(sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1)) )  #endif
612       &          *_maskS(i,j,k,bi,bj)*maskp1       &                 )*SlopeX(i,j)*taperFct(i,j)
613             ENDDO
614            ENDDO
615            DO j=1,sNy
616             DO i=1,sNx+1
617    C-        Vertical gradients interpolated to U points
618              dTdz = (
619         &     +recip_drC(k)*
620         &       ( maskC(i-1,j,k,bi,bj)*
621         &           (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
622         &        +maskC( i ,j,k,bi,bj)*
623         &           (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
624         &       )
625         &     +recip_drC(kp1)*
626         &       ( maskC(i-1,j,kp1,bi,bj)*
627         &           (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
628         &        +maskC( i ,j,kp1,bi,bj)*
629         &           (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
630         &       )      ) * 0.25 _d 0
631               tmp1k(i,j) = dyG(i,j,bi,bj)*drF(k)
632         &                * _hFacW(i,j,k,bi,bj)
633         &                * tmp1k(i,j) * dTdz
634             ENDDO
635            ENDDO
636            CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KuzTz', k,1,2,bi,bj,myThid)
637           ENDIF
638    #endif /* ALLOW_DIAGNOSTICS */
639    
640    C-- end 2nd  loop on vertical level index k
641          ENDDO
642    
643    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
644    C-- 3rd  k loop : compute Tensor Coeff. at V point
645    
646    #ifdef ALLOW_KPP
647          IF ( useKPP ) THEN
648           DO j=2-Oly,sNy+Oly
649            DO i=1-Olx,sNx+Olx
650             locMixLayer(i,j) = ( KPPhbl(i,j-1,bi,bj)
651         &                      + KPPhbl(i, j ,bi,bj) )*op5
652            ENDDO
653           ENDDO
654          ELSE
655    #else
656          IF ( .TRUE. ) THEN
657    #endif
658           DO j=2-Oly,sNy+Oly
659            DO i=1-Olx,sNx+Olx
660             locMixLayer(i,j) = ( hMixLayer(i,j-1,bi,bj)
661         &                      + hMixLayer(i, j ,bi,bj) )*op5
662            ENDDO
663           ENDDO
664          ENDIF
665          DO j=1-Oly,sNy+Oly
666           DO i=1-Olx,sNx+Olx
667             hTransLay(i,j) =  0.
668             baseSlope(i,j) =  0.
669             recipLambda(i,j)= 0.
670         ENDDO         ENDDO
671        ENDDO        ENDDO
672          DO j=2-Oly,sNy+Oly
673           DO i=1-Olx,sNx+Olx
674             hTransLay(i,j) = MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
675           ENDDO
676          ENDDO
677    
678    C     Gradient of Sigma at V points
679          DO k=Nr,1,-1
680           kp1 = MIN(Nr,k+1)
681           maskp1 = 1. _d 0
682           IF (k.GE.Nr) maskp1 = 0. _d 0
683    
684           DO j=1-Oly+1,sNy+Oly-1
685            DO i=1-Olx+1,sNx+Olx-1
686             dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
687         &                       +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k)
688         &                      )*_maskS(i,j,k,bi,bj)
689             dSigmaDy(i,j)=sigmaY(i,j,k)
690         &                       *_maskS(i,j,k,bi,bj)
691             dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k )+sigmaR(i,j, k )
692         &                      +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1
693         &                      )*_maskS(i,j,k,bi,bj)
694            ENDDO
695           ENDDO
696    
697  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
698  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
699  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
700  CADJ STORE dsigmadrreal(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte
701  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
702    
703  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
704        CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
705       U             dSigmadRReal,       O             SlopeX, SlopeY,
      I             rF(K),K,  
      U             SlopeX, SlopeY,  
      U             dSigmaDx, dSigmaDy,  
706       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
707       I             bi, bj, myThid )       U             hTransLay, baseSlope, recipLambda,
708         U             dSigmaDr,
709         I             dSigmaDx, dSigmaDy,
710         I             ldd97_LrhoS, locMixLayer, rC,
711         I             kLow_S,
712         I             k, bi, bj, myTime, myIter, myThid )
713    
714    cph(
715    #ifdef ALLOW_AUTODIFF_TAMC
716    cph(
717    CADJ STORE taperfct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
718    cph)
719    #endif /* ALLOW_AUTODIFF_TAMC */
720    cph)
721    
722  #ifdef GM_NON_UNITY_DIAGONAL  #ifdef GM_NON_UNITY_DIAGONAL
723    c      IF ( GM_nonUnitDiag ) THEN
724          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
725           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
726            Kvy(i,j,k,bi,bj) =            Kvy(i,j,k,bi,bj) =
# Line 368  C     Calculate slopes for use in tensor Line 728  C     Calculate slopes for use in tensor
728  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
729       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
730  #endif  #endif
731       &     )       &     )*taperFct(i,j)
      &     *taperFct(i,j)  
732           ENDDO           ENDDO
733          ENDDO          ENDDO
734  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
735  # ifndef GM_TAPER_ORIG_CLIPPING  # ifdef GM_EXCLUDE_CLIPPING
736  CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
737  # endif  # endif
738  #endif  #endif
# Line 382  CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_b Line 741  CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_b
741            Kvy(i,j,k,bi,bj) = MAX( Kvy(i,j,k,bi,bj), GM_Kmin_horiz )            Kvy(i,j,k,bi,bj) = MAX( Kvy(i,j,k,bi,bj), GM_Kmin_horiz )
742           ENDDO           ENDDO
743          ENDDO          ENDDO
744    c      ENDIF
745  #endif /* GM_NON_UNITY_DIAGONAL */  #endif /* GM_NON_UNITY_DIAGONAL */
746    
747  #ifdef GM_EXTRA_DIAGONAL  #ifdef GM_EXTRA_DIAGONAL
# Line 390  CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_b Line 750  CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_b
750  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
751  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
752  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
753        IF (GM_ExtraDiag) THEN         IF ( GM_ExtraDiag ) THEN
754          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
755           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
756            Kvz(i,j,k,bi,bj) =            Kvz(i,j,k,bi,bj) =
757    #ifdef ALLOW_KAPGM_CONTROL
758         &     ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)
759    #else
760       &     ( GM_isopycK - GM_skewflx*GM_background_K       &     ( GM_isopycK - GM_skewflx*GM_background_K
761    #endif
762  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
763       &     +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
764  #endif  #endif
765       &     )*SlopeY(i,j)*taperFct(i,j)       &     )*SlopeY(i,j)*taperFct(i,j)
766           ENDDO           ENDDO
767          ENDDO          ENDDO
768        ENDIF         ENDIF
769  #endif /* GM_EXTRA_DIAGONAL */  #endif /* GM_EXTRA_DIAGONAL */
770    
771  #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */  #ifdef ALLOW_DIAGNOSTICS
772           IF (doDiagRediFlx) THEN
773    c       km1 = MAX(k-1,1)
774            DO j=1,sNy+1
775             DO i=1,sNx
776    C         store in tmp1k Kvz_Redi
777              tmp1k(i,j) = ( GM_isopycK
778    #ifdef GM_VISBECK_VARIABLE_K
779         &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0
780    #endif
781         &                 )*SlopeY(i,j)*taperFct(i,j)
782             ENDDO
783            ENDDO
784            DO j=1,sNy+1
785             DO i=1,sNx
786    C-        Vertical gradients interpolated to U points
787              dTdz = (
788         &     +recip_drC(k)*
789         &       ( maskC(i,j-1,k,bi,bj)*
790         &           (theta(i,j-1,km1,bi,bj)-theta(i,j-1,k,bi,bj))
791         &        +maskC(i, j ,k,bi,bj)*
792         &           (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))
793         &       )
794         &     +recip_drC(kp1)*
795         &       ( maskC(i,j-1,kp1,bi,bj)*
796         &           (theta(i,j-1,k,bi,bj)-theta(i,j-1,kp1,bi,bj))
797         &        +maskC(i, j ,kp1,bi,bj)*
798         &           (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))
799         &       )      ) * 0.25 _d 0
800               tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)
801         &                * _hFacS(i,j,k,bi,bj)
802         &                * tmp1k(i,j) * dTdz
803             ENDDO
804            ENDDO
805            CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KvzTz', k,1,2,bi,bj,myThid)
806           ENDIF
807    #endif /* ALLOW_DIAGNOSTICS */
808    
809  #ifdef ALLOW_TIMEAVE  C-- end 3rd  loop on vertical level index k
 C--   Time-average  
       DO j=1-Oly+1,sNy+Oly-1  
        DO i=1-Olx+1,sNx+Olx-1  
         GM_Kwx_T(i,j,k,bi,bj)=GM_Kwx_T(i,j,k,bi,bj)  
      &                       +Kwx(i,j,k,bi,bj)*deltaTclock  
         GM_Kwy_T(i,j,k,bi,bj)=GM_Kwy_T(i,j,k,bi,bj)  
      &                       +Kwy(i,j,k,bi,bj)*deltaTclock  
         GM_Kwz_T(i,j,k,bi,bj)=GM_Kwz_T(i,j,k,bi,bj)  
      &                       +Kwz(i,j,k,bi,bj)*deltaTclock  
        ENDDO  
810        ENDDO        ENDDO
       GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock  
 #endif /* ALLOW_TIMEAVE */  
811    
812  C-- end 2nd loop on vertical level index k  #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
       ENDDO  
813    
814    
815  #ifdef GM_BOLUS_ADVEC  #ifdef GM_BOLUS_ADVEC
816        IF (GM_AdvForm) THEN        IF (GM_AdvForm) THEN
817          CALL GMREDI_CALC_PSI_B(         CALL GMREDI_CALC_PSI_B(
818       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
819       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
820       I             myThid )       I             ldd97_LrhoW, ldd97_LrhoS,
821         I             myThid )
822        ENDIF        ENDIF
823  #endif  #endif
824    
825    #ifdef ALLOW_TIMEAVE
826    C--   Time-average
827          IF ( taveFreq.GT.0. ) THEN
828    
829             CALL TIMEAVE_CUMULATE( GM_Kwx_T, Kwx, Nr,
830         &                          deltaTclock, bi, bj, myThid )
831             CALL TIMEAVE_CUMULATE( GM_Kwy_T, Kwy, Nr,
832         &                          deltaTclock, bi, bj, myThid )
833             CALL TIMEAVE_CUMULATE( GM_Kwz_T, Kwz, Nr,
834         &                          deltaTclock, bi, bj, myThid )
835    #ifdef GM_VISBECK_VARIABLE_K
836           IF ( GM_Visbeck_alpha.NE.0. ) THEN
837             CALL TIMEAVE_CUMULATE( Visbeck_K_T, VisbeckK, 1,
838         &                          deltaTclock, bi, bj, myThid )
839           ENDIF
840    #endif
841    #ifdef GM_BOLUS_ADVEC
842           IF ( GM_AdvForm ) THEN
843             CALL TIMEAVE_CUMULATE( GM_PsiXtave, GM_PsiX, Nr,
844         &                          deltaTclock, bi, bj, myThid )
845             CALL TIMEAVE_CUMULATE( GM_PsiYtave, GM_PsiY, Nr,
846         &                          deltaTclock, bi, bj, myThid )
847           ENDIF
848    #endif
849           DO k=1,Nr
850             GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock
851           ENDDO
852    
853          ENDIF
854    #endif /* ALLOW_TIMEAVE */
855    
856    #ifdef ALLOW_DIAGNOSTICS
857          IF ( useDiagnostics ) THEN
858            CALL GMREDI_DIAGNOSTICS_FILL(bi,bj,myThid)
859          ENDIF
860    #endif /* ALLOW_DIAGNOSTICS */
861    
862  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
863    
864        RETURN        RETURN
# Line 441  C-- end 2nd loop on vertical level index Line 866  C-- end 2nd loop on vertical level index
866    
867    
868        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
869       I             bi, bj, iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
870       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
871       I             myThid )       I             bi, bj, myTime, myIter, myThid )
872  C     /==========================================================\  C     /==========================================================\
873  C     | SUBROUTINE GMREDI_CALC_TENSOR                            |  C     | SUBROUTINE GMREDI_CALC_TENSOR                            |
874  C     | o Calculate tensor elements for GM/Redi tensor.          |  C     | o Calculate tensor elements for GM/Redi tensor.          |
# Line 453  C     \================================= Line 878  C     \=================================
878    
879  C     == Global variables ==  C     == Global variables ==
880  #include "SIZE.h"  #include "SIZE.h"
 #include "GRID.h"  
 #include "DYNVARS.h"  
881  #include "EEPARAMS.h"  #include "EEPARAMS.h"
 #include "PARAMS.h"  
882  #include "GMREDI.h"  #include "GMREDI.h"
883    
884  C     == Routine arguments ==  C     == Routine arguments ==
# Line 464  C Line 886  C
886        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
887        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
888        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
889        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
890          INTEGER bi, bj
891          _RL     myTime
892          INTEGER myIter
893        INTEGER myThid        INTEGER myThid
894  CEndOfInterface  CEndOfInterface
895    
       INTEGER i, j, k  
   
896  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
897    
898          INTEGER i, j, k
899    
900        DO k=1,Nr        DO k=1,Nr
901         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
902          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
Legend:
Removed from v.1.14  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22