/[MITgcm]/MITgcm_contrib/submesoscale/code/gmredi_calc_tensor.F
ViewVC logotype

Annotation of /MITgcm_contrib/submesoscale/code/gmredi_calc_tensor.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.7 - (hide annotations) (download)
Fri Mar 19 19:29:47 2010 UTC (15 years, 7 months ago) by zhc
Branch: MAIN
CVS Tags: HEAD
Changes since 1.6: +3 -3 lines 
*** empty log message ***
1 zhc 1.7 C $Header: /u/gcmpack/MITgcm_contrib/submesoscale/code/gmredi_calc_tensor.F,v 1.6 2010/03/19 19:17:45 zhc Exp $
2 dimitri 1.1 C $Name: $
3    
4     #include "GMREDI_OPTIONS.h"
5     #ifdef ALLOW_KPP
6     # include "KPP_OPTIONS.h"
7     #endif
8    
9     CBOP
10     C !ROUTINE: GMREDI_CALC_TENSOR
11     C !INTERFACE:
12     SUBROUTINE GMREDI_CALC_TENSOR(
13     I iMin, iMax, jMin, jMax,
14     I sigmaX, sigmaY, sigmaR,
15     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     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     #include "GMREDI_TAVE.h"
36     #ifdef ALLOW_KPP
37     # include "KPP.h"
38     #endif
39    
40     #ifdef ALLOW_AUTODIFF_TAMC
41     #include "tamc.h"
42     #include "tamc_keys.h"
43     #endif /* ALLOW_AUTODIFF_TAMC */
44    
45     C !INPUT/OUTPUT PARAMETERS:
46     C == Routine arguments ==
47     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     C
52     INTEGER iMin,iMax,jMin,jMax
53     _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     INTEGER bi, bj
57     _RL myTime
58     INTEGER myIter
59     INTEGER myThid
60     CEOP
61    
62     #ifdef ALLOW_GMREDI
63    
64     C !LOCAL VARIABLES:
65     C == Local variables ==
66 zhc 1.5 INTEGER i,j,k
67 dimitri 1.1 _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
68     _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
69     _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
70     _RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
71     _RL dSigmaDr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
72     _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
73     _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
74 zhc 1.5 _RL Kgm_tmp
75 dimitri 1.1 _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    
80     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 zhc 1.5 c#if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
87     INTEGER kp1
88     _RL maskp1
89     c#endif
90 dimitri 1.1
91 dimitri 1.2 #ifdef GM_SUBMESO
92     _RL dBdxAV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
93     _RL dBdyAV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
94     _RL SM_Lf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
95     _RL SM_PsiX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
96     _RL SM_PsiY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
97     _RL SM_PsiXm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
98     _RL SM_PsiYm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
99     _RL hsqmu, hml, recip_hml, qfac, dS, mlmax
100     #endif
101 zhc 1.5
102 dimitri 1.2
103 zhc 1.5 c#ifdef GM_VISBECK_VARIABLE_K
104 dimitri 1.1 #ifdef OLD_VISBECK_CALC
105     _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
106     #else
107 zhc 1.5 _RL dSigmaH, dSigmaR
108     _RL Sloc, M2loc
109 dimitri 1.1 #endif
110 zhc 1.5 _RL recipMaxSlope
111     _RL deltaH, integrDepth
112     _RL N2loc, SNloc
113     c#endif /* GM_VISBECK_VARIABLE_K */
114 dimitri 1.1
115     #ifdef ALLOW_DIAGNOSTICS
116     LOGICAL doDiagRediFlx
117     LOGICAL DIAGNOSTICS_IS_ON
118     EXTERNAL DIAGNOSTICS_IS_ON
119 zhc 1.5 c#if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
120 dimitri 1.1 INTEGER km1
121 dimitri 1.2 _RL dTdz, dTdx, dTdy
122 dimitri 1.1 _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123 zhc 1.5 c#endif
124     #endif /* ALLOW_DIAGNOSTICS */
125 dimitri 1.1
126     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
127    
128     #ifdef ALLOW_AUTODIFF_TAMC
129     act1 = bi - myBxLo(myThid)
130     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
131     act2 = bj - myByLo(myThid)
132     max2 = myByHi(myThid) - myByLo(myThid) + 1
133     act3 = myThid - 1
134     max3 = nTx*nTy
135     act4 = ikey_dynamics - 1
136     igmkey = (act1 + 1) + act2*max1
137     & + act3*max1*max2
138     & + act4*max1*max2*max3
139     #endif /* ALLOW_AUTODIFF_TAMC */
140    
141     #ifdef ALLOW_DIAGNOSTICS
142     doDiagRediFlx = .FALSE.
143     IF ( useDiagnostics ) THEN
144     doDiagRediFlx = DIAGNOSTICS_IS_ON('GM_KuzTz', myThid )
145     doDiagRediFlx = doDiagRediFlx .OR.
146     & DIAGNOSTICS_IS_ON('GM_KvzTz', myThid )
147     ENDIF
148     #endif
149    
150     #ifdef GM_VISBECK_VARIABLE_K
151 zhc 1.5 recipMaxSlope = 0. _d 0
152     IF ( GM_Visbeck_maxSlope.GT.0. _d 0 ) THEN
153     recipMaxSlope = 1. _d 0 / GM_Visbeck_maxSlope
154     ENDIF
155 dimitri 1.1 DO j=1-Oly,sNy+Oly
156     DO i=1-Olx,sNx+Olx
157     VisbeckK(i,j,bi,bj) = 0. _d 0
158     ENDDO
159     ENDDO
160     #endif
161    
162     C-- set ldd97_Lrho (for tapering scheme ldd97):
163     IF ( GM_taper_scheme.EQ.'ldd97' .OR.
164     & GM_taper_scheme.EQ.'fm07' ) THEN
165     Cspd = 2. _d 0
166     LrhoInf = 15. _d 3
167     LrhoSup = 100. _d 3
168     C- Tracer point location (center):
169     DO j=1-Oly,sNy+Oly
170     DO i=1-Olx,sNx+Olx
171     IF (fCori(i,j,bi,bj).NE.0.) THEN
172     ldd97_LrhoC(i,j) = Cspd/ABS(fCori(i,j,bi,bj))
173     ELSE
174     ldd97_LrhoC(i,j) = LrhoSup
175     ENDIF
176     ldd97_LrhoC(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoC(i,j),LrhoSup))
177     ENDDO
178     ENDDO
179     C- U point location (West):
180     DO j=1-Oly,sNy+Oly
181     kLow_W(1-Olx,j) = 0
182     ldd97_LrhoW(1-Olx,j) = LrhoSup
183     DO i=1-Olx+1,sNx+Olx
184     kLow_W(i,j) = MIN(kLowC(i-1,j,bi,bj),kLowC(i,j,bi,bj))
185     fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))
186     IF (fCoriLoc.NE.0.) THEN
187     ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)
188     ELSE
189     ldd97_LrhoW(i,j) = LrhoSup
190     ENDIF
191     ldd97_LrhoW(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoW(i,j),LrhoSup))
192     ENDDO
193     ENDDO
194     C- V point location (South):
195     DO i=1-Olx+1,sNx+Olx
196     kLow_S(i,1-Oly) = 0
197     ldd97_LrhoS(i,1-Oly) = LrhoSup
198     ENDDO
199     DO j=1-Oly+1,sNy+Oly
200     DO i=1-Olx,sNx+Olx
201     kLow_S(i,j) = MIN(kLowC(i,j-1,bi,bj),kLowC(i,j,bi,bj))
202     fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))
203     IF (fCoriLoc.NE.0.) THEN
204     ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)
205     ELSE
206     ldd97_LrhoS(i,j) = LrhoSup
207     ENDIF
208     ldd97_LrhoS(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoS(i,j),LrhoSup))
209     ENDDO
210     ENDDO
211     ELSE
212     C- Just initialize to zero (not use anyway)
213     DO j=1-Oly,sNy+Oly
214     DO i=1-Olx,sNx+Olx
215     ldd97_LrhoC(i,j) = 0. _d 0
216     ldd97_LrhoW(i,j) = 0. _d 0
217     ldd97_LrhoS(i,j) = 0. _d 0
218     ENDDO
219     ENDDO
220     ENDIF
221    
222     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
223     C-- 1rst loop on k : compute Tensor Coeff. at W points.
224    
225     DO j=1-Oly,sNy+Oly
226     DO i=1-Olx,sNx+Olx
227     hTransLay(i,j) = R_low(i,j,bi,bj)
228     baseSlope(i,j) = 0. _d 0
229     recipLambda(i,j) = 0. _d 0
230     locMixLayer(i,j) = 0. _d 0
231     ENDDO
232     ENDDO
233 zhc 1.5 C SM(1)
234     mlmax=0. _d 0
235 dimitri 1.1 #ifdef ALLOW_KPP
236     IF ( useKPP ) THEN
237     DO j=1-Oly,sNy+Oly
238     DO i=1-Olx,sNx+Olx
239     locMixLayer(i,j) = KPPhbl(i,j,bi,bj)
240 zhc 1.5 C SM(1)
241     mlmax=max(mlmax,locMixLayer(i,j))
242 dimitri 1.1 ENDDO
243     ENDDO
244     ELSE
245     #else
246     IF ( .TRUE. ) THEN
247     #endif
248     DO j=1-Oly,sNy+Oly
249     DO i=1-Olx,sNx+Olx
250     locMixLayer(i,j) = hMixLayer(i,j,bi,bj)
251 zhc 1.5 C SM(1)
252     mlmax=max(mlmax,locMixLayer(i,j))
253 dimitri 1.1 ENDDO
254     ENDDO
255     ENDIF
256    
257 dimitri 1.2 #ifdef GM_SUBMESO
258 zhc 1.5 DO j=1-Oly,sNy+Oly
259 dimitri 1.2 DO i=1-Olx,sNx+Olx
260     dBdxAV(i,j) = 0. _d 0
261     dBdyAV(i,j) = 0. _d 0
262     SM_Lf(i,j) = 0. _d 0
263     SM_PsiX(i,j) = 0. _d 0
264     SM_PsiY(i,j) = 0. _d 0
265     SM_PsiXm1(i,j) = 0. _d 0
266 zhc 1.5 SM_PsiYm1(i,j) = 0. _d 0
267 dimitri 1.2 ENDDO
268     ENDDO
269     #endif
270 zhc 1.5
271 dimitri 1.2
272 dimitri 1.1 DO k=Nr,2,-1
273    
274     #ifdef ALLOW_AUTODIFF_TAMC
275     kkey = (igmkey-1)*Nr + k
276     DO j=1-Oly,sNy+Oly
277     DO i=1-Olx,sNx+Olx
278     SlopeX(i,j) = 0. _d 0
279     SlopeY(i,j) = 0. _d 0
280     dSigmaDx(i,j) = 0. _d 0
281     dSigmaDy(i,j) = 0. _d 0
282     dSigmaDr(i,j) = 0. _d 0
283     SlopeSqr(i,j) = 0. _d 0
284     taperFct(i,j) = 0. _d 0
285     Kwx(i,j,k,bi,bj) = 0. _d 0
286     Kwy(i,j,k,bi,bj) = 0. _d 0
287     Kwz(i,j,k,bi,bj) = 0. _d 0
288     # ifdef GM_NON_UNITY_DIAGONAL
289     Kux(i,j,k,bi,bj) = 0. _d 0
290     Kvy(i,j,k,bi,bj) = 0. _d 0
291     # endif
292     # ifdef GM_EXTRA_DIAGONAL
293     Kuz(i,j,k,bi,bj) = 0. _d 0
294     Kvz(i,j,k,bi,bj) = 0. _d 0
295     # endif
296     # ifdef GM_BOLUS_ADVEC
297     GM_PsiX(i,j,k,bi,bj) = 0. _d 0
298     GM_PsiY(i,j,k,bi,bj) = 0. _d 0
299     # endif
300     ENDDO
301     ENDDO
302 zhc 1.5 #endif /* ALLOW_AUTODIFF_TAMC */
303 dimitri 1.1
304     DO j=1-Oly+1,sNy+Oly-1
305     DO i=1-Olx+1,sNx+Olx-1
306     C Gradient of Sigma at rVel points
307     dSigmaDx(i,j)=op25*( sigmaX(i+1,j,k-1)+sigmaX(i,j,k-1)
308     & +sigmaX(i+1,j, k )+sigmaX(i,j, k )
309     & )*maskC(i,j,k,bi,bj)
310     dSigmaDy(i,j)=op25*( sigmaY(i,j+1,k-1)+sigmaY(i,j,k-1)
311     & +sigmaY(i,j+1, k )+sigmaY(i,j, k )
312     & )*maskC(i,j,k,bi,bj)
313     dSigmaDr(i,j)=sigmaR(i,j,k)
314 dimitri 1.2 #ifdef GM_SUBMESO
315     #ifdef GM_SUBMESO_VARYLf
316     C-- Depth average of SigmaR at W points
317     C compute depth average from surface down to the MixLayer depth
318     IF (-rC(k-1).LT.locMixLayer(i,j) ) THEN
319 zhc 1.5 IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
320     integrDepth = -rC( k )
321 dimitri 1.2 C- in 2 steps to avoid mix of RS & RL type in min fct. arguments
322     integrDepth = MIN( integrDepth, locMixLayer(i,j) )
323     C Distance between level center above and the integration depth
324     deltaH = integrDepth + rC(k-1)
325     C If negative then we are below the integration level
326     C (cannot be the case with 2 conditions on maskC & -rC(k-1))
327     C If positive we limit this to the distance from center above
328     deltaH = MIN( deltaH, drC(k) )
329     C Now we convert deltaH to a non-dimensional fraction
330     deltaH = deltaH/( integrDepth+rC(1) )
331     C-- Store db/dr in SM_Lf for now.
332     SM_Lf(i,j) = SM_Lf(i,j)
333     & -gravity*recip_rhoConst*dSigmaDr(i,j)*deltaH
334     ENDIF
335     ENDIF
336     #endif
337     #endif
338 dimitri 1.1 ENDDO
339     ENDDO
340    
341     #ifdef GM_VISBECK_VARIABLE_K
342     #ifndef OLD_VISBECK_CALC
343     IF ( GM_Visbeck_alpha.GT.0. .AND.
344     & -rC(k-1).LT.GM_Visbeck_depth ) THEN
345    
346 zhc 1.5 DO j=1-Oly,sNy+Oly
347     DO i=1-Olx,sNx+Olx
348     dSigmaDr(i,j) = MIN( sigmaR(i,j,k), 0. _d 0 )
349     ENDDO
350     ENDDO
351    
352 dimitri 1.1 C-- Depth average of f/sqrt(Ri) = M^2/N^2 * N
353     C M^2 and N^2 are horizontal & vertical gradient of buoyancy.
354    
355     C Calculate terms for mean Richardson number which is used
356     C in the "variable K" parameterisaton:
357     C compute depth average from surface down to the bottom or
358     C GM_Visbeck_depth, whatever is the shallower.
359    
360     DO j=1-Oly+1,sNy+Oly-1
361     DO i=1-Olx+1,sNx+Olx-1
362     IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
363     integrDepth = -rC( kLowC(i,j,bi,bj) )
364     C- in 2 steps to avoid mix of RS & RL type in min fct. arguments
365     integrDepth = MIN( integrDepth, GM_Visbeck_depth )
366 zhc 1.5 C- to recover "old-visbeck" form with Visbeck_minDepth = Visbeck_depth
367     integrDepth = MAX( integrDepth, GM_Visbeck_minDepth )
368 dimitri 1.1 C Distance between level center above and the integration depth
369     deltaH = integrDepth + rC(k-1)
370     C If negative then we are below the integration level
371     C (cannot be the case with 2 conditions on maskC & -rC(k-1))
372     C If positive we limit this to the distance from center above
373     deltaH = MIN( deltaH, drC(k) )
374     C Now we convert deltaH to a non-dimensional fraction
375     deltaH = deltaH/( integrDepth+rC(1) )
376    
377 zhc 1.5 C-- compute: ( M^2 * S )^1/2 (= S*N since S=M^2/N^2 )
378     C a 5 points average gives a more "homogeneous" formulation
379     C (same stencil and same weights as for dSigmaH calculation)
380     dSigmaR = ( dSigmaDr(i,j)*4. _d 0
381     & + dSigmaDr(i-1,j)
382     & + dSigmaDr(i+1,j)
383     & + dSigmaDr(i,j-1)
384     & + dSigmaDr(i,j+1)
385     & )/( 4. _d 0
386     & + maskC(i-1,j,k,bi,bj)
387     & + maskC(i+1,j,k,bi,bj)
388     & + maskC(i,j-1,k,bi,bj)
389     & + maskC(i,j+1,k,bi,bj)
390     & )
391 dimitri 1.1 dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)
392     & + dSigmaDy(i,j)*dSigmaDy(i,j)
393     IF ( dSigmaH .GT. 0. _d 0 ) THEN
394     dSigmaH = SQRT( dSigmaH )
395 zhc 1.5 C- compute slope, limited by GM_Visbeck_maxSlope:
396     IF ( -dSigmaR.GT.dSigmaH*recipMaxSlope ) THEN
397     Sloc = dSigmaH / ( -dSigmaR )
398 dimitri 1.1 ELSE
399 zhc 1.5 Sloc = GM_Visbeck_maxSlope
400     ENDIF
401     M2loc = gravity*recip_rhoConst*dSigmaH
402     c SNloc = SQRT( Sloc*M2loc )
403     N2loc = -gravity*recip_rhoConst*dSigmaR
404     c N2loc = -gravity*recip_rhoConst*dSigmaDr(i,j)
405     IF ( N2loc.GT.0. _d 0 ) THEN
406     SNloc = Sloc*SQRT(N2loc)
407     ELSE
408     SNloc = 0. _d 0
409 dimitri 1.1 ENDIF
410     ELSE
411     SNloc = 0. _d 0
412     ENDIF
413     VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
414     & +deltaH*GM_Visbeck_alpha
415     & *GM_Visbeck_length*GM_Visbeck_length*SNloc
416     ENDIF
417     ENDDO
418     ENDDO
419     ENDIF
420     #endif /* ndef OLD_VISBECK_CALC */
421     #endif /* GM_VISBECK_VARIABLE_K */
422 zhc 1.5 DO j=1-Oly,sNy+Oly
423     DO i=1-Olx,sNx+Olx
424     dSigmaDr(i,j)=sigmaR(i,j,k)
425     ENDDO
426     ENDDO
427    
428     #ifdef ALLOW_AUTODIFF_TAMC
429     CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
430     CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
431     CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
432     CADJ STORE baseSlope(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
433     CADJ STORE hTransLay(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
434     CADJ STORE recipLambda(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
435     #endif /* ALLOW_AUTODIFF_TAMC */
436 dimitri 1.1
437     C Calculate slopes for use in tensor, taper and/or clip
438     CALL GMREDI_SLOPE_LIMIT(
439     O SlopeX, SlopeY,
440     O SlopeSqr, taperFct,
441     U hTransLay, baseSlope, recipLambda,
442     U dSigmaDr,
443     I dSigmaDx, dSigmaDy,
444     I ldd97_LrhoC, locMixLayer, rF,
445     I kLowC(1-Olx,1-Oly,bi,bj),
446     I k, bi, bj, myTime, myIter, myThid )
447    
448     DO j=1-Oly+1,sNy+Oly-1
449     DO i=1-Olx+1,sNx+Olx-1
450     C Mask Iso-neutral slopes
451     SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)
452     SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)
453     SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)
454     ENDDO
455     ENDDO
456    
457     #ifdef ALLOW_AUTODIFF_TAMC
458     CADJ STORE SlopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
459     CADJ STORE SlopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
460     CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
461     CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
462     CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
463     #endif /* ALLOW_AUTODIFF_TAMC */
464    
465     C Components of Redi/GM tensor
466     DO j=1-Oly+1,sNy+Oly-1
467     DO i=1-Olx+1,sNx+Olx-1
468     Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
469     Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
470     Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
471     ENDDO
472     ENDDO
473    
474     #ifdef GM_VISBECK_VARIABLE_K
475     #ifdef OLD_VISBECK_CALC
476     DO j=1-Oly+1,sNy+Oly-1
477     DO i=1-Olx+1,sNx+Olx-1
478    
479     C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K
480     C but do not know if *taperFct (or **2 ?) is necessary
481     Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)
482    
483     C-- Depth average of M^2/N^2 * N
484    
485     C Calculate terms for mean Richardson number
486     C which is used in the "variable K" parameterisaton.
487     C Distance between interface above layer and the integration depth
488     deltaH=abs(GM_Visbeck_depth)-abs(rF(k))
489     C If positive we limit this to the layer thickness
490 zhc 1.5 integrDepth = drF(k)
491     deltaH=min(deltaH,integrDepth)
492 dimitri 1.1 C If negative then we are below the integration level
493 zhc 1.5 deltaH=max(deltaH, 0. _d 0)
494 dimitri 1.1 C Now we convert deltaH to a non-dimensional fraction
495     deltaH=deltaH/GM_Visbeck_depth
496    
497     IF ( Ssq(i,j).NE.0. .AND. dSigmaDr(i,j).NE.0. ) THEN
498 zhc 1.5 N2loc = -gravity*recip_rhoConst*dSigmaDr(i,j)
499     SNloc = SQRT(Ssq(i,j)*N2loc )
500     VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
501     & +deltaH*GM_Visbeck_alpha
502     & *GM_Visbeck_length*GM_Visbeck_length*SNloc
503 dimitri 1.1 ENDIF
504    
505     ENDDO
506     ENDDO
507     #endif /* OLD_VISBECK_CALC */
508     #endif /* GM_VISBECK_VARIABLE_K */
509    
510     C-- end 1rst loop on vertical level index k
511     ENDDO
512    
513 zhc 1.5
514 dimitri 1.2 #ifdef GM_SUBMESO
515     CBFK-- Use the dsigmadr average to construct the coefficients of the SM param
516     DO j=1-Oly+1,sNy+Oly-1
517     DO i=1-Olx+1,sNx+Olx-1
518     #ifdef GM_SUBMESO_VARYLf
519    
520     IF (SM_Lf(i,j).gt.0) THEN
521     CBFK ML def. rad. as Lf if available and not too small
522     SM_Lf(i,j)=max(sqrt(SM_Lf(i,j))*locMixLayer(i,j)
523     & /abs(fCori(i,j,bi,bj))
524     & ,GM_SM_Lf)
525     ELSE
526     #else
527     IF (.TRUE.) THEN
528     #endif
529     CBFK Otherwise, store just the fixed number
530     SM_Lf(i,j)=GM_SM_Lf
531     ENDIF
532     CBFK Now do the rest of the coefficient
533     dS=2*dxC(i,j,bi,bj)*dyC(i,j,bi,bj)/
534     & (dxC(i,j,bi,bj)+dyC(i,j,bi,bj))
535     CBFK Scaling only works up to 1 degree.
536     dS=min(dS,GM_SM_Lmax)
537     deltaH=sqrt(fCori(i,j,bi,bj)**2+1 _d 0/(GM_SM_tau**2))
538     SM_Lf(i,j)=GM_SM_Ce*dS/(deltaH*SM_Lf(i,j))
539     ENDDO
540     ENDDO
541     #endif
542 dimitri 1.1
543 zhc 1.5
544 dimitri 1.1 #ifdef GM_VISBECK_VARIABLE_K
545     #ifdef ALLOW_AUTODIFF_TAMC
546     CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
547     #endif
548     IF ( GM_Visbeck_alpha.GT.0. ) THEN
549     C- Limit range that KapGM can take
550     DO j=1-Oly+1,sNy+Oly-1
551     DO i=1-Olx+1,sNx+Olx-1
552     VisbeckK(i,j,bi,bj)=
553 zhc 1.5 & MIN( MAX( VisbeckK(i,j,bi,bj), GM_Visbeck_minVal_K ),
554     & GM_Visbeck_maxVal_K )
555 dimitri 1.1 ENDDO
556     ENDDO
557     ENDIF
558     cph( NEW
559     #ifdef ALLOW_AUTODIFF_TAMC
560     CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
561     #endif
562     cph)
563     #endif /* GM_VISBECK_VARIABLE_K */
564    
565     C- express the Tensor in term of Diffusivity (= m**2 / s )
566     DO k=1,Nr
567     #ifdef ALLOW_AUTODIFF_TAMC
568     kkey = (igmkey-1)*Nr + k
569     # if (defined (GM_NON_UNITY_DIAGONAL) || \
570     defined (GM_VISBECK_VARIABLE_K))
571     CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
572     CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
573     CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
574     # endif
575     #endif
576     DO j=1-Oly+1,sNy+Oly-1
577     DO i=1-Olx+1,sNx+Olx-1
578     #ifdef ALLOW_KAPREDI_CONTROL
579     Kgm_tmp = kapredi(i,j,k,bi,bj)
580     #else
581     Kgm_tmp = GM_isopycK
582     #endif
583 zhc 1.5 #ifdef ALLOW_KAPGM_CONTROL
584 dimitri 1.1 & + GM_skewflx*kapgm(i,j,k,bi,bj)
585     #else
586     & + GM_skewflx*GM_background_K
587     #endif
588     #ifdef GM_VISBECK_VARIABLE_K
589     & + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)
590     #endif
591     Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)
592     Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)
593     #ifdef ALLOW_KAPREDI_CONTROL
594 zhc 1.5 Kwz(i,j,k,bi,bj)= ( kapredi(i,j,k,bi,bj)
595 dimitri 1.1 #else
596     Kwz(i,j,k,bi,bj)= ( GM_isopycK
597     #endif
598     #ifdef GM_VISBECK_VARIABLE_K
599     & + VisbeckK(i,j,bi,bj)
600     #endif
601     & )*Kwz(i,j,k,bi,bj)
602     ENDDO
603     ENDDO
604     ENDDO
605    
606     #ifdef ALLOW_DIAGNOSTICS
607     IF ( useDiagnostics .AND. GM_taper_scheme.EQ.'fm07' ) THEN
608     CALL DIAGNOSTICS_FILL( hTransLay, 'GM_hTrsL', 0,1,2,bi,bj,myThid)
609     CALL DIAGNOSTICS_FILL( baseSlope, 'GM_baseS', 0,1,2,bi,bj,myThid)
610     CALL DIAGNOSTICS_FILL(recipLambda,'GM_rLamb', 0,1,2,bi,bj,myThid)
611     ENDIF
612     #endif /* ALLOW_DIAGNOSTICS */
613    
614 zhc 1.5
615 dimitri 1.1 #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
616     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
617     C-- 2nd k loop : compute Tensor Coeff. at U point
618    
619     #ifdef ALLOW_KPP
620     IF ( useKPP ) THEN
621     DO j=1-Oly,sNy+Oly
622     DO i=2-Olx,sNx+Olx
623     locMixLayer(i,j) = ( KPPhbl(i-1,j,bi,bj)
624     & + KPPhbl( i ,j,bi,bj) )*op5
625     ENDDO
626     ENDDO
627     ELSE
628     #else
629     IF ( .TRUE. ) THEN
630     #endif
631     DO j=1-Oly,sNy+Oly
632     DO i=2-Olx,sNx+Olx
633     locMixLayer(i,j) = ( hMixLayer(i-1,j,bi,bj)
634     & + hMixLayer( i ,j,bi,bj) )*op5
635     ENDDO
636     ENDDO
637     ENDIF
638     DO j=1-Oly,sNy+Oly
639     DO i=1-Olx,sNx+Olx
640     hTransLay(i,j) = 0.
641     baseSlope(i,j) = 0.
642     recipLambda(i,j)= 0.
643     ENDDO
644     DO i=2-Olx,sNx+Olx
645     hTransLay(i,j) = MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
646     ENDDO
647     ENDDO
648    
649     DO k=Nr,1,-1
650     kp1 = MIN(Nr,k+1)
651     maskp1 = 1. _d 0
652     IF (k.GE.Nr) maskp1 = 0. _d 0
653     #ifdef ALLOW_AUTODIFF_TAMC
654     kkey = (igmkey-1)*Nr + k
655     #endif
656    
657     C Gradient of Sigma at U points
658     DO j=1-Oly+1,sNy+Oly-1
659     DO i=1-Olx+1,sNx+Olx-1
660     dSigmaDx(i,j)=sigmaX(i,j,k)
661     & *_maskW(i,j,k,bi,bj)
662     dSigmaDy(i,j)=op25*( sigmaY(i-1,j+1,k)+sigmaY(i,j+1,k)
663     & +sigmaY(i-1, j ,k)+sigmaY(i, j ,k)
664     & )*_maskW(i,j,k,bi,bj)
665     dSigmaDr(i,j)=op25*( sigmaR(i-1,j, k )+sigmaR(i,j, k )
666     & +(sigmaR(i-1,j,kp1)+sigmaR(i,j,kp1))*maskp1
667     & )*_maskW(i,j,k,bi,bj)
668 dimitri 1.2
669     #ifdef GM_SUBMESO
670     C-- Depth average of SigmaX at U points
671     C compute depth average from surface down to the MixLayer depth
672     IF (k.GT.1) THEN
673     IF (-rC(k-1).LT.locMixLayer(i,j) ) THEN
674     IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
675     integrDepth = -rC( k )
676     C- in 2 steps to avoid mix of RS & RL type in min fct. arguments
677     integrDepth = MIN( integrDepth, locMixLayer(i,j) )
678     C Distance between level center above and the integration depth
679     deltaH = integrDepth + rC(k-1)
680     C If negative then we are below the integration level
681     C (cannot be the case with 2 conditions on maskC & -rC(k-1))
682     C If positive we limit this to the distance from center above
683     deltaH = MIN( deltaH, drC(k) )
684     C Now we convert deltaH to a non-dimensional fraction
685     deltaH = deltaH/( integrDepth+rC(1) )
686     C-- compute: ( M^2 * S )^1/2 (= M^2 / N since S=M^2/N^2 )
687     dBdxAV(i,j) = dBdxAV(i,j)
688     & +dSigmaDx(i,j)*deltaH*recip_rhoConst*gravity
689     ENDIF
690     ENDIF
691     ENDIF
692     #endif
693 dimitri 1.1 ENDDO
694     ENDDO
695    
696     #ifdef ALLOW_AUTODIFF_TAMC
697     CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
698     CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
699     CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
700     CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
701     CADJ STORE locMixLayer(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
702     CADJ STORE baseSlope(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
703     CADJ STORE hTransLay(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
704     CADJ STORE recipLambda(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
705     #endif /* ALLOW_AUTODIFF_TAMC */
706    
707     C Calculate slopes for use in tensor, taper and/or clip
708     CALL GMREDI_SLOPE_LIMIT(
709     O SlopeX, SlopeY,
710     O SlopeSqr, taperFct,
711     U hTransLay, baseSlope, recipLambda,
712     U dSigmaDr,
713     I dSigmaDx, dSigmaDy,
714     I ldd97_LrhoW, locMixLayer, rC,
715     I kLow_W,
716     I k, bi, bj, myTime, myIter, myThid )
717    
718     #ifdef ALLOW_AUTODIFF_TAMC
719     CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
720     CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
721     #endif /* ALLOW_AUTODIFF_TAMC */
722    
723     #ifdef GM_NON_UNITY_DIAGONAL
724     c IF ( GM_nonUnitDiag ) THEN
725     DO j=1-Oly+1,sNy+Oly-1
726     DO i=1-Olx+1,sNx+Olx-1
727     Kux(i,j,k,bi,bj) =
728     #ifdef ALLOW_KAPREDI_CONTROL
729     & ( kapredi(i,j,k,bi,bj)
730     #else
731     & ( GM_isopycK
732     #endif
733     #ifdef GM_VISBECK_VARIABLE_K
734     & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
735     #endif
736     & )*taperFct(i,j)
737     ENDDO
738     ENDDO
739     #ifdef ALLOW_AUTODIFF_TAMC
740     # ifdef GM_EXCLUDE_CLIPPING
741     CADJ STORE Kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
742     # endif
743     #endif
744     DO j=1-Oly+1,sNy+Oly-1
745     DO i=1-Olx+1,sNx+Olx-1
746     Kux(i,j,k,bi,bj) = MAX( Kux(i,j,k,bi,bj), GM_Kmin_horiz )
747     ENDDO
748     ENDDO
749     c ENDIF
750     #endif /* GM_NON_UNITY_DIAGONAL */
751    
752     #ifdef GM_EXTRA_DIAGONAL
753    
754     #ifdef ALLOW_AUTODIFF_TAMC
755     CADJ STORE SlopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
756     CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
757     #endif /* ALLOW_AUTODIFF_TAMC */
758     IF ( GM_ExtraDiag ) THEN
759     DO j=1-Oly+1,sNy+Oly-1
760     DO i=1-Olx+1,sNx+Olx-1
761     Kuz(i,j,k,bi,bj) =
762     #ifdef ALLOW_KAPREDI_CONTROL
763     & ( kapredi(i,j,k,bi,bj)
764     #else
765     & ( GM_isopycK
766     #endif
767 zhc 1.5 #ifdef ALLOW_KAPGM_CONTROL
768 dimitri 1.1 & - GM_skewflx*kapgm(i,j,k,bi,bj)
769     #else
770     & - GM_skewflx*GM_background_K
771     #endif
772     #ifdef GM_VISBECK_VARIABLE_K
773     & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*GM_advect
774     #endif
775     & )*SlopeX(i,j)*taperFct(i,j)
776     ENDDO
777     ENDDO
778     ENDIF
779     #endif /* GM_EXTRA_DIAGONAL */
780    
781     #ifdef ALLOW_DIAGNOSTICS
782     IF (doDiagRediFlx) THEN
783     km1 = MAX(k-1,1)
784     DO j=1,sNy
785     DO i=1,sNx+1
786     C store in tmp1k Kuz_Redi
787     #ifdef ALLOW_KAPREDI_CONTROL
788     tmp1k(i,j) = ( kapredi(i,j,k,bi,bj)
789     #else
790     tmp1k(i,j) = ( GM_isopycK
791     #endif
792     #ifdef GM_VISBECK_VARIABLE_K
793     & +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0
794     #endif
795     & )*SlopeX(i,j)*taperFct(i,j)
796     ENDDO
797     ENDDO
798     DO j=1,sNy
799     DO i=1,sNx+1
800     C- Vertical gradients interpolated to U points
801     dTdz = (
802     & +recip_drC(k)*
803     & ( maskC(i-1,j,k,bi,bj)*
804     & (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
805     & +maskC( i ,j,k,bi,bj)*
806     & (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
807     & )
808     & +recip_drC(kp1)*
809     & ( maskC(i-1,j,kp1,bi,bj)*
810     & (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
811     & +maskC( i ,j,kp1,bi,bj)*
812     & (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
813     & ) ) * 0.25 _d 0
814     tmp1k(i,j) = dyG(i,j,bi,bj)*drF(k)
815     & * _hFacW(i,j,k,bi,bj)
816     & * tmp1k(i,j) * dTdz
817     ENDDO
818     ENDDO
819     CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KuzTz', k,1,2,bi,bj,myThid)
820     ENDIF
821     #endif /* ALLOW_DIAGNOSTICS */
822    
823     C-- end 2nd loop on vertical level index k
824     ENDDO
825    
826     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
827     C-- 3rd k loop : compute Tensor Coeff. at V point
828 dimitri 1.3
829 dimitri 1.1 #ifdef ALLOW_KPP
830     IF ( useKPP ) THEN
831     DO j=2-Oly,sNy+Oly
832     DO i=1-Olx,sNx+Olx
833     locMixLayer(i,j) = ( KPPhbl(i,j-1,bi,bj)
834     & + KPPhbl(i, j ,bi,bj) )*op5
835     ENDDO
836     ENDDO
837     ELSE
838     #else
839     IF ( .TRUE. ) THEN
840     #endif
841     DO j=2-Oly,sNy+Oly
842     DO i=1-Olx,sNx+Olx
843     locMixLayer(i,j) = ( hMixLayer(i,j-1,bi,bj)
844     & + hMixLayer(i, j ,bi,bj) )*op5
845     ENDDO
846     ENDDO
847     ENDIF
848     DO j=1-Oly,sNy+Oly
849     DO i=1-Olx,sNx+Olx
850     hTransLay(i,j) = 0.
851     baseSlope(i,j) = 0.
852     recipLambda(i,j)= 0.
853     ENDDO
854     ENDDO
855     DO j=2-Oly,sNy+Oly
856     DO i=1-Olx,sNx+Olx
857     hTransLay(i,j) = MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
858     ENDDO
859     ENDDO
860    
861     C Gradient of Sigma at V points
862     DO k=Nr,1,-1
863     kp1 = MIN(Nr,k+1)
864     maskp1 = 1. _d 0
865     IF (k.GE.Nr) maskp1 = 0. _d 0
866     #ifdef ALLOW_AUTODIFF_TAMC
867     kkey = (igmkey-1)*Nr + k
868     #endif
869    
870     DO j=1-Oly+1,sNy+Oly-1
871     DO i=1-Olx+1,sNx+Olx-1
872     dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
873     & +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k)
874     & )*_maskS(i,j,k,bi,bj)
875     dSigmaDy(i,j)=sigmaY(i,j,k)
876     & *_maskS(i,j,k,bi,bj)
877     dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k )+sigmaR(i,j, k )
878     & +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1
879     & )*_maskS(i,j,k,bi,bj)
880 dimitri 1.2
881     #ifdef GM_SUBMESO
882     C-- Depth average of SigmaY at V points
883     C compute depth average from surface down to the MixLayer depth
884     IF (k.GT.1) THEN
885     IF (-rC(k-1).LT.locMixLayer(i,j) ) THEN
886     IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
887     integrDepth = -rC( k )
888     C- in 2 steps to avoid mix of RS & RL type in min fct. arguments
889     integrDepth = MIN( integrDepth, locMixLayer(i,j) )
890     C Distance between level center above and the integration depth
891     deltaH = integrDepth + rC(k-1)
892     C If negative then we are below the integration level
893     C (cannot be the case with 2 conditions on maskC & -rC(k-1))
894     C If positive we limit this to the distance from center above
895     deltaH = MIN( deltaH, drC(k) )
896     C Now we convert deltaH to a non-dimensional fraction
897     deltaH = deltaH/( integrDepth+rC(1) )
898     dBdyAV(i,j) = dBdyAV(i,j)
899     & +dSigmaDy(i,j)*deltaH*recip_rhoConst*gravity
900     ENDIF
901     ENDIF
902     ENDIF
903     #endif
904 dimitri 1.1 ENDDO
905     ENDDO
906    
907     #ifdef ALLOW_AUTODIFF_TAMC
908     CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
909     CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
910     CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
911     CADJ STORE baseSlope(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
912     CADJ STORE hTransLay(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
913     CADJ STORE recipLambda(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
914     #endif /* ALLOW_AUTODIFF_TAMC */
915    
916     C Calculate slopes for use in tensor, taper and/or clip
917     CALL GMREDI_SLOPE_LIMIT(
918     O SlopeX, SlopeY,
919     O SlopeSqr, taperFct,
920     U hTransLay, baseSlope, recipLambda,
921     U dSigmaDr,
922     I dSigmaDx, dSigmaDy,
923     I ldd97_LrhoS, locMixLayer, rC,
924     I kLow_S,
925     I k, bi, bj, myTime, myIter, myThid )
926    
927     cph(
928     #ifdef ALLOW_AUTODIFF_TAMC
929     cph(
930     CADJ STORE taperfct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
931     cph)
932     #endif /* ALLOW_AUTODIFF_TAMC */
933     cph)
934    
935     #ifdef GM_NON_UNITY_DIAGONAL
936     c IF ( GM_nonUnitDiag ) THEN
937     DO j=1-Oly+1,sNy+Oly-1
938     DO i=1-Olx+1,sNx+Olx-1
939     Kvy(i,j,k,bi,bj) =
940     #ifdef ALLOW_KAPREDI_CONTROL
941     & ( kapredi(i,j,k,bi,bj)
942     #else
943     & ( GM_isopycK
944     #endif
945     #ifdef GM_VISBECK_VARIABLE_K
946     & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
947     #endif
948     & )*taperFct(i,j)
949     ENDDO
950     ENDDO
951     #ifdef ALLOW_AUTODIFF_TAMC
952     # ifdef GM_EXCLUDE_CLIPPING
953     CADJ STORE Kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
954     # endif
955     #endif
956     DO j=1-Oly+1,sNy+Oly-1
957     DO i=1-Olx+1,sNx+Olx-1
958     Kvy(i,j,k,bi,bj) = MAX( Kvy(i,j,k,bi,bj), GM_Kmin_horiz )
959     ENDDO
960     ENDDO
961     c ENDIF
962     #endif /* GM_NON_UNITY_DIAGONAL */
963    
964     #ifdef GM_EXTRA_DIAGONAL
965    
966     #ifdef ALLOW_AUTODIFF_TAMC
967     CADJ STORE SlopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
968     CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
969     #endif /* ALLOW_AUTODIFF_TAMC */
970     IF ( GM_ExtraDiag ) THEN
971     DO j=1-Oly+1,sNy+Oly-1
972     DO i=1-Olx+1,sNx+Olx-1
973     Kvz(i,j,k,bi,bj) =
974     #ifdef ALLOW_KAPREDI_CONTROL
975     & ( kapredi(i,j,k,bi,bj)
976     #else
977 zhc 1.5 & ( GM_isopycK
978 dimitri 1.1 #endif
979 zhc 1.5 #ifdef ALLOW_KAPGM_CONTROL
980 dimitri 1.1 & - GM_skewflx*kapgm(i,j,k,bi,bj)
981     #else
982     & - GM_skewflx*GM_background_K
983     #endif
984     #ifdef GM_VISBECK_VARIABLE_K
985     & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*GM_advect
986     #endif
987     & )*SlopeY(i,j)*taperFct(i,j)
988     ENDDO
989     ENDDO
990     ENDIF
991     #endif /* GM_EXTRA_DIAGONAL */
992    
993     #ifdef ALLOW_DIAGNOSTICS
994     IF (doDiagRediFlx) THEN
995     km1 = MAX(k-1,1)
996     DO j=1,sNy+1
997     DO i=1,sNx
998     C store in tmp1k Kvz_Redi
999     #ifdef ALLOW_KAPREDI_CONTROL
1000     tmp1k(i,j) = ( kapredi(i,j,k,bi,bj)
1001     #else
1002     tmp1k(i,j) = ( GM_isopycK
1003     #endif
1004     #ifdef GM_VISBECK_VARIABLE_K
1005     & +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0
1006     #endif
1007     & )*SlopeY(i,j)*taperFct(i,j)
1008     ENDDO
1009     ENDDO
1010     DO j=1,sNy+1
1011     DO i=1,sNx
1012 zhc 1.5 C- Vertical gradients interpolated to U points
1013     dTdz = (
1014     & +recip_drC(k)*
1015 dimitri 1.1 & ( maskC(i,j-1,k,bi,bj)*
1016     & (theta(i,j-1,km1,bi,bj)-theta(i,j-1,k,bi,bj))
1017     & +maskC(i, j ,k,bi,bj)*
1018     & (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))
1019     & )
1020 zhc 1.5 & +recip_drC(kp1)*
1021 dimitri 1.1 & ( maskC(i,j-1,kp1,bi,bj)*
1022     & (theta(i,j-1,k,bi,bj)-theta(i,j-1,kp1,bi,bj))
1023     & +maskC(i, j ,kp1,bi,bj)*
1024     & (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))
1025 zhc 1.5 & ) ) * 0.25 _d 0
1026 dimitri 1.1 tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)
1027     & * _hFacS(i,j,k,bi,bj)
1028     & * tmp1k(i,j) * dTdz
1029     ENDDO
1030     ENDDO
1031     CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KvzTz', k,1,2,bi,bj,myThid)
1032     ENDIF
1033     #endif /* ALLOW_DIAGNOSTICS */
1034    
1035     C-- end 3rd loop on vertical level index k
1036     ENDDO
1037    
1038     #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
1039    
1040    
1041     #ifdef GM_BOLUS_ADVEC
1042     IF (GM_AdvForm) THEN
1043     CALL GMREDI_CALC_PSI_B(
1044     I bi, bj, iMin, iMax, jMin, jMax,
1045     I sigmaX, sigmaY, sigmaR,
1046     I ldd97_LrhoW, ldd97_LrhoS,
1047     I myThid )
1048     ENDIF
1049     #endif
1050    
1051 dimitri 1.2 #ifdef GM_SUBMESO
1052     CBFK Add the submesoscale contribution, in a 4th k loop
1053     DO k=1,Nr
1054     km1=max(1,k-1)
1055     IF ((k.gt.1).and.(-rF(k-1) .lt. mlmax)) THEN
1056     kp1 = MIN(k+1,Nr)
1057     CBFK Add in the mu vertical structure factor
1058     DO j=1-Oly+1,sNy+Oly-1
1059     DO i=1-Olx+1,sNx+Olx-1
1060     hml=hMixLayer(i,j,bi,bj)
1061     IF (hml.gt.0 _d 0) THEN
1062     recip_hml=1 _d 0/hml
1063     ELSE
1064     recip_hml=0 _d 0
1065     ENDIF
1066     CBFK We calculate the h^2 mu(z) factor only on w points.
1067     CBFK It is possible that we might need to calculate it
1068     CBFK on Psi or u,v points independently to prevent spurious
1069     CBFK entrainment. Unlikely that this will be major
1070     CBFK (it wasnt in offline testing).
1071     qfac=(2*rf(k)*recip_hml+1 _d 0)**2
1072     hsqmu=(1 _d 0-qfac)*(1 _d 0+(5 _d 0)*qfac/21 _d 0)
1073     hsqmu=max(0 _d 0, hsqmu)*hml**2
1074     SM_Lf(i,j)=SM_Lf(i,j)*hsqmu
1075     ENDDO
1076     ENDDO
1077     CBFK Now interpolate to match locations
1078     DO j=1-Oly+1,sNy+Oly-1
1079     DO i=1-Olx+1,sNx+Olx-1
1080     C SM_Lf coefficients are on rVel points
1081     C Psix are on faces above U
1082     SM_PsiX(i,j)=op5*(SM_Lf(i+1,j)+SM_Lf(i,j))*dBdxAV(i,j)
1083     & *_maskW(i,j,k,bi,bj)
1084     C Psiy are on faces above V
1085     SM_PsiY(i,j)=op5*(SM_Lf(i,j+1)+SM_Lf(i,j))*dBdyAV(i,j)
1086     & *_maskS(i,j,k,bi,bj)
1087    
1088 zhc 1.6 c hzhang: clipping here, ref to Baylor paper 3 appendix A MOM
1089     hml=hMixLayer(i,j,bi,bj)
1090     IF (hml .lt. (delR(1)+delR(2)+delR(3)+delR(4))) THEN
1091     SM_PsiX(i,j)=0 _d 0
1092     SM_PsiY(i,j)=0 _d 0
1093     ENDIF
1094    
1095     hml=.5 * delR(k)
1096 zhc 1.7 IF ( abs(SM_PsiX(i,j)) .gt. hml ) THEN
1097 zhc 1.6 SM_PsiX(i,j)=SIGN( hml, SM_PsiX(i,j) )
1098     ENDIF
1099 zhc 1.7 IF ( abs(SM_PsiY(i,j)) .gt. hml ) THEN
1100 zhc 1.6 SM_PsiY(i,j)=SIGN( hml, SM_PsiY(i,j) )
1101     ENDIF
1102     c hzhang: clipping done
1103    
1104    
1105 dimitri 1.2 #ifndef GM_BOLUS_ADVEC
1106     C Kwx,Kwy are on rVel Points
1107     Kwx(i,j,k,bi,bj) = Kwx(i,j,k,bi,bj)
1108     & +op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))
1109     Kwy(i,j,k,bi,bj) = Kwy(i,j,k,bi,bj)
1110     & +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))
1111     #ifdef GM_EXTRA_DIAGONAL
1112     IF (GM_ExtraDiag) THEN
1113     C Kuz,Kvz are on u,v Points
1114     Kuz(i,j,k,bi,bj) = Kuz(i,j,k,bi,bj)
1115     & -op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1116     Kvz(i,j,k,bi,bj) = Kvz(i,j,k,bi,bj)
1117     & -op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1118     ENDIF
1119     #endif
1120     #else
1121     IF (GM_AdvForm) THEN
1122     GM_PsiX(i,j,k,bi,bj)=GM_PsiX(i,j,k,bi,bj)+SM_PsiX(i,j)
1123     GM_PsiY(i,j,k,bi,bj)=GM_PsiY(i,j,k,bi,bj)+SM_PsiY(i,j)
1124     ENDIF
1125     #endif
1126     ENDDO
1127     ENDDO
1128     ELSE
1129     DO j=1-Oly+1,sNy+Oly-1
1130     DO i=1-Olx+1,sNx+Olx-1
1131     SM_PsiX(i,j)=0. _d 0
1132     SM_PsiY(i,j)=0. _d 0
1133     ENDDO
1134     ENDDO
1135     ENDIF
1136    
1137     #ifdef ALLOW_DIAGNOSTICS
1138     IF ( useDiagnostics ) THEN
1139     IF ( DIAGNOSTICS_IS_ON('SM_PsiX ',myThid) ) THEN
1140     CALL DIAGNOSTICS_FILL(SM_PsiX,'SM_PsiX ',k,1,2,bi,bj,myThid)
1141     ENDIF
1142     IF ( DIAGNOSTICS_IS_ON('SM_PsiY ',myThid) ) THEN
1143     CALL DIAGNOSTICS_FILL(SM_PsiY,'SM_PsiY ',k,1,2,bi,bj,myThid)
1144     ENDIF
1145    
1146     CBFK Note: for comparision, you can diagnose the bolus form
1147     CBFK or the Kappa form in the same simulation, regardless of other
1148     CBFK settings
1149     IF ( DIAGNOSTICS_IS_ON('SM_ubT ',myThid) ) THEN
1150     DO j=jMin,jMax
1151     DO i=iMin,iMax
1152     tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiX(i,j)
1153     & -SM_PsiXm1(i,j) )
1154     & *maskW(i,j,km1,bi,bj)
1155     & *op5*(Theta(i,j,km1,bi,bj)+Theta(i-1,j,km1,bi,bj))
1156     ENDDO
1157     ENDDO
1158     CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT ', km1,1,2,bi,bj,myThid)
1159     ENDIF
1160    
1161     IF ( DIAGNOSTICS_IS_ON('SM_vbT ',myThid) ) THEN
1162     DO j=jMin,jMax
1163     DO i=iMin,iMax
1164     tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiY(i,j)
1165     & -SM_PsiYm1(i,j) )
1166     & *maskS(i,j,km1,bi,bj)
1167     & *op5*(Theta(i,j,km1,bi,bj)+Theta(i,j-1,km1,bi,bj))
1168     ENDDO
1169     ENDDO
1170     CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT ', km1,1,2,bi,bj,myThid)
1171     ENDIF
1172    
1173     IF ( DIAGNOSTICS_IS_ON('SM_wbT ',myThid) ) THEN
1174     DO j=jMin,jMax
1175     DO i=iMin,iMax
1176     tmp1k(i,j) =
1177     & (dyG(i+1,j,bi,bj)*SM_PsiX(i+1,j)
1178     & -dyG( i ,j,bi,bj)*SM_PsiX( i ,j)
1179     & +dxG(i,j+1,bi,bj)*SM_PsiY(i,j+1)
1180     & -dxG(i, j ,bi,bj)*SM_PsiY(i, j ))
1181     & *op5*(Theta(i,j,k,bi,bj)+Theta(i,j,km1,bi,bj))
1182     ENDDO
1183     ENDDO
1184     CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT ', k,1,2,bi,bj,myThid)
1185     C print *,'SM_wbT',k,tmp1k
1186     ENDIF
1187    
1188     IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1189     DO j=1,sNy
1190     DO i=1,sNx+1
1191     C- Vertical gradients interpolated to U points
1192     dTdz = (
1193     & +recip_drC(k)*
1194     & ( maskC(i-1,j,k,bi,bj)*
1195     & (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
1196     & +maskC( i ,j,k,bi,bj)*
1197     & (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
1198     & )
1199     & +recip_drC(kp1)*
1200     & ( maskC(i-1,j,kp1,bi,bj)*
1201     & (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
1202     & +maskC( i ,j,kp1,bi,bj)*
1203     & (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
1204     & ) ) * 0.25 _d 0
1205     tmp1k(i,j) = - dyG(i,j,bi,bj)*drF(k)
1206     & * _hFacW(i,j,k,bi,bj)
1207     & *op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1208     & * dTdz
1209     ENDDO
1210     ENDDO
1211     CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KuzTz', k,1,2,bi,bj,myThid)
1212     C print *,'SM_KuzTz',k,tmp1k
1213     ENDIF
1214    
1215     IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1216     DO j=1,sNy+1
1217     DO i=1,sNx
1218     C- Vertical gradients interpolated to V points
1219     dTdz = op5*(
1220     & +op5*recip_drC(k)*
1221     & ( maskC(i,j-1,k,bi,bj)*
1222     & (Theta(i,j-1,km1,bi,bj)-Theta(i,j-1,k,bi,bj))
1223     & +maskC(i, j ,k,bi,bj)*
1224     & (Theta(i, j ,km1,bi,bj)-Theta(i, j ,k,bi,bj))
1225     & )
1226     & +op5*recip_drC(kp1)*
1227     & ( maskC(i,j-1,kp1,bi,bj)*
1228     & (Theta(i,j-1,k,bi,bj)-Theta(i,j-1,kp1,bi,bj))
1229     & +maskC(i, j ,kp1,bi,bj)*
1230     & (Theta(i, j ,k,bi,bj)-Theta(i, j ,kp1,bi,bj))
1231     & ) )
1232     tmp1k(i,j) = - dxG(i,j,bi,bj)*drF(k)
1233     & * _hFacS(i,j,k,bi,bj)
1234     & *op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1235     & * dTdz
1236     ENDDO
1237     ENDDO
1238     CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KvzTz', k,1,2,bi,bj,myThid)
1239     C print *,'SM_KvzTz',k,tmp1k
1240     ENDIF
1241    
1242     IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1243     DO j=jMin,jMax
1244     DO i=iMin,iMax
1245     C- Horizontal gradients interpolated to W points
1246     dTdx = op5*(
1247     & +op5*(_maskW(i+1,j,k,bi,bj)
1248     & *_recip_dxC(i+1,j,bi,bj)*
1249     & (Theta(i+1,j,k,bi,bj)-Theta(i,j,k,bi,bj))
1250     & +_maskW(i,j,k,bi,bj)
1251     & *_recip_dxC(i,j,bi,bj)*
1252     & (Theta(i,j,k,bi,bj)-Theta(i-1,j,k,bi,bj)))
1253     & +op5*(_maskW(i+1,j,k-1,bi,bj)
1254     & *_recip_dxC(i+1,j,bi,bj)*
1255     & (Theta(i+1,j,k-1,bi,bj)-Theta(i,j,k-1,bi,bj))
1256     & +_maskW(i,j,k-1,bi,bj)
1257     & *_recip_dxC(i,j,bi,bj)*
1258     & (Theta(i,j,k-1,bi,bj)-Theta(i-1,j,k-1,bi,bj)))
1259     & )
1260    
1261     dTdy = op5*(
1262     & +op5*(_maskS(i,j,k,bi,bj)
1263     & *_recip_dyC(i,j,bi,bj)*
1264     & (Theta(i,j,k,bi,bj)-Theta(i,j-1,k,bi,bj))
1265     & +_maskS(i,j+1,k,bi,bj)
1266     & *_recip_dyC(i,j+1,bi,bj)*
1267     & (Theta(i,j+1,k,bi,bj)-Theta(i,j,k,bi,bj)))
1268     & +op5*(_maskS(i,j,k-1,bi,bj)
1269     & *_recip_dyC(i,j,bi,bj)*
1270     & (Theta(i,j,k-1,bi,bj)-Theta(i,j-1,k-1,bi,bj))
1271     & +_maskS(i,j+1,k-1,bi,bj)
1272     & *_recip_dyC(i,j+1,bi,bj)*
1273     & (Theta(i,j+1,k-1,bi,bj)-Theta(i,j,k-1,bi,bj)))
1274     & )
1275    
1276     tmp1k(i,j) = - _rA(i,j,bi,bj)
1277     & *(op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))*dTdx
1278     & +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))*dTdy)
1279     ENDDO
1280     ENDDO
1281     CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', k,1,2,bi,bj,myThid)
1282     C print *,'SM_KrddT',k,tmp1k
1283     ENDIF
1284     ENDIF
1285     #endif
1286     DO j=1-Oly+1,sNy+Oly-1
1287     DO i=1-Olx+1,sNx+Olx-1
1288     SM_PsiXm1(i,j)=SM_PsiX(i,j)
1289     SM_PsiYm1(i,j)=SM_PsiY(i,j)
1290     tmp1k(i,j)=0 _d 0
1291     ENDDO
1292     ENDDO
1293     ENDDO
1294    
1295     CBFK Always Zero at the bottom.
1296     IF ( DIAGNOSTICS_IS_ON('SM_ubT ',myThid) ) THEN
1297     CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT ', Nr,1,2,bi,bj,myThid)
1298     ENDIF
1299     IF ( DIAGNOSTICS_IS_ON('SM_vbT ',myThid) ) THEN
1300     CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT ', Nr,1,2,bi,bj,myThid)
1301     ENDIF
1302     IF ( DIAGNOSTICS_IS_ON('SM_wbT ',myThid) ) THEN
1303     CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT ', Nr,1,2,bi,bj,myThid)
1304     ENDIF
1305     IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1306     CALL DIAGNOSTICS_FILL(tmp1k,'SM_KuzTz', Nr,1,2,bi,bj,myThid)
1307     ENDIF
1308     IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1309     CALL DIAGNOSTICS_FILL(tmp1k,'SM_KvzTz', Nr,1,2,bi,bj,myThid)
1310     ENDIF
1311     IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1312     CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', Nr,1,2,bi,bj,myThid)
1313     ENDIF
1314     #endif
1315    
1316 dimitri 1.1 #ifdef ALLOW_TIMEAVE
1317     C-- Time-average
1318     IF ( taveFreq.GT.0. ) THEN
1319    
1320     CALL TIMEAVE_CUMULATE( GM_Kwx_T, Kwx, Nr,
1321     & deltaTclock, bi, bj, myThid )
1322     CALL TIMEAVE_CUMULATE( GM_Kwy_T, Kwy, Nr,
1323     & deltaTclock, bi, bj, myThid )
1324     CALL TIMEAVE_CUMULATE( GM_Kwz_T, Kwz, Nr,
1325     & deltaTclock, bi, bj, myThid )
1326     #ifdef GM_VISBECK_VARIABLE_K
1327     IF ( GM_Visbeck_alpha.NE.0. ) THEN
1328     CALL TIMEAVE_CUMULATE( Visbeck_K_T, VisbeckK, 1,
1329     & deltaTclock, bi, bj, myThid )
1330     ENDIF
1331     #endif
1332     #ifdef GM_BOLUS_ADVEC
1333     IF ( GM_AdvForm ) THEN
1334     CALL TIMEAVE_CUMULATE( GM_PsiXtave, GM_PsiX, Nr,
1335     & deltaTclock, bi, bj, myThid )
1336     CALL TIMEAVE_CUMULATE( GM_PsiYtave, GM_PsiY, Nr,
1337     & deltaTclock, bi, bj, myThid )
1338     ENDIF
1339     #endif
1340 zhc 1.5 GM_timeAve(bi,bj) = GM_timeAve(bi,bj)+deltaTclock
1341 dimitri 1.1
1342     ENDIF
1343     #endif /* ALLOW_TIMEAVE */
1344    
1345     #ifdef ALLOW_DIAGNOSTICS
1346     IF ( useDiagnostics ) THEN
1347     CALL GMREDI_DIAGNOSTICS_FILL(bi,bj,myThid)
1348     ENDIF
1349     #endif /* ALLOW_DIAGNOSTICS */
1350    
1351     #endif /* ALLOW_GMREDI */
1352    
1353     RETURN
1354     END
1355    
1356 zhc 1.5 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1357 dimitri 1.1
1358     SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
1359     I iMin, iMax, jMin, jMax,
1360     I sigmaX, sigmaY, sigmaR,
1361     I bi, bj, myTime, myIter, myThid )
1362     C /==========================================================\
1363     C | SUBROUTINE GMREDI_CALC_TENSOR |
1364     C | o Calculate tensor elements for GM/Redi tensor. |
1365     C |==========================================================|
1366     C \==========================================================/
1367     IMPLICIT NONE
1368    
1369     C == Global variables ==
1370     #include "SIZE.h"
1371     #include "EEPARAMS.h"
1372     #include "GMREDI.h"
1373    
1374     C == Routine arguments ==
1375     C
1376     _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
1377     _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
1378     _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
1379     INTEGER iMin,iMax,jMin,jMax
1380     INTEGER bi, bj
1381     _RL myTime
1382     INTEGER myIter
1383     INTEGER myThid
1384     CEndOfInterface
1385    
1386     #ifdef ALLOW_GMREDI
1387    
1388     INTEGER i, j, k
1389    
1390     DO k=1,Nr
1391     DO j=1-Oly+1,sNy+Oly-1
1392     DO i=1-Olx+1,sNx+Olx-1
1393     Kwx(i,j,k,bi,bj) = 0.0
1394     Kwy(i,j,k,bi,bj) = 0.0
1395     Kwz(i,j,k,bi,bj) = 0.0
1396     ENDDO
1397     ENDDO
1398     ENDDO
1399     #endif /* ALLOW_GMREDI */
1400    
1401     RETURN
1402     END
  ViewVC Help
Powered by ViewVC 1.1.22