/[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.13 by heimbach, Thu Nov 28 17:30:34 2002 UTC revision 1.31 by gforget, Sat Feb 2 02:35:53 2008 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  #ifdef GM_VISBECK_VARIABLE_K        _RL ldd97_LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
78        _RS deltaH,zero_rs        _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
79        PARAMETER(zero_rs=0.)        _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
89    #ifdef OLD_VISBECK_CALC
90          _RL deltaH,zero_rs
91          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        DO j=1-Oly,sNy+Oly
206         IF (k.LE.1) maskm1 = 0. _d 0         DO i=1-Olx,sNx+Olx
207             hTransLay(i,j) = R_low(i,j,bi,bj)
208             baseSlope(i,j) =  0. _d 0
209             recipLambda(i,j) = 0. _d 0
210             locMixLayer(i,j) = 0. _d 0
211           ENDDO
212          ENDDO
213    #ifdef ALLOW_KPP
214          IF ( useKPP ) THEN
215           DO j=1-Oly,sNy+Oly
216            DO i=1-Olx,sNx+Olx
217             locMixLayer(i,j) = KPPhbl(i,j,bi,bj)
218            ENDDO
219           ENDDO
220          ELSE
221    #else
222          IF ( .TRUE. ) THEN
223    #endif
224           DO j=1-Oly,sNy+Oly
225            DO i=1-Olx,sNx+Olx
226             locMixLayer(i,j) = hMixLayer(i,j,bi,bj)
227            ENDDO
228           ENDDO
229          ENDIF
230    
231          DO k=Nr,2,-1
232    
233  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
234         kkey = (igmkey-1)*Nr + k         kkey = (igmkey-1)*Nr + k
# Line 93  C-- 1rst loop on k : compute Tensor Coef Line 238  C-- 1rst loop on k : compute Tensor Coef
238           SlopeY(i,j)       = 0. _d 0           SlopeY(i,j)       = 0. _d 0
239           dSigmaDx(i,j)     = 0. _d 0           dSigmaDx(i,j)     = 0. _d 0
240           dSigmaDy(i,j)     = 0. _d 0           dSigmaDy(i,j)     = 0. _d 0
241           dSigmaDrReal(i,j) = 0. _d 0           dSigmaDr(i,j)     = 0. _d 0
242           SlopeSqr(i,j)     = 0. _d 0           SlopeSqr(i,j)     = 0. _d 0
243           taperFct(i,j)     = 0. _d 0           taperFct(i,j)     = 0. _d 0
244           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 260  C-- 1rst loop on k : compute Tensor Coef
260         ENDDO         ENDDO
261  #endif  #endif
262    
263        DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
264         DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
265  C      Gradient of Sigma at rVel points  C      Gradient of Sigma at rVel points
266          dSigmaDx(i,j)=0.25*( sigmaX(i+1, j ,km1) +sigmaX(i,j,km1)           dSigmaDx(i,j)=op25*( sigmaX(i+1,j,k-1)+sigmaX(i,j,k-1)
267       &                    +sigmaX(i+1, j , k ) +sigmaX(i,j, k ) )       &                       +sigmaX(i+1,j, k )+sigmaX(i,j, k )
268       &                  *maskC(i,j,k,bi,bj)*maskm1       &                      )*maskC(i,j,k,bi,bj)
269          dSigmaDy(i,j)=0.25*( sigmaY( i ,j+1,km1) +sigmaY(i,j,km1)           dSigmaDy(i,j)=op25*( sigmaY(i,j+1,k-1)+sigmaY(i,j,k-1)
270       &                    +sigmaY( i ,j+1, k ) +sigmaY(i,j, k ) )       &                       +sigmaY(i,j+1, k )+sigmaY(i,j, k )
271       &                  *maskC(i,j,k,bi,bj)*maskm1       &                      )*maskC(i,j,k,bi,bj)
272          dSigmaDrReal(i,j)=sigmaR(i,j,k)*maskm1           dSigmaDr(i,j)=sigmaR(i,j,k)
273            ENDDO
274         ENDDO         ENDDO
       ENDDO  
275    
276  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
277  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
278  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
279  CADJ STORE dsigmadrreal(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
280    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
281    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
282    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
283  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
284    
285  C     Calculate slopes for use in tensor, taper and/or clip  #ifdef GM_VISBECK_VARIABLE_K
286        CALL GMREDI_SLOPE_LIMIT(  #ifndef OLD_VISBECK_CALC
287       U             dSigmadRReal,         IF ( GM_Visbeck_alpha.GT.0. .AND.
288       I             rF(K),K,       &      -rC(k-1).LT.GM_Visbeck_depth ) THEN
289       U             SlopeX, SlopeY,  
290       U             dSigmaDx, dSigmaDy,  C--     Depth average of f/sqrt(Ri) = M^2/N^2 * N
291       O             SlopeSqr, taperFct,  C       M^2 and N^2 are horizontal & vertical gradient of buoyancy.
292       I             bi, bj, myThid )  
293    C       Calculate terms for mean Richardson number which is used
294    C       in the "variable K" parameterisaton:
295    C       compute depth average from surface down to the bottom or
296    C       GM_Visbeck_depth, whatever is the shallower.
297    
298        DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
299         DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
300              IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
301               integrDepth = -rC( kLowC(i,j,bi,bj) )
302    C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments
303               integrDepth = MIN( integrDepth, GM_Visbeck_depth )
304    C       Distance between level center above and the integration depth
305               deltaH = integrDepth + rC(k-1)
306    C       If negative then we are below the integration level
307    C       (cannot be the case with 2 conditions on maskC & -rC(k-1))
308    C       If positive we limit this to the distance from center above
309               deltaH = MIN( deltaH, drC(k) )
310    C       Now we convert deltaH to a non-dimensional fraction
311               deltaH = deltaH/( integrDepth+rC(1) )
312    
313  C       Mask Iso-neutral slopes  C--      compute: ( M^2 * S )^1/2   (= M^2 / N since S=M^2/N^2 )
314          SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)*maskm1             dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)
315          SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)*maskm1       &             + dSigmaDy(i,j)*dSigmaDy(i,j)
316          SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)*maskm1             IF ( dSigmaH .GT. 0. _d 0 ) THEN
317                 dSigmaH = SQRT( dSigmaH )
318    C-       compute slope, limited by GM_maxSlope:
319                 IF ( -dSigmaDr(i,j).GT.dSigmaH*GM_rMaxSlope ) THEN
320                  Sloc = dSigmaH / ( -dSigmaDr(i,j) )
321                 ELSE
322                  Sloc = GM_maxSlope
323                 ENDIF
324                 M2loc = Gravity*recip_RhoConst*dSigmaH
325                 SNloc = SQRT( Sloc*M2loc )
326               ELSE
327                 SNloc = 0. _d 0
328               ENDIF
329               VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
330         &       +deltaH*GM_Visbeck_alpha
331         &              *GM_Visbeck_length*GM_Visbeck_length*SNloc
332              ENDIF
333             ENDDO
334            ENDDO
335           ENDIF
336    #endif /* ndef OLD_VISBECK_CALC */
337    #endif /* GM_VISBECK_VARIABLE_K */
338    
339    C     Calculate slopes for use in tensor, taper and/or clip
340           CALL GMREDI_SLOPE_LIMIT(
341         O             SlopeX, SlopeY,
342         O             SlopeSqr, taperFct,
343         U             hTransLay, baseSlope, recipLambda,
344         U             dSigmaDr,
345         I             dSigmaDx, dSigmaDy,
346         I             ldd97_LrhoC, locMixLayer, rF,
347         I             kLowC(1-Olx,1-Oly,bi,bj),
348         I             k, bi, bj, myTime, myIter, myThid )
349    
350           DO j=1-Oly+1,sNy+Oly-1
351            DO i=1-Olx+1,sNx+Olx-1
352    C      Mask Iso-neutral slopes
353             SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)
354             SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)
355             SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)
356            ENDDO
357         ENDDO         ENDDO
       ENDDO  
