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