/[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.4 - (hide annotations) (download)
Fri May 30 23:18:08 2008 UTC (17 years, 5 months ago) by dimitri
Branch: MAIN
Changes since 1.3: +1 -5 lines 
Replacing KPPhbl with hMixLayer for submesoscale parameterization.
1 dimitri 1.4 C $Header: /u/gcmpack/MITgcm_contrib/submesoscale/code/gmredi_calc_tensor.F,v 1.3 2008/05/30 22:24:25 dimitri 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     hml=hMixLayer(i,j,bi,bj)
1014     IF (hml.gt.0 _d 0) THEN
1015     recip_hml=1 _d 0/hml
1016     ELSE
1017     recip_hml=0 _d 0
1018     ENDIF
1019     CBFK We calculate the h^2 mu(z) factor only on w points.
1020     CBFK It is possible that we might need to calculate it
1021     CBFK on Psi or u,v points independently to prevent spurious
1022     CBFK entrainment. Unlikely that this will be major
1023     CBFK (it wasnt in offline testing).
1024     qfac=(2*rf(k)*recip_hml+1 _d 0)**2
1025     hsqmu=(1 _d 0-qfac)*(1 _d 0+(5 _d 0)*qfac/21 _d 0)
1026     hsqmu=max(0 _d 0, hsqmu)*hml**2
1027     SM_Lf(i,j)=SM_Lf(i,j)*hsqmu
1028     ENDDO
1029     ENDDO
1030     CBFK Now interpolate to match locations
1031     DO j=1-Oly+1,sNy+Oly-1
1032     DO i=1-Olx+1,sNx+Olx-1
1033     C SM_Lf coefficients are on rVel points
1034     C Psix are on faces above U
1035     SM_PsiX(i,j)=op5*(SM_Lf(i+1,j)+SM_Lf(i,j))*dBdxAV(i,j)
1036     & *_maskW(i,j,k,bi,bj)
1037     C Psiy are on faces above V
1038     SM_PsiY(i,j)=op5*(SM_Lf(i,j+1)+SM_Lf(i,j))*dBdyAV(i,j)
1039     & *_maskS(i,j,k,bi,bj)
1040    
1041     #ifndef GM_BOLUS_ADVEC
1042     C Kwx,Kwy are on rVel Points
1043     Kwx(i,j,k,bi,bj) = Kwx(i,j,k,bi,bj)
1044     & +op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))
1045     Kwy(i,j,k,bi,bj) = Kwy(i,j,k,bi,bj)
1046     & +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))
1047     #ifdef GM_EXTRA_DIAGONAL
1048     IF (GM_ExtraDiag) THEN
1049     C Kuz,Kvz are on u,v Points
1050     Kuz(i,j,k,bi,bj) = Kuz(i,j,k,bi,bj)
1051     & -op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1052     Kvz(i,j,k,bi,bj) = Kvz(i,j,k,bi,bj)
1053     & -op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1054     ENDIF
1055     #endif
1056     #else
1057     IF (GM_AdvForm) THEN
1058     GM_PsiX(i,j,k,bi,bj)=GM_PsiX(i,j,k,bi,bj)+SM_PsiX(i,j)
1059     GM_PsiY(i,j,k,bi,bj)=GM_PsiY(i,j,k,bi,bj)+SM_PsiY(i,j)
1060     ENDIF
1061     #endif
1062     ENDDO
1063     ENDDO
1064     ELSE
1065     DO j=1-Oly+1,sNy+Oly-1
1066     DO i=1-Olx+1,sNx+Olx-1
1067     SM_PsiX(i,j)=0. _d 0
1068     SM_PsiY(i,j)=0. _d 0
1069     ENDDO
1070     ENDDO
1071     ENDIF
1072    
1073     #ifdef ALLOW_DIAGNOSTICS
1074     IF ( useDiagnostics ) THEN
1075     IF ( DIAGNOSTICS_IS_ON('SM_PsiX ',myThid) ) THEN
1076     CALL DIAGNOSTICS_FILL(SM_PsiX,'SM_PsiX ',k,1,2,bi,bj,myThid)
1077     ENDIF
1078     IF ( DIAGNOSTICS_IS_ON('SM_PsiY ',myThid) ) THEN
1079     CALL DIAGNOSTICS_FILL(SM_PsiY,'SM_PsiY ',k,1,2,bi,bj,myThid)
1080     ENDIF
1081    
1082     CBFK Note: for comparision, you can diagnose the bolus form
1083     CBFK or the Kappa form in the same simulation, regardless of other
1084     CBFK settings
1085     IF ( DIAGNOSTICS_IS_ON('SM_ubT ',myThid) ) THEN
1086     DO j=jMin,jMax
1087     DO i=iMin,iMax
1088     tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiX(i,j)
1089     & -SM_PsiXm1(i,j) )
1090     & *maskW(i,j,km1,bi,bj)
1091     & *op5*(Theta(i,j,km1,bi,bj)+Theta(i-1,j,km1,bi,bj))
1092     ENDDO
1093     ENDDO
1094     CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT ', km1,1,2,bi,bj,myThid)
1095     ENDIF
1096    
1097     IF ( DIAGNOSTICS_IS_ON('SM_vbT ',myThid) ) THEN
1098     DO j=jMin,jMax
1099     DO i=iMin,iMax
1100     tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiY(i,j)
1101     & -SM_PsiYm1(i,j) )
1102     & *maskS(i,j,km1,bi,bj)
1103     & *op5*(Theta(i,j,km1,bi,bj)+Theta(i,j-1,km1,bi,bj))
1104     ENDDO
1105     ENDDO
1106     CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT ', km1,1,2,bi,bj,myThid)
1107     ENDIF
1108    
1109     IF ( DIAGNOSTICS_IS_ON('SM_wbT ',myThid) ) THEN
1110     DO j=jMin,jMax
1111     DO i=iMin,iMax
1112     tmp1k(i,j) =
1113     & (dyG(i+1,j,bi,bj)*SM_PsiX(i+1,j)
1114     & -dyG( i ,j,bi,bj)*SM_PsiX( i ,j)
1115     & +dxG(i,j+1,bi,bj)*SM_PsiY(i,j+1)
1116     & -dxG(i, j ,bi,bj)*SM_PsiY(i, j ))
1117     & *op5*(Theta(i,j,k,bi,bj)+Theta(i,j,km1,bi,bj))
1118     ENDDO
1119     ENDDO
1120     CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT ', k,1,2,bi,bj,myThid)
1121     C print *,'SM_wbT',k,tmp1k
1122     ENDIF
1123    
1124     IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1125     DO j=1,sNy
1126     DO i=1,sNx+1
1127     C- Vertical gradients interpolated to U points
1128     dTdz = (
1129     & +recip_drC(k)*
1130     & ( maskC(i-1,j,k,bi,bj)*
1131     & (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
1132     & +maskC( i ,j,k,bi,bj)*
1133     & (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
1134     & )
1135     & +recip_drC(kp1)*
1136     & ( maskC(i-1,j,kp1,bi,bj)*
1137     & (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
1138     & +maskC( i ,j,kp1,bi,bj)*
1139     & (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
1140     & ) ) * 0.25 _d 0
1141     tmp1k(i,j) = - dyG(i,j,bi,bj)*drF(k)
1142     & * _hFacW(i,j,k,bi,bj)
1143     & *op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1144     & * dTdz
1145     ENDDO
1146     ENDDO
1147     CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KuzTz', k,1,2,bi,bj,myThid)
1148     C print *,'SM_KuzTz',k,tmp1k
1149     ENDIF
1150    
1151     IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1152     DO j=1,sNy+1
1153     DO i=1,sNx
1154     C- Vertical gradients interpolated to V points
1155     dTdz = op5*(
1156     & +op5*recip_drC(k)*
1157     & ( maskC(i,j-1,k,bi,bj)*
1158     & (Theta(i,j-1,km1,bi,bj)-Theta(i,j-1,k,bi,bj))
1159     & +maskC(i, j ,k,bi,bj)*
1160     & (Theta(i, j ,km1,bi,bj)-Theta(i, j ,k,bi,bj))
1161     & )
1162     & +op5*recip_drC(kp1)*
1163     & ( maskC(i,j-1,kp1,bi,bj)*
1164     & (Theta(i,j-1,k,bi,bj)-Theta(i,j-1,kp1,bi,bj))
1165     & +maskC(i, j ,kp1,bi,bj)*
1166     & (Theta(i, j ,k,bi,bj)-Theta(i, j ,kp1,bi,bj))
1167     & ) )
1168     tmp1k(i,j) = - dxG(i,j,bi,bj)*drF(k)
1169     & * _hFacS(i,j,k,bi,bj)
1170     & *op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1171     & * dTdz
1172     ENDDO
1173     ENDDO
1174     CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KvzTz', k,1,2,bi,bj,myThid)
1175     C print *,'SM_KvzTz',k,tmp1k
1176     ENDIF
1177    
1178     IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1179     DO j=jMin,jMax
1180     DO i=iMin,iMax
1181     C- Horizontal gradients interpolated to W points
1182     dTdx = op5*(
1183     & +op5*(_maskW(i+1,j,k,bi,bj)
1184     & *_recip_dxC(i+1,j,bi,bj)*
1185     & (Theta(i+1,j,k,bi,bj)-Theta(i,j,k,bi,bj))
1186     & +_maskW(i,j,k,bi,bj)
1187     & *_recip_dxC(i,j,bi,bj)*
1188     & (Theta(i,j,k,bi,bj)-Theta(i-1,j,k,bi,bj)))
1189     & +op5*(_maskW(i+1,j,k-1,bi,bj)
1190     & *_recip_dxC(i+1,j,bi,bj)*
1191     & (Theta(i+1,j,k-1,bi,bj)-Theta(i,j,k-1,bi,bj))
1192     & +_maskW(i,j,k-1,bi,bj)
1193     & *_recip_dxC(i,j,bi,bj)*
1194     & (Theta(i,j,k-1,bi,bj)-Theta(i-1,j,k-1,bi,bj)))
1195     & )
1196    
1197     dTdy = op5*(
1198     & +op5*(_maskS(i,j,k,bi,bj)
1199     & *_recip_dyC(i,j,bi,bj)*
1200     & (Theta(i,j,k,bi,bj)-Theta(i,j-1,k,bi,bj))
1201     & +_maskS(i,j+1,k,bi,bj)
1202     & *_recip_dyC(i,j+1,bi,bj)*
1203     & (Theta(i,j+1,k,bi,bj)-Theta(i,j,k,bi,bj)))
1204     & +op5*(_maskS(i,j,k-1,bi,bj)
1205     & *_recip_dyC(i,j,bi,bj)*
1206     & (Theta(i,j,k-1,bi,bj)-Theta(i,j-1,k-1,bi,bj))
1207     & +_maskS(i,j+1,k-1,bi,bj)
1208     & *_recip_dyC(i,j+1,bi,bj)*
1209     & (Theta(i,j+1,k-1,bi,bj)-Theta(i,j,k-1,bi,bj)))
1210     & )
1211    
1212     tmp1k(i,j) = - _rA(i,j,bi,bj)
1213     & *(op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))*dTdx
1214     & +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))*dTdy)
1215     ENDDO
1216     ENDDO
1217     CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', k,1,2,bi,bj,myThid)
1218     C print *,'SM_KrddT',k,tmp1k
1219     ENDIF
1220     ENDIF
1221     #endif
1222     DO j=1-Oly+1,sNy+Oly-1
1223     DO i=1-Olx+1,sNx+Olx-1
1224     SM_PsiXm1(i,j)=SM_PsiX(i,j)
1225     SM_PsiYm1(i,j)=SM_PsiY(i,j)
1226     tmp1k(i,j)=0 _d 0
1227     ENDDO
1228     ENDDO
1229     ENDDO
1230    
1231     CBFK Always Zero at the bottom.
1232     IF ( DIAGNOSTICS_IS_ON('SM_ubT ',myThid) ) THEN
1233     CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT ', Nr,1,2,bi,bj,myThid)
1234     ENDIF
1235     IF ( DIAGNOSTICS_IS_ON('SM_vbT ',myThid) ) THEN
1236     CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT ', Nr,1,2,bi,bj,myThid)
1237     ENDIF
1238     IF ( DIAGNOSTICS_IS_ON('SM_wbT ',myThid) ) THEN
1239     CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT ', Nr,1,2,bi,bj,myThid)
1240     ENDIF
1241     IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1242     CALL DIAGNOSTICS_FILL(tmp1k,'SM_KuzTz', Nr,1,2,bi,bj,myThid)
1243     ENDIF
1244     IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1245     CALL DIAGNOSTICS_FILL(tmp1k,'SM_KvzTz', Nr,1,2,bi,bj,myThid)
1246     ENDIF
1247     IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1248     CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', Nr,1,2,bi,bj,myThid)
1249     ENDIF
1250     #endif
1251    
1252 dimitri 1.1 #ifdef ALLOW_TIMEAVE
1253     C-- Time-average
1254     IF ( taveFreq.GT.0. ) THEN
1255    
1256     CALL TIMEAVE_CUMULATE( GM_Kwx_T, Kwx, Nr,
1257     & deltaTclock, bi, bj, myThid )
1258     CALL TIMEAVE_CUMULATE( GM_Kwy_T, Kwy, Nr,
1259     & deltaTclock, bi, bj, myThid )
1260     CALL TIMEAVE_CUMULATE( GM_Kwz_T, Kwz, Nr,
1261     & deltaTclock, bi, bj, myThid )
1262     #ifdef GM_VISBECK_VARIABLE_K
1263     IF ( GM_Visbeck_alpha.NE.0. ) THEN
1264     CALL TIMEAVE_CUMULATE( Visbeck_K_T, VisbeckK, 1,
1265     & deltaTclock, bi, bj, myThid )
1266     ENDIF
1267     #endif
1268     #ifdef GM_BOLUS_ADVEC
1269     IF ( GM_AdvForm ) THEN
1270     CALL TIMEAVE_CUMULATE( GM_PsiXtave, GM_PsiX, Nr,
1271     & deltaTclock, bi, bj, myThid )
1272     CALL TIMEAVE_CUMULATE( GM_PsiYtave, GM_PsiY, Nr,
1273     & deltaTclock, bi, bj, myThid )
1274     ENDIF
1275     #endif
1276     DO k=1,Nr
1277     GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock
1278     ENDDO
1279    
1280     ENDIF
1281     #endif /* ALLOW_TIMEAVE */
1282    
1283     #ifdef ALLOW_DIAGNOSTICS
1284     IF ( useDiagnostics ) THEN
1285     CALL GMREDI_DIAGNOSTICS_FILL(bi,bj,myThid)
1286     ENDIF
1287     #endif /* ALLOW_DIAGNOSTICS */
1288    
1289     #endif /* ALLOW_GMREDI */
1290    
1291     RETURN
1292     END
1293    
1294    
1295     SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
1296     I iMin, iMax, jMin, jMax,
1297     I sigmaX, sigmaY, sigmaR,
1298     I bi, bj, myTime, myIter, myThid )
1299     C /==========================================================\
1300     C | SUBROUTINE GMREDI_CALC_TENSOR |
1301     C | o Calculate tensor elements for GM/Redi tensor. |
1302     C |==========================================================|
1303     C \==========================================================/
1304     IMPLICIT NONE
1305    
1306     C == Global variables ==
1307     #include "SIZE.h"
1308     #include "EEPARAMS.h"
1309     #include "GMREDI.h"
1310    
1311     C == Routine arguments ==
1312     C
1313     _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
1314     _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
1315     _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
1316     INTEGER iMin,iMax,jMin,jMax
1317     INTEGER bi, bj
1318     _RL myTime
1319     INTEGER myIter
1320     INTEGER myThid
1321     CEndOfInterface
1322    
1323     #ifdef ALLOW_GMREDI
1324    
1325     INTEGER i, j, k
1326    
1327     DO k=1,Nr
1328     DO j=1-Oly+1,sNy+Oly-1
1329     DO i=1-Olx+1,sNx+Olx-1
1330     Kwx(i,j,k,bi,bj) = 0.0
1331     Kwy(i,j,k,bi,bj) = 0.0
1332     Kwz(i,j,k,bi,bj) = 0.0
1333     ENDDO
1334     ENDDO
1335     ENDDO
1336     #endif /* ALLOW_GMREDI */
1337    
1338     RETURN
1339     END
  ViewVC Help
Powered by ViewVC 1.1.22