/[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.24 by jmc, Tue Jun 20 22:55:08 2006 UTC revision 1.33 by jmc, Tue Oct 21 22:10:55 2008 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 23  C     == Global variables == Line 33  C     == Global variables ==
33  #include "PARAMS.h"  #include "PARAMS.h"
34  #include "GMREDI.h"  #include "GMREDI.h"
35  #include "GMREDI_TAVE.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)
# Line 55  C     == Local variables == Line 77  C     == Local variables ==
77        _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
78        _RL Cspd, LrhoInf, LrhoSup, fCoriLoc        _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
92          _RL Sloc, M2loc
93    #endif
94          _RL deltaH, integrDepth
95          _RL N2loc, SNloc
96  #endif  #endif
97    
98  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
# Line 90  C---+----1----+----2----+----3----+----4 Line 123  C---+----1----+----2----+----3----+----4
123        doDiagRediFlx = .FALSE.        doDiagRediFlx = .FALSE.
124        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
125          doDiagRediFlx = DIAGNOSTICS_IS_ON('GM_KuzTz', myThid )          doDiagRediFlx = DIAGNOSTICS_IS_ON('GM_KuzTz', myThid )
126          doDiagRediFlx = doDiagRediFlx .OR.          doDiagRediFlx = doDiagRediFlx .OR.
127       &                  DIAGNOSTICS_IS_ON('GM_KvzTz', myThid )       &                  DIAGNOSTICS_IS_ON('GM_KvzTz', myThid )
128        ENDIF        ENDIF
129  #endif  #endif
130        
131  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
132        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
133         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
# Line 104  C---+----1----+----2----+----3----+----4 Line 137  C---+----1----+----2----+----3----+----4
137  #endif  #endif
138    
139  C--   set ldd97_Lrho (for tapering scheme ldd97):  C--   set ldd97_Lrho (for tapering scheme ldd97):
140        IF (GM_taper_scheme.EQ.'ldd97') THEN        IF ( GM_taper_scheme.EQ.'ldd97' .OR.
141         &     GM_taper_scheme.EQ.'fm07' ) THEN
142         Cspd = 2. _d 0         Cspd = 2. _d 0
143         LrhoInf = 15. _d 3         LrhoInf = 15. _d 3
144         LrhoSup = 100. _d 3         LrhoSup = 100. _d 3
# Line 121  C-     Tracer point location (center): Line 155  C-     Tracer point location (center):
155         ENDDO         ENDDO
156  C-     U point location (West):  C-     U point location (West):
157         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
158            kLow_W(1-Olx,j) = 0
159          ldd97_LrhoW(1-Olx,j) = LrhoSup          ldd97_LrhoW(1-Olx,j) = LrhoSup
160          DO i=1-Olx+1,sNx+Olx          DO i=1-Olx+1,sNx+Olx
161             kLow_W(i,j) = MIN(kLowC(i-1,j,bi,bj),kLowC(i,j,bi,bj))
162           fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))           fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))
163           IF (fCoriLoc.NE.0.) THEN           IF (fCoriLoc.NE.0.) THEN
164             ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)             ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)
# Line 134  C-     U point location (West): Line 170  C-     U point location (West):
170         ENDDO         ENDDO
171  C-     V point location (South):  C-     V point location (South):
172         DO i=1-Olx+1,sNx+Olx         DO i=1-Olx+1,sNx+Olx
173             kLow_S(i,1-Oly) = 0
174           ldd97_LrhoS(i,1-Oly) = LrhoSup           ldd97_LrhoS(i,1-Oly) = LrhoSup
175         ENDDO         ENDDO
176         DO j=1-Oly+1,sNy+Oly         DO j=1-Oly+1,sNy+Oly
177          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
178             kLow_S(i,j) = MIN(kLowC(i,j-1,bi,bj),kLowC(i,j,bi,bj))
179           fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))           fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))
180           IF (fCoriLoc.NE.0.) THEN           IF (fCoriLoc.NE.0.) THEN
181             ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)             ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)
# Line 157  C-    Just initialize to zero (not use a Line 195  C-    Just initialize to zero (not use a
195          ENDDO          ENDDO
196         ENDDO         ENDDO
197        ENDIF        ENDIF
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
198    
199        DO k=2,Nr  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
200  C-- 1rst loop on k : compute Tensor Coeff. at W points.  C-- 1rst loop on k : compute Tensor Coeff. at W points.
201    
202          DO j=1-Oly,sNy+Oly
203           DO i=1-Olx,sNx+Olx
204             hTransLay(i,j) = R_low(i,j,bi,bj)
205             baseSlope(i,j) =  0. _d 0
206             recipLambda(i,j) = 0. _d 0
207             locMixLayer(i,j) = 0. _d 0
208           ENDDO
209          ENDDO
210    #ifdef ALLOW_KPP
211          IF ( useKPP ) THEN
212           DO j=1-Oly,sNy+Oly
213            DO i=1-Olx,sNx+Olx
214             locMixLayer(i,j) = KPPhbl(i,j,bi,bj)
215            ENDDO
216           ENDDO
217          ELSE
218    #else
219          IF ( .TRUE. ) THEN
220    #endif
221           DO j=1-Oly,sNy+Oly
222            DO i=1-Olx,sNx+Olx
223             locMixLayer(i,j) = hMixLayer(i,j,bi,bj)
224            ENDDO
225           ENDDO
226          ENDIF
227    
228          DO k=Nr,2,-1
229    
230  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
231         kkey = (igmkey-1)*Nr + k         kkey = (igmkey-1)*Nr + k
232         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
# Line 192  C-- 1rst loop on k : compute Tensor Coef Line 257  C-- 1rst loop on k : compute Tensor Coef
257         ENDDO         ENDDO
258  #endif  #endif
259    
260        DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
261         DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
262  C      Gradient of Sigma at rVel points  C      Gradient of Sigma at rVel points
263          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)
264       &                    +sigmaX(i+1, j , k ) +sigmaX(i,j, k ) )       &                       +sigmaX(i+1,j, k )+sigmaX(i,j, k )
265       &                  *maskC(i,j,k,bi,bj)       &                      )*maskC(i,j,k,bi,bj)
266          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)
267       &                    +sigmaY( i ,j+1, k ) +sigmaY(i,j, k ) )       &                       +sigmaY(i,j+1, k )+sigmaY(i,j, k )
268       &                  *maskC(i,j,k,bi,bj)       &                      )*maskC(i,j,k,bi,bj)
269          dSigmaDr(i,j)=sigmaR(i,j,k)           dSigmaDr(i,j)=sigmaR(i,j,k)
270            ENDDO
271         ENDDO         ENDDO
       ENDDO  
