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