/[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.26 by jmc, Thu May 24 22:34:38 2007 UTC revision 1.40 by jmc, Wed Jul 13 22:59:53 2011 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  #undef OLD_VISBECK_CALC  #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 24  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
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)
# Line 50  C     == Local variables == Line 71  C     == Local variables ==
71        _RL dSigmaDr(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)
       _RL maskp1, Kgm_tmp  
74        _RL ldd97_LrhoC(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL ldd97_LrhoC(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
75        _RL ldd97_LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL ldd97_LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
76        _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
77        _RL Cspd, LrhoInf, LrhoSup, fCoriLoc        _RL Cspd, LrhoInf, LrhoSup, fCoriLoc
78          _RL Kgm_tmp, isopycK, bolus_K
79    
80          INTEGER kLow_W (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
81          INTEGER kLow_S (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
82          _RL locMixLayer(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
83          _RL baseSlope  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
84          _RL hTransLay  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
85          _RL recipLambda(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
86          INTEGER  km1
87    #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
88          INTEGER kp1
89          _RL maskp1
90    #endif
91    
92  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
93  #ifdef OLD_VISBECK_CALC  #ifdef OLD_VISBECK_CALC
       _RL deltaH,zero_rs  
       PARAMETER(zero_rs=0.D0)  
       _RL N2,SN  
94        _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
95  #else  #else
96        _RL dSigmaH        _RL dSigmaH, dSigmaR
97        _RL deltaH, integrDepth        _RL Sloc, M2loc
       _RL Sloc, M2loc, SNloc  
 #endif  
98  #endif  #endif
99          _RL recipMaxSlope
100          _RL deltaH, integrDepth
101          _RL N2loc, SNloc
102    #endif /* GM_VISBECK_VARIABLE_K */
103    
104  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
105        LOGICAL  doDiagRediFlx        LOGICAL  doDiagRediFlx
106        LOGICAL  DIAGNOSTICS_IS_ON        LOGICAL  DIAGNOSTICS_IS_ON
107        EXTERNAL DIAGNOSTICS_IS_ON        EXTERNAL DIAGNOSTICS_IS_ON
108        INTEGER  km1  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
109        _RL dTdz        _RL dTdz
110        _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111  #endif  #endif
112    #endif /* ALLOW_DIAGNOSTICS */
113    
114  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
115    
# Line 103  C---+----1----+----2----+----3----+----4 Line 136  C---+----1----+----2----+----3----+----4
136  #endif  #endif
137    
138  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
139          recipMaxSlope = 0. _d 0
140          IF ( GM_Visbeck_maxSlope.GT.0. _d 0 ) THEN
141            recipMaxSlope = 1. _d 0 / GM_Visbeck_maxSlope
142          ENDIF
143        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
144         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
145          VisbeckK(i,j,bi,bj) = 0. _d 0          VisbeckK(i,j,bi,bj) = 0. _d 0
# Line 111  C---+----1----+----2----+----3----+----4 Line 148  C---+----1----+----2----+----3----+----4
148  #endif  #endif
149    
150  C--   set ldd97_Lrho (for tapering scheme ldd97):  C--   set ldd97_Lrho (for tapering scheme ldd97):
151        IF (GM_taper_scheme.EQ.'ldd97') THEN        IF ( GM_taper_scheme.EQ.'ldd97' .OR.
152         &     GM_taper_scheme.EQ.'fm07' ) THEN
153         Cspd = 2. _d 0         Cspd = 2. _d 0
154         LrhoInf = 15. _d 3         LrhoInf = 15. _d 3
155         LrhoSup = 100. _d 3         LrhoSup = 100. _d 3
# Line 128  C-     Tracer point location (center): Line 166  C-     Tracer point location (center):
166         ENDDO         ENDDO
167  C-     U point location (West):  C-     U point location (West):
168         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
169            kLow_W(1-Olx,j) = 0
170          ldd97_LrhoW(1-Olx,j) = LrhoSup          ldd97_LrhoW(1-Olx,j) = LrhoSup
171          DO i=1-Olx+1,sNx+Olx          DO i=1-Olx+1,sNx+Olx
172             kLow_W(i,j) = MIN(kLowC(i-1,j,bi,bj),kLowC(i,j,bi,bj))
173           fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))           fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))
174           IF (fCoriLoc.NE.0.) THEN           IF (fCoriLoc.NE.0.) THEN
175             ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)             ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)
# Line 141  C-     U point location (West): Line 181  C-     U point location (West):
181         ENDDO         ENDDO
182  C-     V point location (South):  C-     V point location (South):
183         DO i=1-Olx+1,sNx+Olx         DO i=1-Olx+1,sNx+Olx
184             kLow_S(i,1-Oly) = 0
185           ldd97_LrhoS(i,1-Oly) = LrhoSup           ldd97_LrhoS(i,1-Oly) = LrhoSup
186         ENDDO         ENDDO
187         DO j=1-Oly+1,sNy+Oly         DO j=1-Oly+1,sNy+Oly
188          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
189             kLow_S(i,j) = MIN(kLowC(i,j-1,bi,bj),kLowC(i,j,bi,bj))
190           fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))           fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))
191           IF (fCoriLoc.NE.0.) THEN           IF (fCoriLoc.NE.0.) THEN
192             ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)             ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)
# Line 164  C-    Just initialize to zero (not use a Line 206  C-    Just initialize to zero (not use a
206          ENDDO          ENDDO
207         ENDDO         ENDDO
208        ENDIF        ENDIF
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
   
       DO k=2,Nr  
 C-- 1rst loop on k : compute Tensor Coeff. at W points.  
