/[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.4 by adcroft, Fri Feb 2 21:36:29 2001 UTC revision 1.19 by jmc, Sun Nov 21 15:55:38 2004 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "GMREDI_OPTIONS.h"  #include "GMREDI_OPTIONS.h"
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 23  C     == Global variables == Line 24  C     == Global variables ==
24  #include "GMREDI.h"  #include "GMREDI.h"
25  #include "GMREDI_DIAGS.h"  #include "GMREDI_DIAGS.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 dRdSigmaLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49        _RL Ssq        _RL dSigmaDr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50          _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51          _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(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
63  #endif  #endif
64    
65    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
66    
67        km1=max(1,K-1)  #ifdef ALLOW_AUTODIFF_TAMC
68        kp1=min(Nr,K)            act1 = bi - myBxLo(myThid)
69              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
70              act2 = bj - myByLo(myThid)
71              max2 = myByHi(myThid) - myByLo(myThid) + 1
72              act3 = myThid - 1
73              max3 = nTx*nTy
74              act4 = ikey_dynamics - 1
75              igmkey = (act1 + 1) + act2*max1
76         &                      + act3*max1*max2
77         &                      + act4*max1*max2*max3
78    #endif /* ALLOW_AUTODIFF_TAMC */
79    
80    #ifdef GM_VISBECK_VARIABLE_K
81          DO j=1-Oly,sNy+Oly
82           DO i=1-Olx,sNx+Olx
83            VisbeckK(i,j,bi,bj) = 0. _d 0
84           ENDDO
85          ENDDO
86    #endif
87    
88    C--   set ldd97_Lrho (for tapering scheme ldd97):
89          IF (GM_taper_scheme.EQ.'ldd97') THEN
90           Cspd = 2. _d 0
91           LrhoInf = 15. _d 3
92           LrhoSup = 100. _d 3
93    C-     Tracer point location (center):
94           DO j=1-Oly,sNy+Oly
95            DO i=1-Olx,sNx+Olx
96             IF (fCori(i,j,bi,bj).NE.0.) THEN
97               ldd97_LrhoC(i,j) = Cspd/ABS(fCori(i,j,bi,bj))
98             ELSE
99               ldd97_LrhoC(i,j) = LrhoSup
100             ENDIF
101             ldd97_LrhoC(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoC(i,j),LrhoSup))
102            ENDDO
103           ENDDO
104    C-     U point location (West):
105           DO j=1-Oly,sNy+Oly
106            ldd97_LrhoW(1-Olx,j) = LrhoSup
107            DO i=1-Olx+1,sNx+Olx
108             fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))
109             IF (fCoriLoc.NE.0.) THEN
110               ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)
111             ELSE
112               ldd97_LrhoW(i,j) = LrhoSup
113             ENDIF
114             ldd97_LrhoW(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoW(i,j),LrhoSup))
115            ENDDO
116           ENDDO
117    C-     V point location (South):
118           DO i=1-Olx+1,sNx+Olx
119             ldd97_LrhoS(i,1-Oly) = LrhoSup
120           ENDDO
121           DO j=1-Oly+1,sNy+Oly
122            DO i=1-Olx,sNx+Olx
123             fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))
124             IF (fCoriLoc.NE.0.) THEN
125               ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)
126             ELSE
127               ldd97_LrhoS(i,j) = LrhoSup
128             ENDIF
129             ldd97_LrhoS(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoS(i,j),LrhoSup))
130            ENDDO
131           ENDDO
132          ELSE
133    C-    Just initialize to zero (not use anyway)
134           DO j=1-Oly,sNy+Oly
135            DO i=1-Olx,sNx+Olx
136              ldd97_LrhoC(i,j) = 0. _d 0
137              ldd97_LrhoW(i,j) = 0. _d 0
138              ldd97_LrhoS(i,j) = 0. _d 0
139            ENDDO
140           ENDDO
141          ENDIF
142    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
143    
144          DO k=2,Nr
145    C-- 1rst loop on k : compute Tensor Coeff. at W points.
146    
147  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
148  !HPF$ INDEPENDENT         kkey = (igmkey-1)*Nr + k
149           DO j=1-Oly,sNy+Oly
150            DO i=1-Olx,sNx+Olx
151             SlopeX(i,j)       = 0. _d 0
152             SlopeY(i,j)       = 0. _d 0
153             dSigmaDx(i,j)     = 0. _d 0
154             dSigmaDy(i,j)     = 0. _d 0
155             dSigmaDr(i,j)     = 0. _d 0
156             SlopeSqr(i,j)     = 0. _d 0
157             taperFct(i,j)     = 0. _d 0
158             Kwx(i,j,k,bi,bj)  = 0. _d 0
159             Kwy(i,j,k,bi,bj)  = 0. _d 0
160             Kwz(i,j,k,bi,bj)  = 0. _d 0
161    # ifdef GM_NON_UNITY_DIAGONAL
162             Kux(i,j,k,bi,bj)  = 0. _d 0
163             Kvy(i,j,k,bi,bj)  = 0. _d 0
164    # endif
165    # ifdef GM_EXTRA_DIAGONAL
166             Kuz(i,j,k,bi,bj)  = 0. _d 0
167             Kvz(i,j,k,bi,bj)  = 0. _d 0
168    # endif
169    # ifdef GM_BOLUS_ADVEC
170             GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
171             GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
172    # endif
173            ENDDO
174           ENDDO
175  #endif  #endif
176    
177        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
 #ifdef ALLOW_AUTODIFF_TAMC  
 !HPF$ INDEPENDENT  
 #endif  
