/[MITgcm]/MITgcm/pkg/gmredi/gmredi_calc_tensor.F
ViewVC logotype

Annotation of /MITgcm/pkg/gmredi/gmredi_calc_tensor.F

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

Revision 1.33 - (hide annotations) (download)
Tue Oct 21 22:10:55 2008 UTC (15 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint61e
Changes since 1.32: +26 -22 lines 
- fix the "OLD_VISBECK_CALC" option + move the #undef to GMREDI_OPTIONS.h
- change computation of Visbeck-K where Slope > Smax :
  now: N*min(Slope,Smax) (similar to OLD_VISBECK_CALC with gkw91 taperFct)
  previously was: M*sqrt(min(Slope,Smax))
1 jmc 1.33 C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_calc_tensor.F,v 1.32 2008/03/28 18:48:05 heimbach Exp $
2 heimbach 1.13 C $Name: $
3 adcroft 1.1
4     #include "GMREDI_OPTIONS.h"
5 jmc 1.27 #ifdef ALLOW_KPP
6     # include "KPP_OPTIONS.h"
7     #endif
8 adcroft 1.1
9 jmc 1.27 CBOP
10     C !ROUTINE: GMREDI_CALC_TENSOR
11     C !INTERFACE:
12 adcroft 1.1 SUBROUTINE GMREDI_CALC_TENSOR(
13 jmc 1.27 I iMin, iMax, jMin, jMax,
14 adcroft 1.1 I sigmaX, sigmaY, sigmaR,
15 jmc 1.27 I bi, bj, myTime, myIter, myThid )
16    
17     C !DESCRIPTION: \bv
18     C *==========================================================*
19     C | SUBROUTINE GMREDI_CALC_TENSOR
20     C | o Calculate tensor elements for GM/Redi tensor.
21     C *==========================================================*
22     C *==========================================================*
23     C \ev
24    
25     C !USES:
26 adcroft 1.1 IMPLICIT NONE
27    
28     C == Global variables ==
29     #include "SIZE.h"
30     #include "GRID.h"
31     #include "DYNVARS.h"
32     #include "EEPARAMS.h"
33     #include "PARAMS.h"
34     #include "GMREDI.h"
35 jmc 1.20 #include "GMREDI_TAVE.h"
36 jmc 1.27 #ifdef ALLOW_KPP
37     # include "KPP.h"
38     #endif
39 adcroft 1.1
40 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
41     #include "tamc.h"
42     #include "tamc_keys.h"
43     #endif /* ALLOW_AUTODIFF_TAMC */
44    
45 jmc 1.27 C !INPUT/OUTPUT PARAMETERS:
46 adcroft 1.1 C == Routine arguments ==
47 jmc 1.27 C bi, bj :: tile indices
48     C myTime :: Current time in simulation
49     C myIter :: Current iteration number in simulation
50     C myThid :: My Thread Id. number
51 adcroft 1.1 C
52 jmc 1.27 INTEGER iMin,iMax,jMin,jMax
53 adcroft 1.1 _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
54     _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
55     _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
56 jmc 1.27 INTEGER bi, bj
57     _RL myTime
58     INTEGER myIter
59 adcroft 1.1 INTEGER myThid
60 jmc 1.27 CEOP
61 adcroft 1.1
62     #ifdef ALLOW_GMREDI
63    
64 jmc 1.27 C !LOCAL VARIABLES:
65 adcroft 1.1 C == Local variables ==
66 jmc 1.15 INTEGER i,j,k,kp1
67 adcroft 1.1 _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
68     _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
69 heimbach 1.12 _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
70     _RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
71 jmc 1.19 _RL dSigmaDr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
72 jmc 1.8 _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
73     _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
74 jmc 1.15 _RL maskp1, Kgm_tmp
75 jmc 1.19 _RL ldd97_LrhoC(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
76     _RL ldd97_LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
77     _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
78     _RL Cspd, LrhoInf, LrhoSup, fCoriLoc
79 adcroft 1.1
80 jmc 1.27 INTEGER kLow_W (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
81     INTEGER kLow_S (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
82     _RL locMixLayer(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
83     _RL baseSlope (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
84     _RL hTransLay (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
85     _RL recipLambda(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
86    
87 adcroft 1.1 #ifdef GM_VISBECK_VARIABLE_K
88 jmc 1.26 #ifdef OLD_VISBECK_CALC
89 heimbach 1.10 _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
90 jmc 1.26 #else
91     _RL dSigmaH
92 jmc 1.33 _RL Sloc, M2loc
93     #endif
94 jmc 1.26 _RL deltaH, integrDepth
95 jmc 1.33 _RL N2loc, SNloc
96 adcroft 1.1 #endif
97    
98 jmc 1.22 #ifdef ALLOW_DIAGNOSTICS
99     LOGICAL doDiagRediFlx
100     LOGICAL DIAGNOSTICS_IS_ON
101     EXTERNAL DIAGNOSTICS_IS_ON
102     INTEGER km1
103     _RL dTdz
104     _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105     #endif
106    
107 jmc 1.19 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
108    
109 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
110     act1 = bi - myBxLo(myThid)
111     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
112     act2 = bj - myByLo(myThid)
113     max2 = myByHi(myThid) - myByLo(myThid) + 1
114     act3 = myThid - 1
115     max3 = nTx*nTy
116     act4 = ikey_dynamics - 1
117 heimbach 1.12 igmkey = (act1 + 1) + act2*max1
118 heimbach 1.10 & + act3*max1*max2
119     & + act4*max1*max2*max3
120     #endif /* ALLOW_AUTODIFF_TAMC */
121    
122 jmc 1.22 #ifdef ALLOW_DIAGNOSTICS
123     doDiagRediFlx = .FALSE.
124     IF ( useDiagnostics ) THEN
125     doDiagRediFlx = DIAGNOSTICS_IS_ON('GM_KuzTz', myThid )
126 jmc 1.26 doDiagRediFlx = doDiagRediFlx .OR.
127 jmc 1.22 & DIAGNOSTICS_IS_ON('GM_KvzTz', myThid )
128     ENDIF
129     #endif
130 jmc 1.26
131 heimbach 1.12 #ifdef GM_VISBECK_VARIABLE_K
132     DO j=1-Oly,sNy+Oly
133     DO i=1-Olx,sNx+Olx
134     VisbeckK(i,j,bi,bj) = 0. _d 0
135     ENDDO
136     ENDDO
137     #endif
138    
139 jmc 1.19 C-- set ldd97_Lrho (for tapering scheme ldd97):
140 jmc 1.27 IF ( GM_taper_scheme.EQ.'ldd97' .OR.
141     & GM_taper_scheme.EQ.'fm07' ) THEN
142 jmc 1.19 Cspd = 2. _d 0
143     LrhoInf = 15. _d 3
144     LrhoSup = 100. _d 3
145     C- Tracer point location (center):
146     DO j=1-Oly,sNy+Oly
147     DO i=1-Olx,sNx+Olx
148     IF (fCori(i,j,bi,bj).NE.0.) THEN
149     ldd97_LrhoC(i,j) = Cspd/ABS(fCori(i,j,bi,bj))
150     ELSE
151     ldd97_LrhoC(i,j) = LrhoSup
152     ENDIF
153     ldd97_LrhoC(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoC(i,j),LrhoSup))
154     ENDDO
155     ENDDO
156     C- U point location (West):
157     DO j=1-Oly,sNy+Oly
158 jmc 1.27 kLow_W(1-Olx,j) = 0
159 jmc 1.19 ldd97_LrhoW(1-Olx,j) = LrhoSup
160     DO i=1-Olx+1,sNx+Olx
161 jmc 1.27 kLow_W(i,j) = MIN(kLowC(i-1,j,bi,bj),kLowC(i,j,bi,bj))
162 jmc 1.19 fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))
163     IF (fCoriLoc.NE.0.) THEN
164     ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)
165     ELSE
166     ldd97_LrhoW(i,j) = LrhoSup
167     ENDIF
168     ldd97_LrhoW(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoW(i,j),LrhoSup))
169     ENDDO
170     ENDDO
171     C- V point location (South):
172     DO i=1-Olx+1,sNx+Olx
173 jmc 1.27 kLow_S(i,1-Oly) = 0
174 jmc 1.19 ldd97_LrhoS(i,1-Oly) = LrhoSup
175     ENDDO
176     DO j=1-Oly+1,sNy+Oly
177     DO i=1-Olx,sNx+Olx
178 jmc 1.27 kLow_S(i,j) = MIN(kLowC(i,j-1,bi,bj),kLowC(i,j,bi,bj))
179 jmc 1.19 fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))
180     IF (fCoriLoc.NE.0.) THEN
181     ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)
182     ELSE
183     ldd97_LrhoS(i,j) = LrhoSup
184     ENDIF
185     ldd97_LrhoS(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoS(i,j),LrhoSup))
186     ENDDO
187     ENDDO
188     ELSE
189     C- Just initialize to zero (not use anyway)
190     DO j=1-Oly,sNy+Oly
191     DO i=1-Olx,sNx+Olx
192     ldd97_LrhoC(i,j) = 0. _d 0
193     ldd97_LrhoW(i,j) = 0. _d 0
194     ldd97_LrhoS(i,j) = 0. _d 0
195     ENDDO
196     ENDDO
197     ENDIF
198 jmc 1.27
199 jmc 1.19 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
200 jmc 1.27 C-- 1rst loop on k : compute Tensor Coeff. at W points.
201    
202 heimbach 1.28 DO j=1-Oly,sNy+Oly
203     DO i=1-Olx,sNx+Olx
204     hTransLay(i,j) = R_low(i,j,bi,bj)
205     baseSlope(i,j) = 0. _d 0
206     recipLambda(i,j) = 0. _d 0
207     locMixLayer(i,j) = 0. _d 0
208     ENDDO
209     ENDDO
210 jmc 1.27 #ifdef ALLOW_KPP
211     IF ( useKPP ) THEN
212     DO j=1-Oly,sNy+Oly
213     DO i=1-Olx,sNx+Olx
214     locMixLayer(i,j) = KPPhbl(i,j,bi,bj)
215     ENDDO
216     ENDDO
217     ELSE
218     #else
219     IF ( .TRUE. ) THEN
220     #endif
221     DO j=1-Oly,sNy+Oly
222     DO i=1-Olx,sNx+Olx
223     locMixLayer(i,j) = hMixLayer(i,j,bi,bj)
224     ENDDO
225     ENDDO
226     ENDIF
227 jmc 1.19
228 jmc 1.27 DO k=Nr,2,-1
229 adcroft 1.1
230     #ifdef ALLOW_AUTODIFF_TAMC
231 heimbach 1.12 kkey = (igmkey-1)*Nr + k
232 heimbach 1.10 DO j=1-Oly,sNy+Oly
233     DO i=1-Olx,sNx+Olx
234     SlopeX(i,j) = 0. _d 0
235     SlopeY(i,j) = 0. _d 0
236 heimbach 1.12 dSigmaDx(i,j) = 0. _d 0
237     dSigmaDy(i,j) = 0. _d 0
238 jmc 1.19 dSigmaDr(i,j) = 0. _d 0
239 heimbach 1.10 SlopeSqr(i,j) = 0. _d 0
240     taperFct(i,j) = 0. _d 0
241     Kwx(i,j,k,bi,bj) = 0. _d 0
242     Kwy(i,j,k,bi,bj) = 0. _d 0
243     Kwz(i,j,k,bi,bj) = 0. _d 0
244 heimbach 1.12 # ifdef GM_NON_UNITY_DIAGONAL
245     Kux(i,j,k,bi,bj) = 0. _d 0
246     Kvy(i,j,k,bi,bj) = 0. _d 0
247     # endif
248     # ifdef GM_EXTRA_DIAGONAL
249     Kuz(i,j,k,bi,bj) = 0. _d 0
250     Kvz(i,j,k,bi,bj) = 0. _d 0
251     # endif
252     # ifdef GM_BOLUS_ADVEC
253     GM_PsiX(i,j,k,bi,bj) = 0. _d 0
254     GM_PsiY(i,j,k,bi,bj) = 0. _d 0
255     # endif
256 heimbach 1.10 ENDDO
257     ENDDO
258 adcroft 1.1 #endif
259 heimbach 1.10
260 jmc 1.26 DO j=1-Oly+1,sNy+Oly-1
261     DO i=1-Olx+1,sNx+Olx-1
262 adcroft 1.1 C Gradient of Sigma at rVel points
263 jmc 1.26 dSigmaDx(i,j)=op25*( sigmaX(i+1,j,k-1)+sigmaX(i,j,k-1)
264     & +sigmaX(i+1,j, k )+sigmaX(i,j, k )
265     & )*maskC(i,j,k,bi,bj)
266     dSigmaDy(i,j)=op25*( sigmaY(i,j+1,k-1)+sigmaY(i,j,k-1)
267     & +sigmaY(i,j+1, k )+sigmaY(i,j, k )
268     & )*maskC(i,j,k,bi,bj)
269     dSigmaDr(i,j)=sigmaR(i,j,k)
270     ENDDO
271 adcroft 1.1 ENDDO
272    
273 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
274 heimbach 1.12 CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
275     CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
276 jmc 1.19 CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
277 heimbach 1.28 CADJ STORE baseSlope(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
278     CADJ STORE hTransLay(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
279     CADJ STORE recipLambda(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
280 heimbach 1.10 #endif /* ALLOW_AUTODIFF_TAMC */
281    
282 jmc 1.26 #ifdef GM_VISBECK_VARIABLE_K
283     #ifndef OLD_VISBECK_CALC
284     IF ( GM_Visbeck_alpha.GT.0. .AND.
285     & -rC(k-1).LT.GM_Visbeck_depth ) THEN
286    
287     C-- Depth average of f/sqrt(Ri) = M^2/N^2 * N
288     C M^2 and N^2 are horizontal & vertical gradient of buoyancy.
289    
290     C Calculate terms for mean Richardson number which is used
291     C in the "variable K" parameterisaton:
292     C compute depth average from surface down to the bottom or
293     C GM_Visbeck_depth, whatever is the shallower.
294    
295     DO j=1-Oly+1,sNy+Oly-1
296     DO i=1-Olx+1,sNx+Olx-1
297     IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
298     integrDepth = -rC( kLowC(i,j,bi,bj) )
299     C- in 2 steps to avoid mix of RS & RL type in min fct. arguments
300     integrDepth = MIN( integrDepth, GM_Visbeck_depth )
301     C Distance between level center above and the integration depth
302     deltaH = integrDepth + rC(k-1)
303     C If negative then we are below the integration level
304     C (cannot be the case with 2 conditions on maskC & -rC(k-1))
305     C If positive we limit this to the distance from center above
306     deltaH = MIN( deltaH, drC(k) )
307     C Now we convert deltaH to a non-dimensional fraction
308     deltaH = deltaH/( integrDepth+rC(1) )
309    
310 jmc 1.33 C-- compute: ( M^2 * S )^1/2 (= S*N since S=M^2/N^2 )
311 jmc 1.26 dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)
312     & + dSigmaDy(i,j)*dSigmaDy(i,j)
313     IF ( dSigmaH .GT. 0. _d 0 ) THEN
314     dSigmaH = SQRT( dSigmaH )
315     C- compute slope, limited by GM_maxSlope:
316     IF ( -dSigmaDr(i,j).GT.dSigmaH*GM_rMaxSlope ) THEN
317     Sloc = dSigmaH / ( -dSigmaDr(i,j) )
318     ELSE
319     Sloc = GM_maxSlope
320     ENDIF
321 jmc 1.33 M2loc = gravity*recip_rhoConst*dSigmaH
322     c SNloc = SQRT( Sloc*M2loc )
323     N2loc = -gravity*recip_rhoConst*dSigmaDr(i,j)
324     IF ( N2loc.GT.0. _d 0 ) THEN
325     SNloc = Sloc*SQRT(N2loc)
326     ELSE
327     SNloc = 0. _d 0
328     ENDIF
329 jmc 1.26 ELSE
330     SNloc = 0. _d 0
331     ENDIF
332     VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
333     & +deltaH*GM_Visbeck_alpha
334     & *GM_Visbeck_length*GM_Visbeck_length*SNloc
335     ENDIF
336     ENDDO
337     ENDDO
338     ENDIF
339     #endif /* ndef OLD_VISBECK_CALC */
340     #endif /* GM_VISBECK_VARIABLE_K */
341    
342 adcroft 1.1 C Calculate slopes for use in tensor, taper and/or clip
343 jmc 1.26 CALL GMREDI_SLOPE_LIMIT(
344 jmc 1.19 O SlopeX, SlopeY,
345 jmc 1.8 O SlopeSqr, taperFct,
346 jmc 1.27 U hTransLay, baseSlope, recipLambda,
347 jmc 1.19 U dSigmaDr,
348     I dSigmaDx, dSigmaDy,
349 jmc 1.27 I ldd97_LrhoC, locMixLayer, rF,
350     I kLowC(1-Olx,1-Oly,bi,bj),
351     I k, bi, bj, myTime, myIter, myThid )
352 adcroft 1.1
353 jmc 1.26 DO j=1-Oly+1,sNy+Oly-1
354     DO i=1-Olx+1,sNx+Olx-1
355     C Mask Iso-neutral slopes
356     SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)
357     SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)
358     SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)
359     ENDDO
360 heimbach 1.10 ENDDO
361    
362     #ifdef ALLOW_AUTODIFF_TAMC
363 heimbach 1.16 CADJ STORE SlopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
364     CADJ STORE SlopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
365 heimbach 1.14 CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
366 jmc 1.19 CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
367     CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
368 heimbach 1.10 #endif /* ALLOW_AUTODIFF_TAMC */
369    
370 jmc 1.27 C Components of Redi/GM tensor
371 jmc 1.26 DO j=1-Oly+1,sNy+Oly-1
372     DO i=1-Olx+1,sNx+Olx-1
373 jmc 1.27 Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
374     Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
375     Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
376     ENDDO
377     ENDDO
378 adcroft 1.1
379     #ifdef GM_VISBECK_VARIABLE_K
380 jmc 1.26 #ifdef OLD_VISBECK_CALC
381 jmc 1.27 DO j=1-Oly+1,sNy+Oly-1
382     DO i=1-Olx+1,sNx+Olx-1
383 jmc 1.8
384 jmc 1.27 C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K
385     C but do not know if *taperFct (or **2 ?) is necessary
386 heimbach 1.10 Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)
387 jmc 1.8
388 adcroft 1.1 C-- Depth average of M^2/N^2 * N
389    
390     C Calculate terms for mean Richardson number
391     C which is used in the "variable K" parameterisaton.
392     C Distance between interface above layer and the integration depth
393     deltaH=abs(GM_Visbeck_depth)-abs(rF(k))
394     C If positive we limit this to the layer thickness
395 jmc 1.33 integrDepth = drF(k)
396     deltaH=min(deltaH,integrDepth)
397 adcroft 1.1 C If negative then we are below the integration level
398 jmc 1.33 deltaH=max(deltaH, 0. _d 0)
399 adcroft 1.1 C Now we convert deltaH to a non-dimensional fraction
400     deltaH=deltaH/GM_Visbeck_depth
401    
402 jmc 1.19 IF ( Ssq(i,j).NE.0. .AND. dSigmaDr(i,j).NE.0. ) THEN
403 jmc 1.33 N2loc = -gravity*recip_rhoConst*dSigmaDr(i,j)
404     SNloc = SQRT(Ssq(i,j)*N2loc )
405     VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
406     & +deltaH*GM_Visbeck_alpha
407     & *GM_Visbeck_length*GM_Visbeck_length*SNloc
408 jmc 1.8 ENDIF
409 adcroft 1.1
410 jmc 1.27 ENDDO
411     ENDDO
412 jmc 1.26 #endif /* OLD_VISBECK_CALC */
413 jmc 1.9 #endif /* GM_VISBECK_VARIABLE_K */
414    
415     C-- end 1rst loop on vertical level index k
416     ENDDO
417    
418 adcroft 1.1
419 jmc 1.9 #ifdef GM_VISBECK_VARIABLE_K
420 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
421     CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
422     #endif
423 jmc 1.27 IF ( GM_Visbeck_alpha.GT.0. ) THEN
424 jmc 1.9 C- Limit range that KapGM can take
425     DO j=1-Oly+1,sNy+Oly-1
426     DO i=1-Olx+1,sNx+Olx-1
427     VisbeckK(i,j,bi,bj)=
428     & MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)
429     ENDDO
430     ENDDO
431     ENDIF
432 heimbach 1.16 cph( NEW
433     #ifdef ALLOW_AUTODIFF_TAMC
434     CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
435     #endif
436     cph)
437 adcroft 1.1 #endif /* GM_VISBECK_VARIABLE_K */
438    
439 heimbach 1.29 C- express the Tensor in term of Diffusivity (= m**2 / s )
440     DO k=1,Nr
441 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
442     kkey = (igmkey-1)*Nr + k
443 heimbach 1.29 # if (defined (GM_NON_UNITY_DIAGONAL) || \
444 heimbach 1.16 defined (GM_VISBECK_VARIABLE_K))
445 heimbach 1.14 CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
446     CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
447     CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
448 heimbach 1.29 # endif
449 heimbach 1.12 #endif
450 jmc 1.26 DO j=1-Oly+1,sNy+Oly-1
451     DO i=1-Olx+1,sNx+Olx-1
452 gforget 1.31 #ifdef ALLOW_KAPREDI_CONTROL
453     Kgm_tmp = kapredi(i,j,k,bi,bj)
454     #else
455     Kgm_tmp = GM_isopycK
456     #endif
457 jmc 1.33 #ifdef ALLOW_KAPGM_CONTROL
458 gforget 1.31 & + GM_skewflx*kapgm(i,j,k,bi,bj)
459 heimbach 1.25 #else
460 gforget 1.31 & + GM_skewflx*GM_background_K
461 heimbach 1.25 #endif
462 jmc 1.9 #ifdef GM_VISBECK_VARIABLE_K
463 jmc 1.26 & + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)
464 jmc 1.9 #endif
465 jmc 1.26 Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)
466     Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)
467 gforget 1.31 #ifdef ALLOW_KAPREDI_CONTROL
468 jmc 1.33 Kwz(i,j,k,bi,bj)= ( kapredi(i,j,k,bi,bj)
469 gforget 1.31 #else
470 jmc 1.26 Kwz(i,j,k,bi,bj)= ( GM_isopycK
471 gforget 1.31 #endif
472 adcroft 1.1 #ifdef GM_VISBECK_VARIABLE_K
473 jmc 1.26 & + VisbeckK(i,j,bi,bj)
474 adcroft 1.1 #endif
475 jmc 1.26 & )*Kwz(i,j,k,bi,bj)
476     ENDDO
477 adcroft 1.1 ENDDO
478 jmc 1.27 ENDDO
479    
480     #ifdef ALLOW_DIAGNOSTICS
481     IF ( useDiagnostics .AND. GM_taper_scheme.EQ.'fm07' ) THEN
482     CALL DIAGNOSTICS_FILL( hTransLay, 'GM_hTrsL', 0,1,2,bi,bj,myThid)
483     CALL DIAGNOSTICS_FILL( baseSlope, 'GM_baseS', 0,1,2,bi,bj,myThid)
484     CALL DIAGNOSTICS_FILL(recipLambda,'GM_rLamb', 0,1,2,bi,bj,myThid)
485     ENDIF
486     #endif /* ALLOW_DIAGNOSTICS */
487    
488 adcroft 1.4
489 jmc 1.9 #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
490 jmc 1.27 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
491     C-- 2nd k loop : compute Tensor Coeff. at U point
492    
493     #ifdef ALLOW_KPP
494     IF ( useKPP ) THEN
495     DO j=1-Oly,sNy+Oly
496     DO i=2-Olx,sNx+Olx
497     locMixLayer(i,j) = ( KPPhbl(i-1,j,bi,bj)
498     & + KPPhbl( i ,j,bi,bj) )*op5
499     ENDDO
500     ENDDO
501     ELSE
502     #else
503     IF ( .TRUE. ) THEN
504     #endif
505     DO j=1-Oly,sNy+Oly
506     DO i=2-Olx,sNx+Olx
507     locMixLayer(i,j) = ( hMixLayer(i-1,j,bi,bj)
508     & + hMixLayer( i ,j,bi,bj) )*op5
509     ENDDO
510     ENDDO
511     ENDIF
512     DO j=1-Oly,sNy+Oly
513     DO i=1-Olx,sNx+Olx
514     hTransLay(i,j) = 0.
515     baseSlope(i,j) = 0.
516     recipLambda(i,j)= 0.
517     ENDDO
518     DO i=2-Olx,sNx+Olx
519     hTransLay(i,j) = MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
520     ENDDO
521     ENDDO
522    
523     DO k=Nr,1,-1
524     kp1 = MIN(Nr,k+1)
525     maskp1 = 1. _d 0
526     IF (k.GE.Nr) maskp1 = 0. _d 0
527 heimbach 1.29 #ifdef ALLOW_AUTODIFF_TAMC
528     kkey = (igmkey-1)*Nr + k
529     #endif
530 adcroft 1.1
531     C Gradient of Sigma at U points
532 jmc 1.26 DO j=1-Oly+1,sNy+Oly-1
533     DO i=1-Olx+1,sNx+Olx-1
534     dSigmaDx(i,j)=sigmaX(i,j,k)
535     & *_maskW(i,j,k,bi,bj)
536     dSigmaDy(i,j)=op25*( sigmaY(i-1,j+1,k)+sigmaY(i,j+1,k)
537     & +sigmaY(i-1, j ,k)+sigmaY(i, j ,k)
538     & )*_maskW(i,j,k,bi,bj)
539     dSigmaDr(i,j)=op25*( sigmaR(i-1,j, k )+sigmaR(i,j, k )
540     & +(sigmaR(i-1,j,kp1)+sigmaR(i,j,kp1))*maskp1
541     & )*_maskW(i,j,k,bi,bj)
542     ENDDO
543 adcroft 1.1 ENDDO
544    
545 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
546 heimbach 1.17 CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
547 heimbach 1.12 CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
548     CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
549 heimbach 1.28 CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
550     CADJ STORE locMixLayer(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
551     CADJ STORE baseSlope(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
552     CADJ STORE hTransLay(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
553     CADJ STORE recipLambda(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
554 heimbach 1.12 #endif /* ALLOW_AUTODIFF_TAMC */
555    
556 adcroft 1.1 C Calculate slopes for use in tensor, taper and/or clip
557 jmc 1.26 CALL GMREDI_SLOPE_LIMIT(
558 jmc 1.19 O SlopeX, SlopeY,
559 jmc 1.8 O SlopeSqr, taperFct,
560 jmc 1.27 U hTransLay, baseSlope, recipLambda,
561 jmc 1.19 U dSigmaDr,
562     I dSigmaDx, dSigmaDy,
563 jmc 1.27 I ldd97_LrhoW, locMixLayer, rC,
564     I kLow_W,
565     I k, bi, bj, myTime, myIter, myThid )
566 adcroft 1.1
567 heimbach 1.16 #ifdef ALLOW_AUTODIFF_TAMC
568 heimbach 1.17 CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
569 jmc 1.19 CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
570 heimbach 1.16 #endif /* ALLOW_AUTODIFF_TAMC */
571    
572 jmc 1.9 #ifdef GM_NON_UNITY_DIAGONAL
573 jmc 1.26 c IF ( GM_nonUnitDiag ) THEN
574 jmc 1.9 DO j=1-Oly+1,sNy+Oly-1
575     DO i=1-Olx+1,sNx+Olx-1
576     Kux(i,j,k,bi,bj) =
577 gforget 1.31 #ifdef ALLOW_KAPREDI_CONTROL
578     & ( kapredi(i,j,k,bi,bj)
579     #else
580 jmc 1.9 & ( GM_isopycK
581 gforget 1.31 #endif
582 jmc 1.9 #ifdef GM_VISBECK_VARIABLE_K
583 heimbach 1.14 & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
584 jmc 1.9 #endif
585 jmc 1.27 & )*taperFct(i,j)
586 heimbach 1.10 ENDDO
587     ENDDO
588 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
589 heimbach 1.16 # ifdef GM_EXCLUDE_CLIPPING
590 heimbach 1.12 CADJ STORE Kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
591     # endif
592     #endif
593 heimbach 1.10 DO j=1-Oly+1,sNy+Oly-1
594     DO i=1-Olx+1,sNx+Olx-1
595 jmc 1.9 Kux(i,j,k,bi,bj) = MAX( Kux(i,j,k,bi,bj), GM_Kmin_horiz )
596     ENDDO
597     ENDDO
598 jmc 1.26 c ENDIF
599 jmc 1.9 #endif /* GM_NON_UNITY_DIAGONAL */
600    
601     #ifdef GM_EXTRA_DIAGONAL
602 heimbach 1.12
603     #ifdef ALLOW_AUTODIFF_TAMC
604     CADJ STORE SlopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
605     CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
606     #endif /* ALLOW_AUTODIFF_TAMC */
607 jmc 1.27 IF ( GM_ExtraDiag ) THEN
608 jmc 1.9 DO j=1-Oly+1,sNy+Oly-1
609     DO i=1-Olx+1,sNx+Olx-1
610     Kuz(i,j,k,bi,bj) =
611 gforget 1.31 #ifdef ALLOW_KAPREDI_CONTROL
612     & ( kapredi(i,j,k,bi,bj)
613     #else
614     & ( GM_isopycK
615     #endif
616 jmc 1.33 #ifdef ALLOW_KAPGM_CONTROL
617 gforget 1.31 & - GM_skewflx*kapgm(i,j,k,bi,bj)
618 heimbach 1.25 #else
619 gforget 1.31 & - GM_skewflx*GM_background_K
620 heimbach 1.25 #endif
621 jmc 1.9 #ifdef GM_VISBECK_VARIABLE_K
622 heimbach 1.14 & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*GM_advect
623 jmc 1.9 #endif
624     & )*SlopeX(i,j)*taperFct(i,j)
625     ENDDO
626     ENDDO
627 jmc 1.26 ENDIF
628 jmc 1.9 #endif /* GM_EXTRA_DIAGONAL */
629 adcroft 1.1
630 jmc 1.22 #ifdef ALLOW_DIAGNOSTICS
631 jmc 1.26 IF (doDiagRediFlx) THEN
632 jmc 1.22 km1 = MAX(k-1,1)
633     DO j=1,sNy
634     DO i=1,sNx+1
635     C store in tmp1k Kuz_Redi
636 gforget 1.31 #ifdef ALLOW_KAPREDI_CONTROL
637     tmp1k(i,j) = ( kapredi(i,j,k,bi,bj)
638     #else
639 jmc 1.22 tmp1k(i,j) = ( GM_isopycK
640 gforget 1.31 #endif
641 jmc 1.22 #ifdef GM_VISBECK_VARIABLE_K
642     & +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0
643     #endif
644     & )*SlopeX(i,j)*taperFct(i,j)
645     ENDDO
646     ENDDO
647     DO j=1,sNy
648     DO i=1,sNx+1
649     C- Vertical gradients interpolated to U points
650     dTdz = (
651     & +recip_drC(k)*
652     & ( maskC(i-1,j,k,bi,bj)*
653     & (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
654     & +maskC( i ,j,k,bi,bj)*
655     & (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
656     & )
657     & +recip_drC(kp1)*
658     & ( maskC(i-1,j,kp1,bi,bj)*
659     & (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
660     & +maskC( i ,j,kp1,bi,bj)*
661     & (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
662     & ) ) * 0.25 _d 0
663 heimbach 1.23 tmp1k(i,j) = dyG(i,j,bi,bj)*drF(k)
664     & * _hFacW(i,j,k,bi,bj)
665 jmc 1.22 & * tmp1k(i,j) * dTdz
666     ENDDO
667     ENDDO
668     CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KuzTz', k,1,2,bi,bj,myThid)
669 jmc 1.26 ENDIF
670 jmc 1.22 #endif /* ALLOW_DIAGNOSTICS */
671    
672 jmc 1.27 C-- end 2nd loop on vertical level index k
673     ENDDO
674    
675     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
676     C-- 3rd k loop : compute Tensor Coeff. at V point
677    
678     #ifdef ALLOW_KPP
679     IF ( useKPP ) THEN
680     DO j=2-Oly,sNy+Oly
681     DO i=1-Olx,sNx+Olx
682     locMixLayer(i,j) = ( KPPhbl(i,j-1,bi,bj)
683     & + KPPhbl(i, j ,bi,bj) )*op5
684     ENDDO
685     ENDDO
686     ELSE
687     #else
688     IF ( .TRUE. ) THEN
689     #endif
690     DO j=2-Oly,sNy+Oly
691     DO i=1-Olx,sNx+Olx
692     locMixLayer(i,j) = ( hMixLayer(i,j-1,bi,bj)
693     & + hMixLayer(i, j ,bi,bj) )*op5
694     ENDDO
695     ENDDO
696     ENDIF
697     DO j=1-Oly,sNy+Oly
698     DO i=1-Olx,sNx+Olx
699     hTransLay(i,j) = 0.
700     baseSlope(i,j) = 0.
701     recipLambda(i,j)= 0.
702     ENDDO
703     ENDDO
704     DO j=2-Oly,sNy+Oly
705     DO i=1-Olx,sNx+Olx
706     hTransLay(i,j) = MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
707     ENDDO
708     ENDDO
709    
710 adcroft 1.1 C Gradient of Sigma at V points
711 jmc 1.27 DO k=Nr,1,-1
712     kp1 = MIN(Nr,k+1)
713     maskp1 = 1. _d 0
714     IF (k.GE.Nr) maskp1 = 0. _d 0
715 heimbach 1.29 #ifdef ALLOW_AUTODIFF_TAMC
716     kkey = (igmkey-1)*Nr + k
717     #endif
718 jmc 1.27
719 jmc 1.26 DO j=1-Oly+1,sNy+Oly-1
720     DO i=1-Olx+1,sNx+Olx-1
721     dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
722     & +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k)
723     & )*_maskS(i,j,k,bi,bj)
724     dSigmaDy(i,j)=sigmaY(i,j,k)
725     & *_maskS(i,j,k,bi,bj)
726     dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k )+sigmaR(i,j, k )
727     & +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1
728     & )*_maskS(i,j,k,bi,bj)
729     ENDDO
730 adcroft 1.1 ENDDO
731    
732 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
733     CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
734     CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
735 heimbach 1.28 CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
736     CADJ STORE baseSlope(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
737     CADJ STORE hTransLay(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
738     CADJ STORE recipLambda(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
739 heimbach 1.12 #endif /* ALLOW_AUTODIFF_TAMC */
740    
741 adcroft 1.1 C Calculate slopes for use in tensor, taper and/or clip
742 jmc 1.26 CALL GMREDI_SLOPE_LIMIT(
743 jmc 1.19 O SlopeX, SlopeY,
744 jmc 1.8 O SlopeSqr, taperFct,
745 jmc 1.27 U hTransLay, baseSlope, recipLambda,
746 jmc 1.19 U dSigmaDr,
747     I dSigmaDx, dSigmaDy,
748 jmc 1.27 I ldd97_LrhoS, locMixLayer, rC,
749     I kLow_S,
750     I k, bi, bj, myTime, myIter, myThid )
751 adcroft 1.1
752 heimbach 1.16 cph(
753     #ifdef ALLOW_AUTODIFF_TAMC
754     cph(
755     CADJ STORE taperfct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
756     cph)
757     #endif /* ALLOW_AUTODIFF_TAMC */
758     cph)
759    
760 jmc 1.9 #ifdef GM_NON_UNITY_DIAGONAL
761 jmc 1.26 c IF ( GM_nonUnitDiag ) THEN
762 jmc 1.9 DO j=1-Oly+1,sNy+Oly-1
763     DO i=1-Olx+1,sNx+Olx-1
764     Kvy(i,j,k,bi,bj) =
765 gforget 1.31 #ifdef ALLOW_KAPREDI_CONTROL
766     & ( kapredi(i,j,k,bi,bj)
767     #else
768 jmc 1.9 & ( GM_isopycK
769 gforget 1.31 #endif
770 jmc 1.9 #ifdef GM_VISBECK_VARIABLE_K
771 heimbach 1.14 & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
772 jmc 1.9 #endif
773 jmc 1.27 & )*taperFct(i,j)
774 heimbach 1.10 ENDDO
775     ENDDO
776 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
777 heimbach 1.16 # ifdef GM_EXCLUDE_CLIPPING
778 heimbach 1.12 CADJ STORE Kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
779     # endif
780     #endif
781 heimbach 1.10 DO j=1-Oly+1,sNy+Oly-1
782     DO i=1-Olx+1,sNx+Olx-1
783 jmc 1.9 Kvy(i,j,k,bi,bj) = MAX( Kvy(i,j,k,bi,bj), GM_Kmin_horiz )
784     ENDDO
785     ENDDO
786 jmc 1.26 c ENDIF
787 jmc 1.9 #endif /* GM_NON_UNITY_DIAGONAL */
788    
789     #ifdef GM_EXTRA_DIAGONAL
790 heimbach 1.12
791     #ifdef ALLOW_AUTODIFF_TAMC
792     CADJ STORE SlopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
793     CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
794     #endif /* ALLOW_AUTODIFF_TAMC */
795 jmc 1.27 IF ( GM_ExtraDiag ) THEN
796 jmc 1.9 DO j=1-Oly+1,sNy+Oly-1
797     DO i=1-Olx+1,sNx+Olx-1
798     Kvz(i,j,k,bi,bj) =
799 gforget 1.31 #ifdef ALLOW_KAPREDI_CONTROL
800     & ( kapredi(i,j,k,bi,bj)
801     #else
802 jmc 1.33 & ( GM_isopycK
803 gforget 1.31 #endif
804 jmc 1.33 #ifdef ALLOW_KAPGM_CONTROL
805 gforget 1.31 & - GM_skewflx*kapgm(i,j,k,bi,bj)
806 heimbach 1.25 #else
807 gforget 1.31 & - GM_skewflx*GM_background_K
808 heimbach 1.25 #endif
809 jmc 1.9 #ifdef GM_VISBECK_VARIABLE_K
810 heimbach 1.14 & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*GM_advect
811 jmc 1.9 #endif
812     & )*SlopeY(i,j)*taperFct(i,j)
813     ENDDO
814     ENDDO
815 jmc 1.26 ENDIF
816 jmc 1.9 #endif /* GM_EXTRA_DIAGONAL */
817    
818 jmc 1.22 #ifdef ALLOW_DIAGNOSTICS
819 jmc 1.26 IF (doDiagRediFlx) THEN
820 dfer 1.30 km1 = MAX(k-1,1)
821 jmc 1.22 DO j=1,sNy+1
822     DO i=1,sNx
823     C store in tmp1k Kvz_Redi
824 gforget 1.31 #ifdef ALLOW_KAPREDI_CONTROL
825     tmp1k(i,j) = ( kapredi(i,j,k,bi,bj)
826     #else
827 jmc 1.22 tmp1k(i,j) = ( GM_isopycK
828 gforget 1.31 #endif
829 jmc 1.22 #ifdef GM_VISBECK_VARIABLE_K
830     & +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0
831     #endif
832     & )*SlopeY(i,j)*taperFct(i,j)
833     ENDDO
834     ENDDO
835     DO j=1,sNy+1
836     DO i=1,sNx
837     C- Vertical gradients interpolated to U points
838     dTdz = (
839     & +recip_drC(k)*
840     & ( maskC(i,j-1,k,bi,bj)*
841     & (theta(i,j-1,km1,bi,bj)-theta(i,j-1,k,bi,bj))
842     & +maskC(i, j ,k,bi,bj)*
843     & (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))
844     & )
845     & +recip_drC(kp1)*
846     & ( maskC(i,j-1,kp1,bi,bj)*
847     & (theta(i,j-1,k,bi,bj)-theta(i,j-1,kp1,bi,bj))
848     & +maskC(i, j ,kp1,bi,bj)*
849     & (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))
850     & ) ) * 0.25 _d 0
851 heimbach 1.23 tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)
852     & * _hFacS(i,j,k,bi,bj)
853 jmc 1.22 & * tmp1k(i,j) * dTdz
854     ENDDO
855     ENDDO
856     CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KvzTz', k,1,2,bi,bj,myThid)
857 jmc 1.26 ENDIF
858 jmc 1.22 #endif /* ALLOW_DIAGNOSTICS */
859    
860 jmc 1.27 C-- end 3rd loop on vertical level index k
861     ENDDO
862    
863 jmc 1.9 #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
864    
865 adcroft 1.1
866 jmc 1.9 #ifdef GM_BOLUS_ADVEC
867     IF (GM_AdvForm) THEN
868 jmc 1.26 CALL GMREDI_CALC_PSI_B(
869 jmc 1.9 I bi, bj, iMin, iMax, jMin, jMax,
870     I sigmaX, sigmaY, sigmaR,
871 jmc 1.19 I ldd97_LrhoW, ldd97_LrhoS,
872 jmc 1.26 I myThid )
873 jmc 1.9 ENDIF
874     #endif
875 adcroft 1.1
876 jmc 1.19 #ifdef ALLOW_TIMEAVE
877     C-- Time-average
878     IF ( taveFreq.GT.0. ) THEN
879    
880     CALL TIMEAVE_CUMULATE( GM_Kwx_T, Kwx, Nr,
881     & deltaTclock, bi, bj, myThid )
882     CALL TIMEAVE_CUMULATE( GM_Kwy_T, Kwy, Nr,
883     & deltaTclock, bi, bj, myThid )
884     CALL TIMEAVE_CUMULATE( GM_Kwz_T, Kwz, Nr,
885     & deltaTclock, bi, bj, myThid )
886     #ifdef GM_VISBECK_VARIABLE_K
887     IF ( GM_Visbeck_alpha.NE.0. ) THEN
888     CALL TIMEAVE_CUMULATE( Visbeck_K_T, VisbeckK, 1,
889     & deltaTclock, bi, bj, myThid )
890     ENDIF
891     #endif
892     #ifdef GM_BOLUS_ADVEC
893     IF ( GM_AdvForm ) THEN
894     CALL TIMEAVE_CUMULATE( GM_PsiXtave, GM_PsiX, Nr,
895     & deltaTclock, bi, bj, myThid )
896     CALL TIMEAVE_CUMULATE( GM_PsiYtave, GM_PsiY, Nr,
897     & deltaTclock, bi, bj, myThid )
898     ENDIF
899     #endif
900     DO k=1,Nr
901     GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock
902     ENDDO
903    
904     ENDIF
905     #endif /* ALLOW_TIMEAVE */
906    
907 jmc 1.20 #ifdef ALLOW_DIAGNOSTICS
908     IF ( useDiagnostics ) THEN
909 jmc 1.24 CALL GMREDI_DIAGNOSTICS_FILL(bi,bj,myThid)
910 jmc 1.20 ENDIF
911     #endif /* ALLOW_DIAGNOSTICS */
912    
913 adcroft 1.1 #endif /* ALLOW_GMREDI */
914    
915     RETURN
916     END
917 heimbach 1.2
918    
919     SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
920 jmc 1.27 I iMin, iMax, jMin, jMax,
921 heimbach 1.2 I sigmaX, sigmaY, sigmaR,
922 jmc 1.27 I bi, bj, myTime, myIter, myThid )
923 heimbach 1.2 C /==========================================================\
924     C | SUBROUTINE GMREDI_CALC_TENSOR |
925     C | o Calculate tensor elements for GM/Redi tensor. |
926     C |==========================================================|
927     C \==========================================================/
928     IMPLICIT NONE
929    
930     C == Global variables ==
931     #include "SIZE.h"
932     #include "EEPARAMS.h"
933     #include "GMREDI.h"
934    
935     C == Routine arguments ==
936     C
937     _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
938     _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
939     _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
940 jmc 1.27 INTEGER iMin,iMax,jMin,jMax
941     INTEGER bi, bj
942     _RL myTime
943     INTEGER myIter
944 heimbach 1.2 INTEGER myThid
945     CEndOfInterface
946    
947 jmc 1.27 #ifdef ALLOW_GMREDI
948    
949 jmc 1.9 INTEGER i, j, k
950 heimbach 1.2
951 jmc 1.9 DO k=1,Nr
952     DO j=1-Oly+1,sNy+Oly-1
953     DO i=1-Olx+1,sNx+Olx-1
954     Kwx(i,j,k,bi,bj) = 0.0
955     Kwy(i,j,k,bi,bj) = 0.0
956     Kwz(i,j,k,bi,bj) = 0.0
957     ENDDO
958 heimbach 1.2 ENDDO
959     ENDDO
960     #endif /* ALLOW_GMREDI */
961    
962 jmc 1.9 RETURN
963     END
  ViewVC Help
Powered by ViewVC 1.1.22