272    
273  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
274  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
275  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
276  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
277    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
278    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
279    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
280  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
281    
282    #ifdef GM_VISBECK_VARIABLE_K
283    #ifndef OLD_VISBECK_CALC
284           IF ( GM_Visbeck_alpha.GT.0. .AND.
285         &      -rC(k-1).LT.GM_Visbeck_depth ) THEN
286    
287    C--     Depth average of f/sqrt(Ri) = M^2/N^2 * N
288    C       M^2 and N^2 are horizontal & vertical gradient of buoyancy.
289    
290    C       Calculate terms for mean Richardson number which is used
291    C       in the "variable K" parameterisaton:
292    C       compute depth average from surface down to the bottom or
293    C       GM_Visbeck_depth, whatever is the shallower.
294    
295            DO j=1-Oly+1,sNy+Oly-1
296             DO i=1-Olx+1,sNx+Olx-1
297              IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
298               integrDepth = -rC( kLowC(i,j,bi,bj) )
299    C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments
300               integrDepth = MIN( integrDepth, GM_Visbeck_depth )
301    C       Distance between level center above and the integration depth
302               deltaH = integrDepth + rC(k-1)
303    C       If negative then we are below the integration level
304    C       (cannot be the case with 2 conditions on maskC & -rC(k-1))
305    C       If positive we limit this to the distance from center above
306               deltaH = MIN( deltaH, drC(k) )
307    C       Now we convert deltaH to a non-dimensional fraction
308               deltaH = deltaH/( integrDepth+rC(1) )
309    
310    C--      compute: ( M^2 * S )^1/2   (= S*N since S=M^2/N^2 )
311               dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)
312         &             + dSigmaDy(i,j)*dSigmaDy(i,j)
313               IF ( dSigmaH .GT. 0. _d 0 ) THEN
314                 dSigmaH = SQRT( dSigmaH )
315    C-       compute slope, limited by GM_maxSlope:
316                 IF ( -dSigmaDr(i,j).GT.dSigmaH*GM_rMaxSlope ) THEN
317                  Sloc = dSigmaH / ( -dSigmaDr(i,j) )
318                 ELSE
319                  Sloc = GM_maxSlope
320                 ENDIF
321                 M2loc = gravity*recip_rhoConst*dSigmaH
322    c            SNloc = SQRT( Sloc*M2loc )
323                 N2loc = -gravity*recip_rhoConst*dSigmaDr(i,j)
324                 IF ( N2loc.GT.0. _d 0 ) THEN
325                   SNloc = Sloc*SQRT(N2loc)
326                 ELSE
327                   SNloc = 0. _d 0
328                 ENDIF
329               ELSE
330                 SNloc = 0. _d 0
331               ENDIF
332               VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
333         &       +deltaH*GM_Visbeck_alpha
334         &              *GM_Visbeck_length*GM_Visbeck_length*SNloc
335              ENDIF
336             ENDDO
337            ENDDO
338           ENDIF
339    #endif /* ndef OLD_VISBECK_CALC */
340    #endif /* GM_VISBECK_VARIABLE_K */
341    
342  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
343        CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
344       O             SlopeX, SlopeY,       O             SlopeX, SlopeY,
345       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
346         U             hTransLay, baseSlope, recipLambda,
347       U             dSigmaDr,       U             dSigmaDr,
348       I             dSigmaDx, dSigmaDy,       I             dSigmaDx, dSigmaDy,
349       I             ldd97_LrhoC,rF(k),k,       I             ldd97_LrhoC, locMixLayer, rF,
350       I             bi, bj, myThid )       I             kLowC(1-Olx,1-Oly,bi,bj),
351         I             k, bi, bj, myTime, myIter, myThid )
       DO j=1-Oly+1,sNy+Oly-1  
        DO i=1-Olx+1,sNx+Olx-1  
   
 C       Mask Iso-neutral slopes  
         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)  