178         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
   
179  C      Gradient of Sigma at rVel points  C      Gradient of Sigma at rVel points
180          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)
181       &                    +sigmaX(i+1, j , k ) +sigmaX(i,j, k ) )       &                    +sigmaX(i+1, j , k ) +sigmaX(i,j, k ) )
182          SlopeY(i,j)=0.25*( sigmaY( i ,j+1,km1) +sigmaY(i,j,km1)       &                  *maskC(i,j,k,bi,bj)
183            dSigmaDy(i,j)=op25*( sigmaY( i ,j+1,k-1) +sigmaY(i,j,k-1)
184       &                    +sigmaY( i ,j+1, k ) +sigmaY(i,j, k ) )       &                    +sigmaY( i ,j+1, k ) +sigmaY(i,j, k ) )
185          dSigmaDrReal(i,j)=sigmaR(i,j,k)       &                  *maskC(i,j,k,bi,bj)
186            dSigmaDr(i,j)=sigmaR(i,j,k)
         if (hFacC(i,j,k,bi,bj).eq.0.) then  
          SlopeX(i,j)=0.  
          SlopeY(i,j)=0.  
         endif  
   
187         ENDDO         ENDDO
188        ENDDO        ENDDO
189    
190    #ifdef ALLOW_AUTODIFF_TAMC
191    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
192    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
193    CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
194    #endif /* ALLOW_AUTODIFF_TAMC */
195    
196  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
197        CALL GMREDI_SLOPE_LIMIT(        CALL GMREDI_SLOPE_LIMIT(
198       I             dSigmadRReal,       O             SlopeX, SlopeY,
199       I             rF(K),       O             SlopeSqr, taperFct,
200       U             SlopeX, SlopeY,       U             dSigmaDr,
201       O             dRdSigmaLtd,       I             dSigmaDx, dSigmaDy,
202         I             ldd97_LrhoC,rF(k),k,
203       I             bi, bj, myThid )       I             bi, bj, myThid )
204    
205        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
206         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
207    
208  C       Mask Iso-neutral slopes  C       Mask Iso-neutral slopes
209          if (hFacC(i,j,k,bi,bj).eq.0.) then          SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)
210           SlopeX(i,j)=0.          SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)
211           SlopeY(i,j)=0.          SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)
212          endif  
213          Ssq=SlopeX(i,j)*SlopeX(i,j)+SlopeY(i,j)*SlopeY(i,j)         ENDDO
214          ENDDO
215  C       Components of Redi/GM tensor  
216          Kwx(i,j,k,bi,bj)=2.*SlopeX(i,j)  #ifdef ALLOW_AUTODIFF_TAMC
217          Kwy(i,j,k,bi,bj)=2.*SlopeY(i,j)  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
218          Kwz(i,j,k,bi,bj)=Ssq  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
219    CADJ STORE SlopeSqr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
220    CADJ STORE dSigmaDr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
221    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
222    #endif /* ALLOW_AUTODIFF_TAMC */
223    
224          DO j=1-Oly+1,sNy+Oly-1
225           DO i=1-Olx+1,sNx+Olx-1
226    
227    C       Components of Redi/GM tensor
228            Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
229            Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
230            Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
231    
232  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
233    
234    C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K
235    C           but do not know if *taperFct (or **2 ?) is necessary
236            Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)
237    
238  C--     Depth average of M^2/N^2 * N  C--     Depth average of M^2/N^2 * N
239    
240  C       Calculate terms for mean Richardson number  C       Calculate terms for mean Richardson number
# Line 114  C       If negative then we are below th Line 248  C       If negative then we are below th
248  C       Now we convert deltaH to a non-dimensional fraction  C       Now we convert deltaH to a non-dimensional fraction
249          deltaH=deltaH/GM_Visbeck_depth          deltaH=deltaH/GM_Visbeck_depth
250    
251          if (K.eq.2) VisbeckK(i,j,bi,bj)=0.          IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.
252  Calt?   if (dSigmaDrReal(i,j).NE.0.) then          IF ( Ssq(i,j).NE.0. .AND. dSigmaDr(i,j).NE.0. ) THEN
253  Calt?    N2=(-Gravity*recip_Rhonil)*dSigmaDrReal(i,j)           N2= -Gravity*recip_RhoConst*dSigmaDr(i,j)
254          if (dRdSigmaLtd(i,j).NE.0.) then           SN=sqrt(Ssq(i,j)*N2)
          N2=(-Gravity*recip_Rhonil)/dRdSigmaLtd(i,j)  
          SN=sqrt(Ssq*N2)  
