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