/[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.3 by heimbach, Mon Jan 29 20:07:39 2001 UTC revision 1.28 by heimbach, Tue Jun 26 15:34:15 2007 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    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, K,       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 21  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
42    #include "tamc.h"
43    #include "tamc_keys.h"
44    #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,K        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,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 dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
71        _RL dRdSigmaLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
72        _RL Ssq        _RL dSigmaDr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
73          _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
74          _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
75          _RL maskp1, Kgm_tmp
76          _RL ldd97_LrhoC(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
77          _RL ldd97_LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
78          _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
79          _RL Cspd, LrhoInf, LrhoSup, fCoriLoc
80    
81          INTEGER kLow_W (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
82          INTEGER kLow_S (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
83          _RL locMixLayer(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
84          _RL baseSlope  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
85          _RL hTransLay  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
86          _RL recipLambda(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
87    
88  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
89        _RS deltaH,zero_rs  #ifdef OLD_VISBECK_CALC
90        PARAMETER(zero_rs=0.)        _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)
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        km1=max(1,K-1)  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
       kp1=min(Nr,K)  
   
111    
112  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
113  !HPF$ INDEPENDENT            act1 = bi - myBxLo(myThid)
114              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
115              act2 = bj - myByLo(myThid)
116              max2 = myByHi(myThid) - myByLo(myThid) + 1
117              act3 = myThid - 1
118              max3 = nTx*nTy
119              act4 = ikey_dynamics - 1
120              igmkey = (act1 + 1) + act2*max1
121         &                      + act3*max1*max2
122         &                      + act4*max1*max2*max3
123    #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  #endif
133        DO j=1-Oly+1,sNy+Oly-1  
134  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef GM_VISBECK_VARIABLE_K
135  !HPF$ INDEPENDENT        DO j=1-Oly,sNy+Oly
136           DO i=1-Olx,sNx+Olx
137            VisbeckK(i,j,bi,bj) = 0. _d 0
138           ENDDO
139          ENDDO
140  #endif  #endif
        DO i=1-Olx+1,sNx+Olx-1  
141    
142  C      Gradient of Sigma at rVel points  C--   set ldd97_Lrho (for tapering scheme ldd97):
143          SlopeX(i,j)=0.25*( sigmaX(i+1, j ,km1) +sigmaX(i,j,km1)        IF ( GM_taper_scheme.EQ.'ldd97' .OR.
144       &                    +sigmaX(i+1, j , k ) +sigmaX(i,j, k ) )       &     GM_taper_scheme.EQ.'fm07' ) THEN
145          SlopeY(i,j)=0.25*( sigmaY( i ,j+1,km1) +sigmaY(i,j,km1)         Cspd = 2. _d 0
146       &                    +sigmaY( i ,j+1, k ) +sigmaY(i,j, k ) )         LrhoInf = 15. _d 3
147          dSigmaDrReal(i,j)=sigmaR(i,j,k)         LrhoSup = 100. _d 3
148    C-     Tracer point location (center):
149          if (hFacC(i,j,k,bi,bj).eq.0.) then         DO j=1-Oly,sNy+Oly
150           SlopeX(i,j)=0.          DO i=1-Olx,sNx+Olx
151           SlopeY(i,j)=0.           IF (fCori(i,j,bi,bj).NE.0.) THEN
152          endif             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.
204    
205          DO j=1-Oly,sNy+Oly
206           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         ENDDO
212        ENDDO        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
234           kkey = (igmkey-1)*Nr + k
235           DO j=1-Oly,sNy+Oly
236            DO i=1-Olx,sNx+Olx
237             SlopeX(i,j)       = 0. _d 0
238             SlopeY(i,j)       = 0. _d 0
239             dSigmaDx(i,j)     = 0. _d 0
240             dSigmaDy(i,j)     = 0. _d 0
241             dSigmaDr(i,j)     = 0. _d 0
242             SlopeSqr(i,j)     = 0. _d 0
243             taperFct(i,j)     = 0. _d 0
244             Kwx(i,j,k,bi,bj)  = 0. _d 0
245             Kwy(i,j,k,bi,bj)  = 0. _d 0
246             Kwz(i,j,k,bi,bj)  = 0. _d 0
247    # ifdef GM_NON_UNITY_DIAGONAL
248             Kux(i,j,k,bi,bj)  = 0. _d 0
249             Kvy(i,j,k,bi,bj)  = 0. _d 0
250    # endif
251    # ifdef GM_EXTRA_DIAGONAL
252             Kuz(i,j,k,bi,bj)  = 0. _d 0
253             Kvz(i,j,k,bi,bj)  = 0. _d 0
254    # endif
255    # ifdef GM_BOLUS_ADVEC
256             GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
257             GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
258    # endif
259            ENDDO
260           ENDDO
261    #endif
262    
263           DO j=1-Oly+1,sNy+Oly-1
264            DO i=1-Olx+1,sNx+Olx-1
265    C      Gradient of Sigma at rVel points
266             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 )
268         &                      )*maskC(i,j,k,bi,bj)
269             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 )
271         &                      )*maskC(i,j,k,bi,bj)
272             dSigmaDr(i,j)=sigmaR(i,j,k)
273            ENDDO
274           ENDDO
275    
276    #ifdef ALLOW_AUTODIFF_TAMC
277    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
278    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
279    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 */
284    
285    #ifdef GM_VISBECK_VARIABLE_K
286    #ifndef OLD_VISBECK_CALC
287           IF ( GM_Visbeck_alpha.GT.0. .AND.
288         &      -rC(k-1).LT.GM_Visbeck_depth ) THEN
289    
290    C--     Depth average of f/sqrt(Ri) = M^2/N^2 * N
291    C       M^2 and N^2 are horizontal & vertical gradient of buoyancy.
292    
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
299             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--      compute: ( M^2 * S )^1/2   (= M^2 / N since S=M^2/N^2 )
314               dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)
315         &             + dSigmaDy(i,j)*dSigmaDy(i,j)
316               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  C     Calculate slopes for use in tensor, taper and/or clip
340        CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
341       I             dSigmadRReal,       O             SlopeX, SlopeY,
342       I             rF(K),       O             SlopeSqr, taperFct,
343       U             SlopeX, SlopeY,       U             hTransLay, baseSlope, recipLambda,
344       O             dRdSigmaLtd,       U             dSigmaDr,
345       I             bi, bj, myThid )       I             dSigmaDx, dSigmaDy,
346         I             ldd97_LrhoC, locMixLayer, rF,
347        DO j=1-Oly+1,sNy+Oly-1       I             kLowC(1-Olx,1-Oly,bi,bj),
348         DO i=1-Olx+1,sNx+Olx-1       I             k, bi, bj, myTime, myIter, myThid )
349    
350  C       Mask Iso-neutral slopes         DO j=1-Oly+1,sNy+Oly-1
351          if (hFacC(i,j,k,bi,bj).eq.0.) then          DO i=1-Olx+1,sNx+Olx-1
352           SlopeX(i,j)=0.  C      Mask Iso-neutral slopes
353           SlopeY(i,j)=0.           SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)
354          endif           SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)
355          Ssq=SlopeX(i,j)*SlopeX(i,j)+SlopeY(i,j)*SlopeY(i,j)           SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)
356            ENDDO
357  C       Components of Redi/GM tensor         ENDDO
358          Kwx(i,j,k,bi,bj)=2.*SlopeX(i,j)  
359          Kwy(i,j,k,bi,bj)=2.*SlopeY(i,j)  #ifdef ALLOW_AUTODIFF_TAMC
360          Kwz(i,j,k,bi,bj)=Ssq  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
361    CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
362    CADJ STORE SlopeSqr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
363    CADJ STORE dSigmaDr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
364    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
365    #endif /* ALLOW_AUTODIFF_TAMC */
366    
367    C      Components of Redi/GM tensor
368           DO j=1-Oly+1,sNy+Oly-1
369            DO i=1-Olx+1,sNx+Olx-1
370              Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
371              Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
372              Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
373            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
382    C           but do not know if *taperFct (or **2 ?) is necessary
383            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
386    
387  C       Calculate terms for mean Richardson number  C       Calculate terms for mean Richardson number
# Line 114  C       If negative then we are below th Line 395  C       If negative then we are below th
395  C       Now we convert deltaH to a non-dimensional fraction  C       Now we convert deltaH to a non-dimensional fraction
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  Calt?   if (dSigmaDrReal(i,j).NE.0.) then          IF ( Ssq(i,j).NE.0. .AND. dSigmaDr(i,j).NE.0. ) THEN
400  Calt?    N2=(-Gravity*recip_Rhonil)*dSigmaDrReal(i,j)           N2= -Gravity*recip_RhoConst*dSigmaDr(i,j)
401          if (dRdSigmaLtd(i,j).NE.0.) then           SN=sqrt(Ssq(i,j)*N2)
          N2=(-Gravity*recip_Rhonil)/dRdSigmaLtd(i,j)  
          SN=sqrt(Ssq*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  C       Limit range that KapGM can take          ENDDO
407          VisbeckK(i,j,bi,bj)=         ENDDO
408       &     min(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)  #endif /* OLD_VISBECK_CALC */
409    #endif /* GM_VISBECK_VARIABLE_K */
410    
411    C-- end 1rst loop on vertical level index k
412          ENDDO
413    
414    
415    #ifdef GM_VISBECK_VARIABLE_K
416    #ifdef ALLOW_AUTODIFF_TAMC
417    CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
418    #endif
419          IF ( GM_Visbeck_alpha.GT.0. ) THEN
420    C-     Limit range that KapGM can take
421           DO j=1-Oly+1,sNy+Oly-1
422            DO i=1-Olx+1,sNx+Olx-1
423             VisbeckK(i,j,bi,bj)=
424         &       MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)
425            ENDDO
426           ENDDO
427          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    #ifdef ALLOW_AUTODIFF_TAMC
436           kkey = (igmkey-1)*Nr + k
437    #if (defined (GM_NON_UNITY_DIAGONAL) || \
438         defined (GM_VISBECK_VARIABLE_K))
439    CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
440    CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
441    CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
442    #endif
443    #endif
444    
445  #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE  C-    express the Tensor in term of Diffusivity (= m**2 / s )
446  C--     Time-average        DO k=1,Nr
447          GM_Kwx_T(i,j,k,bi,bj)=GM_Kwx_T(i,j,k,bi,bj)         DO j=1-Oly+1,sNy+Oly-1
448       &                       +Kwx(i,j,k,bi,bj)*deltaTclock          DO i=1-Olx+1,sNx+Olx-1
449          GM_Kwy_T(i,j,k,bi,bj)=GM_Kwy_T(i,j,k,bi,bj)  #ifdef ALLOW_KAPGM_CONTROL
450       &                       +Kwy(i,j,k,bi,bj)*deltaTclock           Kgm_tmp = GM_isopycK + GM_skewflx*kapgm(i,j,k,bi,bj)
451          GM_Kwz_T(i,j,k,bi,bj)=GM_Kwz_T(i,j,k,bi,bj)  #else
452       &                       +Kwz(i,j,k,bi,bj)*deltaTclock           Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K
453  #ifdef GM_VISBECK_VARIABLE_K  #endif
454          IF (K.EQ.Nr)  #ifdef GM_VISBECK_VARIABLE_K
455       &  Visbeck_K_T(i,j,bi,bj)=Visbeck_K_T(i,j,bi,bj)       &           + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)
456       &                       +VisbeckK(i,j,bi,bj)*deltaTclock  #endif
457             Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)
458             Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)
459             Kwz(i,j,k,bi,bj)= ( GM_isopycK
460    #ifdef GM_VISBECK_VARIABLE_K
461         &                     + VisbeckK(i,j,bi,bj)
462  #endif  #endif
463         &                     )*Kwz(i,j,k,bi,bj)
464            ENDDO
465         ENDDO         ENDDO
466        ENDDO        ENDDO
       GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock  
 #endif /* INCLUDE_DIAGNOSTICS_INTERFACE_CODE */  