209    
210    #ifdef GM_BOLUS_ADVEC
211          DO k=1,Nr
212           DO j=1-Oly,sNy+Oly
213            DO i=1-Olx,sNx+Olx
214             GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
215             GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
216            ENDDO
217           ENDDO
218          ENDDO
219    #endif /* GM_BOLUS_ADVEC */
220  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
221         kkey = (igmkey-1)*Nr + k        DO k=1,Nr
222         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
223          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
          SlopeX(i,j)       = 0. _d 0  
          SlopeY(i,j)       = 0. _d 0  
          dSigmaDx(i,j)     = 0. _d 0  
          dSigmaDy(i,j)     = 0. _d 0  
          dSigmaDr(i,j)     = 0. _d 0  
          SlopeSqr(i,j)     = 0. _d 0  
          taperFct(i,j)     = 0. _d 0  
224           Kwx(i,j,k,bi,bj)  = 0. _d 0           Kwx(i,j,k,bi,bj)  = 0. _d 0
225           Kwy(i,j,k,bi,bj)  = 0. _d 0           Kwy(i,j,k,bi,bj)  = 0. _d 0
226           Kwz(i,j,k,bi,bj)  = 0. _d 0           Kwz(i,j,k,bi,bj)  = 0. _d 0
# Line 191  C-- 1rst loop on k : compute Tensor Coef Line 232  C-- 1rst loop on k : compute Tensor Coef
232           Kuz(i,j,k,bi,bj)  = 0. _d 0           Kuz(i,j,k,bi,bj)  = 0. _d 0
233           Kvz(i,j,k,bi,bj)  = 0. _d 0           Kvz(i,j,k,bi,bj)  = 0. _d 0
234  # endif  # endif
 # ifdef GM_BOLUS_ADVEC  
          GM_PsiX(i,j,k,bi,bj)  = 0. _d 0  
          GM_PsiY(i,j,k,bi,bj)  = 0. _d 0  
 # endif  
