/[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.10 by heimbach, Sun Mar 24 02:33:16 2002 UTC revision 1.25 by heimbach, Wed Feb 7 00:01:15 2007 UTC
# Line 22  C     == Global variables == Line 22  C     == Global variables ==
22  #include "EEPARAMS.h"  #include "EEPARAMS.h"
23  #include "PARAMS.h"  #include "PARAMS.h"
24  #include "GMREDI.h"  #include "GMREDI.h"
25  #include "GMREDI_DIAGS.h"  #include "GMREDI_TAVE.h"
26    
27  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
28  #include "tamc.h"  #include "tamc.h"
# Line 41  CEndOfInterface Line 41  CEndOfInterface
41  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
42    
43  C     == Local variables ==  C     == Local variables ==
44        INTEGER i,j,k,km1,kp1        INTEGER i,j,k,kp1
45        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
46        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47        _RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48          _RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49          _RL dSigmaDr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51        _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
52        _RL maskp1, maskm1, Kgm_tmp        _RL maskp1, Kgm_tmp
53          _RL ldd97_LrhoC(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
54          _RL ldd97_LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
55          _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56          _RL Cspd, LrhoInf, LrhoSup, fCoriLoc
57    
58  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
59        _RS deltaH,zero_rs        _RL deltaH,zero_rs
60        PARAMETER(zero_rs=0.)        PARAMETER(zero_rs=0.D0)
61        _RL N2,SN        _RL N2,SN
62        _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
63  #endif  #endif
64    
65    #ifdef ALLOW_DIAGNOSTICS
66          LOGICAL  doDiagRediFlx
67          LOGICAL  DIAGNOSTICS_IS_ON
68          EXTERNAL DIAGNOSTICS_IS_ON
69          INTEGER  km1
70          _RL dTdz
71          _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72    #endif
73    
74    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
75    
76  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
77            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
78            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
# Line 64  C     == Local variables == Line 81  C     == Local variables ==
81            act3 = myThid - 1            act3 = myThid - 1
82            max3 = nTx*nTy            max3 = nTx*nTy
83            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
84            ikey = (act1 + 1) + act2*max1            igmkey = (act1 + 1) + act2*max1
85       &                      + act3*max1*max2       &                      + act3*max1*max2
86       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
87  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
88    
89    #ifdef ALLOW_DIAGNOSTICS
90          doDiagRediFlx = .FALSE.
91          IF ( useDiagnostics ) THEN
92            doDiagRediFlx = DIAGNOSTICS_IS_ON('GM_KuzTz', myThid )
93            doDiagRediFlx = doDiagRediFlx .OR.
94         &                  DIAGNOSTICS_IS_ON('GM_KvzTz', myThid )
95          ENDIF
96    #endif
97        
98    #ifdef GM_VISBECK_VARIABLE_K
99          DO j=1-Oly,sNy+Oly
100           DO i=1-Olx,sNx+Olx
101            VisbeckK(i,j,bi,bj) = 0. _d 0
102           ENDDO
103          ENDDO
104    #endif
105    
106    C--   set ldd97_Lrho (for tapering scheme ldd97):
107          IF (GM_taper_scheme.EQ.'ldd97') THEN
108           Cspd = 2. _d 0
109           LrhoInf = 15. _d 3
110           LrhoSup = 100. _d 3
111    C-     Tracer point location (center):
112           DO j=1-Oly,sNy+Oly
113            DO i=1-Olx,sNx+Olx
114             IF (fCori(i,j,bi,bj).NE.0.) THEN
115               ldd97_LrhoC(i,j) = Cspd/ABS(fCori(i,j,bi,bj))
116             ELSE
117               ldd97_LrhoC(i,j) = LrhoSup
118             ENDIF
119             ldd97_LrhoC(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoC(i,j),LrhoSup))
120            ENDDO
121           ENDDO
122    C-     U point location (West):
123           DO j=1-Oly,sNy+Oly
124            ldd97_LrhoW(1-Olx,j) = LrhoSup
125            DO i=1-Olx+1,sNx+Olx
126             fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))
127             IF (fCoriLoc.NE.0.) THEN
128               ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)
129             ELSE
130               ldd97_LrhoW(i,j) = LrhoSup
131             ENDIF
132             ldd97_LrhoW(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoW(i,j),LrhoSup))
133            ENDDO
134           ENDDO
135    C-     V point location (South):
136           DO i=1-Olx+1,sNx+Olx
137             ldd97_LrhoS(i,1-Oly) = LrhoSup
138           ENDDO
139           DO j=1-Oly+1,sNy+Oly
140            DO i=1-Olx,sNx+Olx
141             fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))
142             IF (fCoriLoc.NE.0.) THEN
143               ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)
144             ELSE
145               ldd97_LrhoS(i,j) = LrhoSup
146             ENDIF
147             ldd97_LrhoS(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoS(i,j),LrhoSup))
148            ENDDO
149           ENDDO
150          ELSE
151    C-    Just initialize to zero (not use anyway)
152           DO j=1-Oly,sNy+Oly
153            DO i=1-Olx,sNx+Olx
154              ldd97_LrhoC(i,j) = 0. _d 0
155              ldd97_LrhoW(i,j) = 0. _d 0
156              ldd97_LrhoS(i,j) = 0. _d 0
157            ENDDO
158           ENDDO
159          ENDIF
160    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
161    
162        DO k=2,Nr        DO k=2,Nr
163  C-- 1rst loop on k : compute Tensor Coeff. at W points.  C-- 1rst loop on k : compute Tensor Coeff. at W points.
        km1 = MAX(1,k-1)  
        maskm1 = 1. _d 0  
        IF (k.LE.1) maskm1 = 0. _d 0  