467    
468    #ifdef ALLOW_DIAGNOSTICS
469          IF ( useDiagnostics .AND. GM_taper_scheme.EQ.'fm07' ) THEN
470           CALL DIAGNOSTICS_FILL( hTransLay, 'GM_hTrsL', 0,1,2,bi,bj,myThid)
471           CALL DIAGNOSTICS_FILL( baseSlope, 'GM_baseS', 0,1,2,bi,bj,myThid)
472           CALL DIAGNOSTICS_FILL(recipLambda,'GM_rLamb', 0,1,2,bi,bj,myThid)
473          ENDIF
474    #endif /* ALLOW_DIAGNOSTICS */
475    
476    
477    #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
478    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
479    C-- 2nd  k loop : compute Tensor Coeff. at U point
480    
481    #ifdef ALLOW_KPP
482          IF ( useKPP ) THEN
483           DO j=1-Oly,sNy+Oly
484            DO i=2-Olx,sNx+Olx
485             locMixLayer(i,j) = ( KPPhbl(i-1,j,bi,bj)
486         &                      + KPPhbl( i ,j,bi,bj) )*op5
487            ENDDO
488           ENDDO
489          ELSE
490    #else
491          IF ( .TRUE. ) THEN
492    #endif
493           DO j=1-Oly,sNy+Oly
494            DO i=2-Olx,sNx+Olx
495             locMixLayer(i,j) = ( hMixLayer(i-1,j,bi,bj)
496         &                      + hMixLayer( i ,j,bi,bj) )*op5
497            ENDDO
498           ENDDO
499          ENDIF
500          DO j=1-Oly,sNy+Oly
501           DO i=1-Olx,sNx+Olx
502             hTransLay(i,j) =  0.
503             baseSlope(i,j) =  0.
504             recipLambda(i,j)= 0.
505           ENDDO
506           DO i=2-Olx,sNx+Olx
507             hTransLay(i,j) = MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
508           ENDDO
509          ENDDO
510    
511          DO k=Nr,1,-1
512           kp1 = MIN(Nr,k+1)
513           maskp1 = 1. _d 0
514           IF (k.GE.Nr) maskp1 = 0. _d 0
515    
 #ifdef GM_NON_UNITY_DIAGONAL  