235          ENDDO          ENDDO
236         ENDDO         ENDDO
237          ENDDO
238    #endif /* ALLOW_AUTODIFF_TAMC */
239    
240    C--   Initialise Mixed Layer related array:
241          DO j=1-Oly,sNy+Oly
242           DO i=1-Olx,sNx+Olx
243             hTransLay(i,j) = R_low(i,j,bi,bj)
244             baseSlope(i,j) =  0. _d 0
245             recipLambda(i,j) = 0. _d 0
246             locMixLayer(i,j) = 0. _d 0
247           ENDDO
248          ENDDO
249    #ifdef ALLOW_KPP
250          IF ( useKPP ) THEN
251           DO j=1-Oly,sNy+Oly
252            DO i=1-Olx,sNx+Olx
253             locMixLayer(i,j) = KPPhbl(i,j,bi,bj)
254            ENDDO
255           ENDDO
256          ELSE
257    #else
258          IF ( .TRUE. ) THEN
259  #endif  #endif
260           DO j=1-Oly,sNy+Oly
261            DO i=1-Olx,sNx+Olx
262             locMixLayer(i,j) = hMixLayer(i,j,bi,bj)
263            ENDDO
264           ENDDO
265          ENDIF
266    
267    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
268    C-- 1rst loop on k : compute Tensor Coeff. at W points.
269    
270          DO k=Nr,2,-1
271    
272    #ifdef ALLOW_AUTODIFF_TAMC
273           kkey = (igmkey-1)*Nr + k
274           DO j=1-Oly,sNy+Oly
275            DO i=1-Olx,sNx+Olx
276             SlopeX(i,j)       = 0. _d 0
277             SlopeY(i,j)       = 0. _d 0
278             dSigmaDx(i,j)     = 0. _d 0
279             dSigmaDy(i,j)     = 0. _d 0
280             dSigmaDr(i,j)     = 0. _d 0
281             SlopeSqr(i,j)     = 0. _d 0
282             taperFct(i,j)     = 0. _d 0
283            ENDDO
284           ENDDO
285    #endif /* ALLOW_AUTODIFF_TAMC */
286    
287         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
288          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
# Line 208  C      Gradient of Sigma at rVel points Line 293  C      Gradient of Sigma at rVel points
293           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)
294       &                       +sigmaY(i,j+1, k )+sigmaY(i,j, k )       &                       +sigmaY(i,j+1, k )+sigmaY(i,j, k )
295       &                      )*maskC(i,j,k,bi,bj)       &                      )*maskC(i,j,k,bi,bj)
296           dSigmaDr(i,j)=sigmaR(i,j,k)  c        dSigmaDr(i,j)=sigmaR(i,j,k)
297          ENDDO          ENDDO
298         ENDDO         ENDDO
299    
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
300  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
301  #ifndef OLD_VISBECK_CALC  #ifndef OLD_VISBECK_CALC
302         IF ( GM_Visbeck_alpha.GT.0. .AND.         IF ( GM_Visbeck_alpha.GT.0. .AND.
303       &      -rC(k-1).LT.GM_Visbeck_depth ) THEN       &      -rC(k-1).LT.GM_Visbeck_depth ) THEN
304    
305            DO j=1-Oly,sNy+Oly
306             DO i=1-Olx,sNx+Olx
307              dSigmaDr(i,j) = MIN( sigmaR(i,j,k), 0. _d 0 )
308             ENDDO
309            ENDDO
310    
311  C--     Depth average of f/sqrt(Ri) = M^2/N^2 * N  C--     Depth average of f/sqrt(Ri) = M^2/N^2 * N
312  C       M^2 and N^2 are horizontal & vertical gradient of buoyancy.  C       M^2 and N^2 are horizontal & vertical gradient of buoyancy.
313    
# Line 237  C       GM_Visbeck_depth, whatever is th Line 322  C       GM_Visbeck_depth, whatever is th
322             integrDepth = -rC( kLowC(i,j,bi,bj) )             integrDepth = -rC( kLowC(i,j,bi,bj) )
323  C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments  C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments
324             integrDepth = MIN( integrDepth, GM_Visbeck_depth )             integrDepth = MIN( integrDepth, GM_Visbeck_depth )
325    C-      to recover "old-visbeck" form with Visbeck_minDepth = Visbeck_depth
326               integrDepth = MAX( integrDepth, GM_Visbeck_minDepth )
327  C       Distance between level center above and the integration depth  C       Distance between level center above and the integration depth
328             deltaH = integrDepth + rC(k-1)             deltaH = integrDepth + rC(k-1)
329  C       If negative then we are below the integration level  C       If negative then we are below the integration level
# Line 246  C       If positive we limit this to the Line 333  C       If positive we limit this to the
333  C       Now we convert deltaH to a non-dimensional fraction  C       Now we convert deltaH to a non-dimensional fraction
334             deltaH = deltaH/( integrDepth+rC(1) )             deltaH = deltaH/( integrDepth+rC(1) )
335    
336  C--      compute: ( M^2 * S )^1/2   (= M^2 / N since S=M^2/N^2 )  C--      compute: ( M^2 * S )^1/2   (= S*N since S=M^2/N^2 )
337    C        a 5 points average gives a more "homogeneous" formulation
338    C        (same stencil and same weights as for dSigmaH calculation)
339               dSigmaR = ( dSigmaDr(i,j)*4. _d 0
340         &               + dSigmaDr(i-1,j)
341         &               + dSigmaDr(i+1,j)
342         &               + dSigmaDr(i,j-1)
343         &               + dSigmaDr(i,j+1)
344         &               )/( 4. _d 0
345         &                 + maskC(i-1,j,k,bi,bj)
346         &                 + maskC(i+1,j,k,bi,bj)
347         &                 + maskC(i,j-1,k,bi,bj)
348         &                 + maskC(i,j+1,k,bi,bj)
349         &                 )
350             dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)             dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)
351       &             + dSigmaDy(i,j)*dSigmaDy(i,j)       &             + dSigmaDy(i,j)*dSigmaDy(i,j)
352             IF ( dSigmaH .GT. 0. _d 0 ) THEN             IF ( dSigmaH .GT. 0. _d 0 ) THEN
353               dSigmaH = SQRT( dSigmaH )               dSigmaH = SQRT( dSigmaH )
354  C-       compute slope, limited by GM_maxSlope:  C-       compute slope, limited by GM_Visbeck_maxSlope:
355               IF ( -dSigmaDr(i,j).GT.dSigmaH*GM_rMaxSlope ) THEN               IF ( -dSigmaR.GT.dSigmaH*recipMaxSlope ) THEN
356                Sloc = dSigmaH / ( -dSigmaDr(i,j) )                Sloc = dSigmaH / ( -dSigmaR )
357               ELSE               ELSE
358                Sloc = GM_maxSlope                Sloc = GM_Visbeck_maxSlope
359                 ENDIF
360                 M2loc = gravity*recip_rhoConst*dSigmaH
361    c            SNloc = SQRT( Sloc*M2loc )
362                 N2loc = -gravity*recip_rhoConst*dSigmaR
363    c            N2loc = -gravity*recip_rhoConst*dSigmaDr(i,j)
364                 IF ( N2loc.GT.0. _d 0 ) THEN
365                   SNloc = Sloc*SQRT(N2loc)
366                 ELSE
367                   SNloc = 0. _d 0
368               ENDIF               ENDIF
              M2loc = Gravity*recip_RhoConst*dSigmaH  
              SNloc = SQRT( Sloc*M2loc )  
