/[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.27 - (show annotations) (download)
Thu Jun 21 01:33:01 2007 UTC (16 years, 11 months ago) by jmc
Branch: MAIN
Changes since 1.26: +205 -52 lines
add option (GM_taper_scheme='fm07') for Ferrari & McWilliams 2007 Scheme.

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

  ViewVC Help
Powered by ViewVC 1.1.22