/[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.2 - (hide annotations) (download)
Fri May 30 22:13:42 2008 UTC (17 years, 1 month ago) by dimitri
Branch: MAIN
Changes since 1.1: +393 -10 lines 
Initial code submitted by Baylor on May 19, 2008.  See:
http://forge.csail.mit.edu/pipermail/mitgcm-devel/2008-May/003392.html
1 dimitri 1.1 C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_calc_tensor.F,v 1.31 2008/02/02 02:35:53 gforget Exp $
2     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     #ifdef ALLOW_KAPGM_CONTROL
535     & + 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     #ifdef ALLOW_KAPGM_CONTROL
719     & - 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     #ifdef ALLOW_KPP
780     IF ( useKPP ) THEN
781     DO j=2-Oly,sNy+Oly
782     DO i=1-Olx,sNx+Olx
783     locMixLayer(i,j) = ( KPPhbl(i,j-1,bi,bj)
784     & + KPPhbl(i, j ,bi,bj) )*op5
785     ENDDO
786     ENDDO
787     ELSE
788     #else
789     IF ( .TRUE. ) THEN
790     #endif
791     DO j=2-Oly,sNy+Oly
792     DO i=1-Olx,sNx+Olx
793     locMixLayer(i,j) = ( hMixLayer(i,j-1,bi,bj)
794     & + hMixLayer(i, j ,bi,bj) )*op5
795     ENDDO
796     ENDDO
797     ENDIF
798     DO j=1-Oly,sNy+Oly
799     DO i=1-Olx,sNx+Olx
800     hTransLay(i,j) = 0.
801     baseSlope(i,j) = 0.
802     recipLambda(i,j)= 0.
803     ENDDO
804     ENDDO
805     DO j=2-Oly,sNy+Oly
806     DO i=1-Olx,sNx+Olx
807     hTransLay(i,j) = MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
808     ENDDO
809     ENDDO
810    
811     C Gradient of Sigma at V points
812     DO k=Nr,1,-1
813     kp1 = MIN(Nr,k+1)
814     maskp1 = 1. _d 0
815     IF (k.GE.Nr) maskp1 = 0. _d 0
816     #ifdef ALLOW_AUTODIFF_TAMC
817     kkey = (igmkey-1)*Nr + k
818     #endif
819    
820     DO j=1-Oly+1,sNy+Oly-1
821     DO i=1-Olx+1,sNx+Olx-1
822     dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
823     & +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k)
824     & )*_maskS(i,j,k,bi,bj)
825     dSigmaDy(i,j)=sigmaY(i,j,k)
826     & *_maskS(i,j,k,bi,bj)
827     dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k )+sigmaR(i,j, k )
828     & +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1
829     & )*_maskS(i,j,k,bi,bj)
830 dimitri 1.2
831     #ifdef GM_SUBMESO
832     C-- Depth average of SigmaY at V points
833     C compute depth average from surface down to the MixLayer depth
834     IF (k.GT.1) THEN
835     IF (-rC(k-1).LT.locMixLayer(i,j) ) THEN
836     IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
837     integrDepth = -rC( k )
838     C- in 2 steps to avoid mix of RS & RL type in min fct. arguments
839     integrDepth = MIN( integrDepth, locMixLayer(i,j) )
840     C Distance between level center above and the integration depth
841     deltaH = integrDepth + rC(k-1)
842     C If negative then we are below the integration level
843     C (cannot be the case with 2 conditions on maskC & -rC(k-1))
844     C If positive we limit this to the distance from center above
845     deltaH = MIN( deltaH, drC(k) )
846     C Now we convert deltaH to a non-dimensional fraction
847     deltaH = deltaH/( integrDepth+rC(1) )
848     dBdyAV(i,j) = dBdyAV(i,j)
849     & +dSigmaDy(i,j)*deltaH*recip_rhoConst*gravity
850     ENDIF
851     ENDIF
852     ENDIF
853     #endif
854    
855 dimitri 1.1 ENDDO
856     ENDDO
857    
858     #ifdef ALLOW_AUTODIFF_TAMC
859     CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
860     CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
861     CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
862     CADJ STORE baseSlope(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
863     CADJ STORE hTransLay(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
864     CADJ STORE recipLambda(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
865     #endif /* ALLOW_AUTODIFF_TAMC */
866    
867     C Calculate slopes for use in tensor, taper and/or clip
868     CALL GMREDI_SLOPE_LIMIT(
869     O SlopeX, SlopeY,
870     O SlopeSqr, taperFct,
871     U hTransLay, baseSlope, recipLambda,
872     U dSigmaDr,
873     I dSigmaDx, dSigmaDy,
874     I ldd97_LrhoS, locMixLayer, rC,
875     I kLow_S,
876     I k, bi, bj, myTime, myIter, myThid )
877    
878     cph(
879     #ifdef ALLOW_AUTODIFF_TAMC
880     cph(
881     CADJ STORE taperfct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
882     cph)
883     #endif /* ALLOW_AUTODIFF_TAMC */
884     cph)
885    
886     #ifdef GM_NON_UNITY_DIAGONAL
887     c IF ( GM_nonUnitDiag ) THEN
888     DO j=1-Oly+1,sNy+Oly-1
889     DO i=1-Olx+1,sNx+Olx-1
890     Kvy(i,j,k,bi,bj) =
891     #ifdef ALLOW_KAPREDI_CONTROL
892     & ( kapredi(i,j,k,bi,bj)
893     #else
894     & ( GM_isopycK
895     #endif
896     #ifdef GM_VISBECK_VARIABLE_K
897     & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
898     #endif
899     & )*taperFct(i,j)
900     ENDDO
901     ENDDO
902     #ifdef ALLOW_AUTODIFF_TAMC
903     # ifdef GM_EXCLUDE_CLIPPING
904     CADJ STORE Kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
905     # endif
906     #endif
907     DO j=1-Oly+1,sNy+Oly-1
908     DO i=1-Olx+1,sNx+Olx-1
909     Kvy(i,j,k,bi,bj) = MAX( Kvy(i,j,k,bi,bj), GM_Kmin_horiz )
910     ENDDO
911     ENDDO
912     c ENDIF
913     #endif /* GM_NON_UNITY_DIAGONAL */
914    
915     #ifdef GM_EXTRA_DIAGONAL
916    
917     #ifdef ALLOW_AUTODIFF_TAMC
918     CADJ STORE SlopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
919     CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
920     #endif /* ALLOW_AUTODIFF_TAMC */
921     IF ( GM_ExtraDiag ) THEN
922     DO j=1-Oly+1,sNy+Oly-1
923     DO i=1-Olx+1,sNx+Olx-1
924     Kvz(i,j,k,bi,bj) =
925     #ifdef ALLOW_KAPREDI_CONTROL
926     & ( kapredi(i,j,k,bi,bj)
927     #else
928     & ( GM_isopycK
929     #endif
930     #ifdef ALLOW_KAPGM_CONTROL
931     & - GM_skewflx*kapgm(i,j,k,bi,bj)
932     #else
933     & - GM_skewflx*GM_background_K
934     #endif
935     #ifdef GM_VISBECK_VARIABLE_K
936     & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*GM_advect
937     #endif
938     & )*SlopeY(i,j)*taperFct(i,j)
939     ENDDO
940     ENDDO
941     ENDIF
942     #endif /* GM_EXTRA_DIAGONAL */
943    
944     #ifdef ALLOW_DIAGNOSTICS
945     IF (doDiagRediFlx) THEN
946     km1 = MAX(k-1,1)
947     DO j=1,sNy+1
948     DO i=1,sNx
949     C store in tmp1k Kvz_Redi
950     #ifdef ALLOW_KAPREDI_CONTROL
951     tmp1k(i,j) = ( kapredi(i,j,k,bi,bj)
952     #else
953     tmp1k(i,j) = ( GM_isopycK
954     #endif
955     #ifdef GM_VISBECK_VARIABLE_K
956     & +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0
957     #endif
958     & )*SlopeY(i,j)*taperFct(i,j)
959     ENDDO
960     ENDDO
961     DO j=1,sNy+1
962     DO i=1,sNx
963 dimitri 1.2 C- Vertical gradients interpolated to V points
964     dTdz = op5*(
965     & +op5*recip_drC(k)*
966 dimitri 1.1 & ( maskC(i,j-1,k,bi,bj)*
967     & (theta(i,j-1,km1,bi,bj)-theta(i,j-1,k,bi,bj))
968     & +maskC(i, j ,k,bi,bj)*
969     & (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))
970     & )
971 dimitri 1.2 & +op5*recip_drC(kp1)*
972 dimitri 1.1 & ( maskC(i,j-1,kp1,bi,bj)*
973     & (theta(i,j-1,k,bi,bj)-theta(i,j-1,kp1,bi,bj))
974     & +maskC(i, j ,kp1,bi,bj)*
975     & (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))
976 dimitri 1.2 & ) )
977 dimitri 1.1 tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)
978     & * _hFacS(i,j,k,bi,bj)
979     & * tmp1k(i,j) * dTdz
980     ENDDO
981     ENDDO
982     CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KvzTz', k,1,2,bi,bj,myThid)
983     ENDIF
984     #endif /* ALLOW_DIAGNOSTICS */
985    
986     C-- end 3rd loop on vertical level index k
987     ENDDO
988    
989     #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
990    
991    
992     #ifdef GM_BOLUS_ADVEC
993     IF (GM_AdvForm) THEN
994     CALL GMREDI_CALC_PSI_B(
995     I bi, bj, iMin, iMax, jMin, jMax,
996     I sigmaX, sigmaY, sigmaR,
997     I ldd97_LrhoW, ldd97_LrhoS,
998     I myThid )
999     ENDIF
1000     #endif
1001    
1002 dimitri 1.2
1003     #ifdef GM_SUBMESO
1004     CBFK Add the submesoscale contribution, in a 4th k loop
1005     DO k=1,Nr
1006     km1=max(1,k-1)
1007     IF ((k.gt.1).and.(-rF(k-1) .lt. mlmax)) THEN
1008     kp1 = MIN(k+1,Nr)
1009     CBFK Add in the mu vertical structure factor
1010     DO j=1-Oly+1,sNy+Oly-1
1011     DO i=1-Olx+1,sNx+Olx-1
1012     #ifdef ALLOW_KPP
1013     hml=KPPhbl(i,j,bi,bj)
1014     #else
1015     hml=hMixLayer(i,j,bi,bj)
1016     #endif
1017     IF (hml.gt.0 _d 0) THEN
1018     recip_hml=1 _d 0/hml
1019     ELSE
1020     recip_hml=0 _d 0
1021     ENDIF
1022     CBFK We calculate the h^2 mu(z) factor only on w points.
1023     CBFK It is possible that we might need to calculate it
1024     CBFK on Psi or u,v points independently to prevent spurious
1025     CBFK entrainment. Unlikely that this will be major
1026     CBFK (it wasnt in offline testing).
1027     qfac=(2*rf(k)*recip_hml+1 _d 0)**2
1028     hsqmu=(1 _d 0-qfac)*(1 _d 0+(5 _d 0)*qfac/21 _d 0)
1029     hsqmu=max(0 _d 0, hsqmu)*hml**2
1030     SM_Lf(i,j)=SM_Lf(i,j)*hsqmu
1031     ENDDO
1032     ENDDO
1033     CBFK Now interpolate to match locations
1034     DO j=1-Oly+1,sNy+Oly-1
1035     DO i=1-Olx+1,sNx+Olx-1
1036     C SM_Lf coefficients are on rVel points
1037     C Psix are on faces above U
1038     SM_PsiX(i,j)=op5*(SM_Lf(i+1,j)+SM_Lf(i,j))*dBdxAV(i,j)
1039     & *_maskW(i,j,k,bi,bj)
1040     C Psiy are on faces above V
1041     SM_PsiY(i,j)=op5*(SM_Lf(i,j+1)+SM_Lf(i,j))*dBdyAV(i,j)
1042     & *_maskS(i,j,k,bi,bj)
1043    
1044     #ifndef GM_BOLUS_ADVEC
1045     C Kwx,Kwy are on rVel Points
1046     Kwx(i,j,k,bi,bj) = Kwx(i,j,k,bi,bj)
1047     & +op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))
1048     Kwy(i,j,k,bi,bj) = Kwy(i,j,k,bi,bj)
1049     & +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))
1050     #ifdef GM_EXTRA_DIAGONAL
1051     IF (GM_ExtraDiag) THEN
1052     C Kuz,Kvz are on u,v Points
1053     Kuz(i,j,k,bi,bj) = Kuz(i,j,k,bi,bj)
1054     & -op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1055     Kvz(i,j,k,bi,bj) = Kvz(i,j,k,bi,bj)
1056     & -op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1057     ENDIF
1058     #endif
1059     #else
1060     IF (GM_AdvForm) THEN
1061     GM_PsiX(i,j,k,bi,bj)=GM_PsiX(i,j,k,bi,bj)+SM_PsiX(i,j)
1062     GM_PsiY(i,j,k,bi,bj)=GM_PsiY(i,j,k,bi,bj)+SM_PsiY(i,j)
1063     ENDIF
1064     #endif
1065     ENDDO
1066     ENDDO
1067     ELSE
1068     DO j=1-Oly+1,sNy+Oly-1
1069     DO i=1-Olx+1,sNx+Olx-1
1070     SM_PsiX(i,j)=0. _d 0
1071     SM_PsiY(i,j)=0. _d 0
1072     ENDDO
1073     ENDDO
1074     ENDIF
1075    
1076     #ifdef ALLOW_DIAGNOSTICS
1077     IF ( useDiagnostics ) THEN
1078     IF ( DIAGNOSTICS_IS_ON('SM_PsiX ',myThid) ) THEN
1079     CALL DIAGNOSTICS_FILL(SM_PsiX,'SM_PsiX ',k,1,2,bi,bj,myThid)
1080     ENDIF
1081     IF ( DIAGNOSTICS_IS_ON('SM_PsiY ',myThid) ) THEN
1082     CALL DIAGNOSTICS_FILL(SM_PsiY,'SM_PsiY ',k,1,2,bi,bj,myThid)
1083     ENDIF
1084    
1085     CBFK Note: for comparision, you can diagnose the bolus form
1086     CBFK or the Kappa form in the same simulation, regardless of other
1087     CBFK settings
1088     IF ( DIAGNOSTICS_IS_ON('SM_ubT ',myThid) ) THEN
1089     DO j=jMin,jMax
1090     DO i=iMin,iMax
1091     tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiX(i,j)
1092     & -SM_PsiXm1(i,j) )
1093     & *maskW(i,j,km1,bi,bj)
1094     & *op5*(Theta(i,j,km1,bi,bj)+Theta(i-1,j,km1,bi,bj))
1095     ENDDO
1096     ENDDO
1097     CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT ', km1,1,2,bi,bj,myThid)
1098     ENDIF
1099    
1100     IF ( DIAGNOSTICS_IS_ON('SM_vbT ',myThid) ) THEN
1101     DO j=jMin,jMax
1102     DO i=iMin,iMax
1103     tmp1k(i,j) = dyG(i,j,bi,bj)*( SM_PsiY(i,j)
1104     & -SM_PsiYm1(i,j) )
1105     & *maskS(i,j,km1,bi,bj)
1106     & *op5*(Theta(i,j,km1,bi,bj)+Theta(i,j-1,km1,bi,bj))
1107     ENDDO
1108     ENDDO
1109     CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT ', km1,1,2,bi,bj,myThid)
1110     ENDIF
1111    
1112     IF ( DIAGNOSTICS_IS_ON('SM_wbT ',myThid) ) THEN
1113     DO j=jMin,jMax
1114     DO i=iMin,iMax
1115     tmp1k(i,j) =
1116     & (dyG(i+1,j,bi,bj)*SM_PsiX(i+1,j)
1117     & -dyG( i ,j,bi,bj)*SM_PsiX( i ,j)
1118     & +dxG(i,j+1,bi,bj)*SM_PsiY(i,j+1)
1119     & -dxG(i, j ,bi,bj)*SM_PsiY(i, j ))
1120     & *op5*(Theta(i,j,k,bi,bj)+Theta(i,j,km1,bi,bj))
1121     ENDDO
1122     ENDDO
1123     CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT ', k,1,2,bi,bj,myThid)
1124     C print *,'SM_wbT',k,tmp1k
1125     ENDIF
1126    
1127     IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1128     DO j=1,sNy
1129     DO i=1,sNx+1
1130     C- Vertical gradients interpolated to U points
1131     dTdz = (
1132     & +recip_drC(k)*
1133     & ( maskC(i-1,j,k,bi,bj)*
1134     & (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
1135     & +maskC( i ,j,k,bi,bj)*
1136     & (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
1137     & )
1138     & +recip_drC(kp1)*
1139     & ( maskC(i-1,j,kp1,bi,bj)*
1140     & (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
1141     & +maskC( i ,j,kp1,bi,bj)*
1142     & (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
1143     & ) ) * 0.25 _d 0
1144     tmp1k(i,j) = - dyG(i,j,bi,bj)*drF(k)
1145     & * _hFacW(i,j,k,bi,bj)
1146     & *op5*(SM_PsiX(i,j)+SM_PsiXm1(i+1,j))
1147     & * dTdz
1148     ENDDO
1149     ENDDO
1150     CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KuzTz', k,1,2,bi,bj,myThid)
1151     C print *,'SM_KuzTz',k,tmp1k
1152     ENDIF
1153    
1154     IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1155     DO j=1,sNy+1
1156     DO i=1,sNx
1157     C- Vertical gradients interpolated to V points
1158     dTdz = op5*(
1159     & +op5*recip_drC(k)*
1160     & ( maskC(i,j-1,k,bi,bj)*
1161     & (Theta(i,j-1,km1,bi,bj)-Theta(i,j-1,k,bi,bj))
1162     & +maskC(i, j ,k,bi,bj)*
1163     & (Theta(i, j ,km1,bi,bj)-Theta(i, j ,k,bi,bj))
1164     & )
1165     & +op5*recip_drC(kp1)*
1166     & ( maskC(i,j-1,kp1,bi,bj)*
1167     & (Theta(i,j-1,k,bi,bj)-Theta(i,j-1,kp1,bi,bj))
1168     & +maskC(i, j ,kp1,bi,bj)*
1169     & (Theta(i, j ,k,bi,bj)-Theta(i, j ,kp1,bi,bj))
1170     & ) )
1171     tmp1k(i,j) = - dxG(i,j,bi,bj)*drF(k)
1172     & * _hFacS(i,j,k,bi,bj)
1173     & *op5*(SM_PsiY(i,j)+SM_PsiYm1(i+1,j))
1174     & * dTdz
1175     ENDDO
1176     ENDDO
1177     CALL DIAGNOSTICS_FILL(tmp1k, 'SM_KvzTz', k,1,2,bi,bj,myThid)
1178     C print *,'SM_KvzTz',k,tmp1k
1179     ENDIF
1180    
1181     IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1182     DO j=jMin,jMax
1183     DO i=iMin,iMax
1184     C- Horizontal gradients interpolated to W points
1185     dTdx = op5*(
1186     & +op5*(_maskW(i+1,j,k,bi,bj)
1187     & *_recip_dxC(i+1,j,bi,bj)*
1188     & (Theta(i+1,j,k,bi,bj)-Theta(i,j,k,bi,bj))
1189     & +_maskW(i,j,k,bi,bj)
1190     & *_recip_dxC(i,j,bi,bj)*
1191     & (Theta(i,j,k,bi,bj)-Theta(i-1,j,k,bi,bj)))
1192     & +op5*(_maskW(i+1,j,k-1,bi,bj)
1193     & *_recip_dxC(i+1,j,bi,bj)*
1194     & (Theta(i+1,j,k-1,bi,bj)-Theta(i,j,k-1,bi,bj))
1195     & +_maskW(i,j,k-1,bi,bj)
1196     & *_recip_dxC(i,j,bi,bj)*
1197     & (Theta(i,j,k-1,bi,bj)-Theta(i-1,j,k-1,bi,bj)))
1198     & )
1199    
1200     dTdy = op5*(
1201     & +op5*(_maskS(i,j,k,bi,bj)
1202     & *_recip_dyC(i,j,bi,bj)*
1203     & (Theta(i,j,k,bi,bj)-Theta(i,j-1,k,bi,bj))
1204     & +_maskS(i,j+1,k,bi,bj)
1205     & *_recip_dyC(i,j+1,bi,bj)*
1206     & (Theta(i,j+1,k,bi,bj)-Theta(i,j,k,bi,bj)))
1207     & +op5*(_maskS(i,j,k-1,bi,bj)
1208     & *_recip_dyC(i,j,bi,bj)*
1209     & (Theta(i,j,k-1,bi,bj)-Theta(i,j-1,k-1,bi,bj))
1210     & +_maskS(i,j+1,k-1,bi,bj)
1211     & *_recip_dyC(i,j+1,bi,bj)*
1212     & (Theta(i,j+1,k-1,bi,bj)-Theta(i,j,k-1,bi,bj)))
1213     & )
1214    
1215     tmp1k(i,j) = - _rA(i,j,bi,bj)
1216     & *(op5*(SM_PsiX(i,j)+SM_PsiX(i+1,j))*dTdx
1217     & +op5*(SM_PsiX(i,j+1)+SM_PsiX(i,j))*dTdy)
1218     ENDDO
1219     ENDDO
1220     CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', k,1,2,bi,bj,myThid)
1221     C print *,'SM_KrddT',k,tmp1k
1222     ENDIF
1223     ENDIF
1224     #endif
1225     DO j=1-Oly+1,sNy+Oly-1
1226     DO i=1-Olx+1,sNx+Olx-1
1227     SM_PsiXm1(i,j)=SM_PsiX(i,j)
1228     SM_PsiYm1(i,j)=SM_PsiY(i,j)
1229     tmp1k(i,j)=0 _d 0
1230     ENDDO
1231     ENDDO
1232     ENDDO
1233    
1234     CBFK Always Zero at the bottom.
1235     IF ( DIAGNOSTICS_IS_ON('SM_ubT ',myThid) ) THEN
1236     CALL DIAGNOSTICS_FILL(tmp1k,'SM_ubT ', Nr,1,2,bi,bj,myThid)
1237     ENDIF
1238     IF ( DIAGNOSTICS_IS_ON('SM_vbT ',myThid) ) THEN
1239     CALL DIAGNOSTICS_FILL(tmp1k,'SM_vbT ', Nr,1,2,bi,bj,myThid)
1240     ENDIF
1241     IF ( DIAGNOSTICS_IS_ON('SM_wbT ',myThid) ) THEN
1242     CALL DIAGNOSTICS_FILL(tmp1k,'SM_wbT ', Nr,1,2,bi,bj,myThid)
1243     ENDIF
1244     IF ( DIAGNOSTICS_IS_ON('SM_KuzTz',myThid) ) THEN
1245     CALL DIAGNOSTICS_FILL(tmp1k,'SM_KuzTz', Nr,1,2,bi,bj,myThid)
1246     ENDIF
1247     IF ( DIAGNOSTICS_IS_ON('SM_KvzTz',myThid) ) THEN
1248     CALL DIAGNOSTICS_FILL(tmp1k,'SM_KvzTz', Nr,1,2,bi,bj,myThid)
1249     ENDIF
1250     IF ( DIAGNOSTICS_IS_ON('SM_KrddT',myThid) ) THEN
1251     CALL DIAGNOSTICS_FILL(tmp1k,'SM_KrddT', Nr,1,2,bi,bj,myThid)
1252     ENDIF
1253     #endif
1254    
1255 dimitri 1.1 #ifdef ALLOW_TIMEAVE
1256     C-- Time-average
1257     IF ( taveFreq.GT.0. ) THEN
1258    
1259     CALL TIMEAVE_CUMULATE( GM_Kwx_T, Kwx, Nr,
1260     & deltaTclock, bi, bj, myThid )
1261     CALL TIMEAVE_CUMULATE( GM_Kwy_T, Kwy, Nr,
1262     & deltaTclock, bi, bj, myThid )
1263     CALL TIMEAVE_CUMULATE( GM_Kwz_T, Kwz, Nr,
1264     & deltaTclock, bi, bj, myThid )
1265     #ifdef GM_VISBECK_VARIABLE_K
1266     IF ( GM_Visbeck_alpha.NE.0. ) THEN
1267     CALL TIMEAVE_CUMULATE( Visbeck_K_T, VisbeckK, 1,
1268     & deltaTclock, bi, bj, myThid )
1269     ENDIF
1270     #endif
1271     #ifdef GM_BOLUS_ADVEC
1272     IF ( GM_AdvForm ) THEN
1273     CALL TIMEAVE_CUMULATE( GM_PsiXtave, GM_PsiX, Nr,
1274     & deltaTclock, bi, bj, myThid )
1275     CALL TIMEAVE_CUMULATE( GM_PsiYtave, GM_PsiY, Nr,
1276     & deltaTclock, bi, bj, myThid )
1277     ENDIF
1278     #endif
1279     DO k=1,Nr
1280     GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock
1281     ENDDO
1282    
1283     ENDIF
1284     #endif /* ALLOW_TIMEAVE */
1285    
1286     #ifdef ALLOW_DIAGNOSTICS
1287     IF ( useDiagnostics ) THEN
1288     CALL GMREDI_DIAGNOSTICS_FILL(bi,bj,myThid)
1289     ENDIF
1290     #endif /* ALLOW_DIAGNOSTICS */
1291    
1292     #endif /* ALLOW_GMREDI */
1293    
1294     RETURN
1295     END
1296    
1297    
1298     SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
1299     I iMin, iMax, jMin, jMax,
1300     I sigmaX, sigmaY, sigmaR,
1301     I bi, bj, myTime, myIter, myThid )
1302     C /==========================================================\
1303     C | SUBROUTINE GMREDI_CALC_TENSOR |
1304     C | o Calculate tensor elements for GM/Redi tensor. |
1305     C |==========================================================|
1306     C \==========================================================/
1307     IMPLICIT NONE
1308    
1309     C == Global variables ==
1310     #include "SIZE.h"
1311     #include "EEPARAMS.h"
1312     #include "GMREDI.h"
1313    
1314     C == Routine arguments ==
1315     C
1316     _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
1317     _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
1318     _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
1319     INTEGER iMin,iMax,jMin,jMax
1320     INTEGER bi, bj
1321     _RL myTime
1322     INTEGER myIter
1323     INTEGER myThid
1324     CEndOfInterface
1325    
1326     #ifdef ALLOW_GMREDI
1327    
1328     INTEGER i, j, k
1329    
1330     DO k=1,Nr
1331     DO j=1-Oly+1,sNy+Oly-1
1332     DO i=1-Olx+1,sNx+Olx-1
1333     Kwx(i,j,k,bi,bj) = 0.0
1334     Kwy(i,j,k,bi,bj) = 0.0
1335     Kwz(i,j,k,bi,bj) = 0.0
1336     ENDDO
1337     ENDDO
1338     ENDDO
1339     #endif /* ALLOW_GMREDI */
1340    
1341     RETURN
1342     END
  ViewVC Help
Powered by ViewVC 1.1.22