352    
353           DO j=1-Oly+1,sNy+Oly-1
354            DO i=1-Olx+1,sNx+Olx-1
355    C      Mask Iso-neutral slopes
356             SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)
357             SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)
358             SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)
359            ENDDO
360         ENDDO         ENDDO
       ENDDO  
361    
362  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
363  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
# Line 239  CADJ STORE dSigmaDr(:,:)     = comlev1_b Line 367  CADJ STORE dSigmaDr(:,:)     = comlev1_b
367  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
368  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
369    
370        DO j=1-Oly+1,sNy+Oly-1  C      Components of Redi/GM tensor
371         DO i=1-Olx+1,sNx+Olx-1         DO j=1-Oly+1,sNy+Oly-1
372            DO i=1-Olx+1,sNx+Olx-1
373  C       Components of Redi/GM tensor            Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
374          Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)            Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
375          Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)            Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
376          Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)          ENDDO
377           ENDDO
378    
379  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
380    #ifdef OLD_VISBECK_CALC
381           DO j=1-Oly+1,sNy+Oly-1
382            DO i=1-Olx+1,sNx+Olx-1
383    
384  C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K  C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K
385  C           but do not know if *taperFct (or **2 ?) is necessary  C           but do not know if *taperFct (or **2 ?) is necessary
386          Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)          Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)
387    
388  C--     Depth average of M^2/N^2 * N  C--     Depth average of M^2/N^2 * N
# Line 260  C       which is used in the "variable K Line 392  C       which is used in the "variable K
392  C       Distance between interface above layer and the integration depth  C       Distance between interface above layer and the integration depth
393          deltaH=abs(GM_Visbeck_depth)-abs(rF(k))          deltaH=abs(GM_Visbeck_depth)-abs(rF(k))
394  C       If positive we limit this to the layer thickness  C       If positive we limit this to the layer thickness
395          deltaH=min(deltaH,drF(k))          integrDepth = drF(k)
396            deltaH=min(deltaH,integrDepth)
397  C       If negative then we are below the integration level  C       If negative then we are below the integration level
398          deltaH=max(deltaH,zero_rs)          deltaH=max(deltaH, 0. _d 0)
399  C       Now we convert deltaH to a non-dimensional fraction  C       Now we convert deltaH to a non-dimensional fraction
400          deltaH=deltaH/GM_Visbeck_depth          deltaH=deltaH/GM_Visbeck_depth
401    
         IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.  