164    
165  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
166         kkey = (ikey-1)*Nr + k         kkey = (igmkey-1)*Nr + k
167         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
168          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
169           SlopeX(i,j)       = 0. _d 0           SlopeX(i,j)       = 0. _d 0
170           SlopeY(i,j)       = 0. _d 0           SlopeY(i,j)       = 0. _d 0
171           dSigmaDrReal(i,j) = 0. _d 0           dSigmaDx(i,j)     = 0. _d 0
172             dSigmaDy(i,j)     = 0. _d 0
173             dSigmaDr(i,j)     = 0. _d 0
174           SlopeSqr(i,j)     = 0. _d 0           SlopeSqr(i,j)     = 0. _d 0
175           taperFct(i,j)     = 0. _d 0           taperFct(i,j)     = 0. _d 0
176           Kwx(i,j,k,bi,bj)  = 0. _d 0           Kwx(i,j,k,bi,bj)  = 0. _d 0
177           Kwy(i,j,k,bi,bj)  = 0. _d 0           Kwy(i,j,k,bi,bj)  = 0. _d 0
178           Kwz(i,j,k,bi,bj)  = 0. _d 0           Kwz(i,j,k,bi,bj)  = 0. _d 0
179    # ifdef GM_NON_UNITY_DIAGONAL
180             Kux(i,j,k,bi,bj)  = 0. _d 0
181             Kvy(i,j,k,bi,bj)  = 0. _d 0
182    # endif
183    # ifdef GM_EXTRA_DIAGONAL
184             Kuz(i,j,k,bi,bj)  = 0. _d 0
185             Kvz(i,j,k,bi,bj)  = 0. _d 0
186    # endif
187    # ifdef GM_BOLUS_ADVEC
188             GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
189             GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
190    # endif
191          ENDDO          ENDDO
192         ENDDO         ENDDO
193  #endif  #endif
194    
195        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
196         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
   
197  C      Gradient of Sigma at rVel points  C      Gradient of Sigma at rVel points
198          SlopeX(i,j)=0.25*( sigmaX(i+1, j ,km1) +sigmaX(i,j,km1)          dSigmaDx(i,j)=op25*( sigmaX(i+1, j ,k-1) +sigmaX(i,j,k-1)
199       &                    +sigmaX(i+1, j , k ) +sigmaX(i,j, k ) )       &                    +sigmaX(i+1, j , k ) +sigmaX(i,j, k ) )
200       &                  *maskC(i,j,k,bi,bj)*maskm1       &                  *maskC(i,j,k,bi,bj)
201          SlopeY(i,j)=0.25*( sigmaY( i ,j+1,km1) +sigmaY(i,j,km1)          dSigmaDy(i,j)=op25*( sigmaY( i ,j+1,k-1) +sigmaY(i,j,k-1)
202       &                    +sigmaY( i ,j+1, k ) +sigmaY(i,j, k ) )       &                    +sigmaY( i ,j+1, k ) +sigmaY(i,j, k ) )
203       &                  *maskC(i,j,k,bi,bj)*maskm1       &                  *maskC(i,j,k,bi,bj)
204          dSigmaDrReal(i,j)=sigmaR(i,j,k)*maskm1          dSigmaDr(i,j)=sigmaR(i,j,k)
   
