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