402          IF ( Ssq(i,j).NE.0. .AND. dSigmaDr(i,j).NE.0. ) THEN          IF ( Ssq(i,j).NE.0. .AND. dSigmaDr(i,j).NE.0. ) THEN
403           N2= -Gravity*recip_RhoConst*dSigmaDr(i,j)           N2loc = -gravity*recip_rhoConst*dSigmaDr(i,j)
404           SN=sqrt(Ssq(i,j)*N2)           SNloc = SQRT(Ssq(i,j)*N2loc )
405           VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH           VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
406       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN       &       +deltaH*GM_Visbeck_alpha
407         &              *GM_Visbeck_length*GM_Visbeck_length*SNloc
408          ENDIF          ENDIF
409    
410  #endif /* GM_VISBECK_VARIABLE_K */          ENDDO
   
411         ENDDO         ENDDO
412        ENDDO  #endif /* OLD_VISBECK_CALC */
413    #endif /* GM_VISBECK_VARIABLE_K */
414    
415  C-- end 1rst loop on vertical level index k  C-- end 1rst loop on vertical level index k
416        ENDDO        ENDDO
# Line 287  C-- end 1rst loop on vertical level inde Line 420  C-- end 1rst loop on vertical level inde
420  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
421  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
422  #endif  #endif
423        IF ( GM_Visbeck_alpha.NE.0. ) THEN        IF ( GM_Visbeck_alpha.GT.0. ) THEN
424  C-     Limit range that KapGM can take  C-     Limit range that KapGM can take
425         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
426          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
# Line 303  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1 Line 436  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1
436  cph)  cph)
437  #endif /* GM_VISBECK_VARIABLE_K */  #endif /* GM_VISBECK_VARIABLE_K */
438    
439    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.  
440        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  
   