205         ENDDO         ENDDO
206        ENDDO        ENDDO
207    
208  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
209  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
210  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
211  CADJ STORE dsigmadrreal(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
212  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
213    
214  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
215        CALL GMREDI_SLOPE_LIMIT(        CALL GMREDI_SLOPE_LIMIT(
216       U             dSigmadRReal,       O             SlopeX, SlopeY,
      I             rF(K),  
      U             SlopeX, SlopeY,  
217       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
218         U             dSigmaDr,
219         I             dSigmaDx, dSigmaDy,
220         I             ldd97_LrhoC,rF(k),k,
221       I             bi, bj, myThid )       I             bi, bj, myThid )
222    
223        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
224         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
225    
226  C       Mask Iso-neutral slopes  C       Mask Iso-neutral slopes
227          SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)*maskm1          SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)
228          SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)*maskm1          SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)
229          SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)*maskm1          SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)
230    
231         ENDDO         ENDDO
232        ENDDO        ENDDO
# Line 135  C       Mask Iso-neutral slopes Line 235  C       Mask Iso-neutral slopes
235  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
236  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
237  CADJ STORE SlopeSqr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeSqr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
238    CADJ STORE dSigmaDr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
239    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
240  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
241    
242        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
# Line 148  C       Components of Redi/GM tensor Line 250  C       Components of Redi/GM tensor
250  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
251    
252  C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K  C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K
253  C           but don't know if *taperFct (or **2 ?) is necessary  C           but do not know if *taperFct (or **2 ?) is necessary
254          Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)          Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)
255    
256  C--     Depth average of M^2/N^2 * N  C--     Depth average of M^2/N^2 * N
# Line 165  C       Now we convert deltaH to a non-d Line 267  C       Now we convert deltaH to a non-d
267          deltaH=deltaH/GM_Visbeck_depth          deltaH=deltaH/GM_Visbeck_depth
268    
269          IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.          IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.
270          IF (Ssq(i,j).NE.0.) THEN          IF ( Ssq(i,j).NE.0. .AND. dSigmaDr(i,j).NE.0. ) THEN
271           N2= -Gravity*recip_Rhonil*dSigmaDrReal(i,j)           N2= -Gravity*recip_RhoConst*dSigmaDr(i,j)
272           SN=sqrt(Ssq(i,j)*N2)           SN=sqrt(Ssq(i,j)*N2)
273           VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH           VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH
274       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN
# Line 182  C-- end 1rst loop on vertical level inde Line 284  C-- end 1rst loop on vertical level inde
284    
285    
286  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
287    #ifdef ALLOW_AUTODIFF_TAMC
288    CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
289    #endif
290        IF ( GM_Visbeck_alpha.NE.0. ) THEN        IF ( GM_Visbeck_alpha.NE.0. ) THEN
291  C-     Limit range that KapGM can take  C-     Limit range that KapGM can take
292         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
293          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
294           VisbeckK(i,j,bi,bj)=           VisbeckK(i,j,bi,bj)=
295       &       MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)       &       MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)
 #ifdef ALLOW_TIMEAVE  
          Visbeck_K_T(i,j,bi,bj)=Visbeck_K_T(i,j,bi,bj)  
      &                         +VisbeckK(i,j,bi,bj)*deltaTclock  
 #endif  
