/[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.35 - (hide annotations) (download)
Sat Jan 2 23:10:47 2010 UTC (14 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62a
Changes since 1.34: +3 -4 lines
time-ave: use simpler (no level index) cumulative-time counter: GM_timeAve(bi,bj)

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

  ViewVC Help
Powered by ViewVC 1.1.22