/[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.5 - (hide annotations) (download)
Fri Mar 12 18:31:00 2010 UTC (14 years, 3 months ago) by zhc
Branch: MAIN
Changes since 1.4: +108 -62 lines 
updated with checkpoint62c
1 zhc 1.5 C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_calc_tensor.F,v 1.36 2010/01/20 01:20:29 jmc 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     #ifndef GM_BOLUS_ADVEC
1089     C Kwx,Kwy are on rVel Points
1090     Kwx(i,j,k,bi,bj) = Kwx(i,j,k,bi,bj)
1091     & +op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))
1092     Kwy(i,j,k,bi,bj) = Kwy(i,j,k,bi,bj)
1093     & +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))
1094     #ifdef GM_EXTRA_DIAGONAL
1095     IF (GM_ExtraDiag) THEN
1096     C Kuz,Kvz are on u,v Points
1097     Kuz(i,j,k,bi,bj) = Kuz(i,j,k,bi,bj)
1098     & -op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1099     Kvz(i,j,k,bi,bj) = Kvz(i,j,k,bi,bj)
1100     & -op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1101     ENDIF
1102     #endif
1103     #else
1104     IF (GM_AdvForm) THEN
1105     GM_PsiX(i,j,k,bi,bj)=GM_PsiX(i,j,k,bi,bj)+SM_PsiX(i,j)
1106     GM_PsiY(i,j,k,bi,bj)=GM_PsiY(i,j,k,bi,bj)+SM_PsiY(i,j)
1107     ENDIF
1108     #endif
1109     ENDDO
1110     ENDDO
1111     ELSE
1112     DO j=1-Oly+1,sNy+Oly-1
1113     DO i=1-Olx+1,sNx+Olx-1
1114     SM_PsiX(i,j)=0. _d 0
1115     SM_PsiY(i,j)=0. _d 0
1116     ENDDO
1117     ENDDO
1118     ENDIF
1119    
1120     #ifdef ALLOW_DIAGNOSTICS
1121     IF ( useDiagnostics ) THEN
1122     IF ( DIAGNOSTICS_IS_ON('SM_PsiX ',myThid) ) THEN
1123     CALL DIAGNOSTICS_FILL(SM_PsiX,'SM_PsiX ',k,1,2,bi,bj,myThid)
1124     ENDIF
1125     IF ( DIAGNOSTICS_IS_ON('SM_PsiY ',myThid) ) THEN
1126     CALL DIAGNOSTICS_FILL(SM_PsiY,'SM_PsiY ',k,1,2,bi,bj,myThid)
1127     ENDIF
1128    
1129     CBFK Note: for comparision, you can diagnose the bolus form
1130     CBFK or the Kappa form in the same simulation, regardless of other
1131     CBFK settings
1132     IF ( DIAGNOSTICS_IS_ON('SM_ubT ',myThid) ) THEN
1133     DO j=jMin,jMax
1134     DO i=iMin,iMax
1135     tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiX(i,j)
1136     & -SM_PsiXm1(i,j) )
1137     & *maskW(i,j,km1,bi,bj)
1138     & *op5*(Theta(i,j,km1,bi,bj)+Theta(i-1,j,km1,bi,bj))
1139     ENDDO
1140     ENDDO
1141     CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT ', km1,1,2,bi,bj,myThid)
1142     ENDIF
1143    
1144     IF ( DIAGNOSTICS_IS_ON('SM_vbT ',myThid) ) THEN
1145     DO j=jMin,jMax
1146     DO i=iMin,iMax
1147     tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiY(i,j)
1148     & -SM_PsiYm1(i,j) )
1149     & *maskS(i,j,km1,bi,bj)
1150     & *op5*(Theta(i,j,km1,bi,bj)+Theta(i,j-1,km1,bi,bj))
1151     ENDDO
1152     ENDDO
1153     CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT ', km1,1,2,bi,bj,myThid)
1154     ENDIF
1155    
1156     IF ( DIAGNOSTICS_IS_ON('SM_wbT ',myThid) ) THEN
1157     DO j=jMin,jMax
1158     DO i=iMin,iMax
1159     tmp1k(i,j) =
1160     & (dyG(i+1,j,bi,bj)*SM_PsiX(i+1,j)
1161     & -dyG( i ,j,bi,bj)*SM_PsiX( i ,j)
1162     & +dxG(i,j+1,bi,bj)*SM_PsiY(i,j+1)
1163     & -dxG(i, j ,bi,bj)*SM_PsiY(i, j ))
1164     & *op5*(Theta(i,j,k,bi,bj)+Theta(i,j,km1,bi,bj))
1165     ENDDO
1166     ENDDO
1167     CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT ', k,1,2,bi,bj,myThid)
1168     C print *,'SM_wbT',k,tmp1k
1169     ENDIF
1170    
1171     IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1172     DO j=1,sNy
1173     DO i=1,sNx+1
1174     C- Vertical gradients interpolated to U points
1175     dTdz = (
1176     & +recip_drC(k)*
1177     & ( maskC(i-1,j,k,bi,bj)*
1178     & (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
1179     & +maskC( i ,j,k,bi,bj)*
1180     & (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
1181     & )
1182     & +recip_drC(kp1)*
1183     & ( maskC(i-1,j,kp1,bi,bj)*
1184     & (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
1185     & +maskC( i ,j,kp1,bi,bj)*
1186     & (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
1187     & ) ) * 0.25 _d 0
1188     tmp1k(i,j) = - dyG(i,j,bi,bj)*drF(k)
1189     & * _hFacW(i,j,k,bi,bj)
1190     & *op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1191     & * dTdz
1192     ENDDO
1193     ENDDO
1194     CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KuzTz', k,1,2,bi,bj,myThid)
1195     C print *,'SM_KuzTz',k,tmp1k
1196     ENDIF
1197    
1198     IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1199     DO j=1,sNy+1
1200     DO i=1,sNx
1201     C- Vertical gradients interpolated to V points
1202     dTdz = op5*(
1203     & +op5*recip_drC(k)*
1204     & ( maskC(i,j-1,k,bi,bj)*
1205     & (Theta(i,j-1,km1,bi,bj)-Theta(i,j-1,k,bi,bj))
1206     & +maskC(i, j ,k,bi,bj)*
1207     & (Theta(i, j ,km1,bi,bj)-Theta(i, j ,k,bi,bj))
1208     & )
1209     & +op5*recip_drC(kp1)*
1210     & ( maskC(i,j-1,kp1,bi,bj)*
1211     & (Theta(i,j-1,k,bi,bj)-Theta(i,j-1,kp1,bi,bj))
1212     & +maskC(i, j ,kp1,bi,bj)*
1213     & (Theta(i, j ,k,bi,bj)-Theta(i, j ,kp1,bi,bj))
1214     & ) )
1215     tmp1k(i,j) = - dxG(i,j,bi,bj)*drF(k)
1216     & * _hFacS(i,j,k,bi,bj)
1217     & *op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1218     & * dTdz
1219     ENDDO
1220     ENDDO
1221     CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KvzTz', k,1,2,bi,bj,myThid)
1222     C print *,'SM_KvzTz',k,tmp1k
1223     ENDIF
1224    
1225     IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1226     DO j=jMin,jMax
1227     DO i=iMin,iMax
1228     C- Horizontal gradients interpolated to W points
1229     dTdx = op5*(
1230     & +op5*(_maskW(i+1,j,k,bi,bj)
1231     & *_recip_dxC(i+1,j,bi,bj)*
1232     & (Theta(i+1,j,k,bi,bj)-Theta(i,j,k,bi,bj))
1233     & +_maskW(i,j,k,bi,bj)
1234     & *_recip_dxC(i,j,bi,bj)*
1235     & (Theta(i,j,k,bi,bj)-Theta(i-1,j,k,bi,bj)))
1236     & +op5*(_maskW(i+1,j,k-1,bi,bj)
1237     & *_recip_dxC(i+1,j,bi,bj)*
1238     & (Theta(i+1,j,k-1,bi,bj)-Theta(i,j,k-1,bi,bj))
1239     & +_maskW(i,j,k-1,bi,bj)
1240     & *_recip_dxC(i,j,bi,bj)*
1241     & (Theta(i,j,k-1,bi,bj)-Theta(i-1,j,k-1,bi,bj)))
1242     & )
1243    
1244     dTdy = op5*(
1245     & +op5*(_maskS(i,j,k,bi,bj)
1246     & *_recip_dyC(i,j,bi,bj)*
1247     & (Theta(i,j,k,bi,bj)-Theta(i,j-1,k,bi,bj))
1248     & +_maskS(i,j+1,k,bi,bj)
1249     & *_recip_dyC(i,j+1,bi,bj)*
1250     & (Theta(i,j+1,k,bi,bj)-Theta(i,j,k,bi,bj)))
1251     & +op5*(_maskS(i,j,k-1,bi,bj)
1252     & *_recip_dyC(i,j,bi,bj)*
1253     & (Theta(i,j,k-1,bi,bj)-Theta(i,j-1,k-1,bi,bj))
1254     & +_maskS(i,j+1,k-1,bi,bj)
1255     & *_recip_dyC(i,j+1,bi,bj)*
1256     & (Theta(i,j+1,k-1,bi,bj)-Theta(i,j,k-1,bi,bj)))
1257     & )
1258    
1259     tmp1k(i,j) = - _rA(i,j,bi,bj)
1260     & *(op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))*dTdx
1261     & +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))*dTdy)
1262     ENDDO
1263     ENDDO
1264     CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', k,1,2,bi,bj,myThid)
1265     C print *,'SM_KrddT',k,tmp1k
1266     ENDIF
1267     ENDIF
1268     #endif
1269     DO j=1-Oly+1,sNy+Oly-1
1270     DO i=1-Olx+1,sNx+Olx-1
1271     SM_PsiXm1(i,j)=SM_PsiX(i,j)
1272     SM_PsiYm1(i,j)=SM_PsiY(i,j)
1273     tmp1k(i,j)=0 _d 0
1274     ENDDO
1275     ENDDO
1276     ENDDO
1277    
1278     CBFK Always Zero at the bottom.
1279     IF ( DIAGNOSTICS_IS_ON('SM_ubT ',myThid) ) THEN
1280     CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT ', Nr,1,2,bi,bj,myThid)
1281     ENDIF
1282     IF ( DIAGNOSTICS_IS_ON('SM_vbT ',myThid) ) THEN
1283     CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT ', Nr,1,2,bi,bj,myThid)
1284     ENDIF
1285     IF ( DIAGNOSTICS_IS_ON('SM_wbT ',myThid) ) THEN
1286     CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT ', Nr,1,2,bi,bj,myThid)
1287     ENDIF
1288     IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1289     CALL DIAGNOSTICS_FILL(tmp1k,'SM_KuzTz', Nr,1,2,bi,bj,myThid)
1290     ENDIF
1291     IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1292     CALL DIAGNOSTICS_FILL(tmp1k,'SM_KvzTz', Nr,1,2,bi,bj,myThid)
1293     ENDIF
1294     IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1295     CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', Nr,1,2,bi,bj,myThid)
1296     ENDIF
1297     #endif
1298    
1299 dimitri 1.1 #ifdef ALLOW_TIMEAVE
1300     C-- Time-average
1301     IF ( taveFreq.GT.0. ) THEN
1302    
1303     CALL TIMEAVE_CUMULATE( GM_Kwx_T, Kwx, Nr,
1304     & deltaTclock, bi, bj, myThid )
1305     CALL TIMEAVE_CUMULATE( GM_Kwy_T, Kwy, Nr,
1306     & deltaTclock, bi, bj, myThid )
1307     CALL TIMEAVE_CUMULATE( GM_Kwz_T, Kwz, Nr,
1308     & deltaTclock, bi, bj, myThid )
1309     #ifdef GM_VISBECK_VARIABLE_K
1310     IF ( GM_Visbeck_alpha.NE.0. ) THEN
1311     CALL TIMEAVE_CUMULATE( Visbeck_K_T, VisbeckK, 1,
1312     & deltaTclock, bi, bj, myThid )
1313     ENDIF
1314     #endif
1315     #ifdef GM_BOLUS_ADVEC
1316     IF ( GM_AdvForm ) THEN
1317     CALL TIMEAVE_CUMULATE( GM_PsiXtave, GM_PsiX, Nr,
1318     & deltaTclock, bi, bj, myThid )
1319     CALL TIMEAVE_CUMULATE( GM_PsiYtave, GM_PsiY, Nr,
1320     & deltaTclock, bi, bj, myThid )
1321     ENDIF
1322     #endif
1323 zhc 1.5 GM_timeAve(bi,bj) = GM_timeAve(bi,bj)+deltaTclock
1324 dimitri 1.1
1325     ENDIF
1326     #endif /* ALLOW_TIMEAVE */
1327    
1328     #ifdef ALLOW_DIAGNOSTICS
1329     IF ( useDiagnostics ) THEN
1330     CALL GMREDI_DIAGNOSTICS_FILL(bi,bj,myThid)
1331     ENDIF
1332     #endif /* ALLOW_DIAGNOSTICS */
1333    
1334     #endif /* ALLOW_GMREDI */
1335    
1336     RETURN
1337     END
1338    
1339 zhc 1.5 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1340 dimitri 1.1
1341     SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
1342     I iMin, iMax, jMin, jMax,
1343     I sigmaX, sigmaY, sigmaR,
1344     I bi, bj, myTime, myIter, myThid )
1345     C /==========================================================\
1346     C | SUBROUTINE GMREDI_CALC_TENSOR |
1347     C | o Calculate tensor elements for GM/Redi tensor. |
1348     C |==========================================================|
1349     C \==========================================================/
1350     IMPLICIT NONE
1351    
1352     C == Global variables ==
1353     #include "SIZE.h"
1354     #include "EEPARAMS.h"
1355     #include "GMREDI.h"
1356    
1357     C == Routine arguments ==
1358     C
1359     _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
1360     _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
1361     _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
1362     INTEGER iMin,iMax,jMin,jMax
1363     INTEGER bi, bj
1364     _RL myTime
1365     INTEGER myIter
1366     INTEGER myThid
1367     CEndOfInterface
1368    
1369     #ifdef ALLOW_GMREDI
1370    
1371     INTEGER i, j, k
1372    
1373     DO k=1,Nr
1374     DO j=1-Oly+1,sNy+Oly-1
1375     DO i=1-Olx+1,sNx+Olx-1
1376     Kwx(i,j,k,bi,bj) = 0.0
1377     Kwy(i,j,k,bi,bj) = 0.0
1378     Kwz(i,j,k,bi,bj) = 0.0
1379     ENDDO
1380     ENDDO
1381     ENDDO
1382     #endif /* ALLOW_GMREDI */
1383    
1384     RETURN
1385     END
  ViewVC Help
Powered by ViewVC 1.1.22