516  C     Gradient of Sigma at U points  C     Gradient of Sigma at U points
517        DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
518         DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
519          SlopeX(i,j)=sigmaX(i,j,km1)           dSigmaDx(i,j)=sigmaX(i,j,k)
520       &          *_maskW(i,j,k,bi,bj)       &                       *_maskW(i,j,k,bi,bj)
521          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)
522       &                    +sigmaY(i-1, j ,k) +sigmaY(i, j ,k) )       &                       +sigmaY(i-1, j ,k)+sigmaY(i, j ,k)
523       &          *_maskW(i,j,k,bi,bj)       &                      )*_maskW(i,j,k,bi,bj)
524          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 )
525       &                          +sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1) )       &                      +(sigmaR(i-1,j,kp1)+sigmaR(i,j,kp1))*maskp1
526       &          *_maskW(i,j,k,bi,bj)       &                      )*_maskW(i,j,k,bi,bj)
527            ENDDO
528         ENDDO         ENDDO
529        ENDDO  
530    #ifdef ALLOW_AUTODIFF_TAMC
531    CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
532    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
533    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
534    CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
535    CADJ STORE locMixLayer(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
536    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
537    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
538    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
539    #endif /* ALLOW_AUTODIFF_TAMC */
540    
541  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
542        CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
543       I             dSigmadRReal,       O             SlopeX, SlopeY,
544       I             rF(K),       O             SlopeSqr, taperFct,
545       U             SlopeX, SlopeY,       U             hTransLay, baseSlope, recipLambda,
546       O             dRdSigmaLtd,       U             dSigmaDr,
547       I             bi, bj, myThid )       I             dSigmaDx, dSigmaDy,
548         I             ldd97_LrhoW, locMixLayer, rC,
549        DO j=1-Oly+1,sNy+Oly-1       I             kLow_W,
550         DO i=1-Olx+1,sNx+Olx-1       I             k, bi, bj, myTime, myIter, myThid )
551          Kux(i,j,k,bi,bj)=(dSigmaDrReal(i,j)*dRdSigmaLtd(i,j))**2  
552    cph( NEW
553    #ifdef ALLOW_AUTODIFF_TAMC
554    cph(
555    CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
556    CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
557    cph)
558    #endif /* ALLOW_AUTODIFF_TAMC */
559    cph)
560    
561    #ifdef GM_NON_UNITY_DIAGONAL
562    c      IF ( GM_nonUnitDiag ) THEN
563            DO j=1-Oly+1,sNy+Oly-1
564             DO i=1-Olx+1,sNx+Olx-1
565              Kux(i,j,k,bi,bj) =
566         &     ( GM_isopycK
567    #ifdef GM_VISBECK_VARIABLE_K
568         &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
569    #endif
570         &     )*taperFct(i,j)
571             ENDDO
572            ENDDO
573    #ifdef ALLOW_AUTODIFF_TAMC
574    # ifdef GM_EXCLUDE_CLIPPING
575    CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
576    # endif
577    #endif
578            DO j=1-Oly+1,sNy+Oly-1
579             DO i=1-Olx+1,sNx+Olx-1
580              Kux(i,j,k,bi,bj) = MAX( Kux(i,j,k,bi,bj), GM_Kmin_horiz )
581             ENDDO
582            ENDDO
583    c      ENDIF
584    #endif /* GM_NON_UNITY_DIAGONAL */
585    
586    #ifdef GM_EXTRA_DIAGONAL
587    
588    #ifdef ALLOW_AUTODIFF_TAMC
589    CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
590    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
591    #endif /* ALLOW_AUTODIFF_TAMC */
592           IF ( GM_ExtraDiag ) THEN
593            DO j=1-Oly+1,sNy+Oly-1
594             DO i=1-Olx+1,sNx+Olx-1
595              Kuz(i,j,k,bi,bj) =
596    #ifdef ALLOW_KAPGM_CONTROL
597         &     ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)
598    #else
599         &     ( GM_isopycK - GM_skewflx*GM_background_K
600    #endif
601    #ifdef GM_VISBECK_VARIABLE_K
602         &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*GM_advect
603    #endif
604         &     )*SlopeX(i,j)*taperFct(i,j)
605             ENDDO
606            ENDDO
607           ENDIF
608    #endif /* GM_EXTRA_DIAGONAL */
609    
610    #ifdef ALLOW_DIAGNOSTICS
611           IF (doDiagRediFlx) THEN
612            km1 = MAX(k-1,1)
613            DO j=1,sNy
614             DO i=1,sNx+1
615    C         store in tmp1k Kuz_Redi
616              tmp1k(i,j) = ( GM_isopycK
617    #ifdef GM_VISBECK_VARIABLE_K
618         &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0
619    #endif
620         &                 )*SlopeX(i,j)*taperFct(i,j)
621             ENDDO
622            ENDDO
623            DO j=1,sNy
624             DO i=1,sNx+1
625    C-        Vertical gradients interpolated to U points
626              dTdz = (
627         &     +recip_drC(k)*
628         &       ( maskC(i-1,j,k,bi,bj)*
629         &           (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
630         &        +maskC( i ,j,k,bi,bj)*
631         &           (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
632         &       )
633         &     +recip_drC(kp1)*
634         &       ( maskC(i-1,j,kp1,bi,bj)*
635         &           (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
636         &        +maskC( i ,j,kp1,bi,bj)*
637         &           (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
638         &       )      ) * 0.25 _d 0
639               tmp1k(i,j) = dyG(i,j,bi,bj)*drF(k)
640         &                * _hFacW(i,j,k,bi,bj)
641         &                * tmp1k(i,j) * dTdz
642             ENDDO
643            ENDDO
644            CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KuzTz', k,1,2,bi,bj,myThid)
645           ENDIF
646    #endif /* ALLOW_DIAGNOSTICS */
647    
648    C-- end 2nd  loop on vertical level index k
649          ENDDO
650    
651    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
652    C-- 3rd  k loop : compute Tensor Coeff. at V point
653    
654    #ifdef ALLOW_KPP
655          IF ( useKPP ) THEN
656           DO j=2-Oly,sNy+Oly
657            DO i=1-Olx,sNx+Olx
658             locMixLayer(i,j) = ( KPPhbl(i,j-1,bi,bj)
659         &                      + KPPhbl(i, j ,bi,bj) )*op5
660            ENDDO
661           ENDDO
662          ELSE
663    #else
664          IF ( .TRUE. ) THEN
665    #endif
666           DO j=2-Oly,sNy+Oly
667            DO i=1-Olx,sNx+Olx
668             locMixLayer(i,j) = ( hMixLayer(i,j-1,bi,bj)
669         &                      + hMixLayer(i, j ,bi,bj) )*op5
670            ENDDO
671           ENDDO
672          ENDIF
673          DO j=1-Oly,sNy+Oly
674           DO i=1-Olx,sNx+Olx
675             hTransLay(i,j) =  0.
676             baseSlope(i,j) =  0.
677             recipLambda(i,j)= 0.
678           ENDDO
679          ENDDO
680          DO j=2-Oly,sNy+Oly
681           DO i=1-Olx,sNx+Olx
682             hTransLay(i,j) = MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
683         ENDDO         ENDDO
684        ENDDO        ENDDO
685    
686  C     Gradient of Sigma at V points  C     Gradient of Sigma at V points
687        DO j=1-Oly+1,sNy+Oly-1        DO k=Nr,1,-1
688         DO i=1-Olx+1,sNx+Olx-1         kp1 = MIN(Nr,k+1)
689          SlopeX(i,j)=0.25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)         maskp1 = 1. _d 0
690       &                    +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k) )         IF (k.GE.Nr) maskp1 = 0. _d 0
691       &          *_maskS(i,j,k,bi,bj)  
692          SlopeY(i,j)=sigmaY(i,j,km1)         DO j=1-Oly+1,sNy+Oly-1
693       &          *_maskS(i,j,k,bi,bj)          DO i=1-Olx+1,sNx+Olx-1
694          dSigmaDrReal(i,j)=0.25*( sigmaR(i,j-1, k ) +sigmaR(i,j, k )           dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
695       &                          +sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1) )       &                       +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k)
696       &          *_maskS(i,j,k,bi,bj)       &                      )*_maskS(i,j,k,bi,bj)
697             dSigmaDy(i,j)=sigmaY(i,j,k)
698         &                       *_maskS(i,j,k,bi,bj)
699             dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k )+sigmaR(i,j, k )
700         &                      +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1
701         &                      )*_maskS(i,j,k,bi,bj)
702            ENDDO
703         ENDDO         ENDDO
704        ENDDO  
705    #ifdef ALLOW_AUTODIFF_TAMC
706    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
707    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
708    CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
709    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
710    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
711    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
712    #endif /* ALLOW_AUTODIFF_TAMC */
713    
714  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
715        CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
716       I             dSigmadRReal,       O             SlopeX, SlopeY,
717       I             rF(K),       O             SlopeSqr, taperFct,
718       U             SlopeX, SlopeY,       U             hTransLay, baseSlope, recipLambda,
719       O             dRdSigmaLtd,       U             dSigmaDr,
720       I             bi, bj, myThid )       I             dSigmaDx, dSigmaDy,
721         I             ldd97_LrhoS, locMixLayer, rC,
722        DO j=1-Oly+1,sNy+Oly-1       I             kLow_S,
723         DO i=1-Olx+1,sNx+Olx-1       I             k, bi, bj, myTime, myIter, myThid )
         Kvy(i,j,k,bi,bj)=(dSigmaDrReal(i,j)*dRdSigmaLtd(i,j))**2  
        ENDDO  
       ENDDO  