358    
359  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
360  CADJ STORE SlopeX(:,:)          = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
361  CADJ STORE SlopeY(:,:)          = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
362  CADJ STORE SlopeSqr(:,:)        = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeSqr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
363  CADJ STORE taperFct(:,:)        = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
364  #ifdef GM_VISBECK_VARIABLE_K  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
 CADJ STORE dSigmaDrReal(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif  
365  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
366    
367        DO j=1-Oly+1,sNy+Oly-1  C      Components of Redi/GM tensor
368         DO i=1-Olx+1,sNx+Olx-1         DO j=1-Oly+1,sNy+Oly-1
369            DO i=1-Olx+1,sNx+Olx-1
370  C       Components of Redi/GM tensor            Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
371          Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)            Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
372          Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)            Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
373          Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)          ENDDO
374           ENDDO
375    
376  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
377    #ifdef OLD_VISBECK_CALC
378           DO j=1-Oly+1,sNy+Oly-1
379            DO i=1-Olx+1,sNx+Olx-1
380    
381  C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K  C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K
382  C           but don't know if *taperFct (or **2 ?) is necessary  C           but do not know if *taperFct (or **2 ?) is necessary
383          Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)          Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)
384    
385  C--     Depth average of M^2/N^2 * N  C--     Depth average of M^2/N^2 * N
# Line 192  C       Now we convert deltaH to a non-d Line 396  C       Now we convert deltaH to a non-d
396          deltaH=deltaH/GM_Visbeck_depth          deltaH=deltaH/GM_Visbeck_depth
397    
398          IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.          IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.
399          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
400           N2= -Gravity*recip_RhoConst*dSigmaDrReal(i,j)           N2= -Gravity*recip_RhoConst*dSigmaDr(i,j)
401           SN=sqrt(Ssq(i,j)*N2)           SN=sqrt(Ssq(i,j)*N2)
402           VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH           VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH
403       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN
404          ENDIF          ENDIF
405    
406  #endif /* GM_VISBECK_VARIABLE_K */          ENDDO
   