296          ENDDO          ENDDO
297         ENDDO         ENDDO
298        ENDIF        ENDIF
299    cph( NEW
300    #ifdef ALLOW_AUTODIFF_TAMC
301    CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
302    #endif
303    cph)
304  #endif /* GM_VISBECK_VARIABLE_K */  #endif /* GM_VISBECK_VARIABLE_K */
305    
306    
# Line 206  C-- 2nd loop on k : compute Tensor Coeff Line 312  C-- 2nd loop on k : compute Tensor Coeff
312         maskp1 = 1. _d 0         maskp1 = 1. _d 0
313         IF (k.GE.Nr) maskp1 = 0. _d 0         IF (k.GE.Nr) maskp1 = 0. _d 0
314    
315    #ifdef ALLOW_AUTODIFF_TAMC
316           kkey = (igmkey-1)*Nr + k
317    #if (defined (GM_NON_UNITY_DIAGONAL) || \
318         defined (GM_VISBECK_VARIABLE_K))
319    CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
320    CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
321    CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
322    #endif
323    #endif
324    
325  C-    express the Tensor in term of Diffusivity (= m**2 / s )  C-    express the Tensor in term of Diffusivity (= m**2 / s )
326        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
327         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
328    #ifdef ALLOW_KAPGM_CONTROL
329            Kgm_tmp = GM_isopycK + GM_skewflx*kapgm(i,j,k,bi,bj)
330    #else
331          Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K          Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K
332    #endif
333  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
334       &          + VisbeckK(i,j,bi,bj)*(1.+GM_skewflx)           &          + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)    
335  #endif  #endif
336          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)
337          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)
# Line 228  C-    express the Tensor in term of Diff Line 348  C-    express the Tensor in term of Diff
348  C     Gradient of Sigma at U points  C     Gradient of Sigma at U points
349        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
350         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
351          SlopeX(i,j)=sigmaX(i,j,k)          dSigmaDx(i,j)=sigmaX(i,j,k)
352       &          *_maskW(i,j,k,bi,bj)       &          *_maskW(i,j,k,bi,bj)
353          SlopeY(i,j)=0.25*( sigmaY(i-1,j+1,k) +sigmaY(i,j+1,k)          dSigmaDy(i,j)=op25*( sigmaY(i-1,j+1,k) +sigmaY(i,j+1,k)
354       &                    +sigmaY(i-1, j ,k) +sigmaY(i, j ,k) )       &                      +sigmaY(i-1, j ,k) +sigmaY(i, j ,k) )
355       &          *_maskW(i,j,k,bi,bj)       &          *_maskW(i,j,k,bi,bj)
356          dSigmaDrReal(i,j)=0.25*( sigmaR(i-1,j, k ) +sigmaR(i,j, k )          dSigmaDr(i,j)=op25*( sigmaR(i-1,j, k ) +sigmaR(i,j, k )
357       &                  +maskp1*(sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1)) )       &                  +maskp1*(sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1)) )
358       &          *_maskW(i,j,k,bi,bj)       &          *_maskW(i,j,k,bi,bj)
359         ENDDO         ENDDO
360        ENDDO        ENDDO
361    
362    #ifdef ALLOW_AUTODIFF_TAMC
363    CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
364    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
365    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
366    CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte
367    #endif /* ALLOW_AUTODIFF_TAMC */
368    
369  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
370        CALL GMREDI_SLOPE_LIMIT(        CALL GMREDI_SLOPE_LIMIT(
371       U             dSigmadRReal,       O             SlopeX, SlopeY,
      I             rF(K),  
      U             SlopeX, SlopeY,  
372       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
373         U             dSigmaDr,
374         I             dSigmaDx, dSigmaDy,
375         I             ldd97_LrhoW,rC(k),k,
376       I             bi, bj, myThid )       I             bi, bj, myThid )
377    
378    cph( NEW
379    #ifdef ALLOW_AUTODIFF_TAMC
380    cph(
381    CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
382    CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
383    cph)
384    #endif /* ALLOW_AUTODIFF_TAMC */
385    cph)
386    
387  #ifdef GM_NON_UNITY_DIAGONAL  #ifdef GM_NON_UNITY_DIAGONAL
388          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
389           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
390            Kux(i,j,k,bi,bj) =            Kux(i,j,k,bi,bj) =
391       &     ( GM_isopycK       &     ( GM_isopycK
392  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
393       &     +0.5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
394  #endif  #endif
395       &     )       &     )
396       &     *taperFct(i,j)       &     *taperFct(i,j)
397           ENDDO           ENDDO
398          ENDDO          ENDDO
399    #ifdef ALLOW_AUTODIFF_TAMC
400    # ifdef GM_EXCLUDE_CLIPPING
401    CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
402    # endif
403    #endif
404          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
405           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
406            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 )
# Line 267  C     Calculate slopes for use in tensor Line 409  C     Calculate slopes for use in tensor
409  #endif /* GM_NON_UNITY_DIAGONAL */  #endif /* GM_NON_UNITY_DIAGONAL */
410    
411  #ifdef GM_EXTRA_DIAGONAL  #ifdef GM_EXTRA_DIAGONAL
412    
413    #ifdef ALLOW_AUTODIFF_TAMC
414    CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
415    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
416    #endif /* ALLOW_AUTODIFF_TAMC */
417        IF (GM_ExtraDiag) THEN        IF (GM_ExtraDiag) THEN
418          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
419           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
420            Kuz(i,j,k,bi,bj) =            Kuz(i,j,k,bi,bj) =
421    #ifdef ALLOW_KAPGM_CONTROL
422         &     ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)
423    #else
424       &     ( GM_isopycK - GM_skewflx*GM_background_K       &     ( GM_isopycK - GM_skewflx*GM_background_K
425    #endif
426  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
427       &     +0.5*(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
428  #endif  #endif
429       &     )*SlopeX(i,j)*taperFct(i,j)       &     )*SlopeX(i,j)*taperFct(i,j)
430           ENDDO           ENDDO
# Line 281  C     Calculate slopes for use in tensor Line 432  C     Calculate slopes for use in tensor
432        ENDIF        ENDIF
433  #endif /* GM_EXTRA_DIAGONAL */  #endif /* GM_EXTRA_DIAGONAL */
434    
435    #ifdef ALLOW_DIAGNOSTICS
436          IF (doDiagRediFlx) THEN
437            km1 = MAX(k-1,1)
438            DO j=1,sNy
439             DO i=1,sNx+1
440    C         store in tmp1k Kuz_Redi
441              tmp1k(i,j) = ( GM_isopycK
442    #ifdef GM_VISBECK_VARIABLE_K
443         &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0
444    #endif
445         &                 )*SlopeX(i,j)*taperFct(i,j)
446             ENDDO
447            ENDDO
448            DO j=1,sNy
449             DO i=1,sNx+1
450    C-        Vertical gradients interpolated to U points
451              dTdz = (
452         &     +recip_drC(k)*
453         &       ( maskC(i-1,j,k,bi,bj)*
454         &           (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
455         &        +maskC( i ,j,k,bi,bj)*
456         &           (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
457         &       )
458         &     +recip_drC(kp1)*
459         &       ( maskC(i-1,j,kp1,bi,bj)*
460         &           (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
461         &        +maskC( i ,j,kp1,bi,bj)*
462         &           (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
463         &       )      ) * 0.25 _d 0
464               tmp1k(i,j) = dyG(i,j,bi,bj)*drF(k)
465         &                * _hFacW(i,j,k,bi,bj)
466         &                * tmp1k(i,j) * dTdz
467             ENDDO
468            ENDDO
469            CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KuzTz', k,1,2,bi,bj,myThid)
470          ENDIF
471    #endif /* ALLOW_DIAGNOSTICS */
472    
473  C     Gradient of Sigma at V points  C     Gradient of Sigma at V points
474        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
475         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
476          SlopeX(i,j)=0.25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)          dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
477       &                    +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k) )       &                    +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k) )
478       &          *_maskS(i,j,k,bi,bj)       &          *_maskS(i,j,k,bi,bj)
479          SlopeY(i,j)=sigmaY(i,j,k)          dSigmaDy(i,j)=sigmaY(i,j,k)
480       &          *_maskS(i,j,k,bi,bj)       &          *_maskS(i,j,k,bi,bj)
481          dSigmaDrReal(i,j)=0.25*( sigmaR(i,j-1, k ) +sigmaR(i,j, k )          dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k ) +sigmaR(i,j, k )
482       &                  +maskp1*(sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1)) )       &                  +maskp1*(sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1)) )
483       &          *_maskS(i,j,k,bi,bj)       &          *_maskS(i,j,k,bi,bj)
484         ENDDO         ENDDO
485        ENDDO        ENDDO
486    
487    #ifdef ALLOW_AUTODIFF_TAMC
488    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
489    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
490    CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte
491    #endif /* ALLOW_AUTODIFF_TAMC */
492    
493  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
494        CALL GMREDI_SLOPE_LIMIT(        CALL GMREDI_SLOPE_LIMIT(
495       U             dSigmadRReal,       O             SlopeX, SlopeY,
      I             rF(K),  
      U             SlopeX, SlopeY,  
496       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
497         U             dSigmaDr,
498         I             dSigmaDx, dSigmaDy,
499         I             ldd97_LrhoS,rC(k),k,
500       I             bi, bj, myThid )       I             bi, bj, myThid )
501    
502    cph(
503    #ifdef ALLOW_AUTODIFF_TAMC
504    cph(
505    CADJ STORE taperfct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
506    cph)
507    #endif /* ALLOW_AUTODIFF_TAMC */
508    cph)
509    
510  #ifdef GM_NON_UNITY_DIAGONAL  #ifdef GM_NON_UNITY_DIAGONAL
511          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
512           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
513            Kvy(i,j,k,bi,bj) =            Kvy(i,j,k,bi,bj) =
514       &     ( GM_isopycK       &     ( GM_isopycK
515  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
516       &     +0.5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
517  #endif  #endif
518       &     )       &     )
519       &     *taperFct(i,j)       &     *taperFct(i,j)
520           ENDDO           ENDDO
521          ENDDO          ENDDO
522    #ifdef ALLOW_AUTODIFF_TAMC
523    # ifdef GM_EXCLUDE_CLIPPING
524    CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
525    # endif
526    #endif
527          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
528           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
529            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 )
# Line 323  C     Calculate slopes for use in tensor Line 532  C     Calculate slopes for use in tensor
532  #endif /* GM_NON_UNITY_DIAGONAL */  #endif /* GM_NON_UNITY_DIAGONAL */
533    
534  #ifdef GM_EXTRA_DIAGONAL  #ifdef GM_EXTRA_DIAGONAL
535    
536    #ifdef ALLOW_AUTODIFF_TAMC
537    CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
538    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
539    #endif /* ALLOW_AUTODIFF_TAMC */
540        IF (GM_ExtraDiag) THEN        IF (GM_ExtraDiag) THEN
541          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
542           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
543            Kvz(i,j,k,bi,bj) =            Kvz(i,j,k,bi,bj) =
544    #ifdef ALLOW_KAPGM_CONTROL
545         &     ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)
546    #else
547       &     ( GM_isopycK - GM_skewflx*GM_background_K       &     ( GM_isopycK - GM_skewflx*GM_background_K
548    #endif
549  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
550       &     +0.5*(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
551  #endif  #endif
552       &     )*SlopeY(i,j)*taperFct(i,j)       &     )*SlopeY(i,j)*taperFct(i,j)
553           ENDDO           ENDDO
# Line 337  C     Calculate slopes for use in tensor Line 555  C     Calculate slopes for use in tensor
555        ENDIF        ENDIF
556  #endif /* GM_EXTRA_DIAGONAL */  #endif /* GM_EXTRA_DIAGONAL */
557    
558  #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */  #ifdef ALLOW_DIAGNOSTICS
559          IF (doDiagRediFlx) THEN
560    c       km1 = MAX(k-1,1)
561            DO j=1,sNy+1
562             DO i=1,sNx
563    C         store in tmp1k Kvz_Redi
564              tmp1k(i,j) = ( GM_isopycK
565    #ifdef GM_VISBECK_VARIABLE_K
566         &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0
567    #endif
568         &                 )*SlopeY(i,j)*taperFct(i,j)
569             ENDDO
570            ENDDO
571            DO j=1,sNy+1
572             DO i=1,sNx
573    C-        Vertical gradients interpolated to U points
574              dTdz = (
575         &     +recip_drC(k)*
576         &       ( maskC(i,j-1,k,bi,bj)*
577         &           (theta(i,j-1,km1,bi,bj)-theta(i,j-1,k,bi,bj))
578         &        +maskC(i, j ,k,bi,bj)*
579         &           (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))
580         &       )
581         &     +recip_drC(kp1)*
582         &       ( maskC(i,j-1,kp1,bi,bj)*
583         &           (theta(i,j-1,k,bi,bj)-theta(i,j-1,kp1,bi,bj))
584         &        +maskC(i, j ,kp1,bi,bj)*
585         &           (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))
586         &       )      ) * 0.25 _d 0
587               tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)
588         &                * _hFacS(i,j,k,bi,bj)
589         &                * tmp1k(i,j) * dTdz
590             ENDDO
591            ENDDO
592            CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KvzTz', k,1,2,bi,bj,myThid)
593          ENDIF
594    #endif /* ALLOW_DIAGNOSTICS */
595    
596  #ifdef ALLOW_TIMEAVE  #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
 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  
       ENDDO  
       GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock  
 #endif /* ALLOW_TIMEAVE */  