369             ELSE             ELSE
370               SNloc = 0. _d 0               SNloc = 0. _d 0
371             ENDIF             ENDIF
# Line 271  C-       compute slope, limited by GM_ma Line 378  C-       compute slope, limited by GM_ma
378         ENDIF         ENDIF
379  #endif /* ndef OLD_VISBECK_CALC */  #endif /* ndef OLD_VISBECK_CALC */
380  #endif /* GM_VISBECK_VARIABLE_K */  #endif /* GM_VISBECK_VARIABLE_K */
381           DO j=1-Oly,sNy+Oly
382            DO i=1-Olx,sNx+Olx
383             dSigmaDr(i,j)=sigmaR(i,j,k)
384            ENDDO
385           ENDDO
386    
387    #ifdef ALLOW_AUTODIFF_TAMC
388    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
389    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
390    CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
391    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
392    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
393    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
394    #endif /* ALLOW_AUTODIFF_TAMC */
395    
396  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
397         CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
398       O             SlopeX, SlopeY,       O             SlopeX, SlopeY,
399       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
400         U             hTransLay, baseSlope, recipLambda,
401       U             dSigmaDr,       U             dSigmaDr,
402       I             dSigmaDx, dSigmaDy,       I             dSigmaDx, dSigmaDy,
403       I             ldd97_LrhoC,rF(k),k,       I             ldd97_LrhoC, locMixLayer, rF,
404       I             bi, bj, myThid )       I             kLowC(1-Olx,1-Oly,bi,bj),
405         I             k, bi, bj, myTime, myIter, myThid )
406    
407         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
408          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
# Line 298  CADJ STORE dSigmaDr(:,:)     = comlev1_b Line 421  CADJ STORE dSigmaDr(:,:)     = comlev1_b
421  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
422  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
423    
424    C      Components of Redi/GM tensor
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
427  C      Components of Redi/GM tensor            Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
428           Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)            Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
429           Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)            Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
430           Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)          ENDDO
431           ENDDO
432    
433  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
434  #ifdef OLD_VISBECK_CALC  #ifdef OLD_VISBECK_CALC
435           DO j=1-Oly+1,sNy+Oly-1
436            DO i=1-Olx+1,sNx+Olx-1
437    
438  C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K  C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K
439  C           but do not know if *taperFct (or **2 ?) is necessary  C           but do not know if *taperFct (or **2 ?) is necessary
440          Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)          Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)
441    
442  C--     Depth average of M^2/N^2 * N  C--     Depth average of M^2/N^2 * N
# Line 319  C       which is used in the "variable K Line 446  C       which is used in the "variable K
446  C       Distance between interface above layer and the integration depth  C       Distance between interface above layer and the integration depth
447          deltaH=abs(GM_Visbeck_depth)-abs(rF(k))          deltaH=abs(GM_Visbeck_depth)-abs(rF(k))
448  C       If positive we limit this to the layer thickness  C       If positive we limit this to the layer thickness
449          deltaH=min(deltaH,drF(k))          integrDepth = drF(k)
450            deltaH=min(deltaH,integrDepth)
451  C       If negative then we are below the integration level  C       If negative then we are below the integration level
452          deltaH=max(deltaH,zero_rs)          deltaH=max(deltaH, 0. _d 0)
453  C       Now we convert deltaH to a non-dimensional fraction  C       Now we convert deltaH to a non-dimensional fraction
454          deltaH=deltaH/GM_Visbeck_depth          deltaH=deltaH/GM_Visbeck_depth
455    
         IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.  
456          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
457           N2= -Gravity*recip_RhoConst*dSigmaDr(i,j)           N2loc = -gravity*recip_rhoConst*dSigmaDr(i,j)
458           SN=sqrt(Ssq(i,j)*N2)           SNloc = SQRT(Ssq(i,j)*N2loc )
459           VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH           VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
460       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN       &       +deltaH*GM_Visbeck_alpha
461         &              *GM_Visbeck_length*GM_Visbeck_length*SNloc
462          ENDIF          ENDIF
463    
 #endif /* OLD_VISBECK_CALC */  
 #endif /* GM_VISBECK_VARIABLE_K */  
464          ENDDO          ENDDO
465         ENDDO         ENDDO
466    #endif /* OLD_VISBECK_CALC */
467    #endif /* GM_VISBECK_VARIABLE_K */
468    
469  C-- end 1rst loop on vertical level index k  C-- end 1rst loop on vertical level index k
470        ENDDO        ENDDO
# Line 346  C-- end 1rst loop on vertical level inde Line 474  C-- end 1rst loop on vertical level inde
474  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
475  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
476  #endif  #endif
477        IF ( GM_Visbeck_alpha.NE.0. ) THEN        IF ( GM_Visbeck_alpha.GT.0. ) THEN
478  C-     Limit range that KapGM can take  C-     Limit range that KapGM can take
479         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
480          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
481           VisbeckK(i,j,bi,bj)=           VisbeckK(i,j,bi,bj)=
482       &       MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)       &       MIN( MAX( VisbeckK(i,j,bi,bj), GM_Visbeck_minVal_K ),
483         &            GM_Visbeck_maxVal_K )
484          ENDDO          ENDDO
485         ENDDO         ENDDO
486        ENDIF        ENDIF
# Line 362  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1 Line 491  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1
491  cph)  cph)
492  #endif /* GM_VISBECK_VARIABLE_K */  #endif /* GM_VISBECK_VARIABLE_K */
493    
494    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.  
495        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  
   