441  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
442         kkey = (igmkey-1)*Nr + k         kkey = (igmkey-1)*Nr + k
443  #if (defined (GM_NON_UNITY_DIAGONAL) || \  # if (defined (GM_NON_UNITY_DIAGONAL) || \
444       defined (GM_VISBECK_VARIABLE_K))       defined (GM_VISBECK_VARIABLE_K))
445  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
446  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
447  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
448    # endif
449  #endif  #endif
450           DO j=1-Oly+1,sNy+Oly-1
451            DO i=1-Olx+1,sNx+Olx-1
452    #ifdef ALLOW_KAPREDI_CONTROL
453             Kgm_tmp = kapredi(i,j,k,bi,bj)
454    #else
455             Kgm_tmp = GM_isopycK
456    #endif
457    #ifdef ALLOW_KAPGM_CONTROL
458         &           + GM_skewflx*kapgm(i,j,k,bi,bj)
459    #else
460         &           + GM_skewflx*GM_background_K
461  #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  
462  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
463       &          + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)           &           + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)
464    #endif
465             Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)
466             Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)
467    #ifdef ALLOW_KAPREDI_CONTROL
468             Kwz(i,j,k,bi,bj)= ( kapredi(i,j,k,bi,bj)
469    #else
470             Kwz(i,j,k,bi,bj)= ( GM_isopycK
471  #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  
472  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
473       &                    + VisbeckK(i,j,bi,bj)       &                     + VisbeckK(i,j,bi,bj)
474  #endif  #endif
475       &                    )*Kwz(i,j,k,bi,bj)       &                     )*Kwz(i,j,k,bi,bj)
476            ENDDO
477         ENDDO         ENDDO
478        ENDDO        ENDDO
479    
480    #ifdef ALLOW_DIAGNOSTICS
481          IF ( useDiagnostics .AND. GM_taper_scheme.EQ.'fm07' ) THEN
482           CALL DIAGNOSTICS_FILL( hTransLay, 'GM_hTrsL', 0,1,2,bi,bj,myThid)
483           CALL DIAGNOSTICS_FILL( baseSlope, 'GM_baseS', 0,1,2,bi,bj,myThid)
484           CALL DIAGNOSTICS_FILL(recipLambda,'GM_rLamb', 0,1,2,bi,bj,myThid)
485          ENDIF
486    #endif /* ALLOW_DIAGNOSTICS */
487    
488    
489  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
490    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
491    C-- 2nd  k loop : compute Tensor Coeff. at U point
492    
493  C     Gradient of Sigma at U points  #ifdef ALLOW_KPP
494        DO j=1-Oly+1,sNy+Oly-1        IF ( useKPP ) THEN
495         DO i=1-Olx+1,sNx+Olx-1         DO j=1-Oly,sNy+Oly
496          dSigmaDx(i,j)=sigmaX(i,j,k)          DO i=2-Olx,sNx+Olx
497       &          *_maskW(i,j,k,bi,bj)           locMixLayer(i,j) = ( KPPhbl(i-1,j,bi,bj)
498          dSigmaDy(i,j)=op25*( sigmaY(i-1,j+1,k) +sigmaY(i,j+1,k)       &                      + KPPhbl( i ,j,bi,bj) )*op5
499       &                      +sigmaY(i-1, j ,k) +sigmaY(i, j ,k) )          ENDDO
500       &          *_maskW(i,j,k,bi,bj)         ENDDO
501          dSigmaDr(i,j)=op25*( sigmaR(i-1,j, k ) +sigmaR(i,j, k )        ELSE
502       &                  +maskp1*(sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1)) )  #else
503       &          *_maskW(i,j,k,bi,bj)        IF ( .TRUE. ) THEN
504    #endif
505           DO j=1-Oly,sNy+Oly
506            DO i=2-Olx,sNx+Olx
507             locMixLayer(i,j) = ( hMixLayer(i-1,j,bi,bj)
508         &                      + hMixLayer( i ,j,bi,bj) )*op5
509            ENDDO
510           ENDDO
511          ENDIF
512          DO j=1-Oly,sNy+Oly
513           DO i=1-Olx,sNx+Olx
514             hTransLay(i,j) =  0.
515             baseSlope(i,j) =  0.
516             recipLambda(i,j)= 0.
517           ENDDO
518           DO i=2-Olx,sNx+Olx
519             hTransLay(i,j) = MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
520         ENDDO         ENDDO
521        ENDDO        ENDDO
522    
523          DO k=Nr,1,-1
524           kp1 = MIN(Nr,k+1)
525           maskp1 = 1. _d 0
526           IF (k.GE.Nr) maskp1 = 0. _d 0
527    #ifdef ALLOW_AUTODIFF_TAMC
528           kkey = (igmkey-1)*Nr + k
529    #endif
530    
531    C     Gradient of Sigma at U points
532           DO j=1-Oly+1,sNy+Oly-1
533            DO i=1-Olx+1,sNx+Olx-1
534             dSigmaDx(i,j)=sigmaX(i,j,k)
535         &                       *_maskW(i,j,k,bi,bj)
536             dSigmaDy(i,j)=op25*( sigmaY(i-1,j+1,k)+sigmaY(i,j+1,k)
537         &                       +sigmaY(i-1, j ,k)+sigmaY(i, j ,k)
538         &                      )*_maskW(i,j,k,bi,bj)
539             dSigmaDr(i,j)=op25*( sigmaR(i-1,j, k )+sigmaR(i,j, k )
540         &                      +(sigmaR(i-1,j,kp1)+sigmaR(i,j,kp1))*maskp1
541         &                      )*_maskW(i,j,k,bi,bj)
542            ENDDO
543           ENDDO
544    
545  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
546  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
547  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
548  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
549  CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
550    CADJ STORE locMixLayer(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
551    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
552    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
553    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
554  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
555    
556  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
557        CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
558       O             SlopeX, SlopeY,       O             SlopeX, SlopeY,
559       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
560         U             hTransLay, baseSlope, recipLambda,
561       U             dSigmaDr,       U             dSigmaDr,
562       I             dSigmaDx, dSigmaDy,       I             dSigmaDx, dSigmaDy,
563       I             ldd97_LrhoW,rC(k),k,       I             ldd97_LrhoW, locMixLayer, rC,
564       I             bi, bj, myThid )       I             kLow_W,
565         I             k, bi, bj, myTime, myIter, myThid )
566    
 cph( NEW  
567  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
 cph(  
568  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
569  CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
 cph)  
570  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
 cph)  
