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