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