/[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.8 by jmc, Tue Dec 4 15:09:34 2001 UTC revision 1.24 by jmc, Tue Jun 20 22:55:08 2006 UTC
# Line 5  C $Name$ Line 5  C $Name$
5    
6  CStartOfInterface  CStartOfInterface
7        SUBROUTINE GMREDI_CALC_TENSOR(        SUBROUTINE GMREDI_CALC_TENSOR(
8       I             bi, bj, iMin, iMax, jMin, jMax, K,       I             bi, bj, iMin, iMax, jMin, jMax,
9       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
10       I             myThid )       I             myThid )
11  C     /==========================================================\  C     /==========================================================\
# 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
28    #include "tamc.h"
29    #include "tamc_keys.h"
30    #endif /* ALLOW_AUTODIFF_TAMC */
31    
32  C     == Routine arguments ==  C     == Routine arguments ==
33  C  C
34        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
35        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
36        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
37        INTEGER bi,bj,iMin,iMax,jMin,jMax,K        INTEGER bi,bj,iMin,iMax,jMin,jMax
38        INTEGER myThid        INTEGER myThid
39  CEndOfInterface  CEndOfInterface
40    
41  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
42    
43  C     == Local variables ==  C     == Local variables ==
44        INTEGER i,j,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, 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        _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        km1=max(1,K-1)  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
       kp1=min(Nr,K)  
   
75    
76  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
77  !HPF$ INDEPENDENT            act1 = bi - myBxLo(myThid)
78              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
79              act2 = bj - myByLo(myThid)
80              max2 = myByHi(myThid) - myByLo(myThid) + 1
81              act3 = myThid - 1
82              max3 = nTx*nTy
83              act4 = ikey_dynamics - 1
84              igmkey = (act1 + 1) + act2*max1
85         &                      + act3*max1*max2
86         &                      + act4*max1*max2*max3
87    #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  #endif
97        DO j=1-Oly+1,sNy+Oly-1      
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
163    C-- 1rst loop on k : compute Tensor Coeff. at W points.
164    
165  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
166  !HPF$ INDEPENDENT         kkey = (igmkey-1)*Nr + k
167           DO j=1-Oly,sNy+Oly
168            DO i=1-Olx,sNx+Olx
169             SlopeX(i,j)       = 0. _d 0
170             SlopeY(i,j)       = 0. _d 0
171             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
175             taperFct(i,j)     = 0. _d 0
176             Kwx(i,j,k,bi,bj)  = 0. _d 0
177             Kwy(i,j,k,bi,bj)  = 0. _d 0
178             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
192           ENDDO
193  #endif  #endif
        DO i=1-Olx+1,sNx+Olx-1  
194    
195          DO j=1-Oly+1,sNy+Oly-1
196           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)       &                  *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)       &                  *maskC(i,j,k,bi,bj)
204          dSigmaDrReal(i,j)=sigmaR(i,j,k)          dSigmaDr(i,j)=sigmaR(i,j,k)
   
205         ENDDO         ENDDO
206        ENDDO        ENDDO
207    
208    #ifdef ALLOW_AUTODIFF_TAMC
209    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
210    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
211    CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
212    #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       I             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
# Line 91  C       Mask Iso-neutral slopes Line 227  C       Mask Iso-neutral slopes
227          SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)          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)          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)          SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)
 c       Ssq=SlopeX(i,j)*SlopeX(i,j)+SlopeY(i,j)*SlopeY(i,j)  
230    
231  C       Components of Redi/GM tensor         ENDDO
232          Kwx(i,j,k,bi,bj)=2.*SlopeX(i,j)*taperFct(i,j)        ENDDO
233          Kwy(i,j,k,bi,bj)=2.*SlopeY(i,j)*taperFct(i,j)  
234    #ifdef ALLOW_AUTODIFF_TAMC
235    CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
236    CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
237    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 */
241    
242          DO j=1-Oly+1,sNy+Oly-1
243           DO i=1-Olx+1,sNx+Olx-1
244    
245    C       Components of Redi/GM tensor
246            Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
247            Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
248          Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)          Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
249    
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=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
257    
# Line 118  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.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*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
275          ENDIF          ENDIF
276    
 C       Limit range that KapGM can take  
         VisbeckK(i,j,bi,bj)=  
      &     min(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)  
   