407         ENDDO         ENDDO
408        ENDDO  #endif /* OLD_VISBECK_CALC */
409    #endif /* GM_VISBECK_VARIABLE_K */
410    
411  C-- end 1rst loop on vertical level index k  C-- end 1rst loop on vertical level index k
412        ENDDO        ENDDO
# Line 212  C-- end 1rst loop on vertical level inde Line 416  C-- end 1rst loop on vertical level inde
416  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
417  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
418  #endif  #endif
419        IF ( GM_Visbeck_alpha.NE.0. ) THEN        IF ( GM_Visbeck_alpha.GT.0. ) THEN
420  C-     Limit range that KapGM can take  C-     Limit range that KapGM can take
421         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
422          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
423           VisbeckK(i,j,bi,bj)=           VisbeckK(i,j,bi,bj)=
424       &       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  
425          ENDDO          ENDDO
426         ENDDO         ENDDO
427        ENDIF        ENDIF
428    cph( NEW
429    #ifdef ALLOW_AUTODIFF_TAMC
430    CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
431    #endif
432    cph)
433  #endif /* GM_VISBECK_VARIABLE_K */  #endif /* GM_VISBECK_VARIABLE_K */
434    
435    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.  
436        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  
   
437  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
438         kkey = (igmkey-1)*Nr + k         kkey = (igmkey-1)*Nr + k
439  #ifdef GM_VISBECK_VARIABLE_K  # if (defined (GM_NON_UNITY_DIAGONAL) || \
440  CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj, key=kkey, byte=isbyte       defined (GM_VISBECK_VARIABLE_K))
441  CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj, key=kkey, byte=isbyte  CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
442  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj, key=kkey, byte=isbyte  CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
443    CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
444    # endif
445  #endif  #endif
446           DO j=1-Oly+1,sNy+Oly-1
447            DO i=1-Olx+1,sNx+Olx-1
448    #ifdef ALLOW_KAPREDI_CONTROL
449             Kgm_tmp = kapredi(i,j,k,bi,bj)
450    #else
451             Kgm_tmp = GM_isopycK
452    #endif
453    #ifdef ALLOW_KAPGM_CONTROL
454         &           + GM_skewflx*kapgm(i,j,k,bi,bj)
455    #else
456         &           + GM_skewflx*GM_background_K
457  #endif  #endif
   
 C-    express the Tensor in term of Diffusivity (= m**2 / s )  
       DO j=1-Oly+1,sNy+Oly-1  
        DO i=1-Olx+1,sNx+Olx-1  
         Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K  
