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