496  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
497         kkey = (igmkey-1)*Nr + k         kkey = (igmkey-1)*Nr + k
498  #if (defined (GM_NON_UNITY_DIAGONAL) || \  # if (defined (GM_NON_UNITY_DIAGONAL) || \
499       defined (GM_VISBECK_VARIABLE_K))        defined (GM_VISBECK_VARIABLE_K))
500  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
501  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
502  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
503    # endif
504  #endif  #endif
505  #endif         km1 = MAX(k-1,1)
506           isopycK = GM_isopycK
507  C-    express the Tensor in term of Diffusivity (= m**2 / s )       &         *(GM_isoFac1d(km1)+GM_isoFac1d(k))*op5
508           bolus_K = GM_background_K
509         &         *(GM_bolFac1d(km1)+GM_bolFac1d(k))*op5
510         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
511          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
512    #ifdef ALLOW_KAPREDI_CONTROL
513             Kgm_tmp = kapredi(i,j,k,bi,bj)
514    #else
515             Kgm_tmp = isopycK*GM_isoFac2d(i,j,bi,bj)
516    #endif
517  #ifdef ALLOW_KAPGM_CONTROL  #ifdef ALLOW_KAPGM_CONTROL
518           Kgm_tmp = GM_isopycK + GM_skewflx*kapgm(i,j,k,bi,bj)       &           + GM_skewflx*kapgm(i,j,k,bi,bj)
519  #else  #else
520           Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K       &           + GM_skewflx*bolus_K*GM_bolFac2d(i,j,bi,bj)
521  #endif  #endif
522  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
523       &           + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)       &           + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)
524  #endif  #endif
525           Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)           Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)
526           Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)           Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)
527           Kwz(i,j,k,bi,bj)= ( GM_isopycK  #ifdef ALLOW_KAPREDI_CONTROL
528             Kwz(i,j,k,bi,bj)= ( kapredi(i,j,k,bi,bj)
529    #else
530             Kwz(i,j,k,bi,bj)= ( isopycK*GM_isoFac2d(i,j,bi,bj)
531    #endif
532  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
533       &                     + VisbeckK(i,j,bi,bj)       &                     + VisbeckK(i,j,bi,bj)
534  #endif  #endif
535       &                     )*Kwz(i,j,k,bi,bj)       &                     )*Kwz(i,j,k,bi,bj)
536          ENDDO          ENDDO
537         ENDDO         ENDDO
538          ENDDO
539    
540    #ifdef ALLOW_DIAGNOSTICS
541          IF ( useDiagnostics .AND. GM_taper_scheme.EQ.'fm07' ) THEN
542           CALL DIAGNOSTICS_FILL( hTransLay, 'GM_hTrsL', 0,1,2,bi,bj,myThid)
543           CALL DIAGNOSTICS_FILL( baseSlope, 'GM_baseS', 0,1,2,bi,bj,myThid)
544           CALL DIAGNOSTICS_FILL(recipLambda,'GM_rLamb', 0,1,2,bi,bj,myThid)
545          ENDIF
546    #endif /* ALLOW_DIAGNOSTICS */
547    
548    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
549    C--   Calculate Stream-Functions used in Advective Form:
550    
551    #ifdef GM_BOLUS_ADVEC
552          IF (GM_AdvForm) THEN
553    #ifdef GM_BOLUS_BVP
554           IF (GM_UseBVP) THEN
555            CALL GMREDI_CALC_PSI_BVP(
556         I             bi, bj, iMin, iMax, jMin, jMax,
557         I             sigmaX, sigmaY, sigmaR,
558         I             myThid )
559           ELSE
560    #endif
561            CALL GMREDI_CALC_PSI_B(
562         I              bi, bj, iMin, iMax, jMin, jMax,
563         I              sigmaX, sigmaY, sigmaR,
564         I              ldd97_LrhoW, ldd97_LrhoS,
565         I              myThid )
566    #ifdef GM_BOLUS_BVP
567           ENDIF
568    #endif
569          ENDIF
570    #endif
571    
572    #ifndef GM_EXCLUDE_SUBMESO
573          IF ( GM_useSubMeso .AND. GM_AdvForm ) THEN
574            CALL SUBMESO_CALC_PSI(
575         I              bi, bj, iMin, iMax, jMin, jMax,
576         I              sigmaX, sigmaY, sigmaR,
577         I              locMixLayer,
578         I              myIter, myThid )
579          ENDIF
580    #endif /* ndef GM_EXCLUDE_SUBMESO */
581    
582  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
583    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
584    C-- 2nd  k loop : compute Tensor Coeff. at U point
585    
586    #ifdef ALLOW_KPP
587          IF ( useKPP ) THEN
588           DO j=1-Oly,sNy+Oly
589            DO i=2-Olx,sNx+Olx
590             locMixLayer(i,j) = ( KPPhbl(i-1,j,bi,bj)
591         &                      + KPPhbl( i ,j,bi,bj) )*op5
592            ENDDO
593           ENDDO
594          ELSE
595    #else
596          IF ( .TRUE. ) THEN
597    #endif
598           DO j=1-Oly,sNy+Oly
599            DO i=2-Olx,sNx+Olx
600             locMixLayer(i,j) = ( hMixLayer(i-1,j,bi,bj)
601         &                      + hMixLayer( i ,j,bi,bj) )*op5
602            ENDDO
603           ENDDO
604          ENDIF
605          DO j=1-Oly,sNy+Oly
606           DO i=1-Olx,sNx+Olx
607             hTransLay(i,j) =  0.
608             baseSlope(i,j) =  0.
609             recipLambda(i,j)= 0.
610           ENDDO
611           DO i=2-Olx,sNx+Olx
612             hTransLay(i,j) = MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
613           ENDDO
614          ENDDO
615    
616          DO k=Nr,1,-1
617           kp1 = MIN(Nr,k+1)
618           maskp1 = 1. _d 0
619           IF (k.GE.Nr) maskp1 = 0. _d 0
620    #ifdef ALLOW_AUTODIFF_TAMC
621           kkey = (igmkey-1)*Nr + k
622    #endif
623    
624  C     Gradient of Sigma at U points  C     Gradient of Sigma at U points
625         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
# Line 422  C     Gradient of Sigma at U points Line 639  C     Gradient of Sigma at U points
639  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
640  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
641  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
642  CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
643    CADJ STORE locMixLayer(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
644    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
645    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
646    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
647  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
648    
649  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
650         CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
651       O             SlopeX, SlopeY,       O             SlopeX, SlopeY,
652       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
653         U             hTransLay, baseSlope, recipLambda,
654       U             dSigmaDr,       U             dSigmaDr,
655       I             dSigmaDx, dSigmaDy,       I             dSigmaDx, dSigmaDy,
656       I             ldd97_LrhoW,rC(k),k,       I             ldd97_LrhoW, locMixLayer, rC,
657       I             bi, bj, myThid )       I             kLow_W,
658         I             k, bi, bj, myTime, myIter, myThid )
659    
 cph( NEW  
660  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
 cph(  
661  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
662  CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
 cph)  
663  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
 cph)  