724    
725    cph(
726    #ifdef ALLOW_AUTODIFF_TAMC
727    cph(
728    CADJ STORE taperfct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
729    cph)
730    #endif /* ALLOW_AUTODIFF_TAMC */
731    cph)
732    
733    #ifdef GM_NON_UNITY_DIAGONAL
734    c      IF ( GM_nonUnitDiag ) THEN
735            DO j=1-Oly+1,sNy+Oly-1
736             DO i=1-Olx+1,sNx+Olx-1
737              Kvy(i,j,k,bi,bj) =
738         &     ( GM_isopycK
739    #ifdef GM_VISBECK_VARIABLE_K
740         &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
741    #endif
742         &     )*taperFct(i,j)
743             ENDDO
744            ENDDO
745    #ifdef ALLOW_AUTODIFF_TAMC
746    # ifdef GM_EXCLUDE_CLIPPING
747    CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
748    # endif
749    #endif
750            DO j=1-Oly+1,sNy+Oly-1
751             DO i=1-Olx+1,sNx+Olx-1
752              Kvy(i,j,k,bi,bj) = MAX( Kvy(i,j,k,bi,bj), GM_Kmin_horiz )
753             ENDDO
754            ENDDO
755    c      ENDIF
756  #endif /* GM_NON_UNITY_DIAGONAL */  #endif /* GM_NON_UNITY_DIAGONAL */
757    
758    #ifdef GM_EXTRA_DIAGONAL
759    
760    #ifdef ALLOW_AUTODIFF_TAMC
761    CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
762    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
763    #endif /* ALLOW_AUTODIFF_TAMC */
764           IF ( GM_ExtraDiag ) THEN
765            DO j=1-Oly+1,sNy+Oly-1
766             DO i=1-Olx+1,sNx+Olx-1
767              Kvz(i,j,k,bi,bj) =
768    #ifdef ALLOW_KAPGM_CONTROL
769         &     ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)
770    #else
771         &     ( GM_isopycK - GM_skewflx*GM_background_K
772    #endif
773    #ifdef GM_VISBECK_VARIABLE_K
774         &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*GM_advect
775    #endif
776         &     )*SlopeY(i,j)*taperFct(i,j)
777             ENDDO
778            ENDDO
779           ENDIF
780    #endif /* GM_EXTRA_DIAGONAL */
781    
782    #ifdef ALLOW_DIAGNOSTICS
783           IF (doDiagRediFlx) THEN
784    c       km1 = MAX(k-1,1)
785            DO j=1,sNy+1
786             DO i=1,sNx
787    C         store in tmp1k Kvz_Redi
788              tmp1k(i,j) = ( GM_isopycK
789    #ifdef GM_VISBECK_VARIABLE_K
790         &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0
791    #endif
792         &                 )*SlopeY(i,j)*taperFct(i,j)
793             ENDDO
794            ENDDO
795            DO j=1,sNy+1
796             DO i=1,sNx
797    C-        Vertical gradients interpolated to U points
798              dTdz = (
799         &     +recip_drC(k)*
800         &       ( maskC(i,j-1,k,bi,bj)*
801         &           (theta(i,j-1,km1,bi,bj)-theta(i,j-1,k,bi,bj))
802         &        +maskC(i, j ,k,bi,bj)*
803         &           (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))
804         &       )
805         &     +recip_drC(kp1)*
806         &       ( maskC(i,j-1,kp1,bi,bj)*
807         &           (theta(i,j-1,k,bi,bj)-theta(i,j-1,kp1,bi,bj))
808         &        +maskC(i, j ,kp1,bi,bj)*
809         &           (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))
810         &       )      ) * 0.25 _d 0
811               tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)
812         &                * _hFacS(i,j,k,bi,bj)
813         &                * tmp1k(i,j) * dTdz
814             ENDDO
815            ENDDO
816            CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KvzTz', k,1,2,bi,bj,myThid)
817           ENDIF
818    #endif /* ALLOW_DIAGNOSTICS */
819    
820    C-- end 3rd  loop on vertical level index k
821          ENDDO
822    
823    #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
824    
825    
826    #ifdef GM_BOLUS_ADVEC
827          IF (GM_AdvForm) THEN
828           CALL GMREDI_CALC_PSI_B(
829         I             bi, bj, iMin, iMax, jMin, jMax,
830         I             sigmaX, sigmaY, sigmaR,
831         I             ldd97_LrhoW, ldd97_LrhoS,
832         I             myThid )
833          ENDIF
834    #endif
835    
836    #ifdef ALLOW_TIMEAVE
837    C--   Time-average
838          IF ( taveFreq.GT.0. ) THEN
839    
840             CALL TIMEAVE_CUMULATE( GM_Kwx_T, Kwx, Nr,
841         &                          deltaTclock, bi, bj, myThid )
842             CALL TIMEAVE_CUMULATE( GM_Kwy_T, Kwy, Nr,
843         &                          deltaTclock, bi, bj, myThid )
844             CALL TIMEAVE_CUMULATE( GM_Kwz_T, Kwz, Nr,
845         &                          deltaTclock, bi, bj, myThid )
846    #ifdef GM_VISBECK_VARIABLE_K
847           IF ( GM_Visbeck_alpha.NE.0. ) THEN
848             CALL TIMEAVE_CUMULATE( Visbeck_K_T, VisbeckK, 1,
849         &                          deltaTclock, bi, bj, myThid )
850           ENDIF
851    #endif
852    #ifdef GM_BOLUS_ADVEC
853           IF ( GM_AdvForm ) THEN
854             CALL TIMEAVE_CUMULATE( GM_PsiXtave, GM_PsiX, Nr,
855         &                          deltaTclock, bi, bj, myThid )
856             CALL TIMEAVE_CUMULATE( GM_PsiYtave, GM_PsiY, Nr,
857         &                          deltaTclock, bi, bj, myThid )
858           ENDIF
859    #endif
860           DO k=1,Nr
861             GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock
862           ENDDO
863    
864          ENDIF
865    #endif /* ALLOW_TIMEAVE */
866    
867    #ifdef ALLOW_DIAGNOSTICS
868          IF ( useDiagnostics ) THEN
869            CALL GMREDI_DIAGNOSTICS_FILL(bi,bj,myThid)
870          ENDIF
871    #endif /* ALLOW_DIAGNOSTICS */
872    
873  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
874    
# Line 219  C     Calculate slopes for use in tensor Line 877  C     Calculate slopes for use in tensor
877    
878    
879        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
880       I             bi, bj, iMin, iMax, jMin, jMax, K,       I             iMin, iMax, jMin, jMax,
881       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
882       I             myThid )       I             bi, bj, myTime, myIter, myThid )
883  C     /==========================================================\  C     /==========================================================\
884  C     | SUBROUTINE GMREDI_CALC_TENSOR                            |  C     | SUBROUTINE GMREDI_CALC_TENSOR                            |
885  C     | o Calculate tensor elements for GM/Redi tensor.          |  C     | o Calculate tensor elements for GM/Redi tensor.          |
# Line 231  C     \================================= Line 889  C     \=================================
889    
890  C     == Global variables ==  C     == Global variables ==
891  #include "SIZE.h"  #include "SIZE.h"
 #include "GRID.h"  
 #include "DYNVARS.h"  