571    
572  #ifdef GM_NON_UNITY_DIAGONAL  #ifdef GM_NON_UNITY_DIAGONAL
573    c      IF ( GM_nonUnitDiag ) THEN
574          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
575           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
576            Kux(i,j,k,bi,bj) =            Kux(i,j,k,bi,bj) =
577    #ifdef ALLOW_KAPREDI_CONTROL
578         &     ( kapredi(i,j,k,bi,bj)
579    #else
580       &     ( GM_isopycK       &     ( GM_isopycK
581    #endif
582  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
583       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
584  #endif  #endif
585       &     )       &     )*taperFct(i,j)
      &     *taperFct(i,j)  
586           ENDDO           ENDDO
587          ENDDO          ENDDO
588  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 402  CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_b Line 595  CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_b
595            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 )
596           ENDDO           ENDDO
597          ENDDO          ENDDO
598    c      ENDIF
599  #endif /* GM_NON_UNITY_DIAGONAL */  #endif /* GM_NON_UNITY_DIAGONAL */
600    
601  #ifdef GM_EXTRA_DIAGONAL  #ifdef GM_EXTRA_DIAGONAL
# Line 410  CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_b Line 604  CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_b
604  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
605  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
606  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
607        IF (GM_ExtraDiag) THEN         IF ( GM_ExtraDiag ) THEN
608          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
609           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
610            Kuz(i,j,k,bi,bj) =            Kuz(i,j,k,bi,bj) =
611       &     ( GM_isopycK - GM_skewflx*GM_background_K  #ifdef ALLOW_KAPREDI_CONTROL
612         &     ( kapredi(i,j,k,bi,bj)
613    #else
614         &     ( GM_isopycK
615    #endif
616    #ifdef ALLOW_KAPGM_CONTROL
617         &     - GM_skewflx*kapgm(i,j,k,bi,bj)
618    #else
619         &     - GM_skewflx*GM_background_K
620    #endif
621  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
622       &     +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
623  #endif  #endif
624       &     )*SlopeX(i,j)*taperFct(i,j)       &     )*SlopeX(i,j)*taperFct(i,j)
625           ENDDO           ENDDO
626          ENDDO          ENDDO
627        ENDIF         ENDIF
628  #endif /* GM_EXTRA_DIAGONAL */  #endif /* GM_EXTRA_DIAGONAL */
629    
630  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
631        IF (doDiagRediFlx) THEN         IF (doDiagRediFlx) THEN
632          km1 = MAX(k-1,1)          km1 = MAX(k-1,1)
633          DO j=1,sNy          DO j=1,sNy
634           DO i=1,sNx+1           DO i=1,sNx+1
635  C         store in tmp1k Kuz_Redi  C         store in tmp1k Kuz_Redi
636    #ifdef ALLOW_KAPREDI_CONTROL
637              tmp1k(i,j) = ( kapredi(i,j,k,bi,bj)
638    #else
639            tmp1k(i,j) = ( GM_isopycK            tmp1k(i,j) = ( GM_isopycK
640    #endif
641  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
642       &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0       &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0
643  #endif  #endif
# Line 459  C-        Vertical gradients interpolate Line 666  C-        Vertical gradients interpolate
666           ENDDO           ENDDO
667          ENDDO          ENDDO
668          CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KuzTz', k,1,2,bi,bj,myThid)          CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KuzTz', k,1,2,bi,bj,myThid)
669        ENDIF         ENDIF
670  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
671    
672  C     Gradient of Sigma at V points  C-- end 2nd  loop on vertical level index k
673        DO j=1-Oly+1,sNy+Oly-1        ENDDO
674         DO i=1-Olx+1,sNx+Olx-1  
675          dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
676       &                    +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k) )  C-- 3rd  k loop : compute Tensor Coeff. at V point
677       &          *_maskS(i,j,k,bi,bj)  
678          dSigmaDy(i,j)=sigmaY(i,j,k)  #ifdef ALLOW_KPP
679       &          *_maskS(i,j,k,bi,bj)        IF ( useKPP ) THEN
680          dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k ) +sigmaR(i,j, k )         DO j=2-Oly,sNy+Oly
681       &                  +maskp1*(sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1)) )          DO i=1-Olx,sNx+Olx
682       &          *_maskS(i,j,k,bi,bj)           locMixLayer(i,j) = ( KPPhbl(i,j-1,bi,bj)
683         &                      + KPPhbl(i, j ,bi,bj) )*op5
684            ENDDO
685           ENDDO
686          ELSE
687    #else
688          IF ( .TRUE. ) THEN
689    #endif
690           DO j=2-Oly,sNy+Oly
691            DO i=1-Olx,sNx+Olx
692             locMixLayer(i,j) = ( hMixLayer(i,j-1,bi,bj)
693         &                      + hMixLayer(i, j ,bi,bj) )*op5
694            ENDDO
695           ENDDO
696          ENDIF
697          DO j=1-Oly,sNy+Oly
698           DO i=1-Olx,sNx+Olx
699             hTransLay(i,j) =  0.
700             baseSlope(i,j) =  0.
701             recipLambda(i,j)= 0.
702           ENDDO
703          ENDDO
704          DO j=2-Oly,sNy+Oly
705           DO i=1-Olx,sNx+Olx
706             hTransLay(i,j) = MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
707         ENDDO         ENDDO
708        ENDDO        ENDDO
709    
710    C     Gradient of Sigma at V points
711          DO k=Nr,1,-1
712           kp1 = MIN(Nr,k+1)
713           maskp1 = 1. _d 0
714           IF (k.GE.Nr) maskp1 = 0. _d 0
715    #ifdef ALLOW_AUTODIFF_TAMC
716           kkey = (igmkey-1)*Nr + k
717    #endif
718    
719           DO j=1-Oly+1,sNy+Oly-1
720            DO i=1-Olx+1,sNx+Olx-1
721             dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
722         &                       +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k)
723         &                      )*_maskS(i,j,k,bi,bj)
724             dSigmaDy(i,j)=sigmaY(i,j,k)
725         &                       *_maskS(i,j,k,bi,bj)
726             dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k )+sigmaR(i,j, k )
727         &                      +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1
728         &                      )*_maskS(i,j,k,bi,bj)
729            ENDDO
730           ENDDO
731    
732  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
733  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
734  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
735  CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
736    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
737    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
738    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
739  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
740    
741  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
742        CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
743       O             SlopeX, SlopeY,       O             SlopeX, SlopeY,
744       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
745         U             hTransLay, baseSlope, recipLambda,
746       U             dSigmaDr,       U             dSigmaDr,
747       I             dSigmaDx, dSigmaDy,       I             dSigmaDx, dSigmaDy,
748       I             ldd97_LrhoS,rC(k),k,       I             ldd97_LrhoS, locMixLayer, rC,
749       I             bi, bj, myThid )       I             kLow_S,
750         I             k, bi, bj, myTime, myIter, myThid )
751    
752  cph(  cph(
753  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 500  cph) Line 758  cph)
758  cph)  cph)
759    
760  #ifdef GM_NON_UNITY_DIAGONAL  #ifdef GM_NON_UNITY_DIAGONAL
761    c      IF ( GM_nonUnitDiag ) THEN
762          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
763           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
764            Kvy(i,j,k,bi,bj) =            Kvy(i,j,k,bi,bj) =
765    #ifdef ALLOW_KAPREDI_CONTROL
766         &     ( kapredi(i,j,k,bi,bj)
767    #else
768       &     ( GM_isopycK       &     ( GM_isopycK
769    #endif
770  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
771       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
772  #endif  #endif
773       &     )       &     )*taperFct(i,j)
      &     *taperFct(i,j)  
774           ENDDO           ENDDO
775          ENDDO          ENDDO
776  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 521  CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_b Line 783  CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_b
783            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 )
784           ENDDO           ENDDO
785          ENDDO          ENDDO
786    c      ENDIF
787  #endif /* GM_NON_UNITY_DIAGONAL */  #endif /* GM_NON_UNITY_DIAGONAL */
788    
789  #ifdef GM_EXTRA_DIAGONAL  #ifdef GM_EXTRA_DIAGONAL
# Line 529  CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_b Line 792  CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_b
792  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
793  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
794  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
795        IF (GM_ExtraDiag) THEN         IF ( GM_ExtraDiag ) THEN
796          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
797           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
798            Kvz(i,j,k,bi,bj) =            Kvz(i,j,k,bi,bj) =
799       &     ( GM_isopycK - GM_skewflx*GM_background_K  #ifdef ALLOW_KAPREDI_CONTROL
800         &     ( kapredi(i,j,k,bi,bj)
801    #else
802         &     ( GM_isopycK
803    #endif
804    #ifdef ALLOW_KAPGM_CONTROL
805         &     - GM_skewflx*kapgm(i,j,k,bi,bj)
806    #else
807         &     - GM_skewflx*GM_background_K
808    #endif
809  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
810       &     +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
811  #endif  #endif
812       &     )*SlopeY(i,j)*taperFct(i,j)       &     )*SlopeY(i,j)*taperFct(i,j)
813           ENDDO           ENDDO
814          ENDDO          ENDDO
815        ENDIF         ENDIF
816  #endif /* GM_EXTRA_DIAGONAL */  #endif /* GM_EXTRA_DIAGONAL */
817    
818  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
819        IF (doDiagRediFlx) THEN         IF (doDiagRediFlx) THEN
820  c       km1 = MAX(k-1,1)          km1 = MAX(k-1,1)
821          DO j=1,sNy+1          DO j=1,sNy+1
822           DO i=1,sNx           DO i=1,sNx
823  C         store in tmp1k Kvz_Redi  C         store in tmp1k Kvz_Redi
824    #ifdef ALLOW_KAPREDI_CONTROL
825              tmp1k(i,j) = ( kapredi(i,j,k,bi,bj)
826    #else
827            tmp1k(i,j) = ( GM_isopycK            tmp1k(i,j) = ( GM_isopycK
828    #endif
829  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
830       &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0       &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0
831  #endif  #endif
# Line 578  C-        Vertical gradients interpolate Line 854  C-        Vertical gradients interpolate
854           ENDDO           ENDDO
855          ENDDO          ENDDO
856          CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KvzTz', k,1,2,bi,bj,myThid)          CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KvzTz', k,1,2,bi,bj,myThid)
857        ENDIF         ENDIF
858  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
859    
860  #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */  C-- end 3rd  loop on vertical level index k
   
 C-- end 2nd loop on vertical level index k  
861        ENDDO        ENDDO
862    
863    #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
864    
865    
866  #ifdef GM_BOLUS_ADVEC  #ifdef GM_BOLUS_ADVEC
867        IF (GM_AdvForm) THEN        IF (GM_AdvForm) THEN
868          CALL GMREDI_CALC_PSI_B(         CALL GMREDI_CALC_PSI_B(
869       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
870       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
871       I             ldd97_LrhoW, ldd97_LrhoS,       I             ldd97_LrhoW, ldd97_LrhoS,
872       I             myThid )       I             myThid )
873        ENDIF        ENDIF
874  #endif  #endif
875    
# Line 641  C--   Time-average Line 917  C--   Time-average
917    
918    
919        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
920       I             bi, bj, iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
921       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
922       I             myThid )       I             bi, bj, myTime, myIter, myThid )
923  C     /==========================================================\  C     /==========================================================\
924  C     | SUBROUTINE GMREDI_CALC_TENSOR                            |  C     | SUBROUTINE GMREDI_CALC_TENSOR                            |
925  C     | o Calculate tensor elements for GM/Redi tensor.          |  C     | o Calculate tensor elements for GM/Redi tensor.          |
# Line 661  C Line 937  C
937        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
938        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
939        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
940        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
941          INTEGER bi, bj
942          _RL     myTime
943          INTEGER myIter
944        INTEGER myThid        INTEGER myThid
945  CEndOfInterface  CEndOfInterface
946    
       INTEGER i, j, k  
   
947  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
948    
949          INTEGER i, j, k
950    
951        DO k=1,Nr        DO k=1,Nr
952         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
953          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
Legend:
Removed from v.1.24  
changed lines
  Added in v.1.33
  ViewVC Help
Powered by ViewVC 1.1.22