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

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

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

Revision 1.26 - (show annotations) (download)
Thu May 24 22:34:38 2007 UTC (17 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint59d, checkpoint59c
Changes since 1.25: +140 -77 lines 
change calculation of Visbeck K:
 - no longer depend on tapering scheme ; instead, compute directly the slope
   (limited by GM_maxSlope).
 - fix vertical averaging (wrong by 1 level).
1 C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_calc_tensor.F,v 1.25 2007/02/07 00:01:15 heimbach Exp $
2 C $Name: $
3
4 #include "GMREDI_OPTIONS.h"
5 #undef OLD_VISBECK_CALC
6
7 CStartOfInterface
8 SUBROUTINE GMREDI_CALC_TENSOR(
9 I bi, bj, iMin, iMax, jMin, jMax,
10 I sigmaX, sigmaY, sigmaR,
11 I myThid )
12 C /==========================================================\
13 C | SUBROUTINE GMREDI_CALC_TENSOR |
14 C | o Calculate tensor elements for GM/Redi tensor. |
15 C |==========================================================|
16 C \==========================================================/
17 IMPLICIT NONE
18
19 C == Global variables ==
20 #include "SIZE.h"
21 #include "GRID.h"
22 #include "DYNVARS.h"
23 #include "EEPARAMS.h"
24 #include "PARAMS.h"
25 #include "GMREDI.h"
26 #include "GMREDI_TAVE.h"
27
28 #ifdef ALLOW_AUTODIFF_TAMC
29 #include "tamc.h"
30 #include "tamc_keys.h"
31 #endif /* ALLOW_AUTODIFF_TAMC */
32
33 C == Routine arguments ==
34 C
35 _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
36 _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
37 _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
38 INTEGER bi,bj,iMin,iMax,jMin,jMax
39 INTEGER myThid
40 CEndOfInterface
41
42 #ifdef ALLOW_GMREDI
43
44 C == Local variables ==
45 INTEGER i,j,k,kp1
46 _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47 _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48 _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49 _RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50 _RL dSigmaDr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51 _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
52 _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
53 _RL maskp1, Kgm_tmp
54 _RL ldd97_LrhoC(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
55 _RL ldd97_LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56 _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57 _RL Cspd, LrhoInf, LrhoSup, fCoriLoc
58
59 #ifdef GM_VISBECK_VARIABLE_K
60 #ifdef OLD_VISBECK_CALC
61 _RL deltaH,zero_rs
62 PARAMETER(zero_rs=0.D0)
63 _RL N2,SN
64 _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
65 #else
66 _RL dSigmaH
67 _RL deltaH, integrDepth
68 _RL Sloc, M2loc, SNloc
69 #endif
70 #endif
71
72 #ifdef ALLOW_DIAGNOSTICS
73 LOGICAL doDiagRediFlx
74 LOGICAL DIAGNOSTICS_IS_ON
75 EXTERNAL DIAGNOSTICS_IS_ON
76 INTEGER km1
77 _RL dTdz
78 _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79 #endif
80
81 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
82
83 #ifdef ALLOW_AUTODIFF_TAMC
84 act1 = bi - myBxLo(myThid)
85 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
86 act2 = bj - myByLo(myThid)
87 max2 = myByHi(myThid) - myByLo(myThid) + 1
88 act3 = myThid - 1
89 max3 = nTx*nTy
90 act4 = ikey_dynamics - 1
91 igmkey = (act1 + 1) + act2*max1
92 & + act3*max1*max2
93 & + act4*max1*max2*max3
94 #endif /* ALLOW_AUTODIFF_TAMC */
95
96 #ifdef ALLOW_DIAGNOSTICS
97 doDiagRediFlx = .FALSE.
98 IF ( useDiagnostics ) THEN
99 doDiagRediFlx = DIAGNOSTICS_IS_ON('GM_KuzTz', myThid )
100 doDiagRediFlx = doDiagRediFlx .OR.
101 & DIAGNOSTICS_IS_ON('GM_KvzTz', myThid )
102 ENDIF
103 #endif
104
105 #ifdef GM_VISBECK_VARIABLE_K
106 DO j=1-Oly,sNy+Oly
107 DO i=1-Olx,sNx+Olx
108 VisbeckK(i,j,bi,bj) = 0. _d 0
109 ENDDO
110 ENDDO
111 #endif
112
113 C-- set ldd97_Lrho (for tapering scheme ldd97):
114 IF (GM_taper_scheme.EQ.'ldd97') THEN
115 Cspd = 2. _d 0
116 LrhoInf = 15. _d 3
117 LrhoSup = 100. _d 3
118 C- Tracer point location (center):
119 DO j=1-Oly,sNy+Oly
120 DO i=1-Olx,sNx+Olx
121 IF (fCori(i,j,bi,bj).NE.0.) THEN
122 ldd97_LrhoC(i,j) = Cspd/ABS(fCori(i,j,bi,bj))
123 ELSE
124 ldd97_LrhoC(i,j) = LrhoSup
125 ENDIF
126 ldd97_LrhoC(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoC(i,j),LrhoSup))
127 ENDDO
128 ENDDO
129 C- U point location (West):
130 DO j=1-Oly,sNy+Oly
131 ldd97_LrhoW(1-Olx,j) = LrhoSup
132 DO i=1-Olx+1,sNx+Olx
133 fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))
134 IF (fCoriLoc.NE.0.) THEN
135 ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)
136 ELSE
137 ldd97_LrhoW(i,j) = LrhoSup
138 ENDIF
139 ldd97_LrhoW(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoW(i,j),LrhoSup))
140 ENDDO
141 ENDDO
142 C- V point location (South):
143 DO i=1-Olx+1,sNx+Olx
144 ldd97_LrhoS(i,1-Oly) = LrhoSup
145 ENDDO
146 DO j=1-Oly+1,sNy+Oly
147 DO i=1-Olx,sNx+Olx
148 fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))
149 IF (fCoriLoc.NE.0.) THEN
150 ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)
151 ELSE
152 ldd97_LrhoS(i,j) = LrhoSup
153 ENDIF
154 ldd97_LrhoS(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoS(i,j),LrhoSup))
155 ENDDO
156 ENDDO
157 ELSE
158 C- Just initialize to zero (not use anyway)
159 DO j=1-Oly,sNy+Oly
160 DO i=1-Olx,sNx+Olx
161 ldd97_LrhoC(i,j) = 0. _d 0
162 ldd97_LrhoW(i,j) = 0. _d 0
163 ldd97_LrhoS(i,j) = 0. _d 0
164 ENDDO
165 ENDDO
166 ENDIF
167 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
168
169 DO k=2,Nr
170 C-- 1rst loop on k : compute Tensor Coeff. at W points.
171
172 #ifdef ALLOW_AUTODIFF_TAMC
173 kkey = (igmkey-1)*Nr + k
174 DO j=1-Oly,sNy+Oly
175 DO i=1-Olx,sNx+Olx
176 SlopeX(i,j) = 0. _d 0
177 SlopeY(i,j) = 0. _d 0
178 dSigmaDx(i,j) = 0. _d 0
179 dSigmaDy(i,j) = 0. _d 0
180 dSigmaDr(i,j) = 0. _d 0
181 SlopeSqr(i,j) = 0. _d 0
182 taperFct(i,j) = 0. _d 0
183 Kwx(i,j,k,bi,bj) = 0. _d 0
184 Kwy(i,j,k,bi,bj) = 0. _d 0
185 Kwz(i,j,k,bi,bj) = 0. _d 0
186 # ifdef GM_NON_UNITY_DIAGONAL
187 Kux(i,j,k,bi,bj) = 0. _d 0
188 Kvy(i,j,k,bi,bj) = 0. _d 0
189 # endif
190 # ifdef GM_EXTRA_DIAGONAL
191 Kuz(i,j,k,bi,bj) = 0. _d 0
192 Kvz(i,j,k,bi,bj) = 0. _d 0
193 # endif
194 # ifdef GM_BOLUS_ADVEC
195 GM_PsiX(i,j,k,bi,bj) = 0. _d 0
196 GM_PsiY(i,j,k,bi,bj) = 0. _d 0
197 # endif
198 ENDDO
199 ENDDO
200 #endif
201
202 DO j=1-Oly+1,sNy+Oly-1
203 DO i=1-Olx+1,sNx+Olx-1
204 C Gradient of Sigma at rVel points
205 dSigmaDx(i,j)=op25*( sigmaX(i+1,j,k-1)+sigmaX(i,j,k-1)
206 & +sigmaX(i+1,j, k )+sigmaX(i,j, k )
207 & )*maskC(i,j,k,bi,bj)
208 dSigmaDy(i,j)=op25*( sigmaY(i,j+1,k-1)+sigmaY(i,j,k-1)
209 & +sigmaY(i,j+1, k )+sigmaY(i,j, k )
210 & )*maskC(i,j,k,bi,bj)
211 dSigmaDr(i,j)=sigmaR(i,j,k)
212 ENDDO
213 ENDDO
214
215 #ifdef ALLOW_AUTODIFF_TAMC
216 CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
217 CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
218 CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
219 #endif /* ALLOW_AUTODIFF_TAMC */
220
221 #ifdef GM_VISBECK_VARIABLE_K
222 #ifndef OLD_VISBECK_CALC
223 IF ( GM_Visbeck_alpha.GT.0. .AND.
224 & -rC(k-1).LT.GM_Visbeck_depth ) THEN
225
226 C-- Depth average of f/sqrt(Ri) = M^2/N^2 * N
227 C M^2 and N^2 are horizontal & vertical gradient of buoyancy.
228
229 C Calculate terms for mean Richardson number which is used
230 C in the "variable K" parameterisaton:
231 C compute depth average from surface down to the bottom or
232 C GM_Visbeck_depth, whatever is the shallower.
233
234 DO j=1-Oly+1,sNy+Oly-1
235 DO i=1-Olx+1,sNx+Olx-1
236 IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
237 integrDepth = -rC( kLowC(i,j,bi,bj) )
238 C- in 2 steps to avoid mix of RS & RL type in min fct. arguments
239 integrDepth = MIN( integrDepth, GM_Visbeck_depth )
240 C Distance between level center above and the integration depth
241 deltaH = integrDepth + rC(k-1)
242 C If negative then we are below the integration level
243 C (cannot be the case with 2 conditions on maskC & -rC(k-1))
244 C If positive we limit this to the distance from center above
245 deltaH = MIN( deltaH, drC(k) )
246 C Now we convert deltaH to a non-dimensional fraction
247 deltaH = deltaH/( integrDepth+rC(1) )
248
249 C-- compute: ( M^2 * S )^1/2 (= M^2 / N since S=M^2/N^2 )
250 dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)
251 & + dSigmaDy(i,j)*dSigmaDy(i,j)
252 IF ( dSigmaH .GT. 0. _d 0 ) THEN
253 dSigmaH = SQRT( dSigmaH )
254 C- compute slope, limited by GM_maxSlope:
255 IF ( -dSigmaDr(i,j).GT.dSigmaH*GM_rMaxSlope ) THEN
256 Sloc = dSigmaH / ( -dSigmaDr(i,j) )
257 ELSE
258 Sloc = GM_maxSlope
259 ENDIF
260 M2loc = Gravity*recip_RhoConst*dSigmaH
261 SNloc = SQRT( Sloc*M2loc )
262 ELSE
263 SNloc = 0. _d 0
264 ENDIF
265 VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
266 & +deltaH*GM_Visbeck_alpha
267 & *GM_Visbeck_length*GM_Visbeck_length*SNloc
268 ENDIF
269 ENDDO
270 ENDDO
271 ENDIF
272 #endif /* ndef OLD_VISBECK_CALC */
273 #endif /* GM_VISBECK_VARIABLE_K */
274
275 C Calculate slopes for use in tensor, taper and/or clip
276 CALL GMREDI_SLOPE_LIMIT(
277 O SlopeX, SlopeY,
278 O SlopeSqr, taperFct,
279 U dSigmaDr,
280 I dSigmaDx, dSigmaDy,
281 I ldd97_LrhoC,rF(k),k,
282 I bi, bj, myThid )
283
284 DO j=1-Oly+1,sNy+Oly-1
285 DO i=1-Olx+1,sNx+Olx-1
286 C Mask Iso-neutral slopes
287 SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)
288 SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)
289 SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)
290 ENDDO
291 ENDDO
292
293 #ifdef ALLOW_AUTODIFF_TAMC
294 CADJ STORE SlopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
295 CADJ STORE SlopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
296 CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
297 CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
298 CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
299 #endif /* ALLOW_AUTODIFF_TAMC */
300
301 DO j=1-Oly+1,sNy+Oly-1
302 DO i=1-Olx+1,sNx+Olx-1
303 C Components of Redi/GM tensor
304 Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
305 Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
306 Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
307
308 #ifdef GM_VISBECK_VARIABLE_K
309 #ifdef OLD_VISBECK_CALC
310
311 C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K
312 C but do not know if *taperFct (or **2 ?) is necessary
313 Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)
314
315 C-- Depth average of M^2/N^2 * N
316
317 C Calculate terms for mean Richardson number
318 C which is used in the "variable K" parameterisaton.
319 C Distance between interface above layer and the integration depth
320 deltaH=abs(GM_Visbeck_depth)-abs(rF(k))
321 C If positive we limit this to the layer thickness
322 deltaH=min(deltaH,drF(k))
323 C If negative then we are below the integration level
324 deltaH=max(deltaH,zero_rs)
325 C Now we convert deltaH to a non-dimensional fraction
326 deltaH=deltaH/GM_Visbeck_depth
327
328 IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.
329 IF ( Ssq(i,j).NE.0. .AND. dSigmaDr(i,j).NE.0. ) THEN
330 N2= -Gravity*recip_RhoConst*dSigmaDr(i,j)
331 SN=sqrt(Ssq(i,j)*N2)
332 VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH
333 & *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN
334 ENDIF
335
336 #endif /* OLD_VISBECK_CALC */
337 #endif /* GM_VISBECK_VARIABLE_K */
338 ENDDO
339 ENDDO
340
341 C-- end 1rst loop on vertical level index k
342 ENDDO
343
344
345 #ifdef GM_VISBECK_VARIABLE_K
346 #ifdef ALLOW_AUTODIFF_TAMC
347 CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
348 #endif
349 IF ( GM_Visbeck_alpha.NE.0. ) THEN
350 C- Limit range that KapGM can take
351 DO j=1-Oly+1,sNy+Oly-1
352 DO i=1-Olx+1,sNx+Olx-1
353 VisbeckK(i,j,bi,bj)=
354 & MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)
355 ENDDO
356 ENDDO
357 ENDIF
358 cph( NEW
359 #ifdef ALLOW_AUTODIFF_TAMC
360 CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
361 #endif
362 cph)
363 #endif /* GM_VISBECK_VARIABLE_K */
364
365
366 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
367
368 C-- 2nd loop on k : compute Tensor Coeff. at U,V levels.
369 DO k=1,Nr
370 kp1 = MIN(Nr,k+1)
371 maskp1 = 1. _d 0
372 IF (k.GE.Nr) maskp1 = 0. _d 0
373
374 #ifdef ALLOW_AUTODIFF_TAMC
375 kkey = (igmkey-1)*Nr + k
376 #if (defined (GM_NON_UNITY_DIAGONAL) || \
377 defined (GM_VISBECK_VARIABLE_K))
378 CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
379 CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
380 CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
381 #endif
382 #endif
383
384 C- express the Tensor in term of Diffusivity (= m**2 / s )
385 DO j=1-Oly+1,sNy+Oly-1
386 DO i=1-Olx+1,sNx+Olx-1
387 #ifdef ALLOW_KAPGM_CONTROL
388 Kgm_tmp = GM_isopycK + GM_skewflx*kapgm(i,j,k,bi,bj)
389 #else
390 Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K
391 #endif
392 #ifdef GM_VISBECK_VARIABLE_K
393 & + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)
394 #endif
395 Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)
396 Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)
397 Kwz(i,j,k,bi,bj)= ( GM_isopycK
398 #ifdef GM_VISBECK_VARIABLE_K
399 & + VisbeckK(i,j,bi,bj)
400 #endif
401 & )*Kwz(i,j,k,bi,bj)
402 ENDDO
403 ENDDO
404
405 #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
406
407 C Gradient of Sigma at U points
408 DO j=1-Oly+1,sNy+Oly-1
409 DO i=1-Olx+1,sNx+Olx-1
410 dSigmaDx(i,j)=sigmaX(i,j,k)
411 & *_maskW(i,j,k,bi,bj)
412 dSigmaDy(i,j)=op25*( sigmaY(i-1,j+1,k)+sigmaY(i,j+1,k)
413 & +sigmaY(i-1, j ,k)+sigmaY(i, j ,k)
414 & )*_maskW(i,j,k,bi,bj)
415 dSigmaDr(i,j)=op25*( sigmaR(i-1,j, k )+sigmaR(i,j, k )
416 & +(sigmaR(i-1,j,kp1)+sigmaR(i,j,kp1))*maskp1
417 & )*_maskW(i,j,k,bi,bj)
418 ENDDO
419 ENDDO
420
421 #ifdef ALLOW_AUTODIFF_TAMC
422 CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
423 CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
424 CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
425 CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
426 #endif /* ALLOW_AUTODIFF_TAMC */
427
428 C Calculate slopes for use in tensor, taper and/or clip
429 CALL GMREDI_SLOPE_LIMIT(
430 O SlopeX, SlopeY,
431 O SlopeSqr, taperFct,
432 U dSigmaDr,
433 I dSigmaDx, dSigmaDy,
434 I ldd97_LrhoW,rC(k),k,
435 I bi, bj, myThid )
436
437 cph( NEW
438 #ifdef ALLOW_AUTODIFF_TAMC
439 cph(
440 CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
441 CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
442 cph)
443 #endif /* ALLOW_AUTODIFF_TAMC */
444 cph)
445
446 #ifdef GM_NON_UNITY_DIAGONAL
447 c IF ( GM_nonUnitDiag ) THEN
448 DO j=1-Oly+1,sNy+Oly-1
449 DO i=1-Olx+1,sNx+Olx-1
450 Kux(i,j,k,bi,bj) =
451 & ( GM_isopycK
452 #ifdef GM_VISBECK_VARIABLE_K
453 & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
454 #endif
455 & )
456 & *taperFct(i,j)
457 ENDDO
458 ENDDO
459 #ifdef ALLOW_AUTODIFF_TAMC
460 # ifdef GM_EXCLUDE_CLIPPING
461 CADJ STORE Kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
462 # endif
463 #endif
464 DO j=1-Oly+1,sNy+Oly-1
465 DO i=1-Olx+1,sNx+Olx-1
466 Kux(i,j,k,bi,bj) = MAX( Kux(i,j,k,bi,bj), GM_Kmin_horiz )
467 ENDDO
468 ENDDO
469 c ENDIF
470 #endif /* GM_NON_UNITY_DIAGONAL */
471
472 #ifdef GM_EXTRA_DIAGONAL
473
474 #ifdef ALLOW_AUTODIFF_TAMC
475 CADJ STORE SlopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
476 CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
477 #endif /* ALLOW_AUTODIFF_TAMC */
478 IF (GM_ExtraDiag) THEN
479 DO j=1-Oly+1,sNy+Oly-1
480 DO i=1-Olx+1,sNx+Olx-1
481 Kuz(i,j,k,bi,bj) =
482 #ifdef ALLOW_KAPGM_CONTROL
483 & ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)
484 #else
485 & ( GM_isopycK - GM_skewflx*GM_background_K
486 #endif
487 #ifdef GM_VISBECK_VARIABLE_K
488 & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*GM_advect
489 #endif
490 & )*SlopeX(i,j)*taperFct(i,j)
491 ENDDO
492 ENDDO
493 ENDIF
494 #endif /* GM_EXTRA_DIAGONAL */
495
496 #ifdef ALLOW_DIAGNOSTICS
497 IF (doDiagRediFlx) THEN
498 km1 = MAX(k-1,1)
499 DO j=1,sNy
500 DO i=1,sNx+1
501 C store in tmp1k Kuz_Redi
502 tmp1k(i,j) = ( GM_isopycK
503 #ifdef GM_VISBECK_VARIABLE_K
504 & +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0
505 #endif
506 & )*SlopeX(i,j)*taperFct(i,j)
507 ENDDO
508 ENDDO
509 DO j=1,sNy
510 DO i=1,sNx+1
511 C- Vertical gradients interpolated to U points
512 dTdz = (
513 & +recip_drC(k)*
514 & ( maskC(i-1,j,k,bi,bj)*
515 & (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
516 & +maskC( i ,j,k,bi,bj)*
517 & (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
518 & )
519 & +recip_drC(kp1)*
520 & ( maskC(i-1,j,kp1,bi,bj)*
521 & (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
522 & +maskC( i ,j,kp1,bi,bj)*
523 & (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
524 & ) ) * 0.25 _d 0
525 tmp1k(i,j) = dyG(i,j,bi,bj)*drF(k)
526 & * _hFacW(i,j,k,bi,bj)
527 & * tmp1k(i,j) * dTdz
528 ENDDO
529 ENDDO
530 CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KuzTz', k,1,2,bi,bj,myThid)
531 ENDIF
532 #endif /* ALLOW_DIAGNOSTICS */
533
534 C Gradient of Sigma at V points
535 DO j=1-Oly+1,sNy+Oly-1
536 DO i=1-Olx+1,sNx+Olx-1
537 dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
538 & +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k)
539 & )*_maskS(i,j,k,bi,bj)
540 dSigmaDy(i,j)=sigmaY(i,j,k)
541 & *_maskS(i,j,k,bi,bj)
542 dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k )+sigmaR(i,j, k )
543 & +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1
544 & )*_maskS(i,j,k,bi,bj)
545 ENDDO
546 ENDDO
547
548 #ifdef ALLOW_AUTODIFF_TAMC
549 CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
550 CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
551 CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
552 #endif /* ALLOW_AUTODIFF_TAMC */
553
554 C Calculate slopes for use in tensor, taper and/or clip
555 CALL GMREDI_SLOPE_LIMIT(
556 O SlopeX, SlopeY,
557 O SlopeSqr, taperFct,
558 U dSigmaDr,
559 I dSigmaDx, dSigmaDy,
560 I ldd97_LrhoS,rC(k),k,
561 I bi, bj, myThid )
562
563 cph(
564 #ifdef ALLOW_AUTODIFF_TAMC
565 cph(
566 CADJ STORE taperfct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
567 cph)
568 #endif /* ALLOW_AUTODIFF_TAMC */
569 cph)
570
571 #ifdef GM_NON_UNITY_DIAGONAL
572 c IF ( GM_nonUnitDiag ) THEN
573 DO j=1-Oly+1,sNy+Oly-1
574 DO i=1-Olx+1,sNx+Olx-1
575 Kvy(i,j,k,bi,bj) =
576 & ( GM_isopycK
577 #ifdef GM_VISBECK_VARIABLE_K
578 & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
579 #endif
580 & )
581 & *taperFct(i,j)
582 ENDDO
583 ENDDO
584 #ifdef ALLOW_AUTODIFF_TAMC
585 # ifdef GM_EXCLUDE_CLIPPING
586 CADJ STORE Kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
587 # endif
588 #endif
589 DO j=1-Oly+1,sNy+Oly-1
590 DO i=1-Olx+1,sNx+Olx-1
591 Kvy(i,j,k,bi,bj) = MAX( Kvy(i,j,k,bi,bj), GM_Kmin_horiz )
592 ENDDO
593 ENDDO
594 c ENDIF
595 #endif /* GM_NON_UNITY_DIAGONAL */
596
597 #ifdef GM_EXTRA_DIAGONAL
598
599 #ifdef ALLOW_AUTODIFF_TAMC
600 CADJ STORE SlopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
601 CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
602 #endif /* ALLOW_AUTODIFF_TAMC */
603 IF (GM_ExtraDiag) THEN
604 DO j=1-Oly+1,sNy+Oly-1
605 DO i=1-Olx+1,sNx+Olx-1
606 Kvz(i,j,k,bi,bj) =
607 #ifdef ALLOW_KAPGM_CONTROL
608 & ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)
609 #else
610 & ( GM_isopycK - GM_skewflx*GM_background_K
611 #endif
612 #ifdef GM_VISBECK_VARIABLE_K
613 & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*GM_advect
614 #endif
615 & )*SlopeY(i,j)*taperFct(i,j)
616 ENDDO
617 ENDDO
618 ENDIF
619 #endif /* GM_EXTRA_DIAGONAL */
620
621 #ifdef ALLOW_DIAGNOSTICS
622 IF (doDiagRediFlx) THEN
623 c km1 = MAX(k-1,1)
624 DO j=1,sNy+1
625 DO i=1,sNx
626 C store in tmp1k Kvz_Redi
627 tmp1k(i,j) = ( GM_isopycK
628 #ifdef GM_VISBECK_VARIABLE_K
629 & +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0
630 #endif
631 & )*SlopeY(i,j)*taperFct(i,j)
632 ENDDO
633 ENDDO
634 DO j=1,sNy+1
635 DO i=1,sNx
636 C- Vertical gradients interpolated to U points
637 dTdz = (
638 & +recip_drC(k)*
639 & ( maskC(i,j-1,k,bi,bj)*
640 & (theta(i,j-1,km1,bi,bj)-theta(i,j-1,k,bi,bj))
641 & +maskC(i, j ,k,bi,bj)*
642 & (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))
643 & )
644 & +recip_drC(kp1)*
645 & ( maskC(i,j-1,kp1,bi,bj)*
646 & (theta(i,j-1,k,bi,bj)-theta(i,j-1,kp1,bi,bj))
647 & +maskC(i, j ,kp1,bi,bj)*
648 & (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))
649 & ) ) * 0.25 _d 0
650 tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)
651 & * _hFacS(i,j,k,bi,bj)
652 & * tmp1k(i,j) * dTdz
653 ENDDO
654 ENDDO
655 CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KvzTz', k,1,2,bi,bj,myThid)
656 ENDIF
657 #endif /* ALLOW_DIAGNOSTICS */
658
659 #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
660
661 C-- end 2nd loop on vertical level index k
662 ENDDO
663
664
665 #ifdef GM_BOLUS_ADVEC
666 IF (GM_AdvForm) THEN
667 CALL GMREDI_CALC_PSI_B(
668 I bi, bj, iMin, iMax, jMin, jMax,
669 I sigmaX, sigmaY, sigmaR,
670 I ldd97_LrhoW, ldd97_LrhoS,
671 I myThid )
672 ENDIF
673 #endif
674
675 #ifdef ALLOW_TIMEAVE
676 C-- Time-average
677 IF ( taveFreq.GT.0. ) THEN
678
679 CALL TIMEAVE_CUMULATE( GM_Kwx_T, Kwx, Nr,
680 & deltaTclock, bi, bj, myThid )
681 CALL TIMEAVE_CUMULATE( GM_Kwy_T, Kwy, Nr,
682 & deltaTclock, bi, bj, myThid )
683 CALL TIMEAVE_CUMULATE( GM_Kwz_T, Kwz, Nr,
684 & deltaTclock, bi, bj, myThid )
685 #ifdef GM_VISBECK_VARIABLE_K
686 IF ( GM_Visbeck_alpha.NE.0. ) THEN
687 CALL TIMEAVE_CUMULATE( Visbeck_K_T, VisbeckK, 1,
688 & deltaTclock, bi, bj, myThid )
689 ENDIF
690 #endif
691 #ifdef GM_BOLUS_ADVEC
692 IF ( GM_AdvForm ) THEN
693 CALL TIMEAVE_CUMULATE( GM_PsiXtave, GM_PsiX, Nr,
694 & deltaTclock, bi, bj, myThid )
695 CALL TIMEAVE_CUMULATE( GM_PsiYtave, GM_PsiY, Nr,
696 & deltaTclock, bi, bj, myThid )
697 ENDIF
698 #endif
699 DO k=1,Nr
700 GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock
701 ENDDO
702
703 ENDIF
704 #endif /* ALLOW_TIMEAVE */
705
706 #ifdef ALLOW_DIAGNOSTICS
707 IF ( useDiagnostics ) THEN
708 CALL GMREDI_DIAGNOSTICS_FILL(bi,bj,myThid)
709 ENDIF
710 #endif /* ALLOW_DIAGNOSTICS */
711
712 #endif /* ALLOW_GMREDI */
713
714 RETURN
715 END
716
717
718 SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
719 I bi, bj, iMin, iMax, jMin, jMax,
720 I sigmaX, sigmaY, sigmaR,
721 I myThid )
722 C /==========================================================\
723 C | SUBROUTINE GMREDI_CALC_TENSOR |
724 C | o Calculate tensor elements for GM/Redi tensor. |
725 C |==========================================================|
726 C \==========================================================/
727 IMPLICIT NONE
728
729 C == Global variables ==
730 #include "SIZE.h"
731 #include "EEPARAMS.h"
732 #include "GMREDI.h"
733
734 C == Routine arguments ==
735 C
736 _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
737 _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
738 _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
739 INTEGER bi,bj,iMin,iMax,jMin,jMax
740 INTEGER myThid
741 CEndOfInterface
742
743 INTEGER i, j, k
744
745 #ifdef ALLOW_GMREDI
746
747 DO k=1,Nr
748 DO j=1-Oly+1,sNy+Oly-1
749 DO i=1-Olx+1,sNx+Olx-1
750 Kwx(i,j,k,bi,bj) = 0.0
751 Kwy(i,j,k,bi,bj) = 0.0
752 Kwz(i,j,k,bi,bj) = 0.0
753 ENDDO
754 ENDDO
755 ENDDO
756 #endif /* ALLOW_GMREDI */
757
758 RETURN
759 END
  ViewVC Help
Powered by ViewVC 1.1.22