/[MITgcm]/MITgcm/pkg/gmredi/gmredi_calc_tensor.F
ViewVC logotype

Annotation of /MITgcm/pkg/gmredi/gmredi_calc_tensor.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.25 - (hide annotations) (download)
Wed Feb 7 00:01:15 2007 UTC (17 years, 4 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint59, checkpoint58y_post, checkpoint58w_post, checkpoint59a, checkpoint59b, checkpoint58v_post, checkpoint58x_post
Changes since 1.24: +13 -1 lines 
Updating for case ALLOW_KAPGM_CONTROL
1 heimbach 1.25 C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_calc_tensor.F,v 1.24 2006/06/20 22:55:08 jmc Exp $
2 heimbach 1.13 C $Name: $
3 adcroft 1.1
4     #include "GMREDI_OPTIONS.h"
5    
6     CStartOfInterface
7     SUBROUTINE GMREDI_CALC_TENSOR(
8 jmc 1.9 I bi, bj, iMin, iMax, jMin, jMax,
9 adcroft 1.1 I sigmaX, sigmaY, sigmaR,
10     I myThid )
11     C /==========================================================\
12     C | SUBROUTINE GMREDI_CALC_TENSOR |
13     C | o Calculate tensor elements for GM/Redi tensor. |
14     C |==========================================================|
15     C \==========================================================/
16     IMPLICIT NONE
17    
18     C == Global variables ==
19     #include "SIZE.h"
20     #include "GRID.h"
21     #include "DYNVARS.h"
22     #include "EEPARAMS.h"
23     #include "PARAMS.h"
24     #include "GMREDI.h"
25 jmc 1.20 #include "GMREDI_TAVE.h"
26 adcroft 1.1
27 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
28     #include "tamc.h"
29     #include "tamc_keys.h"
30     #endif /* ALLOW_AUTODIFF_TAMC */
31    
32 adcroft 1.1 C == Routine arguments ==
33     C
34     _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
35     _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
36     _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
37 jmc 1.9 INTEGER bi,bj,iMin,iMax,jMin,jMax
38 adcroft 1.1 INTEGER myThid
39     CEndOfInterface
40    
41     #ifdef ALLOW_GMREDI
42    
43     C == Local variables ==
44 jmc 1.15 INTEGER i,j,k,kp1
45 adcroft 1.1 _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
46     _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47 heimbach 1.12 _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48     _RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49 jmc 1.19 _RL dSigmaDr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50 jmc 1.8 _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51     _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
52 jmc 1.15 _RL maskp1, Kgm_tmp
53 jmc 1.19 _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 adcroft 1.1
58     #ifdef GM_VISBECK_VARIABLE_K
59 heimbach 1.14 _RL deltaH,zero_rs
60     PARAMETER(zero_rs=0.D0)
61 adcroft 1.1 _RL N2,SN
62 heimbach 1.10 _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
63 adcroft 1.1 #endif
64    
65 jmc 1.22 #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 jmc 1.19 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
75    
76 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
77     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 heimbach 1.12 igmkey = (act1 + 1) + act2*max1
85 heimbach 1.10 & + act3*max1*max2
86     & + act4*max1*max2*max3
87     #endif /* ALLOW_AUTODIFF_TAMC */
88    
89 jmc 1.22 #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 heimbach 1.12 #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 jmc 1.19 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 jmc 1.9 DO k=2,Nr
163     C-- 1rst loop on k : compute Tensor Coeff. at W points.
164 adcroft 1.1
165     #ifdef ALLOW_AUTODIFF_TAMC
166 heimbach 1.12 kkey = (igmkey-1)*Nr + k
167 heimbach 1.10 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 heimbach 1.12 dSigmaDx(i,j) = 0. _d 0
172     dSigmaDy(i,j) = 0. _d 0
173 jmc 1.19 dSigmaDr(i,j) = 0. _d 0
174 heimbach 1.10 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 heimbach 1.12 # 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 heimbach 1.10 ENDDO
192     ENDDO
193 adcroft 1.1 #endif
194 heimbach 1.10
195 adcroft 1.1 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
198 jmc 1.15 dSigmaDx(i,j)=op25*( sigmaX(i+1, j ,k-1) +sigmaX(i,j,k-1)
199 adcroft 1.1 & +sigmaX(i+1, j , k ) +sigmaX(i,j, k ) )
200 jmc 1.15 & *maskC(i,j,k,bi,bj)
201     dSigmaDy(i,j)=op25*( sigmaY( i ,j+1,k-1) +sigmaY(i,j,k-1)
202 adcroft 1.1 & +sigmaY( i ,j+1, k ) +sigmaY(i,j, k ) )
203 jmc 1.15 & *maskC(i,j,k,bi,bj)
204 jmc 1.19 dSigmaDr(i,j)=sigmaR(i,j,k)
205 adcroft 1.1 ENDDO
206     ENDDO
207    
208 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
209 heimbach 1.12 CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
210     CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
211 jmc 1.19 CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
212 heimbach 1.10 #endif /* ALLOW_AUTODIFF_TAMC */
213    
214 adcroft 1.1 C Calculate slopes for use in tensor, taper and/or clip
215     CALL GMREDI_SLOPE_LIMIT(
216 jmc 1.19 O SlopeX, SlopeY,
217 jmc 1.8 O SlopeSqr, taperFct,
218 jmc 1.19 U dSigmaDr,
219     I dSigmaDx, dSigmaDy,
220     I ldd97_LrhoC,rF(k),k,
221 adcroft 1.1 I bi, bj, myThid )
222    
223     DO j=1-Oly+1,sNy+Oly-1
224     DO i=1-Olx+1,sNx+Olx-1
225    
226     C Mask Iso-neutral slopes
227 jmc 1.15 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)
229     SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)
230 heimbach 1.10
231     ENDDO
232     ENDDO
233    
234     #ifdef ALLOW_AUTODIFF_TAMC
235 heimbach 1.16 CADJ STORE SlopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
236     CADJ STORE SlopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
237 heimbach 1.14 CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
238 jmc 1.19 CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
239     CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
240 heimbach 1.10 #endif /* ALLOW_AUTODIFF_TAMC */
241    
242     DO j=1-Oly+1,sNy+Oly-1
243     DO i=1-Olx+1,sNx+Olx-1
244 adcroft 1.1
245 jmc 1.9 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 jmc 1.8 Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
249 adcroft 1.1
250     #ifdef GM_VISBECK_VARIABLE_K
251 jmc 1.8
252     C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K
253 edhill 1.18 C but do not know if *taperFct (or **2 ?) is necessary
254 heimbach 1.10 Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)
255 jmc 1.8
256 adcroft 1.1 C-- Depth average of M^2/N^2 * N
257    
258     C Calculate terms for mean Richardson number
259     C which is used in the "variable K" parameterisaton.
260     C Distance between interface above layer and the integration depth
261     deltaH=abs(GM_Visbeck_depth)-abs(rF(k))
262     C If positive we limit this to the layer thickness
263     deltaH=min(deltaH,drF(k))
264     C If negative then we are below the integration level
265     deltaH=max(deltaH,zero_rs)
266     C Now we convert deltaH to a non-dimensional fraction
267     deltaH=deltaH/GM_Visbeck_depth
268    
269 jmc 1.8 IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.
270 jmc 1.19 IF ( Ssq(i,j).NE.0. .AND. dSigmaDr(i,j).NE.0. ) THEN
271     N2= -Gravity*recip_RhoConst*dSigmaDr(i,j)
272 heimbach 1.10 SN=sqrt(Ssq(i,j)*N2)
273 heimbach 1.3 VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH
274 adcroft 1.1 & *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN
275 jmc 1.8 ENDIF
276 adcroft 1.1
277 jmc 1.9 #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 adcroft 1.1
286 jmc 1.9 #ifdef GM_VISBECK_VARIABLE_K
287 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
288     CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
289     #endif
290 jmc 1.9 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
298     ENDIF
299 heimbach 1.16 cph( NEW
300     #ifdef ALLOW_AUTODIFF_TAMC
301     CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
302     #endif
303     cph)
304 adcroft 1.1 #endif /* GM_VISBECK_VARIABLE_K */
305    
306    
307 jmc 1.9 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
308 heimbach 1.10
309 jmc 1.9 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 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
316     kkey = (igmkey-1)*Nr + k
317 heimbach 1.16 #if (defined (GM_NON_UNITY_DIAGONAL) || \
318     defined (GM_VISBECK_VARIABLE_K))
319 heimbach 1.14 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 heimbach 1.12 #endif
323     #endif
324    
325 jmc 1.9 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 heimbach 1.25 #ifdef ALLOW_KAPGM_CONTROL
329     Kgm_tmp = GM_isopycK + GM_skewflx*kapgm(i,j,k,bi,bj)
330     #else
331 jmc 1.9 Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K
332 heimbach 1.25 #endif
333 jmc 1.9 #ifdef GM_VISBECK_VARIABLE_K
334 heimbach 1.16 & + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)
335 jmc 1.9 #endif
336     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)
338     Kwz(i,j,k,bi,bj)= ( GM_isopycK
339 adcroft 1.1 #ifdef GM_VISBECK_VARIABLE_K
340 jmc 1.9 & + VisbeckK(i,j,bi,bj)
341 adcroft 1.1 #endif
342 jmc 1.9 & )*Kwz(i,j,k,bi,bj)
343 adcroft 1.1 ENDDO
344     ENDDO
345 adcroft 1.4
346 jmc 1.9 #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
347 adcroft 1.1
348     C Gradient of Sigma at U points
349     DO j=1-Oly+1,sNy+Oly-1
350     DO i=1-Olx+1,sNx+Olx-1
351 heimbach 1.12 dSigmaDx(i,j)=sigmaX(i,j,k)
352 adcroft 1.1 & *_maskW(i,j,k,bi,bj)
353 heimbach 1.14 dSigmaDy(i,j)=op25*( sigmaY(i-1,j+1,k) +sigmaY(i,j+1,k)
354 heimbach 1.12 & +sigmaY(i-1, j ,k) +sigmaY(i, j ,k) )
355 adcroft 1.1 & *_maskW(i,j,k,bi,bj)
356 jmc 1.19 dSigmaDr(i,j)=op25*( sigmaR(i-1,j, k ) +sigmaR(i,j, k )
357 jmc 1.9 & +maskp1*(sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1)) )
358 jmc 1.15 & *_maskW(i,j,k,bi,bj)
359 adcroft 1.1 ENDDO
360     ENDDO
361    
362 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
363 heimbach 1.17 CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
364 heimbach 1.12 CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
365     CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
366 jmc 1.19 CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
367 heimbach 1.12 #endif /* ALLOW_AUTODIFF_TAMC */
368    
369 adcroft 1.1 C Calculate slopes for use in tensor, taper and/or clip
370     CALL GMREDI_SLOPE_LIMIT(
371 jmc 1.19 O SlopeX, SlopeY,
372 jmc 1.8 O SlopeSqr, taperFct,
373 jmc 1.19 U dSigmaDr,
374     I dSigmaDx, dSigmaDy,
375     I ldd97_LrhoW,rC(k),k,
376 adcroft 1.1 I bi, bj, myThid )
377    
378 heimbach 1.16 cph( NEW
379     #ifdef ALLOW_AUTODIFF_TAMC
380     cph(
381 heimbach 1.17 CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
382 jmc 1.19 CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
383 heimbach 1.16 cph)
384     #endif /* ALLOW_AUTODIFF_TAMC */
385     cph)
386    
387 jmc 1.9 #ifdef GM_NON_UNITY_DIAGONAL
388     DO j=1-Oly+1,sNy+Oly-1
389     DO i=1-Olx+1,sNx+Olx-1
390     Kux(i,j,k,bi,bj) =
391     & ( GM_isopycK
392     #ifdef GM_VISBECK_VARIABLE_K
393 heimbach 1.14 & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
394 jmc 1.9 #endif
395 heimbach 1.10 & )
396     & *taperFct(i,j)
397     ENDDO
398     ENDDO
399 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
400 heimbach 1.16 # ifdef GM_EXCLUDE_CLIPPING
401 heimbach 1.12 CADJ STORE Kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
402     # endif
403     #endif
404 heimbach 1.10 DO j=1-Oly+1,sNy+Oly-1
405     DO i=1-Olx+1,sNx+Olx-1
406 jmc 1.9 Kux(i,j,k,bi,bj) = MAX( Kux(i,j,k,bi,bj), GM_Kmin_horiz )
407     ENDDO
408     ENDDO
409     #endif /* GM_NON_UNITY_DIAGONAL */
410    
411     #ifdef GM_EXTRA_DIAGONAL
412 heimbach 1.12
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 jmc 1.9 IF (GM_ExtraDiag) THEN
418     DO j=1-Oly+1,sNy+Oly-1
419     DO i=1-Olx+1,sNx+Olx-1
420     Kuz(i,j,k,bi,bj) =
421 heimbach 1.25 #ifdef ALLOW_KAPGM_CONTROL
422     & ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)
423     #else
424 jmc 1.9 & ( GM_isopycK - GM_skewflx*GM_background_K
425 heimbach 1.25 #endif
426 jmc 1.9 #ifdef GM_VISBECK_VARIABLE_K
427 heimbach 1.14 & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*GM_advect
428 jmc 1.9 #endif
429     & )*SlopeX(i,j)*taperFct(i,j)
430     ENDDO
431     ENDDO
432     ENDIF
433     #endif /* GM_EXTRA_DIAGONAL */
434 adcroft 1.1
435 jmc 1.22 #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 heimbach 1.23 tmp1k(i,j) = dyG(i,j,bi,bj)*drF(k)
465     & * _hFacW(i,j,k,bi,bj)
466 jmc 1.22 & * 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 adcroft 1.1 C Gradient of Sigma at V points
474     DO j=1-Oly+1,sNy+Oly-1
475     DO i=1-Olx+1,sNx+Olx-1
476 heimbach 1.14 dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
477 adcroft 1.1 & +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k) )
478     & *_maskS(i,j,k,bi,bj)
479 heimbach 1.12 dSigmaDy(i,j)=sigmaY(i,j,k)
480 adcroft 1.1 & *_maskS(i,j,k,bi,bj)
481 jmc 1.19 dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k ) +sigmaR(i,j, k )
482 jmc 1.9 & +maskp1*(sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1)) )
483 jmc 1.15 & *_maskS(i,j,k,bi,bj)
484 adcroft 1.1 ENDDO
485     ENDDO
486    
487 heimbach 1.12 #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 jmc 1.19 CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
491 heimbach 1.12 #endif /* ALLOW_AUTODIFF_TAMC */
492    
493 adcroft 1.1 C Calculate slopes for use in tensor, taper and/or clip
494     CALL GMREDI_SLOPE_LIMIT(
495 jmc 1.19 O SlopeX, SlopeY,
496 jmc 1.8 O SlopeSqr, taperFct,
497 jmc 1.19 U dSigmaDr,
498     I dSigmaDx, dSigmaDy,
499     I ldd97_LrhoS,rC(k),k,
500 adcroft 1.1 I bi, bj, myThid )
501    
502 heimbach 1.16 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 jmc 1.9 #ifdef GM_NON_UNITY_DIAGONAL
511     DO j=1-Oly+1,sNy+Oly-1
512     DO i=1-Olx+1,sNx+Olx-1
513     Kvy(i,j,k,bi,bj) =
514     & ( GM_isopycK
515     #ifdef GM_VISBECK_VARIABLE_K
516 heimbach 1.14 & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
517 jmc 1.9 #endif
518 heimbach 1.10 & )
519     & *taperFct(i,j)
520     ENDDO
521     ENDDO
522 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
523 heimbach 1.16 # ifdef GM_EXCLUDE_CLIPPING
524 heimbach 1.12 CADJ STORE Kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
525     # endif
526     #endif
527 heimbach 1.10 DO j=1-Oly+1,sNy+Oly-1
528     DO i=1-Olx+1,sNx+Olx-1
529 jmc 1.9 Kvy(i,j,k,bi,bj) = MAX( Kvy(i,j,k,bi,bj), GM_Kmin_horiz )
530     ENDDO
531     ENDDO
532     #endif /* GM_NON_UNITY_DIAGONAL */
533    
534     #ifdef GM_EXTRA_DIAGONAL
535 heimbach 1.12
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 jmc 1.9 IF (GM_ExtraDiag) THEN
541     DO j=1-Oly+1,sNy+Oly-1
542     DO i=1-Olx+1,sNx+Olx-1
543     Kvz(i,j,k,bi,bj) =
544 heimbach 1.25 #ifdef ALLOW_KAPGM_CONTROL
545     & ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)
546     #else
547 jmc 1.9 & ( GM_isopycK - GM_skewflx*GM_background_K
548 heimbach 1.25 #endif
549 jmc 1.9 #ifdef GM_VISBECK_VARIABLE_K
550 heimbach 1.14 & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*GM_advect
551 jmc 1.9 #endif
552     & )*SlopeY(i,j)*taperFct(i,j)
553     ENDDO
554     ENDDO
555     ENDIF
556     #endif /* GM_EXTRA_DIAGONAL */
557    
558 jmc 1.22 #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 heimbach 1.23 tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)
588     & * _hFacS(i,j,k,bi,bj)
589 jmc 1.22 & * 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 jmc 1.9 #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
597    
598     C-- end 2nd loop on vertical level index k
599     ENDDO
600 adcroft 1.1
601    
602 jmc 1.9 #ifdef GM_BOLUS_ADVEC
603     IF (GM_AdvForm) THEN
604     CALL GMREDI_CALC_PSI_B(
605     I bi, bj, iMin, iMax, jMin, jMax,
606     I sigmaX, sigmaY, sigmaR,
607 jmc 1.19 I ldd97_LrhoW, ldd97_LrhoS,
608 jmc 1.9 I myThid )
609     ENDIF
610     #endif
611 adcroft 1.1
612 jmc 1.19 #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 jmc 1.20 #ifdef ALLOW_DIAGNOSTICS
644     IF ( useDiagnostics ) THEN
645 jmc 1.24 CALL GMREDI_DIAGNOSTICS_FILL(bi,bj,myThid)
646 jmc 1.20 ENDIF
647     #endif /* ALLOW_DIAGNOSTICS */
648    
649 adcroft 1.1 #endif /* ALLOW_GMREDI */
650    
651     RETURN
652     END
653 heimbach 1.2
654    
655     SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
656 jmc 1.9 I bi, bj, iMin, iMax, jMin, jMax,
657 heimbach 1.2 I sigmaX, sigmaY, sigmaR,
658     I myThid )
659     C /==========================================================\
660     C | SUBROUTINE GMREDI_CALC_TENSOR |
661     C | o Calculate tensor elements for GM/Redi tensor. |
662     C |==========================================================|
663     C \==========================================================/
664     IMPLICIT NONE
665    
666     C == Global variables ==
667     #include "SIZE.h"
668     #include "EEPARAMS.h"
669     #include "GMREDI.h"
670    
671     C == Routine arguments ==
672     C
673     _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
674     _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
675     _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
676 jmc 1.9 INTEGER bi,bj,iMin,iMax,jMin,jMax
677 heimbach 1.2 INTEGER myThid
678     CEndOfInterface
679    
680 jmc 1.9 INTEGER i, j, k
681 heimbach 1.2
682     #ifdef ALLOW_GMREDI
683    
684 jmc 1.9 DO k=1,Nr
685     DO j=1-Oly+1,sNy+Oly-1
686     DO i=1-Olx+1,sNx+Olx-1
687     Kwx(i,j,k,bi,bj) = 0.0
688     Kwy(i,j,k,bi,bj) = 0.0
689     Kwz(i,j,k,bi,bj) = 0.0
690     ENDDO
691 heimbach 1.2 ENDDO
692     ENDDO
693     #endif /* ALLOW_GMREDI */
694    
695 jmc 1.9 RETURN
696     END
  ViewVC Help
Powered by ViewVC 1.1.22