597    
598  C-- end 2nd loop on vertical level index k  C-- end 2nd loop on vertical level index k
599        ENDDO        ENDDO
# Line 363  C-- end 2nd loop on vertical level index Line 604  C-- end 2nd loop on vertical level index
604          CALL GMREDI_CALC_PSI_B(          CALL GMREDI_CALC_PSI_B(
605       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
606       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
607         I             ldd97_LrhoW, ldd97_LrhoS,
608       I             myThid )       I             myThid )
609        ENDIF        ENDIF
610  #endif  #endif
611    
612    #ifdef ALLOW_TIMEAVE
613    C--   Time-average
614          IF ( taveFreq.GT.0. ) THEN
615    
616             CALL TIMEAVE_CUMULATE( GM_Kwx_T, Kwx, Nr,
617         &                          deltaTclock, bi, bj, myThid )
618             CALL TIMEAVE_CUMULATE( GM_Kwy_T, Kwy, Nr,
619         &                          deltaTclock, bi, bj, myThid )
620             CALL TIMEAVE_CUMULATE( GM_Kwz_T, Kwz, Nr,
621         &                          deltaTclock, bi, bj, myThid )
622    #ifdef GM_VISBECK_VARIABLE_K
623           IF ( GM_Visbeck_alpha.NE.0. ) THEN
624             CALL TIMEAVE_CUMULATE( Visbeck_K_T, VisbeckK, 1,
625         &                          deltaTclock, bi, bj, myThid )
626           ENDIF
627    #endif
628    #ifdef GM_BOLUS_ADVEC
629           IF ( GM_AdvForm ) THEN
630             CALL TIMEAVE_CUMULATE( GM_PsiXtave, GM_PsiX, Nr,
631         &                          deltaTclock, bi, bj, myThid )
632             CALL TIMEAVE_CUMULATE( GM_PsiYtave, GM_PsiY, Nr,
633         &                          deltaTclock, bi, bj, myThid )
634           ENDIF
635    #endif
636           DO k=1,Nr
637             GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock
638           ENDDO
639    
640          ENDIF
641    #endif /* ALLOW_TIMEAVE */
642    
643    #ifdef ALLOW_DIAGNOSTICS
644          IF ( useDiagnostics ) THEN
645            CALL GMREDI_DIAGNOSTICS_FILL(bi,bj,myThid)
646          ENDIF
647    #endif /* ALLOW_DIAGNOSTICS */
648    
649  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
650    
651        RETURN        RETURN
# Line 386  C     \================================= Line 665  C     \=================================
665    
666  C     == Global variables ==  C     == Global variables ==
667  #include "SIZE.h"  #include "SIZE.h"
 #include "GRID.h"  
 #include "DYNVARS.h"  
668  #include "EEPARAMS.h"  #include "EEPARAMS.h"
 #include "PARAMS.h"  
669  #include "GMREDI.h"  #include "GMREDI.h"
670    
671  C     == Routine arguments ==  C     == Routine arguments ==
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.25
  ViewVC Help
Powered by ViewVC 1.1.22