277  #endif /* GM_VISBECK_VARIABLE_K */  #endif /* GM_VISBECK_VARIABLE_K */
278    
279           ENDDO
280          ENDDO
281    
282    C-- end 1rst loop on vertical level index k
283          ENDDO
284    
285    
286  #ifdef ALLOW_TIMEAVE  #ifdef GM_VISBECK_VARIABLE_K
287  C--     Time-average  #ifdef ALLOW_AUTODIFF_TAMC
288          GM_Kwx_T(i,j,k,bi,bj)=GM_Kwx_T(i,j,k,bi,bj)  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
      &                       +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  
 #ifdef GM_VISBECK_VARIABLE_K  
         IF (K.EQ.Nr)  
      &  Visbeck_K_T(i,j,bi,bj)=Visbeck_K_T(i,j,bi,bj)  
      &                       +VisbeckK(i,j,bi,bj)*deltaTclock  
289  #endif  #endif
290  #endif /* ALLOW_TIMEAVE */        IF ( GM_Visbeck_alpha.NE.0. ) THEN
291    C-     Limit range that KapGM can take
292           DO j=1-Oly+1,sNy+Oly-1
293            DO i=1-Olx+1,sNx+Olx-1
294             VisbeckK(i,j,bi,bj)=
295         &       MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)
296            ENDDO
297         ENDDO         ENDDO
298        ENDDO        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 */
305    
306  #ifdef ALLOW_TIMEAVE  
307        GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
308    
309    C-- 2nd loop on k : compute Tensor Coeff. at U,V levels.
310          DO k=1,Nr
311           kp1 = MIN(Nr,k+1)
312           maskp1 = 1. _d 0
313           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  #endif
324    
325    C-    express the Tensor in term of Diffusivity (= m**2 / s )
326          DO j=1-Oly+1,sNy+Oly-1
327           DO i=1-Olx+1,sNx+Olx-1
328            Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K
329    #ifdef GM_VISBECK_VARIABLE_K
330         &          + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)    
331    #endif
332            Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)
333            Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)
334            Kwz(i,j,k,bi,bj)= ( GM_isopycK
335    #ifdef GM_VISBECK_VARIABLE_K
336         &                    + VisbeckK(i,j,bi,bj)
337    #endif
338         &                    )*Kwz(i,j,k,bi,bj)
339           ENDDO
340          ENDDO
341    
342    #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
343    
 #ifdef GM_NON_UNITY_DIAGONAL  
