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