892  #include "EEPARAMS.h"  #include "EEPARAMS.h"
 #include "PARAMS.h"  
893  #include "GMREDI.h"  #include "GMREDI.h"
894    
895  C     == Routine arguments ==  C     == Routine arguments ==
# Line 242  C Line 897  C
897        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
898        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
899        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
900        INTEGER bi,bj,iMin,iMax,jMin,jMax,K        INTEGER iMin,iMax,jMin,jMax
901          INTEGER bi, bj
902          _RL     myTime
903          INTEGER myIter
904        INTEGER myThid        INTEGER myThid
905  CEndOfInterface  CEndOfInterface
906    
       INTEGER i, j  
   
907  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
908    
909        DO j=1-Oly+1,sNy+Oly-1        INTEGER i, j, k
910         DO i=1-Olx+1,sNx+Olx-1  
911          Kwx(i,j,k,bi,bj) = 0.0        DO k=1,Nr
912          Kwy(i,j,k,bi,bj) = 0.0         DO j=1-Oly+1,sNy+Oly-1
913          Kwz(i,j,k,bi,bj) = 0.0          DO i=1-Olx+1,sNx+Olx-1
914             Kwx(i,j,k,bi,bj) = 0.0
915             Kwy(i,j,k,bi,bj) = 0.0
916             Kwz(i,j,k,bi,bj) = 0.0
917            ENDDO
918         ENDDO         ENDDO
919        ENDDO        ENDDO
920  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
921    
922        end        RETURN
923          END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.28
  ViewVC Help
Powered by ViewVC 1.1.22