664    
665  #ifdef GM_NON_UNITY_DIAGONAL  #ifdef GM_NON_UNITY_DIAGONAL
666  c      IF ( GM_nonUnitDiag ) THEN  c      IF ( GM_nonUnitDiag ) THEN
667          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
668           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
669            Kux(i,j,k,bi,bj) =            Kux(i,j,k,bi,bj) =
670       &     ( GM_isopycK  #ifdef ALLOW_KAPREDI_CONTROL
671         &     ( kapredi(i,j,k,bi,bj)
672    #else
673         &     ( GM_isopycK*GM_isoFac1d(k)
674         &        *op5*(GM_isoFac2d(i-1,j,bi,bj)+GM_isoFac2d(i,j,bi,bj))
675    #endif
676  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
677       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
678  #endif  #endif
679       &     )       &     )*taperFct(i,j)
      &     *taperFct(i,j)  
680           ENDDO           ENDDO
681          ENDDO          ENDDO
682  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 475  c      ENDIF Line 698  c      ENDIF
698  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
699  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
700  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
701         IF (GM_ExtraDiag) THEN         IF ( GM_ExtraDiag ) THEN
702          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
703           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
704            Kuz(i,j,k,bi,bj) =            Kuz(i,j,k,bi,bj) =
705    #ifdef ALLOW_KAPREDI_CONTROL
706         &     ( kapredi(i,j,k,bi,bj)
707    #else
708         &     ( GM_isopycK*GM_isoFac1d(k)
709         &        *op5*(GM_isoFac2d(i-1,j,bi,bj)+GM_isoFac2d(i,j,bi,bj))
710    #endif
711  #ifdef ALLOW_KAPGM_CONTROL  #ifdef ALLOW_KAPGM_CONTROL
712       &     ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)       &     - GM_skewflx*kapgm(i,j,k,bi,bj)
713  #else  #else
714       &     ( GM_isopycK - GM_skewflx*GM_background_K       &     - GM_skewflx*GM_background_K*GM_bolFac1d(k)
715         &        *op5*(GM_bolFac2d(i-1,j,bi,bj)+GM_bolFac2d(i,j,bi,bj))
716  #endif  #endif
717  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
718       &     +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
# Line 499  CADJ STORE taperFct(:,:)     = comlev1_b Line 729  CADJ STORE taperFct(:,:)     = comlev1_b
729          DO j=1,sNy          DO j=1,sNy
730           DO i=1,sNx+1           DO i=1,sNx+1
731  C         store in tmp1k Kuz_Redi  C         store in tmp1k Kuz_Redi
732            tmp1k(i,j) = ( GM_isopycK  #ifdef ALLOW_KAPREDI_CONTROL
733              tmp1k(i,j) = ( kapredi(i,j,k,bi,bj)
734    #else
735              tmp1k(i,j) = ( GM_isopycK*GM_isoFac1d(k)
736         &        *op5*(GM_isoFac2d(i-1,j,bi,bj)+GM_isoFac2d(i,j,bi,bj))
737    #endif
738  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
739       &     +(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
740  #endif  #endif
# Line 531  C-        Vertical gradients interpolate Line 766  C-        Vertical gradients interpolate
766         ENDIF         ENDIF
767  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
768    
769    C-- end 2nd  loop on vertical level index k
770          ENDDO
771    
772    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
773    C-- 3rd  k loop : compute Tensor Coeff. at V point
774    
775    #ifdef ALLOW_KPP
776          IF ( useKPP ) THEN
777           DO j=2-Oly,sNy+Oly
778            DO i=1-Olx,sNx+Olx
779             locMixLayer(i,j) = ( KPPhbl(i,j-1,bi,bj)
780         &                      + KPPhbl(i, j ,bi,bj) )*op5
781            ENDDO
782           ENDDO
783          ELSE
784    #else
785          IF ( .TRUE. ) THEN
786    #endif
787           DO j=2-Oly,sNy+Oly
788            DO i=1-Olx,sNx+Olx
789             locMixLayer(i,j) = ( hMixLayer(i,j-1,bi,bj)
790         &                      + hMixLayer(i, j ,bi,bj) )*op5
791            ENDDO
792           ENDDO
793          ENDIF
794          DO j=1-Oly,sNy+Oly
795           DO i=1-Olx,sNx+Olx
796             hTransLay(i,j) =  0.
797             baseSlope(i,j) =  0.
798             recipLambda(i,j)= 0.
799           ENDDO
800          ENDDO
801          DO j=2-Oly,sNy+Oly
802           DO i=1-Olx,sNx+Olx
803             hTransLay(i,j) = MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
804           ENDDO
805          ENDDO
806    
807  C     Gradient of Sigma at V points  C     Gradient of Sigma at V points
808          DO k=Nr,1,-1
809           kp1 = MIN(Nr,k+1)
810           maskp1 = 1. _d 0
811           IF (k.GE.Nr) maskp1 = 0. _d 0
812    #ifdef ALLOW_AUTODIFF_TAMC
813           kkey = (igmkey-1)*Nr + k
814    #endif
815    
816         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
817          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
818           dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)           dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
# Line 548  C     Gradient of Sigma at V points Line 829  C     Gradient of Sigma at V points
829  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
830  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
831  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
832  CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
833    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
834    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
835    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
836  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
837    
838  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
839         CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
840       O             SlopeX, SlopeY,       O             SlopeX, SlopeY,
841       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
842         U             hTransLay, baseSlope, recipLambda,
843       U             dSigmaDr,       U             dSigmaDr,
844       I             dSigmaDx, dSigmaDy,       I             dSigmaDx, dSigmaDy,
845       I             ldd97_LrhoS,rC(k),k,       I             ldd97_LrhoS, locMixLayer, rC,
846       I             bi, bj, myThid )       I             kLow_S,
847         I             k, bi, bj, myTime, myIter, myThid )
848    
849  cph(  cph(
850  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 573  c      IF ( GM_nonUnitDiag ) THEN Line 859  c      IF ( GM_nonUnitDiag ) THEN
859          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
860           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
861            Kvy(i,j,k,bi,bj) =            Kvy(i,j,k,bi,bj) =
862       &     ( GM_isopycK  #ifdef ALLOW_KAPREDI_CONTROL
863         &     ( kapredi(i,j,k,bi,bj)
864    #else
865         &     ( GM_isopycK*GM_isoFac1d(k)
866         &        *op5*(GM_isoFac2d(i,j-1,bi,bj)+GM_isoFac2d(i,j,bi,bj))
867    #endif
868  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
869       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
870  #endif  #endif
871       &     )       &     )*taperFct(i,j)
      &     *taperFct(i,j)  
872           ENDDO           ENDDO
873          ENDDO          ENDDO
874  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 600  c      ENDIF Line 890  c      ENDIF
890  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
891  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
892  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
893         IF (GM_ExtraDiag) THEN         IF ( GM_ExtraDiag ) THEN
894          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
895           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
896            Kvz(i,j,k,bi,bj) =            Kvz(i,j,k,bi,bj) =
897    #ifdef ALLOW_KAPREDI_CONTROL
898         &     ( kapredi(i,j,k,bi,bj)
899    #else
900         &     ( GM_isopycK*GM_isoFac1d(k)
901         &        *op5*(GM_isoFac2d(i,j-1,bi,bj)+GM_isoFac2d(i,j,bi,bj))
902    #endif
903  #ifdef ALLOW_KAPGM_CONTROL  #ifdef ALLOW_KAPGM_CONTROL
904       &     ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)       &     - GM_skewflx*kapgm(i,j,k,bi,bj)
905  #else  #else
906       &     ( GM_isopycK - GM_skewflx*GM_background_K       &     - GM_skewflx*GM_background_K*GM_bolFac1d(k)
907         &        *op5*(GM_bolFac2d(i,j-1,bi,bj)+GM_bolFac2d(i,j,bi,bj))
908  #endif  #endif
909  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
910       &     +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
# Line 620  CADJ STORE taperFct(:,:)     = comlev1_b Line 917  CADJ STORE taperFct(:,:)     = comlev1_b
917    
918  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
919         IF (doDiagRediFlx) THEN         IF (doDiagRediFlx) THEN
920  c       km1 = MAX(k-1,1)          km1 = MAX(k-1,1)
921          DO j=1,sNy+1          DO j=1,sNy+1
922           DO i=1,sNx           DO i=1,sNx
923  C         store in tmp1k Kvz_Redi  C         store in tmp1k Kvz_Redi
924            tmp1k(i,j) = ( GM_isopycK  #ifdef ALLOW_KAPREDI_CONTROL
925              tmp1k(i,j) = ( kapredi(i,j,k,bi,bj)
926    #else
927              tmp1k(i,j) = ( GM_isopycK*GM_isoFac1d(k)
928         &        *op5*(GM_isoFac2d(i,j-1,bi,bj)+GM_isoFac2d(i,j,bi,bj))
929    #endif
930  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
931       &     +(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
932  #endif  #endif
# Line 656  C-        Vertical gradients interpolate Line 958  C-        Vertical gradients interpolate
958         ENDIF         ENDIF
959  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
960    
961  #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  
962        ENDDO        ENDDO
963    
964    #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
 #ifdef GM_BOLUS_ADVEC  
       IF (GM_AdvForm) THEN  
        CALL GMREDI_CALC_PSI_B(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             sigmaX, sigmaY, sigmaR,  
      I             ldd97_LrhoW, ldd97_LrhoS,  
      I             myThid )  
       ENDIF  
 #endif  
965    
966  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
967  C--   Time-average  C--   Time-average
# Line 696  C--   Time-average Line 987  C--   Time-average
987       &                          deltaTclock, bi, bj, myThid )       &                          deltaTclock, bi, bj, myThid )
988         ENDIF         ENDIF
989  #endif  #endif
990         DO k=1,Nr         GM_timeAve(bi,bj) = GM_timeAve(bi,bj)+deltaTclock
          GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock  
        ENDDO  
991    
992        ENDIF        ENDIF
993  #endif /* ALLOW_TIMEAVE */  #endif /* ALLOW_TIMEAVE */
# Line 714  C--   Time-average Line 1003  C--   Time-average
1003        RETURN        RETURN
1004        END        END
1005    
1006    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1007    
1008    CBOP
1009    C     !ROUTINE: GMREDI_CALC_TENSOR_DUMMY
1010    C     !INTERFACE:
1011        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
1012       I             bi, bj, iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
1013       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
1014       I             myThid )       I             bi, bj, myTime, myIter, myThid )
1015  C     /==========================================================\  
1016  C     | SUBROUTINE GMREDI_CALC_TENSOR                            |  C     !DESCRIPTION: \bv
1017  C     | o Calculate tensor elements for GM/Redi tensor.          |  C     *==========================================================*
1018  C     |==========================================================|  C     | SUBROUTINE GMREDI_CALC_TENSOR_DUMMY
1019  C     \==========================================================/  C     | o Calculate tensor elements for GM/Redi tensor.
1020    C     *==========================================================*
1021    C     \ev
1022    
1023    C     !USES:
1024        IMPLICIT NONE        IMPLICIT NONE
1025    
1026  C     == Global variables ==  C     == Global variables ==
# Line 731  C     == Global variables == Line 1028  C     == Global variables ==
1028  #include "EEPARAMS.h"  #include "EEPARAMS.h"
1029  #include "GMREDI.h"  #include "GMREDI.h"
1030    
1031  C     == Routine arguments ==  C     !INPUT/OUTPUT PARAMETERS:
 C  
1032        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
1033        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
1034        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
1035        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
1036          INTEGER bi, bj
1037          _RL     myTime
1038          INTEGER myIter
1039        INTEGER myThid        INTEGER myThid
1040  CEndOfInterface  CEOP
   
       INTEGER i, j, k  
1041    
1042  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
1043    C     !LOCAL VARIABLES:
1044          INTEGER i, j, k
1045    
1046        DO k=1,Nr        DO k=1,Nr
1047         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
Legend:
Removed from v.1.26  
changed lines
  Added in v.1.40
  ViewVC Help
Powered by ViewVC 1.1.22