/[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.26 - (hide annotations) (download)
Thu May 24 22:34:38 2007 UTC (17 years ago) by jmc
Branch: MAIN
CVS Tags: checkpoint59d, checkpoint59c
Changes since 1.25: +140 -77 lines 
change calculation of Visbeck K:
 - no longer depend on tapering scheme ; instead, compute directly the slope
   (limited by GM_maxSlope).
 - fix vertical averaging (wrong by 1 level).
1 jmc 1.26 C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_calc_tensor.F,v 1.25 2007/02/07 00:01:15 heimbach Exp $
2 heimbach 1.13 C $Name: $
3 adcroft 1.1
4     #include "GMREDI_OPTIONS.h"
5 jmc 1.26 #undef OLD_VISBECK_CALC
6 adcroft 1.1
7     CStartOfInterface
8     SUBROUTINE GMREDI_CALC_TENSOR(
9 jmc 1.9 I bi, bj, iMin, iMax, jMin, jMax,
10 adcroft 1.1 I sigmaX, sigmaY, sigmaR,
11     I myThid )
12     C /==========================================================\
13     C | SUBROUTINE GMREDI_CALC_TENSOR |
14     C | o Calculate tensor elements for GM/Redi tensor. |
15     C |==========================================================|
16     C \==========================================================/
17     IMPLICIT NONE
18    
19     C == Global variables ==
20     #include "SIZE.h"
21     #include "GRID.h"
22     #include "DYNVARS.h"
23     #include "EEPARAMS.h"
24     #include "PARAMS.h"
25     #include "GMREDI.h"
26 jmc 1.20 #include "GMREDI_TAVE.h"
27 adcroft 1.1
28 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
29     #include "tamc.h"
30     #include "tamc_keys.h"
31     #endif /* ALLOW_AUTODIFF_TAMC */
32    
33 adcroft 1.1 C == Routine arguments ==
34     C
35     _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
36     _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
37     _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
38 jmc 1.9 INTEGER bi,bj,iMin,iMax,jMin,jMax
39 adcroft 1.1 INTEGER myThid
40     CEndOfInterface
41    
42     #ifdef ALLOW_GMREDI
43    
44     C == Local variables ==
45 jmc 1.15 INTEGER i,j,k,kp1
46 adcroft 1.1 _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47     _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48 heimbach 1.12 _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49     _RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50 jmc 1.19 _RL dSigmaDr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51 jmc 1.8 _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
52     _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
53 jmc 1.15 _RL maskp1, Kgm_tmp
54 jmc 1.19 _RL ldd97_LrhoC(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
55     _RL ldd97_LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56     _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57     _RL Cspd, LrhoInf, LrhoSup, fCoriLoc
58 adcroft 1.1
59     #ifdef GM_VISBECK_VARIABLE_K
60 jmc 1.26 #ifdef OLD_VISBECK_CALC
61 heimbach 1.14 _RL deltaH,zero_rs
62     PARAMETER(zero_rs=0.D0)
63 adcroft 1.1 _RL N2,SN
64 heimbach 1.10 _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
65 jmc 1.26 #else
66     _RL dSigmaH
67     _RL deltaH, integrDepth
68     _RL Sloc, M2loc, SNloc
69     #endif
70 adcroft 1.1 #endif
71    
72 jmc 1.22 #ifdef ALLOW_DIAGNOSTICS
73     LOGICAL doDiagRediFlx
74     LOGICAL DIAGNOSTICS_IS_ON
75     EXTERNAL DIAGNOSTICS_IS_ON
76     INTEGER km1
77     _RL dTdz
78     _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79     #endif
80    
81 jmc 1.19 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
82    
83 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
84     act1 = bi - myBxLo(myThid)
85     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
86     act2 = bj - myByLo(myThid)
87     max2 = myByHi(myThid) - myByLo(myThid) + 1
88     act3 = myThid - 1
89     max3 = nTx*nTy
90     act4 = ikey_dynamics - 1
91 heimbach 1.12 igmkey = (act1 + 1) + act2*max1
92 heimbach 1.10 & + act3*max1*max2
93     & + act4*max1*max2*max3
94     #endif /* ALLOW_AUTODIFF_TAMC */
95    
96 jmc 1.22 #ifdef ALLOW_DIAGNOSTICS
97     doDiagRediFlx = .FALSE.
98     IF ( useDiagnostics ) THEN
99     doDiagRediFlx = DIAGNOSTICS_IS_ON('GM_KuzTz', myThid )
100 jmc 1.26 doDiagRediFlx = doDiagRediFlx .OR.
101 jmc 1.22 & DIAGNOSTICS_IS_ON('GM_KvzTz', myThid )
102     ENDIF
103     #endif
104 jmc 1.26
105 heimbach 1.12 #ifdef GM_VISBECK_VARIABLE_K
106     DO j=1-Oly,sNy+Oly
107     DO i=1-Olx,sNx+Olx
108     VisbeckK(i,j,bi,bj) = 0. _d 0
109     ENDDO
110     ENDDO
111     #endif
112    
113 jmc 1.19 C-- set ldd97_Lrho (for tapering scheme ldd97):
114     IF (GM_taper_scheme.EQ.'ldd97') THEN
115     Cspd = 2. _d 0
116     LrhoInf = 15. _d 3
117     LrhoSup = 100. _d 3
118     C- Tracer point location (center):
119     DO j=1-Oly,sNy+Oly
120     DO i=1-Olx,sNx+Olx
121     IF (fCori(i,j,bi,bj).NE.0.) THEN
122     ldd97_LrhoC(i,j) = Cspd/ABS(fCori(i,j,bi,bj))
123     ELSE
124     ldd97_LrhoC(i,j) = LrhoSup
125     ENDIF
126     ldd97_LrhoC(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoC(i,j),LrhoSup))
127     ENDDO
128     ENDDO
129     C- U point location (West):
130     DO j=1-Oly,sNy+Oly
131     ldd97_LrhoW(1-Olx,j) = LrhoSup
132     DO i=1-Olx+1,sNx+Olx
133     fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))
134     IF (fCoriLoc.NE.0.) THEN
135     ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)
136     ELSE
137     ldd97_LrhoW(i,j) = LrhoSup
138     ENDIF
139     ldd97_LrhoW(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoW(i,j),LrhoSup))
140     ENDDO
141     ENDDO
142     C- V point location (South):
143     DO i=1-Olx+1,sNx+Olx
144     ldd97_LrhoS(i,1-Oly) = LrhoSup
145     ENDDO
146     DO j=1-Oly+1,sNy+Oly
147     DO i=1-Olx,sNx+Olx
148     fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))
149     IF (fCoriLoc.NE.0.) THEN
150     ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)
151     ELSE
152     ldd97_LrhoS(i,j) = LrhoSup
153     ENDIF
154     ldd97_LrhoS(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoS(i,j),LrhoSup))
155     ENDDO
156     ENDDO
157     ELSE
158     C- Just initialize to zero (not use anyway)
159     DO j=1-Oly,sNy+Oly
160     DO i=1-Olx,sNx+Olx
161     ldd97_LrhoC(i,j) = 0. _d 0
162     ldd97_LrhoW(i,j) = 0. _d 0
163     ldd97_LrhoS(i,j) = 0. _d 0
164     ENDDO
165     ENDDO
166     ENDIF
167     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
168    
169 jmc 1.9 DO k=2,Nr
170     C-- 1rst loop on k : compute Tensor Coeff. at W points.
171 adcroft 1.1
172     #ifdef ALLOW_AUTODIFF_TAMC
173 heimbach 1.12 kkey = (igmkey-1)*Nr + k
174 heimbach 1.10 DO j=1-Oly,sNy+Oly
175     DO i=1-Olx,sNx+Olx
176     SlopeX(i,j) = 0. _d 0
177     SlopeY(i,j) = 0. _d 0
178 heimbach 1.12 dSigmaDx(i,j) = 0. _d 0
179     dSigmaDy(i,j) = 0. _d 0
180 jmc 1.19 dSigmaDr(i,j) = 0. _d 0
181 heimbach 1.10 SlopeSqr(i,j) = 0. _d 0
182     taperFct(i,j) = 0. _d 0
183     Kwx(i,j,k,bi,bj) = 0. _d 0
184     Kwy(i,j,k,bi,bj) = 0. _d 0
185     Kwz(i,j,k,bi,bj) = 0. _d 0
186 heimbach 1.12 # ifdef GM_NON_UNITY_DIAGONAL
187     Kux(i,j,k,bi,bj) = 0. _d 0
188     Kvy(i,j,k,bi,bj) = 0. _d 0
189     # endif
190     # ifdef GM_EXTRA_DIAGONAL
191     Kuz(i,j,k,bi,bj) = 0. _d 0
192     Kvz(i,j,k,bi,bj) = 0. _d 0
193     # endif
194     # ifdef GM_BOLUS_ADVEC
195     GM_PsiX(i,j,k,bi,bj) = 0. _d 0
196     GM_PsiY(i,j,k,bi,bj) = 0. _d 0
197     # endif
198 heimbach 1.10 ENDDO
199     ENDDO
200 adcroft 1.1 #endif
201 heimbach 1.10
202 jmc 1.26 DO j=1-Oly+1,sNy+Oly-1
203     DO i=1-Olx+1,sNx+Olx-1
204 adcroft 1.1 C Gradient of Sigma at rVel points
205 jmc 1.26 dSigmaDx(i,j)=op25*( sigmaX(i+1,j,k-1)+sigmaX(i,j,k-1)
206     & +sigmaX(i+1,j, k )+sigmaX(i,j, k )
207     & )*maskC(i,j,k,bi,bj)
208     dSigmaDy(i,j)=op25*( sigmaY(i,j+1,k-1)+sigmaY(i,j,k-1)
209     & +sigmaY(i,j+1, k )+sigmaY(i,j, k )
210     & )*maskC(i,j,k,bi,bj)
211     dSigmaDr(i,j)=sigmaR(i,j,k)
212     ENDDO
213 adcroft 1.1 ENDDO
214    
215 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
216 heimbach 1.12 CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
217     CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
218 jmc 1.19 CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
219 heimbach 1.10 #endif /* ALLOW_AUTODIFF_TAMC */
220    
221 jmc 1.26 #ifdef GM_VISBECK_VARIABLE_K
222     #ifndef OLD_VISBECK_CALC
223     IF ( GM_Visbeck_alpha.GT.0. .AND.
224     & -rC(k-1).LT.GM_Visbeck_depth ) THEN
225    
226     C-- Depth average of f/sqrt(Ri) = M^2/N^2 * N
227     C M^2 and N^2 are horizontal & vertical gradient of buoyancy.
228    
229     C Calculate terms for mean Richardson number which is used
230     C in the "variable K" parameterisaton:
231     C compute depth average from surface down to the bottom or
232     C GM_Visbeck_depth, whatever is the shallower.
233    
234     DO j=1-Oly+1,sNy+Oly-1
235     DO i=1-Olx+1,sNx+Olx-1
236     IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
237     integrDepth = -rC( kLowC(i,j,bi,bj) )
238     C- in 2 steps to avoid mix of RS & RL type in min fct. arguments
239     integrDepth = MIN( integrDepth, GM_Visbeck_depth )
240     C Distance between level center above and the integration depth
241     deltaH = integrDepth + rC(k-1)
242     C If negative then we are below the integration level
243     C (cannot be the case with 2 conditions on maskC & -rC(k-1))
244     C If positive we limit this to the distance from center above
245     deltaH = MIN( deltaH, drC(k) )
246     C Now we convert deltaH to a non-dimensional fraction
247     deltaH = deltaH/( integrDepth+rC(1) )
248    
249     C-- compute: ( M^2 * S )^1/2 (= M^2 / N since S=M^2/N^2 )
250     dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)
251     & + dSigmaDy(i,j)*dSigmaDy(i,j)
252     IF ( dSigmaH .GT. 0. _d 0 ) THEN
253     dSigmaH = SQRT( dSigmaH )
254     C- compute slope, limited by GM_maxSlope:
255     IF ( -dSigmaDr(i,j).GT.dSigmaH*GM_rMaxSlope ) THEN
256     Sloc = dSigmaH / ( -dSigmaDr(i,j) )
257     ELSE
258     Sloc = GM_maxSlope
259     ENDIF
260     M2loc = Gravity*recip_RhoConst*dSigmaH
261     SNloc = SQRT( Sloc*M2loc )
262     ELSE
263     SNloc = 0. _d 0
264     ENDIF
265     VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
266     & +deltaH*GM_Visbeck_alpha
267     & *GM_Visbeck_length*GM_Visbeck_length*SNloc
268     ENDIF
269     ENDDO
270     ENDDO
271     ENDIF
272     #endif /* ndef OLD_VISBECK_CALC */
273     #endif /* GM_VISBECK_VARIABLE_K */
274    
275 adcroft 1.1 C Calculate slopes for use in tensor, taper and/or clip
276 jmc 1.26 CALL GMREDI_SLOPE_LIMIT(
277 jmc 1.19 O SlopeX, SlopeY,
278 jmc 1.8 O SlopeSqr, taperFct,
279 jmc 1.19 U dSigmaDr,
280     I dSigmaDx, dSigmaDy,
281     I ldd97_LrhoC,rF(k),k,
282 adcroft 1.1 I bi, bj, myThid )
283    
284 jmc 1.26 DO j=1-Oly+1,sNy+Oly-1
285     DO i=1-Olx+1,sNx+Olx-1
286     C Mask Iso-neutral slopes
287     SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)
288     SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)
289     SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)
290     ENDDO
291 heimbach 1.10 ENDDO
292    
293     #ifdef ALLOW_AUTODIFF_TAMC
294 heimbach 1.16 CADJ STORE SlopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
295     CADJ STORE SlopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
296 heimbach 1.14 CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
297 jmc 1.19 CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
298     CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
299 heimbach 1.10 #endif /* ALLOW_AUTODIFF_TAMC */
300    
301 jmc 1.26 DO j=1-Oly+1,sNy+Oly-1
302     DO i=1-Olx+1,sNx+Olx-1
303     C Components of Redi/GM tensor
304     Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
305     Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
306     Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
307 adcroft 1.1
308     #ifdef GM_VISBECK_VARIABLE_K
309 jmc 1.26 #ifdef OLD_VISBECK_CALC
310 jmc 1.8
311     C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K
312 edhill 1.18 C but do not know if *taperFct (or **2 ?) is necessary
313 heimbach 1.10 Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)
314 jmc 1.8
315 adcroft 1.1 C-- Depth average of M^2/N^2 * N
316    
317     C Calculate terms for mean Richardson number
318     C which is used in the "variable K" parameterisaton.
319     C Distance between interface above layer and the integration depth
320     deltaH=abs(GM_Visbeck_depth)-abs(rF(k))
321     C If positive we limit this to the layer thickness
322     deltaH=min(deltaH,drF(k))
323     C If negative then we are below the integration level
324     deltaH=max(deltaH,zero_rs)
325     C Now we convert deltaH to a non-dimensional fraction
326     deltaH=deltaH/GM_Visbeck_depth
327    
328 jmc 1.8 IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.
329 jmc 1.19 IF ( Ssq(i,j).NE.0. .AND. dSigmaDr(i,j).NE.0. ) THEN
330     N2= -Gravity*recip_RhoConst*dSigmaDr(i,j)
331 heimbach 1.10 SN=sqrt(Ssq(i,j)*N2)
332 heimbach 1.3 VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH
333 adcroft 1.1 & *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN
334 jmc 1.8 ENDIF
335 adcroft 1.1
336 jmc 1.26 #endif /* OLD_VISBECK_CALC */
337 jmc 1.9 #endif /* GM_VISBECK_VARIABLE_K */
338 jmc 1.26 ENDDO
339 jmc 1.9 ENDDO
340    
341     C-- end 1rst loop on vertical level index k
342     ENDDO
343    
344 adcroft 1.1
345 jmc 1.9 #ifdef GM_VISBECK_VARIABLE_K
346 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
347     CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
348     #endif
349 jmc 1.9 IF ( GM_Visbeck_alpha.NE.0. ) THEN
350     C- Limit range that KapGM can take
351     DO j=1-Oly+1,sNy+Oly-1
352     DO i=1-Olx+1,sNx+Olx-1
353     VisbeckK(i,j,bi,bj)=
354     & MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)
355     ENDDO
356     ENDDO
357     ENDIF
358 heimbach 1.16 cph( NEW
359     #ifdef ALLOW_AUTODIFF_TAMC
360     CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
361     #endif
362     cph)
363 adcroft 1.1 #endif /* GM_VISBECK_VARIABLE_K */
364    
365    
366 jmc 1.9 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
367 heimbach 1.10
368 jmc 1.9 C-- 2nd loop on k : compute Tensor Coeff. at U,V levels.
369     DO k=1,Nr
370     kp1 = MIN(Nr,k+1)
371     maskp1 = 1. _d 0
372     IF (k.GE.Nr) maskp1 = 0. _d 0
373    
374 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
375     kkey = (igmkey-1)*Nr + k
376 heimbach 1.16 #if (defined (GM_NON_UNITY_DIAGONAL) || \
377     defined (GM_VISBECK_VARIABLE_K))
378 heimbach 1.14 CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
379     CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
380     CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
381 heimbach 1.12 #endif
382     #endif
383    
384 jmc 1.9 C- express the Tensor in term of Diffusivity (= m**2 / s )
385 jmc 1.26 DO j=1-Oly+1,sNy+Oly-1
386     DO i=1-Olx+1,sNx+Olx-1
387 heimbach 1.25 #ifdef ALLOW_KAPGM_CONTROL
388 jmc 1.26 Kgm_tmp = GM_isopycK + GM_skewflx*kapgm(i,j,k,bi,bj)
389 heimbach 1.25 #else
390 jmc 1.26 Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K
391 heimbach 1.25 #endif
392 jmc 1.9 #ifdef GM_VISBECK_VARIABLE_K
393 jmc 1.26 & + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)
394 jmc 1.9 #endif
395 jmc 1.26 Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)
396     Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)
397     Kwz(i,j,k,bi,bj)= ( GM_isopycK
398 adcroft 1.1 #ifdef GM_VISBECK_VARIABLE_K
399 jmc 1.26 & + VisbeckK(i,j,bi,bj)
400 adcroft 1.1 #endif
401 jmc 1.26 & )*Kwz(i,j,k,bi,bj)
402     ENDDO
403 adcroft 1.1 ENDDO
404 adcroft 1.4
405 jmc 1.9 #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
406 adcroft 1.1
407     C Gradient of Sigma at U points
408 jmc 1.26 DO j=1-Oly+1,sNy+Oly-1
409     DO i=1-Olx+1,sNx+Olx-1
410     dSigmaDx(i,j)=sigmaX(i,j,k)
411     & *_maskW(i,j,k,bi,bj)
412     dSigmaDy(i,j)=op25*( sigmaY(i-1,j+1,k)+sigmaY(i,j+1,k)
413     & +sigmaY(i-1, j ,k)+sigmaY(i, j ,k)
414     & )*_maskW(i,j,k,bi,bj)
415     dSigmaDr(i,j)=op25*( sigmaR(i-1,j, k )+sigmaR(i,j, k )
416     & +(sigmaR(i-1,j,kp1)+sigmaR(i,j,kp1))*maskp1
417     & )*_maskW(i,j,k,bi,bj)
418     ENDDO
419 adcroft 1.1 ENDDO
420    
421 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
422 heimbach 1.17 CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
423 heimbach 1.12 CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
424     CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
425 jmc 1.19 CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
426 heimbach 1.12 #endif /* ALLOW_AUTODIFF_TAMC */
427    
428 adcroft 1.1 C Calculate slopes for use in tensor, taper and/or clip
429 jmc 1.26 CALL GMREDI_SLOPE_LIMIT(
430 jmc 1.19 O SlopeX, SlopeY,
431 jmc 1.8 O SlopeSqr, taperFct,
432 jmc 1.19 U dSigmaDr,
433     I dSigmaDx, dSigmaDy,
434     I ldd97_LrhoW,rC(k),k,
435 adcroft 1.1 I bi, bj, myThid )
436    
437 heimbach 1.16 cph( NEW
438     #ifdef ALLOW_AUTODIFF_TAMC
439     cph(
440 heimbach 1.17 CADJ STORE SlopeSqr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
441 jmc 1.19 CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
442 heimbach 1.16 cph)
443     #endif /* ALLOW_AUTODIFF_TAMC */
444     cph)
445    
446 jmc 1.9 #ifdef GM_NON_UNITY_DIAGONAL
447 jmc 1.26 c IF ( GM_nonUnitDiag ) THEN
448 jmc 1.9 DO j=1-Oly+1,sNy+Oly-1
449     DO i=1-Olx+1,sNx+Olx-1
450     Kux(i,j,k,bi,bj) =
451     & ( GM_isopycK
452     #ifdef GM_VISBECK_VARIABLE_K
453 heimbach 1.14 & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
454 jmc 1.9 #endif
455 heimbach 1.10 & )
456     & *taperFct(i,j)
457     ENDDO
458     ENDDO
459 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
460 heimbach 1.16 # ifdef GM_EXCLUDE_CLIPPING
461 heimbach 1.12 CADJ STORE Kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
462     # endif
463     #endif
464 heimbach 1.10 DO j=1-Oly+1,sNy+Oly-1
465     DO i=1-Olx+1,sNx+Olx-1
466 jmc 1.9 Kux(i,j,k,bi,bj) = MAX( Kux(i,j,k,bi,bj), GM_Kmin_horiz )
467     ENDDO
468     ENDDO
469 jmc 1.26 c ENDIF
470 jmc 1.9 #endif /* GM_NON_UNITY_DIAGONAL */
471    
472     #ifdef GM_EXTRA_DIAGONAL
473 heimbach 1.12
474     #ifdef ALLOW_AUTODIFF_TAMC
475     CADJ STORE SlopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
476     CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
477     #endif /* ALLOW_AUTODIFF_TAMC */
478 jmc 1.26 IF (GM_ExtraDiag) THEN
479 jmc 1.9 DO j=1-Oly+1,sNy+Oly-1
480     DO i=1-Olx+1,sNx+Olx-1
481     Kuz(i,j,k,bi,bj) =
482 heimbach 1.25 #ifdef ALLOW_KAPGM_CONTROL
483     & ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)
484     #else
485 jmc 1.9 & ( GM_isopycK - GM_skewflx*GM_background_K
486 heimbach 1.25 #endif
487 jmc 1.9 #ifdef GM_VISBECK_VARIABLE_K
488 heimbach 1.14 & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*GM_advect
489 jmc 1.9 #endif
490     & )*SlopeX(i,j)*taperFct(i,j)
491     ENDDO
492     ENDDO
493 jmc 1.26 ENDIF
494 jmc 1.9 #endif /* GM_EXTRA_DIAGONAL */
495 adcroft 1.1
496 jmc 1.22 #ifdef ALLOW_DIAGNOSTICS
497 jmc 1.26 IF (doDiagRediFlx) THEN
498 jmc 1.22 km1 = MAX(k-1,1)
499     DO j=1,sNy
500     DO i=1,sNx+1
501     C store in tmp1k Kuz_Redi
502     tmp1k(i,j) = ( GM_isopycK
503     #ifdef GM_VISBECK_VARIABLE_K
504     & +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0
505     #endif
506     & )*SlopeX(i,j)*taperFct(i,j)
507     ENDDO
508     ENDDO
509     DO j=1,sNy
510     DO i=1,sNx+1
511     C- Vertical gradients interpolated to U points
512     dTdz = (
513     & +recip_drC(k)*
514     & ( maskC(i-1,j,k,bi,bj)*
515     & (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
516     & +maskC( i ,j,k,bi,bj)*
517     & (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
518     & )
519     & +recip_drC(kp1)*
520     & ( maskC(i-1,j,kp1,bi,bj)*
521     & (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
522     & +maskC( i ,j,kp1,bi,bj)*
523     & (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
524     & ) ) * 0.25 _d 0
525 heimbach 1.23 tmp1k(i,j) = dyG(i,j,bi,bj)*drF(k)
526     & * _hFacW(i,j,k,bi,bj)
527 jmc 1.22 & * tmp1k(i,j) * dTdz
528     ENDDO
529     ENDDO
530     CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KuzTz', k,1,2,bi,bj,myThid)
531 jmc 1.26 ENDIF
532 jmc 1.22 #endif /* ALLOW_DIAGNOSTICS */
533    
534 adcroft 1.1 C Gradient of Sigma at V points
535 jmc 1.26 DO j=1-Oly+1,sNy+Oly-1
536     DO i=1-Olx+1,sNx+Olx-1
537     dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
538     & +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k)
539     & )*_maskS(i,j,k,bi,bj)
540     dSigmaDy(i,j)=sigmaY(i,j,k)
541     & *_maskS(i,j,k,bi,bj)
542     dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k )+sigmaR(i,j, k )
543     & +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1
544     & )*_maskS(i,j,k,bi,bj)
545     ENDDO
546 adcroft 1.1 ENDDO
547    
548 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
549     CADJ STORE dSigmaDx(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
550     CADJ STORE dSigmaDy(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
551 jmc 1.19 CADJ STORE dSigmaDr(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
552 heimbach 1.12 #endif /* ALLOW_AUTODIFF_TAMC */
553    
554 adcroft 1.1 C Calculate slopes for use in tensor, taper and/or clip
555 jmc 1.26 CALL GMREDI_SLOPE_LIMIT(
556 jmc 1.19 O SlopeX, SlopeY,
557 jmc 1.8 O SlopeSqr, taperFct,
558 jmc 1.19 U dSigmaDr,
559     I dSigmaDx, dSigmaDy,
560     I ldd97_LrhoS,rC(k),k,
561 adcroft 1.1 I bi, bj, myThid )
562    
563 heimbach 1.16 cph(
564     #ifdef ALLOW_AUTODIFF_TAMC
565     cph(
566     CADJ STORE taperfct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
567     cph)
568     #endif /* ALLOW_AUTODIFF_TAMC */
569     cph)
570    
571 jmc 1.9 #ifdef GM_NON_UNITY_DIAGONAL
572 jmc 1.26 c IF ( GM_nonUnitDiag ) THEN
573 jmc 1.9 DO j=1-Oly+1,sNy+Oly-1
574     DO i=1-Olx+1,sNx+Olx-1
575     Kvy(i,j,k,bi,bj) =
576     & ( GM_isopycK
577     #ifdef GM_VISBECK_VARIABLE_K
578 heimbach 1.14 & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
579 jmc 1.9 #endif
580 heimbach 1.10 & )
581     & *taperFct(i,j)
582     ENDDO
583     ENDDO
584 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
585 heimbach 1.16 # ifdef GM_EXCLUDE_CLIPPING
586 heimbach 1.12 CADJ STORE Kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
587     # endif
588     #endif
589 heimbach 1.10 DO j=1-Oly+1,sNy+Oly-1
590     DO i=1-Olx+1,sNx+Olx-1
591 jmc 1.9 Kvy(i,j,k,bi,bj) = MAX( Kvy(i,j,k,bi,bj), GM_Kmin_horiz )
592     ENDDO
593     ENDDO
594 jmc 1.26 c ENDIF
595 jmc 1.9 #endif /* GM_NON_UNITY_DIAGONAL */
596    
597     #ifdef GM_EXTRA_DIAGONAL
598 heimbach 1.12
599     #ifdef ALLOW_AUTODIFF_TAMC
600     CADJ STORE SlopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
601     CADJ STORE taperFct(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
602     #endif /* ALLOW_AUTODIFF_TAMC */
603 jmc 1.26 IF (GM_ExtraDiag) THEN
604 jmc 1.9 DO j=1-Oly+1,sNy+Oly-1
605     DO i=1-Olx+1,sNx+Olx-1
606     Kvz(i,j,k,bi,bj) =
607 heimbach 1.25 #ifdef ALLOW_KAPGM_CONTROL
608     & ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)
609     #else
610 jmc 1.9 & ( GM_isopycK - GM_skewflx*GM_background_K
611 heimbach 1.25 #endif
612 jmc 1.9 #ifdef GM_VISBECK_VARIABLE_K
613 heimbach 1.14 & +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*GM_advect
614 jmc 1.9 #endif
615     & )*SlopeY(i,j)*taperFct(i,j)
616     ENDDO
617     ENDDO
618 jmc 1.26 ENDIF
619 jmc 1.9 #endif /* GM_EXTRA_DIAGONAL */
620    
621 jmc 1.22 #ifdef ALLOW_DIAGNOSTICS
622 jmc 1.26 IF (doDiagRediFlx) THEN
623 jmc 1.22 c km1 = MAX(k-1,1)
624     DO j=1,sNy+1
625     DO i=1,sNx
626     C store in tmp1k Kvz_Redi
627     tmp1k(i,j) = ( GM_isopycK
628     #ifdef GM_VISBECK_VARIABLE_K
629     & +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0
630     #endif
631     & )*SlopeY(i,j)*taperFct(i,j)
632     ENDDO
633     ENDDO
634     DO j=1,sNy+1
635     DO i=1,sNx
636     C- Vertical gradients interpolated to U points
637     dTdz = (
638     & +recip_drC(k)*
639     & ( maskC(i,j-1,k,bi,bj)*
640     & (theta(i,j-1,km1,bi,bj)-theta(i,j-1,k,bi,bj))
641     & +maskC(i, j ,k,bi,bj)*
642     & (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))
643     & )
644     & +recip_drC(kp1)*
645     & ( maskC(i,j-1,kp1,bi,bj)*
646     & (theta(i,j-1,k,bi,bj)-theta(i,j-1,kp1,bi,bj))
647     & +maskC(i, j ,kp1,bi,bj)*
648     & (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))
649     & ) ) * 0.25 _d 0
650 heimbach 1.23 tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)
651     & * _hFacS(i,j,k,bi,bj)
652 jmc 1.22 & * tmp1k(i,j) * dTdz
653     ENDDO
654     ENDDO
655     CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KvzTz', k,1,2,bi,bj,myThid)
656 jmc 1.26 ENDIF
657 jmc 1.22 #endif /* ALLOW_DIAGNOSTICS */
658    
659 jmc 1.9 #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
660    
661     C-- end 2nd loop on vertical level index k
662     ENDDO
663 adcroft 1.1
664    
665 jmc 1.9 #ifdef GM_BOLUS_ADVEC
666     IF (GM_AdvForm) THEN
667 jmc 1.26 CALL GMREDI_CALC_PSI_B(
668 jmc 1.9 I bi, bj, iMin, iMax, jMin, jMax,
669     I sigmaX, sigmaY, sigmaR,
670 jmc 1.19 I ldd97_LrhoW, ldd97_LrhoS,
671 jmc 1.26 I myThid )
672 jmc 1.9 ENDIF
673     #endif
674 adcroft 1.1
675 jmc 1.19 #ifdef ALLOW_TIMEAVE
676     C-- Time-average
677     IF ( taveFreq.GT.0. ) THEN
678    
679     CALL TIMEAVE_CUMULATE( GM_Kwx_T, Kwx, Nr,
680     & deltaTclock, bi, bj, myThid )
681     CALL TIMEAVE_CUMULATE( GM_Kwy_T, Kwy, Nr,
682     & deltaTclock, bi, bj, myThid )
683     CALL TIMEAVE_CUMULATE( GM_Kwz_T, Kwz, Nr,
684     & deltaTclock, bi, bj, myThid )
685     #ifdef GM_VISBECK_VARIABLE_K
686     IF ( GM_Visbeck_alpha.NE.0. ) THEN
687     CALL TIMEAVE_CUMULATE( Visbeck_K_T, VisbeckK, 1,
688     & deltaTclock, bi, bj, myThid )
689     ENDIF
690     #endif
691     #ifdef GM_BOLUS_ADVEC
692     IF ( GM_AdvForm ) THEN
693     CALL TIMEAVE_CUMULATE( GM_PsiXtave, GM_PsiX, Nr,
694     & deltaTclock, bi, bj, myThid )
695     CALL TIMEAVE_CUMULATE( GM_PsiYtave, GM_PsiY, Nr,
696     & deltaTclock, bi, bj, myThid )
697     ENDIF
698     #endif
699     DO k=1,Nr
700     GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock
701     ENDDO
702    
703     ENDIF
704     #endif /* ALLOW_TIMEAVE */
705    
706 jmc 1.20 #ifdef ALLOW_DIAGNOSTICS
707     IF ( useDiagnostics ) THEN
708 jmc 1.24 CALL GMREDI_DIAGNOSTICS_FILL(bi,bj,myThid)
709 jmc 1.20 ENDIF
710     #endif /* ALLOW_DIAGNOSTICS */
711    
712 adcroft 1.1 #endif /* ALLOW_GMREDI */
713    
714     RETURN
715     END
716 heimbach 1.2
717    
718     SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
719 jmc 1.9 I bi, bj, iMin, iMax, jMin, jMax,
720 heimbach 1.2 I sigmaX, sigmaY, sigmaR,
721     I myThid )
722     C /==========================================================\
723     C | SUBROUTINE GMREDI_CALC_TENSOR |
724     C | o Calculate tensor elements for GM/Redi tensor. |
725     C |==========================================================|
726     C \==========================================================/
727     IMPLICIT NONE
728    
729     C == Global variables ==
730     #include "SIZE.h"
731     #include "EEPARAMS.h"
732     #include "GMREDI.h"
733    
734     C == Routine arguments ==
735     C
736     _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
737     _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
738     _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
739 jmc 1.9 INTEGER bi,bj,iMin,iMax,jMin,jMax
740 heimbach 1.2 INTEGER myThid
741     CEndOfInterface
742    
743 jmc 1.9 INTEGER i, j, k
744 heimbach 1.2
745     #ifdef ALLOW_GMREDI
746    
747 jmc 1.9 DO k=1,Nr
748     DO j=1-Oly+1,sNy+Oly-1
749     DO i=1-Olx+1,sNx+Olx-1
750     Kwx(i,j,k,bi,bj) = 0.0
751     Kwy(i,j,k,bi,bj) = 0.0
752     Kwz(i,j,k,bi,bj) = 0.0
753     ENDDO
754 heimbach 1.2 ENDDO
755     ENDDO
756     #endif /* ALLOW_GMREDI */
757    
758 jmc 1.9 RETURN
759     END
  ViewVC Help
Powered by ViewVC 1.1.22