344  C     Gradient of Sigma at U points  C     Gradient of Sigma at U points
345        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
346         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
347          SlopeX(i,j)=sigmaX(i,j,km1)          dSigmaDx(i,j)=sigmaX(i,j,k)
348       &          *_maskW(i,j,k,bi,bj)       &          *_maskW(i,j,k,bi,bj)
349          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)
350       &                    +sigmaY(i-1, j ,k) +sigmaY(i, j ,k) )       &                      +sigmaY(i-1, j ,k) +sigmaY(i, j ,k) )
351       &          *_maskW(i,j,k,bi,bj)       &          *_maskW(i,j,k,bi,bj)
352          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 )
353       &                          +sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1) )       &                  +maskp1*(sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1)) )
354       &          *_maskW(i,j,k,bi,bj)       &          *_maskW(i,j,k,bi,bj)
355         ENDDO         ENDDO
356        ENDDO        ENDDO
357    
358    #ifdef ALLOW_AUTODIFF_TAMC
359    CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
360    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
361    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
362    CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte
363    #endif /* ALLOW_AUTODIFF_TAMC */
364    
365  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
366        CALL GMREDI_SLOPE_LIMIT(        CALL GMREDI_SLOPE_LIMIT(
367       I             dSigmadRReal,       O             SlopeX, SlopeY,
      I             rF(K),  
      U             SlopeX, SlopeY,  
368       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
369         U             dSigmaDr,
370         I             dSigmaDx, dSigmaDy,
371         I             ldd97_LrhoW,rC(k),k,
372       I             bi, bj, myThid )       I             bi, bj, myThid )
373    
374        DO j=1-Oly+1,sNy+Oly-1  cph( NEW
375         DO i=1-Olx+1,sNx+Olx-1  #ifdef ALLOW_AUTODIFF_TAMC
376          Kux(i,j,k,bi,bj)=taperFct(i,j)  cph(
377         ENDDO  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
378        ENDDO  CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
379    cph)
380    #endif /* ALLOW_AUTODIFF_TAMC */
381    cph)
382    
383    #ifdef GM_NON_UNITY_DIAGONAL
384            DO j=1-Oly+1,sNy+Oly-1
385             DO i=1-Olx+1,sNx+Olx-1
386              Kux(i,j,k,bi,bj) =
387         &     ( GM_isopycK
388    #ifdef GM_VISBECK_VARIABLE_K
389         &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
390    #endif
391         &     )
392         &     *taperFct(i,j)
393             ENDDO
394            ENDDO
395    #ifdef ALLOW_AUTODIFF_TAMC
396    # ifdef GM_EXCLUDE_CLIPPING
397    CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
398    # endif
399    #endif
400            DO j=1-Oly+1,sNy+Oly-1
401             DO i=1-Olx+1,sNx+Olx-1
402              Kux(i,j,k,bi,bj) = MAX( Kux(i,j,k,bi,bj), GM_Kmin_horiz )
403             ENDDO
404            ENDDO
405    #endif /* GM_NON_UNITY_DIAGONAL */
406    
407    #ifdef GM_EXTRA_DIAGONAL
408    
409    #ifdef ALLOW_AUTODIFF_TAMC
410    CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
411    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
412    #endif /* ALLOW_AUTODIFF_TAMC */
413          IF (GM_ExtraDiag) THEN
414            DO j=1-Oly+1,sNy+Oly-1
415             DO i=1-Olx+1,sNx+Olx-1
416              Kuz(i,j,k,bi,bj) =
417         &     ( GM_isopycK - GM_skewflx*GM_background_K
418    #ifdef GM_VISBECK_VARIABLE_K
419         &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*GM_advect
420    #endif
421         &     )*SlopeX(i,j)*taperFct(i,j)
422             ENDDO
423            ENDDO
424          ENDIF
425    #endif /* GM_EXTRA_DIAGONAL */
426    
427    #ifdef ALLOW_DIAGNOSTICS
428          IF (doDiagRediFlx) THEN
429            km1 = MAX(k-1,1)
430            DO j=1,sNy
431             DO i=1,sNx+1
432    C         store in tmp1k Kuz_Redi
433              tmp1k(i,j) = ( GM_isopycK
434    #ifdef GM_VISBECK_VARIABLE_K
435         &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0
436    #endif
437         &                 )*SlopeX(i,j)*taperFct(i,j)
438             ENDDO
439            ENDDO
440            DO j=1,sNy
441             DO i=1,sNx+1
442    C-        Vertical gradients interpolated to U points
443              dTdz = (
444         &     +recip_drC(k)*
445         &       ( maskC(i-1,j,k,bi,bj)*
446         &           (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
447         &        +maskC( i ,j,k,bi,bj)*
448         &           (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
449         &       )
450         &     +recip_drC(kp1)*
451         &       ( maskC(i-1,j,kp1,bi,bj)*
452         &           (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
453         &        +maskC( i ,j,kp1,bi,bj)*
454         &           (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
455         &       )      ) * 0.25 _d 0
456               tmp1k(i,j) = dyG(i,j,bi,bj)*drF(k)
457         &                * _hFacW(i,j,k,bi,bj)
458         &                * tmp1k(i,j) * dTdz
459             ENDDO
460            ENDDO
461            CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KuzTz', k,1,2,bi,bj,myThid)
462          ENDIF
463    #endif /* ALLOW_DIAGNOSTICS */
464    
465  C     Gradient of Sigma at V points  C     Gradient of Sigma at V points
466        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
467         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
468          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)
469       &                    +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k) )       &                    +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k) )
470       &          *_maskS(i,j,k,bi,bj)       &          *_maskS(i,j,k,bi,bj)
471          SlopeY(i,j)=sigmaY(i,j,km1)          dSigmaDy(i,j)=sigmaY(i,j,k)
472       &          *_maskS(i,j,k,bi,bj)       &          *_maskS(i,j,k,bi,bj)
473          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 )
474       &                          +sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1) )       &                  +maskp1*(sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1)) )
475       &          *_maskS(i,j,k,bi,bj)       &          *_maskS(i,j,k,bi,bj)
476         ENDDO         ENDDO
477        ENDDO        ENDDO
478    
479    #ifdef ALLOW_AUTODIFF_TAMC
480    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
481    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
482    CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte
483    #endif /* ALLOW_AUTODIFF_TAMC */
484    
485  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
486        CALL GMREDI_SLOPE_LIMIT(        CALL GMREDI_SLOPE_LIMIT(
487       I             dSigmadRReal,       O             SlopeX, SlopeY,
      I             rF(K),  
      U             SlopeX, SlopeY,  
488       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
489         U             dSigmaDr,
490         I             dSigmaDx, dSigmaDy,
491         I             ldd97_LrhoS,rC(k),k,
492       I             bi, bj, myThid )       I             bi, bj, myThid )
493    
494        DO j=1-Oly+1,sNy+Oly-1  cph(
495         DO i=1-Olx+1,sNx+Olx-1  #ifdef ALLOW_AUTODIFF_TAMC
496          Kvy(i,j,k,bi,bj)=taperFct(i,j)  cph(
497         ENDDO  CADJ STORE taperfct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
498        ENDDO  cph)
499    #endif /* ALLOW_AUTODIFF_TAMC */
500    cph)
501    
502    #ifdef GM_NON_UNITY_DIAGONAL
503            DO j=1-Oly+1,sNy+Oly-1
504             DO i=1-Olx+1,sNx+Olx-1
505              Kvy(i,j,k,bi,bj) =
506         &     ( GM_isopycK
507    #ifdef GM_VISBECK_VARIABLE_K
508         &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
509    #endif
510         &     )
511         &     *taperFct(i,j)
512             ENDDO
513            ENDDO
514    #ifdef ALLOW_AUTODIFF_TAMC
515    # ifdef GM_EXCLUDE_CLIPPING
516    CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
517    # endif
518    #endif
519            DO j=1-Oly+1,sNy+Oly-1
520             DO i=1-Olx+1,sNx+Olx-1
521              Kvy(i,j,k,bi,bj) = MAX( Kvy(i,j,k,bi,bj), GM_Kmin_horiz )
522             ENDDO
523            ENDDO
524  #endif /* GM_NON_UNITY_DIAGONAL */  #endif /* GM_NON_UNITY_DIAGONAL */
525    
526    #ifdef GM_EXTRA_DIAGONAL
527    
528    #ifdef ALLOW_AUTODIFF_TAMC
529    CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
530    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
531    #endif /* ALLOW_AUTODIFF_TAMC */
532          IF (GM_ExtraDiag) THEN
533            DO j=1-Oly+1,sNy+Oly-1
534             DO i=1-Olx+1,sNx+Olx-1
535              Kvz(i,j,k,bi,bj) =
536         &     ( GM_isopycK - GM_skewflx*GM_background_K
537    #ifdef GM_VISBECK_VARIABLE_K
538         &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*GM_advect
539    #endif
540         &     )*SlopeY(i,j)*taperFct(i,j)
541             ENDDO
542            ENDDO
543          ENDIF
544    #endif /* GM_EXTRA_DIAGONAL */
545    
546    #ifdef ALLOW_DIAGNOSTICS
547          IF (doDiagRediFlx) THEN
548    c       km1 = MAX(k-1,1)
549            DO j=1,sNy+1
550             DO i=1,sNx
551    C         store in tmp1k Kvz_Redi
552              tmp1k(i,j) = ( GM_isopycK
553    #ifdef GM_VISBECK_VARIABLE_K
554         &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0
555    #endif
556         &                 )*SlopeY(i,j)*taperFct(i,j)
557             ENDDO
558            ENDDO
559            DO j=1,sNy+1
560             DO i=1,sNx
561    C-        Vertical gradients interpolated to U points
562              dTdz = (
563         &     +recip_drC(k)*
564         &       ( maskC(i,j-1,k,bi,bj)*
565         &           (theta(i,j-1,km1,bi,bj)-theta(i,j-1,k,bi,bj))
566         &        +maskC(i, j ,k,bi,bj)*
567         &           (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))
568         &       )
569         &     +recip_drC(kp1)*
570         &       ( maskC(i,j-1,kp1,bi,bj)*
571         &           (theta(i,j-1,k,bi,bj)-theta(i,j-1,kp1,bi,bj))
572         &        +maskC(i, j ,kp1,bi,bj)*
573         &           (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))
574         &       )      ) * 0.25 _d 0
575               tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)
576         &                * _hFacS(i,j,k,bi,bj)
577         &                * tmp1k(i,j) * dTdz
578             ENDDO
579            ENDDO
580            CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KvzTz', k,1,2,bi,bj,myThid)
581          ENDIF
582    #endif /* ALLOW_DIAGNOSTICS */
583    
584    #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
585    
586    C-- end 2nd loop on vertical level index k
587          ENDDO
588    
589    
590    #ifdef GM_BOLUS_ADVEC
591          IF (GM_AdvForm) THEN
592            CALL GMREDI_CALC_PSI_B(
593         I             bi, bj, iMin, iMax, jMin, jMax,
594         I             sigmaX, sigmaY, sigmaR,
595         I             ldd97_LrhoW, ldd97_LrhoS,
596         I             myThid )
597          ENDIF
598    #endif
599    
600    #ifdef ALLOW_TIMEAVE
601    C--   Time-average
602          IF ( taveFreq.GT.0. ) THEN
603    
604             CALL TIMEAVE_CUMULATE( GM_Kwx_T, Kwx, Nr,
605         &                          deltaTclock, bi, bj, myThid )
606             CALL TIMEAVE_CUMULATE( GM_Kwy_T, Kwy, Nr,
607         &                          deltaTclock, bi, bj, myThid )
608             CALL TIMEAVE_CUMULATE( GM_Kwz_T, Kwz, Nr,
609         &                          deltaTclock, bi, bj, myThid )
610    #ifdef GM_VISBECK_VARIABLE_K
611           IF ( GM_Visbeck_alpha.NE.0. ) THEN
612             CALL TIMEAVE_CUMULATE( Visbeck_K_T, VisbeckK, 1,
613         &                          deltaTclock, bi, bj, myThid )
614           ENDIF
615    #endif
616    #ifdef GM_BOLUS_ADVEC
617           IF ( GM_AdvForm ) THEN
618             CALL TIMEAVE_CUMULATE( GM_PsiXtave, GM_PsiX, Nr,
619         &                          deltaTclock, bi, bj, myThid )
620             CALL TIMEAVE_CUMULATE( GM_PsiYtave, GM_PsiY, Nr,
621         &                          deltaTclock, bi, bj, myThid )
622           ENDIF
623    #endif
624           DO k=1,Nr
625             GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock
626           ENDDO
627    
628          ENDIF
629    #endif /* ALLOW_TIMEAVE */
630    
631    #ifdef ALLOW_DIAGNOSTICS
632          IF ( useDiagnostics ) THEN
633            CALL GMREDI_DIAGNOSTICS_FILL(bi,bj,myThid)
634          ENDIF
635    #endif /* ALLOW_DIAGNOSTICS */
636    
637  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
638    
# Line 222  C     Calculate slopes for use in tensor Line 641  C     Calculate slopes for use in tensor
641    
642    
643        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
644       I             bi, bj, iMin, iMax, jMin, jMax, K,       I             bi, bj, iMin, iMax, jMin, jMax,
645       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
646       I             myThid )       I             myThid )
647  C     /==========================================================\  C     /==========================================================\
# Line 234  C     \================================= Line 653  C     \=================================
653    
654  C     == Global variables ==  C     == Global variables ==
655  #include "SIZE.h"  #include "SIZE.h"
 #include "GRID.h"  
 #include "DYNVARS.h"  