255           VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH           VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH
256       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN
257          endif          ENDIF
   
 C       Limit range that KapGM can take  
         VisbeckK(i,j,bi,bj)=  
      &     min(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)  
258    
259  #endif /* GM_VISBECK_VARIABLE_K */  #endif /* GM_VISBECK_VARIABLE_K */
260    
261           ENDDO
262          ENDDO
263    
264  #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE  C-- end 1rst loop on vertical level index k
265  C--     Time-average        ENDDO
266          GM_Kwx_T(i,j,k,bi,bj)=GM_Kwx_T(i,j,k,bi,bj)  
267       &                       +Kwx(i,j,k,bi,bj)*deltaTclock  
268          GM_Kwy_T(i,j,k,bi,bj)=GM_Kwy_T(i,j,k,bi,bj)  #ifdef GM_VISBECK_VARIABLE_K
269       &                       +Kwy(i,j,k,bi,bj)*deltaTclock  #ifdef ALLOW_AUTODIFF_TAMC
270          GM_Kwz_T(i,j,k,bi,bj)=GM_Kwz_T(i,j,k,bi,bj)  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
      &                       +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  
271  #endif  #endif
272  #endif /* INCLUDE_DIAGNOSTICS_INTERFACE_CODE */        IF ( GM_Visbeck_alpha.NE.0. ) THEN
273    C-     Limit range that KapGM can take
274           DO j=1-Oly+1,sNy+Oly-1
275            DO i=1-Olx+1,sNx+Olx-1
276             VisbeckK(i,j,bi,bj)=
277         &       MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)
278            ENDDO
279         ENDDO         ENDDO
280        ENDDO        ENDIF
281    cph( NEW
282    #ifdef ALLOW_AUTODIFF_TAMC
283    CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
284    #endif
285    cph)
286    #endif /* GM_VISBECK_VARIABLE_K */
287    
288    
289    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
290    
291  #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE  C-- 2nd loop on k : compute Tensor Coeff. at U,V levels.
292        GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock        DO k=1,Nr
293           kp1 = MIN(Nr,k+1)
294           maskp1 = 1. _d 0
295           IF (k.GE.Nr) maskp1 = 0. _d 0
296    
297    #ifdef ALLOW_AUTODIFF_TAMC
298           kkey = (igmkey-1)*Nr + k
299    #if (defined (GM_NON_UNITY_DIAGONAL) || \
300         defined (GM_VISBECK_VARIABLE_K))
301    CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
302    CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
303    CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
304    #endif
305  #endif  #endif
306    
307    C-    express the Tensor in term of Diffusivity (= m**2 / s )
308          DO j=1-Oly+1,sNy+Oly-1
309           DO i=1-Olx+1,sNx+Olx-1
310            Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K
311    #ifdef GM_VISBECK_VARIABLE_K
312         &          + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)    
313    #endif
314            Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)
315            Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)
316            Kwz(i,j,k,bi,bj)= ( GM_isopycK
317    #ifdef GM_VISBECK_VARIABLE_K
318         &                    + VisbeckK(i,j,bi,bj)
319    #endif
320         &                    )*Kwz(i,j,k,bi,bj)
321           ENDDO
322          ENDDO
323    
324    #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
325    
 #ifdef GM_NON_UNITY_DIAGONAL  
