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