458  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
459       &          + VisbeckK(i,j,bi,bj)*(1.+GM_skewflx)           &           + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)
460    #endif
461             Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)
462             Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)
463    #ifdef ALLOW_KAPREDI_CONTROL
464             Kwz(i,j,k,bi,bj)= ( kapredi(i,j,k,bi,bj)
465    #else
466             Kwz(i,j,k,bi,bj)= ( GM_isopycK
467  #endif  #endif
         Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)  
         Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)  
         Kwz(i,j,k,bi,bj)= ( GM_isopycK  
468  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
469       &                    + VisbeckK(i,j,bi,bj)       &                     + VisbeckK(i,j,bi,bj)
470  #endif  #endif
471       &                    )*Kwz(i,j,k,bi,bj)       &                     )*Kwz(i,j,k,bi,bj)
472            ENDDO
473         ENDDO         ENDDO
474        ENDDO        ENDDO
475  #ifdef ALLOW_AUTODIFF_TAMC  
476  #ifdef GM_VISBECK_VARIABLE_K  #ifdef ALLOW_DIAGNOSTICS
477  CADJ STORE VisbeckK(:,:,bi,bj) =        IF ( useDiagnostics .AND. GM_taper_scheme.EQ.'fm07' ) THEN
478  CADJ &     comlev1_bibj, key=kkey, byte=isbyte         CALL DIAGNOSTICS_FILL( hTransLay, 'GM_hTrsL', 0,1,2,bi,bj,myThid)
479  #endif         CALL DIAGNOSTICS_FILL( baseSlope, 'GM_baseS', 0,1,2,bi,bj,myThid)
480  #endif         CALL DIAGNOSTICS_FILL(recipLambda,'GM_rLamb', 0,1,2,bi,bj,myThid)
481          ENDIF
482    #endif /* ALLOW_DIAGNOSTICS */
483    
484    
485  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
486    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
487    C-- 2nd  k loop : compute Tensor Coeff. at U point
488    
489  C     Gradient of Sigma at U points  #ifdef ALLOW_KPP
490        DO j=1-Oly+1,sNy+Oly-1        IF ( useKPP ) THEN
491         DO i=1-Olx+1,sNx+Olx-1         DO j=1-Oly,sNy+Oly
492          dSigmaDx(i,j)=sigmaX(i,j,k)          DO i=2-Olx,sNx+Olx
493       &          *_maskW(i,j,k,bi,bj)           locMixLayer(i,j) = ( KPPhbl(i-1,j,bi,bj)
494          dSigmaDy(i,j)=0.25*( sigmaY(i-1,j+1,k) +sigmaY(i,j+1,k)       &                      + KPPhbl( i ,j,bi,bj) )*op5
495       &                      +sigmaY(i-1, j ,k) +sigmaY(i, j ,k) )          ENDDO
496       &          *_maskW(i,j,k,bi,bj)         ENDDO
497          dSigmaDrReal(i,j)=0.25*( sigmaR(i-1,j, k ) +sigmaR(i,j, k )        ELSE
498       &                  +maskp1*(sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1)) )  #else
499       &          *_maskW(i,j,k,bi,bj)        IF ( .TRUE. ) THEN
500    #endif
501           DO j=1-Oly,sNy+Oly
502            DO i=2-Olx,sNx+Olx
503             locMixLayer(i,j) = ( hMixLayer(i-1,j,bi,bj)
504         &                      + hMixLayer( i ,j,bi,bj) )*op5
505            ENDDO
506           ENDDO
507          ENDIF
508          DO j=1-Oly,sNy+Oly
509           DO i=1-Olx,sNx+Olx
510             hTransLay(i,j) =  0.
511             baseSlope(i,j) =  0.
512             recipLambda(i,j)= 0.
513           ENDDO
514           DO i=2-Olx,sNx+Olx
515             hTransLay(i,j) = MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
516         ENDDO         ENDDO
517        ENDDO        ENDDO
518    
519          DO k=Nr,1,-1
520           kp1 = MIN(Nr,k+1)
521           maskp1 = 1. _d 0
522           IF (k.GE.Nr) maskp1 = 0. _d 0
523    #ifdef ALLOW_AUTODIFF_TAMC
524           kkey = (igmkey-1)*Nr + k
525    #endif
526    
527    C     Gradient of Sigma at U points
528           DO j=1-Oly+1,sNy+Oly-1
529            DO i=1-Olx+1,sNx+Olx-1
530             dSigmaDx(i,j)=sigmaX(i,j,k)
531         &                       *_maskW(i,j,k,bi,bj)
532             dSigmaDy(i,j)=op25*( sigmaY(i-1,j+1,k)+sigmaY(i,j+1,k)
533         &                       +sigmaY(i-1, j ,k)+sigmaY(i, j ,k)
534         &                      )*_maskW(i,j,k,bi,bj)
535             dSigmaDr(i,j)=op25*( sigmaR(i-1,j, k )+sigmaR(i,j, k )
536         &                      +(sigmaR(i-1,j,kp1)+sigmaR(i,j,kp1))*maskp1
537         &                      )*_maskW(i,j,k,bi,bj)
538            ENDDO
539           ENDDO
540    
541  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
542    CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
543  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
544  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
545  CADJ STORE dsigmadrreal(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
546    CADJ STORE locMixLayer(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
547    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
548    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
549    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
550  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
551    
552  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
553        CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
554       U             dSigmadRReal,       O             SlopeX, SlopeY,
      I             rF(K),K,  
      U             SlopeX, SlopeY,  
      U             dSigmaDx, dSigmaDy,  
555       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
556       I             bi, bj, myThid )       U             hTransLay, baseSlope, recipLambda,
557         U             dSigmaDr,
558         I             dSigmaDx, dSigmaDy,
559         I             ldd97_LrhoW, locMixLayer, rC,
560         I             kLow_W,
561         I             k, bi, bj, myTime, myIter, myThid )
562    
563    #ifdef ALLOW_AUTODIFF_TAMC
564    CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
565    CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
566    #endif /* ALLOW_AUTODIFF_TAMC */
567    
568  #ifdef GM_NON_UNITY_DIAGONAL  #ifdef GM_NON_UNITY_DIAGONAL
569    c      IF ( GM_nonUnitDiag ) THEN
570          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
571           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
572            Kux(i,j,k,bi,bj) =            Kux(i,j,k,bi,bj) =
573    #ifdef ALLOW_KAPREDI_CONTROL
574         &     ( kapredi(i,j,k,bi,bj)
575    #else
576       &     ( GM_isopycK       &     ( GM_isopycK
577    #endif
578  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
579       &     +0.5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
580  #endif  #endif
581       &     )       &     )*taperFct(i,j)
      &     *taperFct(i,j)  
582           ENDDO           ENDDO
583          ENDDO          ENDDO
584  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
585  # ifndef GM_TAPER_ORIG_CLIPPING  # ifdef GM_EXCLUDE_CLIPPING
586  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
587  # endif  # endif
588  #endif  #endif
# Line 321  CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_b Line 591  CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_b
591            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 )
592           ENDDO           ENDDO
593          ENDDO          ENDDO
594    c      ENDIF
595  #endif /* GM_NON_UNITY_DIAGONAL */  #endif /* GM_NON_UNITY_DIAGONAL */
596    
597  #ifdef GM_EXTRA_DIAGONAL  #ifdef GM_EXTRA_DIAGONAL
# Line 329  CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_b Line 600  CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_b
600  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
601  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
602  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
603        IF (GM_ExtraDiag) THEN         IF ( GM_ExtraDiag ) THEN
604          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
605           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
606            Kuz(i,j,k,bi,bj) =            Kuz(i,j,k,bi,bj) =
607       &     ( GM_isopycK - GM_skewflx*GM_background_K  #ifdef ALLOW_KAPREDI_CONTROL
608         &     ( kapredi(i,j,k,bi,bj)
609    #else
610         &     ( GM_isopycK
611    #endif
612    #ifdef ALLOW_KAPGM_CONTROL
613         &     - GM_skewflx*kapgm(i,j,k,bi,bj)
614    #else
615         &     - GM_skewflx*GM_background_K
616    #endif
617  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
618       &     +0.5*(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
619  #endif  #endif
620       &     )*SlopeX(i,j)*taperFct(i,j)       &     )*SlopeX(i,j)*taperFct(i,j)
621           ENDDO           ENDDO
622          ENDDO          ENDDO
623        ENDIF         ENDIF
624  #endif /* GM_EXTRA_DIAGONAL */  #endif /* GM_EXTRA_DIAGONAL */
625    
626  C     Gradient of Sigma at V points  #ifdef ALLOW_DIAGNOSTICS
627        DO j=1-Oly+1,sNy+Oly-1         IF (doDiagRediFlx) THEN
628         DO i=1-Olx+1,sNx+Olx-1          km1 = MAX(k-1,1)
629          dSigmaDx(i,j)=0.25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)          DO j=1,sNy
630       &                    +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k) )           DO i=1,sNx+1
631       &          *_maskS(i,j,k,bi,bj)  C         store in tmp1k Kuz_Redi
632          dSigmaDy(i,j)=sigmaY(i,j,k)  #ifdef ALLOW_KAPREDI_CONTROL
633       &          *_maskS(i,j,k,bi,bj)            tmp1k(i,j) = ( kapredi(i,j,k,bi,bj)
634          dSigmaDrReal(i,j)=0.25*( sigmaR(i,j-1, k ) +sigmaR(i,j, k )  #else
635       &                  +maskp1*(sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1)) )            tmp1k(i,j) = ( GM_isopycK
636       &          *_maskS(i,j,k,bi,bj)  #endif
637    #ifdef GM_VISBECK_VARIABLE_K
638         &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0
639    #endif
640         &                 )*SlopeX(i,j)*taperFct(i,j)
641             ENDDO
642            ENDDO
643            DO j=1,sNy
644             DO i=1,sNx+1
645    C-        Vertical gradients interpolated to U points
646              dTdz = (
647         &     +recip_drC(k)*
648         &       ( maskC(i-1,j,k,bi,bj)*
649         &           (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
650         &        +maskC( i ,j,k,bi,bj)*
651         &           (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
652         &       )
653         &     +recip_drC(kp1)*
654         &       ( maskC(i-1,j,kp1,bi,bj)*
655         &           (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
656         &        +maskC( i ,j,kp1,bi,bj)*
657         &           (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
658         &       )      ) * 0.25 _d 0
659               tmp1k(i,j) = dyG(i,j,bi,bj)*drF(k)
660         &                * _hFacW(i,j,k,bi,bj)
661         &                * tmp1k(i,j) * dTdz
662             ENDDO
663            ENDDO
664            CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KuzTz', k,1,2,bi,bj,myThid)
665           ENDIF
666    #endif /* ALLOW_DIAGNOSTICS */
667    
668    C-- end 2nd  loop on vertical level index k
669          ENDDO
670    
671    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
672    C-- 3rd  k loop : compute Tensor Coeff. at V point
673    
674    #ifdef ALLOW_KPP
675          IF ( useKPP ) THEN
676           DO j=2-Oly,sNy+Oly
677            DO i=1-Olx,sNx+Olx
678             locMixLayer(i,j) = ( KPPhbl(i,j-1,bi,bj)
679         &                      + KPPhbl(i, j ,bi,bj) )*op5
680            ENDDO
681           ENDDO
682          ELSE
683    #else
684          IF ( .TRUE. ) THEN
685    #endif
686           DO j=2-Oly,sNy+Oly
687            DO i=1-Olx,sNx+Olx
688             locMixLayer(i,j) = ( hMixLayer(i,j-1,bi,bj)
689         &                      + hMixLayer(i, j ,bi,bj) )*op5
690            ENDDO
691           ENDDO
692          ENDIF
693          DO j=1-Oly,sNy+Oly
694           DO i=1-Olx,sNx+Olx
695             hTransLay(i,j) =  0.
696             baseSlope(i,j) =  0.
697             recipLambda(i,j)= 0.
698           ENDDO
699          ENDDO
700          DO j=2-Oly,sNy+Oly
701           DO i=1-Olx,sNx+Olx
702             hTransLay(i,j) = MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
703         ENDDO         ENDDO
704        ENDDO        ENDDO
705    
706    C     Gradient of Sigma at V points
707          DO k=Nr,1,-1
708           kp1 = MIN(Nr,k+1)
709           maskp1 = 1. _d 0
710           IF (k.GE.Nr) maskp1 = 0. _d 0
711    #ifdef ALLOW_AUTODIFF_TAMC
712           kkey = (igmkey-1)*Nr + k
713    #endif
714    
715           DO j=1-Oly+1,sNy+Oly-1
716            DO i=1-Olx+1,sNx+Olx-1
717             dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
718         &                       +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k)
719         &                      )*_maskS(i,j,k,bi,bj)
720             dSigmaDy(i,j)=sigmaY(i,j,k)
721         &                       *_maskS(i,j,k,bi,bj)
722             dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k )+sigmaR(i,j, k )
723         &                      +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1
724         &                      )*_maskS(i,j,k,bi,bj)
725            ENDDO
726           ENDDO
727    
728  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
729  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
730  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
731  CADJ STORE dsigmadrreal(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
732    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
733    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
734    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
735  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
736    
737  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
738        CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
739       U             dSigmadRReal,       O             SlopeX, SlopeY,
      I             rF(K),K,  
      U             SlopeX, SlopeY,  
      U             dSigmaDx, dSigmaDy,  
740       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
741       I             bi, bj, myThid )       U             hTransLay, baseSlope, recipLambda,
742         U             dSigmaDr,
743         I             dSigmaDx, dSigmaDy,
744         I             ldd97_LrhoS, locMixLayer, rC,
745         I             kLow_S,
746         I             k, bi, bj, myTime, myIter, myThid )
747    
748    cph(
749    #ifdef ALLOW_AUTODIFF_TAMC
750    cph(
751    CADJ STORE taperfct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
752    cph)
753    #endif /* ALLOW_AUTODIFF_TAMC */
754    cph)
755    
756  #ifdef GM_NON_UNITY_DIAGONAL  #ifdef GM_NON_UNITY_DIAGONAL
757    c      IF ( GM_nonUnitDiag ) THEN
758          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
759           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
760            Kvy(i,j,k,bi,bj) =            Kvy(i,j,k,bi,bj) =
761    #ifdef ALLOW_KAPREDI_CONTROL
762         &     ( kapredi(i,j,k,bi,bj)
763    #else
764       &     ( GM_isopycK       &     ( GM_isopycK
765    #endif
766  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
767       &     +0.5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
768  #endif  #endif
769       &     )       &     )*taperFct(i,j)
      &     *taperFct(i,j)  
770           ENDDO           ENDDO
771          ENDDO          ENDDO
772  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
773  # ifndef GM_TAPER_ORIG_CLIPPING  # ifdef GM_EXCLUDE_CLIPPING
774  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
775  # endif  # endif
776  #endif  #endif
# Line 394  CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_b Line 779  CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_b
779            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 )
780           ENDDO           ENDDO
781          ENDDO          ENDDO
782    c      ENDIF
783  #endif /* GM_NON_UNITY_DIAGONAL */  #endif /* GM_NON_UNITY_DIAGONAL */
784    
785  #ifdef GM_EXTRA_DIAGONAL  #ifdef GM_EXTRA_DIAGONAL
# Line 402  CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_b Line 788  CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_b
788  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
789  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
790  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
791        IF (GM_ExtraDiag) THEN         IF ( GM_ExtraDiag ) THEN
792          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
793           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
794            Kvz(i,j,k,bi,bj) =            Kvz(i,j,k,bi,bj) =
795       &     ( GM_isopycK - GM_skewflx*GM_background_K  #ifdef ALLOW_KAPREDI_CONTROL
796         &     ( kapredi(i,j,k,bi,bj)
797    #else
798         &     ( GM_isopycK
799    #endif
800    #ifdef ALLOW_KAPGM_CONTROL
801         &     - GM_skewflx*kapgm(i,j,k,bi,bj)
802    #else
803         &     - GM_skewflx*GM_background_K
804    #endif
805  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
806       &     +0.5*(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
807  #endif  #endif
808       &     )*SlopeY(i,j)*taperFct(i,j)       &     )*SlopeY(i,j)*taperFct(i,j)
809           ENDDO           ENDDO
810          ENDDO          ENDDO
811        ENDIF         ENDIF
812  #endif /* GM_EXTRA_DIAGONAL */  #endif /* GM_EXTRA_DIAGONAL */
813    
814  #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */  #ifdef ALLOW_DIAGNOSTICS
815           IF (doDiagRediFlx) THEN
816            km1 = MAX(k-1,1)
817            DO j=1,sNy+1
818             DO i=1,sNx
819    C         store in tmp1k Kvz_Redi
820    #ifdef ALLOW_KAPREDI_CONTROL
821              tmp1k(i,j) = ( kapredi(i,j,k,bi,bj)
822    #else
823              tmp1k(i,j) = ( GM_isopycK
824    #endif
825    #ifdef GM_VISBECK_VARIABLE_K
826         &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0
827    #endif
828         &                 )*SlopeY(i,j)*taperFct(i,j)
829             ENDDO
830            ENDDO
831            DO j=1,sNy+1
832             DO i=1,sNx
833    C-        Vertical gradients interpolated to U points
834              dTdz = (
835         &     +recip_drC(k)*
836         &       ( maskC(i,j-1,k,bi,bj)*
837         &           (theta(i,j-1,km1,bi,bj)-theta(i,j-1,k,bi,bj))
838         &        +maskC(i, j ,k,bi,bj)*
839         &           (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))
840         &       )
841         &     +recip_drC(kp1)*
842         &       ( maskC(i,j-1,kp1,bi,bj)*
843         &           (theta(i,j-1,k,bi,bj)-theta(i,j-1,kp1,bi,bj))
844         &        +maskC(i, j ,kp1,bi,bj)*
845         &           (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))
846         &       )      ) * 0.25 _d 0
847               tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)
848         &                * _hFacS(i,j,k,bi,bj)
849         &                * tmp1k(i,j) * dTdz
850             ENDDO
851            ENDDO
852            CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KvzTz', k,1,2,bi,bj,myThid)
853           ENDIF
854    #endif /* ALLOW_DIAGNOSTICS */
855    
856  #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  
857        ENDDO        ENDDO
       GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock  
 #endif /* ALLOW_TIMEAVE */  