326  C     Gradient of Sigma at U points  C     Gradient of Sigma at U points
327        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
328         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
329          SlopeX(i,j)=sigmaX(i,j,km1)          dSigmaDx(i,j)=sigmaX(i,j,k)
330       &          *_maskW(i,j,k,bi,bj)       &          *_maskW(i,j,k,bi,bj)
331          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)
332       &                    +sigmaY(i-1, j ,k) +sigmaY(i, j ,k) )       &                      +sigmaY(i-1, j ,k) +sigmaY(i, j ,k) )
333       &          *_maskW(i,j,k,bi,bj)       &          *_maskW(i,j,k,bi,bj)
334          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 )
335       &                          +sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1) )       &                  +maskp1*(sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1)) )
336       &          *_maskW(i,j,k,bi,bj)       &          *_maskW(i,j,k,bi,bj)
337         ENDDO         ENDDO
338        ENDDO        ENDDO
339    
340    #ifdef ALLOW_AUTODIFF_TAMC
341    CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
342    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
343    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
344    CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte
345    #endif /* ALLOW_AUTODIFF_TAMC */
346    
347  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
348        CALL GMREDI_SLOPE_LIMIT(        CALL GMREDI_SLOPE_LIMIT(
349       I             dSigmadRReal,       O             SlopeX, SlopeY,
350       I             rF(K),       O             SlopeSqr, taperFct,
351       U             SlopeX, SlopeY,       U             dSigmaDr,
352       O             dRdSigmaLtd,       I             dSigmaDx, dSigmaDy,
353         I             ldd97_LrhoW,rC(k),k,
354       I             bi, bj, myThid )       I             bi, bj, myThid )
355    
356        DO j=1-Oly+1,sNy+Oly-1  cph( NEW
357         DO i=1-Olx+1,sNx+Olx-1  #ifdef ALLOW_AUTODIFF_TAMC
358          Kux(i,j,k,bi,bj)=(dSigmaDrReal(i,j)*dRdSigmaLtd(i,j))**2  cph(
359         ENDDO  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
360        ENDDO  CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
361    cph)
362    #endif /* ALLOW_AUTODIFF_TAMC */
363    cph)
364    
365    #ifdef GM_NON_UNITY_DIAGONAL
366            DO j=1-Oly+1,sNy+Oly-1
367             DO i=1-Olx+1,sNx+Olx-1
368              Kux(i,j,k,bi,bj) =
369         &     ( GM_isopycK
370    #ifdef GM_VISBECK_VARIABLE_K
371         &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
372    #endif
373         &     )
374         &     *taperFct(i,j)
375             ENDDO
376            ENDDO
377    #ifdef ALLOW_AUTODIFF_TAMC
378    # ifdef GM_EXCLUDE_CLIPPING
379    CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
380    # endif
381    #endif
382            DO j=1-Oly+1,sNy+Oly-1
383             DO i=1-Olx+1,sNx+Olx-1
384              Kux(i,j,k,bi,bj) = MAX( Kux(i,j,k,bi,bj), GM_Kmin_horiz )
385             ENDDO
386            ENDDO
387    #endif /* GM_NON_UNITY_DIAGONAL */
388    
389    #ifdef GM_EXTRA_DIAGONAL
390    
391    #ifdef ALLOW_AUTODIFF_TAMC
392    CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
393    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
394    #endif /* ALLOW_AUTODIFF_TAMC */
395          IF (GM_ExtraDiag) THEN
396            DO j=1-Oly+1,sNy+Oly-1
397             DO i=1-Olx+1,sNx+Olx-1
398              Kuz(i,j,k,bi,bj) =
399         &     ( GM_isopycK - GM_skewflx*GM_background_K
400    #ifdef GM_VISBECK_VARIABLE_K
401         &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*GM_advect
402    #endif
403         &     )*SlopeX(i,j)*taperFct(i,j)
404             ENDDO
405            ENDDO
406          ENDIF
407    #endif /* GM_EXTRA_DIAGONAL */
408    
409  C     Gradient of Sigma at V points  C     Gradient of Sigma at V points
410        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
411         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
412          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)
413       &                    +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k) )       &                    +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k) )
414       &          *_maskS(i,j,k,bi,bj)       &          *_maskS(i,j,k,bi,bj)
415          SlopeY(i,j)=sigmaY(i,j,km1)          dSigmaDy(i,j)=sigmaY(i,j,k)
416       &          *_maskS(i,j,k,bi,bj)       &          *_maskS(i,j,k,bi,bj)
417          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 )
418       &                          +sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1) )       &                  +maskp1*(sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1)) )
419       &          *_maskS(i,j,k,bi,bj)       &          *_maskS(i,j,k,bi,bj)
420         ENDDO         ENDDO
421        ENDDO        ENDDO
422    
423    #ifdef ALLOW_AUTODIFF_TAMC
424    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
425    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
426    CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte
427    #endif /* ALLOW_AUTODIFF_TAMC */
428    
429  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
430        CALL GMREDI_SLOPE_LIMIT(        CALL GMREDI_SLOPE_LIMIT(
431       I             dSigmadRReal,       O             SlopeX, SlopeY,
432       I             rF(K),       O             SlopeSqr, taperFct,
433       U             SlopeX, SlopeY,       U             dSigmaDr,
434       O             dRdSigmaLtd,       I             dSigmaDx, dSigmaDy,
435         I             ldd97_LrhoS,rC(k),k,
436       I             bi, bj, myThid )       I             bi, bj, myThid )
437    
438        DO j=1-Oly+1,sNy+Oly-1  cph(
439         DO i=1-Olx+1,sNx+Olx-1  #ifdef ALLOW_AUTODIFF_TAMC
440          Kvy(i,j,k,bi,bj)=(dSigmaDrReal(i,j)*dRdSigmaLtd(i,j))**2  cph(
441         ENDDO  CADJ STORE taperfct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
442        ENDDO  cph)
443    #endif /* ALLOW_AUTODIFF_TAMC */
444    cph)
445    
446    #ifdef GM_NON_UNITY_DIAGONAL
447            DO j=1-Oly+1,sNy+Oly-1
448             DO i=1-Olx+1,sNx+Olx-1
449              Kvy(i,j,k,bi,bj) =
450         &     ( GM_isopycK
451    #ifdef GM_VISBECK_VARIABLE_K
452         &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
453    #endif
454         &     )
455         &     *taperFct(i,j)
456             ENDDO
457            ENDDO
458    #ifdef ALLOW_AUTODIFF_TAMC
459    # ifdef GM_EXCLUDE_CLIPPING
460    CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
461    # endif
462    #endif
463            DO j=1-Oly+1,sNy+Oly-1
464             DO i=1-Olx+1,sNx+Olx-1
465              Kvy(i,j,k,bi,bj) = MAX( Kvy(i,j,k,bi,bj), GM_Kmin_horiz )
466             ENDDO
467            ENDDO
468  #endif /* GM_NON_UNITY_DIAGONAL */  #endif /* GM_NON_UNITY_DIAGONAL */
469    
470    #ifdef GM_EXTRA_DIAGONAL
471    
472    #ifdef ALLOW_AUTODIFF_TAMC
473    CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
474    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
475    #endif /* ALLOW_AUTODIFF_TAMC */
476          IF (GM_ExtraDiag) THEN
477            DO j=1-Oly+1,sNy+Oly-1
478             DO i=1-Olx+1,sNx+Olx-1
479              Kvz(i,j,k,bi,bj) =
480         &     ( GM_isopycK - GM_skewflx*GM_background_K
481    #ifdef GM_VISBECK_VARIABLE_K
482         &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*GM_advect
483    #endif
484         &     )*SlopeY(i,j)*taperFct(i,j)
485             ENDDO
486            ENDDO
487          ENDIF
488    #endif /* GM_EXTRA_DIAGONAL */
489    
490    #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
491    
492    C-- end 2nd loop on vertical level index k
493          ENDDO
494    
495    
496    #ifdef GM_BOLUS_ADVEC
497          IF (GM_AdvForm) THEN
498            CALL GMREDI_CALC_PSI_B(
499         I             bi, bj, iMin, iMax, jMin, jMax,
500         I             sigmaX, sigmaY, sigmaR,
501         I             ldd97_LrhoW, ldd97_LrhoS,
502         I             myThid )
503          ENDIF
504    #endif
505    
506    #ifdef ALLOW_TIMEAVE
507    C--   Time-average
508          IF ( taveFreq.GT.0. ) THEN
509    
510             CALL TIMEAVE_CUMULATE( GM_Kwx_T, Kwx, Nr,
511         &                          deltaTclock, bi, bj, myThid )
512             CALL TIMEAVE_CUMULATE( GM_Kwy_T, Kwy, Nr,
513         &                          deltaTclock, bi, bj, myThid )
514             CALL TIMEAVE_CUMULATE( GM_Kwz_T, Kwz, Nr,
515         &                          deltaTclock, bi, bj, myThid )
516    #ifdef GM_VISBECK_VARIABLE_K
517           IF ( GM_Visbeck_alpha.NE.0. ) THEN
518             CALL TIMEAVE_CUMULATE( Visbeck_K_T, VisbeckK, 1,
519         &                          deltaTclock, bi, bj, myThid )
520           ENDIF
521    #endif
522    #ifdef GM_BOLUS_ADVEC
523           IF ( GM_AdvForm ) THEN
524             CALL TIMEAVE_CUMULATE( GM_PsiXtave, GM_PsiX, Nr,
525         &                          deltaTclock, bi, bj, myThid )
526             CALL TIMEAVE_CUMULATE( GM_PsiYtave, GM_PsiY, Nr,
527         &                          deltaTclock, bi, bj, myThid )
528           ENDIF
529    #endif
530           DO k=1,Nr
531             GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock
532           ENDDO
533    
534          ENDIF
535    #endif /* ALLOW_TIMEAVE */
536    
537  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
538    
# Line 221  C     Calculate slopes for use in tensor Line 541  C     Calculate slopes for use in tensor
541    
542    
543        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
544       I             bi, bj, iMin, iMax, jMin, jMax, K,       I             bi, bj, iMin, iMax, jMin, jMax,
545       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
546       I             myThid )       I             myThid )
547  C     /==========================================================\  C     /==========================================================\
# Line 233  C     \================================= Line 553  C     \=================================
553    
554  C     == Global variables ==  C     == Global variables ==
555  #include "SIZE.h"  #include "SIZE.h"
 #include "GRID.h"  
 #include "DYNVARS.h"  