656  #include "EEPARAMS.h"  #include "EEPARAMS.h"
 #include "PARAMS.h"  
657  #include "GMREDI.h"  #include "GMREDI.h"
658    
659  C     == Routine arguments ==  C     == Routine arguments ==
# Line 245  C Line 661  C
661        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
662        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
663        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
664        INTEGER bi,bj,iMin,iMax,jMin,jMax,K        INTEGER bi,bj,iMin,iMax,jMin,jMax
665        INTEGER myThid        INTEGER myThid
666  CEndOfInterface  CEndOfInterface
667    
668        INTEGER i, j        INTEGER i, j, k
669    
670  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
671    
672        DO j=1-Oly+1,sNy+Oly-1        DO k=1,Nr
673         DO i=1-Olx+1,sNx+Olx-1         DO j=1-Oly+1,sNy+Oly-1
674          Kwx(i,j,k,bi,bj) = 0.0          DO i=1-Olx+1,sNx+Olx-1
675          Kwy(i,j,k,bi,bj) = 0.0           Kwx(i,j,k,bi,bj) = 0.0
676          Kwz(i,j,k,bi,bj) = 0.0           Kwy(i,j,k,bi,bj) = 0.0
677             Kwz(i,j,k,bi,bj) = 0.0
678            ENDDO
679         ENDDO         ENDDO
680        ENDDO        ENDDO
681  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
682    
683        end        RETURN
684          END
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22