858    
859  C-- end 2nd loop on vertical level index k  #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
       ENDDO  
860    
861    
862  #ifdef GM_BOLUS_ADVEC  #ifdef GM_BOLUS_ADVEC
863        IF (GM_AdvForm) THEN        IF (GM_AdvForm) THEN
864          CALL GMREDI_CALC_PSI_B(         CALL GMREDI_CALC_PSI_B(
865       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
866       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
867       I             myThid )       I             ldd97_LrhoW, ldd97_LrhoS,
868         I             myThid )
869        ENDIF        ENDIF
870  #endif  #endif
871    
872    #ifdef ALLOW_TIMEAVE
873    C--   Time-average
874          IF ( taveFreq.GT.0. ) THEN
875    
876             CALL TIMEAVE_CUMULATE( GM_Kwx_T, Kwx, Nr,
877         &                          deltaTclock, bi, bj, myThid )
878             CALL TIMEAVE_CUMULATE( GM_Kwy_T, Kwy, Nr,
879         &                          deltaTclock, bi, bj, myThid )
880             CALL TIMEAVE_CUMULATE( GM_Kwz_T, Kwz, Nr,
881         &                          deltaTclock, bi, bj, myThid )
882    #ifdef GM_VISBECK_VARIABLE_K
883           IF ( GM_Visbeck_alpha.NE.0. ) THEN
884             CALL TIMEAVE_CUMULATE( Visbeck_K_T, VisbeckK, 1,
885         &                          deltaTclock, bi, bj, myThid )
886           ENDIF
887    #endif
888    #ifdef GM_BOLUS_ADVEC
889           IF ( GM_AdvForm ) THEN
890             CALL TIMEAVE_CUMULATE( GM_PsiXtave, GM_PsiX, Nr,
891         &                          deltaTclock, bi, bj, myThid )
892             CALL TIMEAVE_CUMULATE( GM_PsiYtave, GM_PsiY, Nr,
893         &                          deltaTclock, bi, bj, myThid )
894           ENDIF
895    #endif
896           DO k=1,Nr
897             GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock
898           ENDDO
899    
900          ENDIF
901    #endif /* ALLOW_TIMEAVE */
902    
903    #ifdef ALLOW_DIAGNOSTICS
904          IF ( useDiagnostics ) THEN
905            CALL GMREDI_DIAGNOSTICS_FILL(bi,bj,myThid)
906          ENDIF
907    #endif /* ALLOW_DIAGNOSTICS */
908    
909  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
910    
911        RETURN        RETURN
# Line 453  C-- end 2nd loop on vertical level index Line 913  C-- end 2nd loop on vertical level index
913    
914    
915        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
916       I             bi, bj, iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
917       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
918       I             myThid )       I             bi, bj, myTime, myIter, myThid )
919  C     /==========================================================\  C     /==========================================================\
920  C     | SUBROUTINE GMREDI_CALC_TENSOR                            |  C     | SUBROUTINE GMREDI_CALC_TENSOR                            |
921  C     | o Calculate tensor elements for GM/Redi tensor.          |  C     | o Calculate tensor elements for GM/Redi tensor.          |
# Line 465  C     \================================= Line 925  C     \=================================
925    
926  C     == Global variables ==  C     == Global variables ==
927  #include "SIZE.h"  #include "SIZE.h"
 #include "GRID.h"  
 #include "DYNVARS.h"  
928  #include "EEPARAMS.h"  #include "EEPARAMS.h"
 #include "PARAMS.h"  
929  #include "GMREDI.h"  #include "GMREDI.h"
930    
931  C     == Routine arguments ==  C     == Routine arguments ==
# Line 476  C Line 933  C
933        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
934        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
935        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
936        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
937          INTEGER bi, bj
938          _RL     myTime
939          INTEGER myIter
940        INTEGER myThid        INTEGER myThid
941  CEndOfInterface  CEndOfInterface
942    
       INTEGER i, j, k  
   
943  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
944    
945          INTEGER i, j, k
946    
947        DO k=1,Nr        DO k=1,Nr
948         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
949          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
Legend:
Removed from v.1.13  
changed lines
  Added in v.1.31
  ViewVC Help
Powered by ViewVC 1.1.22