556  #include "EEPARAMS.h"  #include "EEPARAMS.h"
 #include "PARAMS.h"  
557  #include "GMREDI.h"  #include "GMREDI.h"
558    
559  C     == Routine arguments ==  C     == Routine arguments ==
# Line 244  C Line 561  C
561        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
562        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
563        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
564        INTEGER bi,bj,iMin,iMax,jMin,jMax,K        INTEGER bi,bj,iMin,iMax,jMin,jMax
565        INTEGER myThid        INTEGER myThid
566  CEndOfInterface  CEndOfInterface
567    
568        INTEGER i, j        INTEGER i, j, k
569    
570  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
571    
572        DO j=1-Oly+1,sNy+Oly-1        DO k=1,Nr
573         DO i=1-Olx+1,sNx+Olx-1         DO j=1-Oly+1,sNy+Oly-1
574          Kwx(i,j,k,bi,bj) = 0.0          DO i=1-Olx+1,sNx+Olx-1
575          Kwy(i,j,k,bi,bj) = 0.0           Kwx(i,j,k,bi,bj) = 0.0
576          Kwz(i,j,k,bi,bj) = 0.0           Kwy(i,j,k,bi,bj) = 0.0
577             Kwz(i,j,k,bi,bj) = 0.0
578            ENDDO
579         ENDDO         ENDDO
580        ENDDO        ENDDO
581  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
582    
583        end        RETURN
584          END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22