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

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

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


Revision 1.8 - (hide annotations) (download)
Fri Nov 4 01:31:36 2005 UTC (18 years, 6 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint57y_post, checkpoint57y_pre, checkpoint57x_post
Changes since 1.7: +1 -3 lines
remove unused variables (reduces number of compiler warning)

1 jmc 1.8 C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_slope_psi.F,v 1.7 2004/11/21 15:57:17 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "GMREDI_OPTIONS.h"
5    
6     CStartOfInterface
7 jmc 1.5 SUBROUTINE GMREDI_SLOPE_PSI(
8 jmc 1.7 O taperX, taperY,
9 jmc 1.1 U SlopeX, SlopeY,
10 jmc 1.7 U dSigmaDrW,dSigmaDrS,
11     I LrhoW, LrhoS, depthZ, K,
12 jmc 1.1 I bi,bj, myThid )
13     C /==========================================================\
14 jmc 1.5 C | SUBROUTINE GMREDI_SLOPE_PSI |
15 jmc 1.1 C | o Calculate slopes for use in GM/Redi tensor |
16     C |==========================================================|
17     C | On entry: |
18 jmc 1.7 C | dSigmaDrW,S contains the d/dz Sigma |
19 jmc 1.1 C | SlopeX/Y contains X/Y gradients of sigma |
20 jmc 1.7 C | depthZ contains the height (m) of level |
21 jmc 1.1 C | On exit: |
22 jmc 1.7 C | dSigmaDrW,S contains the effective dSig/dz |
23 jmc 1.1 C | SlopeX/Y contains X/Y slopes |
24     C | taperFct contains tapering funct. value ; |
25     C | = 1 when using no tapering |
26     C \==========================================================/
27     IMPLICIT NONE
28    
29     C == Global variables ==
30     #include "SIZE.h"
31     #include "EEPARAMS.h"
32     #include "GMREDI.h"
33     #include "PARAMS.h"
34    
35 heimbach 1.3 #ifdef ALLOW_AUTODIFF_TAMC
36     #include "tamc.h"
37     #include "tamc_keys.h"
38     #endif /* ALLOW_AUTODIFF_TAMC */
39    
40 jmc 1.1 C == Routine arguments ==
41     C
42 jmc 1.7 _RL taperX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
43     _RL taperY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
44 jmc 1.1 _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
45     _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
46     _RL dSigmaDrW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47     _RL dSigmaDrS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48 jmc 1.7 _RL LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49     _RL LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50 jmc 1.1 _RL depthZ
51 heimbach 1.3 INTEGER K,bi,bj,myThid
52 jmc 1.1 CEndOfInterface
53    
54     #ifdef ALLOW_GMREDI
55 heimbach 1.3 #ifdef GM_BOLUS_ADVEC
56 jmc 1.1
57     C == Local variables ==
58 heimbach 1.3 _RL dSigmaDrLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59 jmc 1.7 _RL f1,Smod,f2,Rnondim
60 jmc 1.1 _RL maxSlopeSqr
61 jmc 1.7 _RL slopeCutoff
62 jmc 1.1 _RL fpi
63     PARAMETER(fpi=3.141592653589793047592d0)
64     INTEGER i,j
65    
66 jmc 1.7 slopeCutoff = SQRT( GM_slopeSqCutoff )
67    
68 heimbach 1.3 #ifdef ALLOW_AUTODIFF_TAMC
69     act1 = bi - myBxLo(myThid)
70     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
71     act2 = bj - myByLo(myThid)
72     max2 = myByHi(myThid) - myByLo(myThid) + 1
73     act3 = myThid - 1
74     max3 = nTx*nTy
75     act4 = ikey_dynamics - 1
76 heimbach 1.6 igmkey = (act1 + 1) + act2*max1
77     & + act3*max1*max2
78     & + act4*max1*max2*max3
79     kkey = (igmkey-1)*Nr + k
80 heimbach 1.3 #endif /* ALLOW_AUTODIFF_TAMC */
81    
82 jmc 1.1 IF (GM_taper_scheme.EQ.'orig' .OR.
83     & GM_taper_scheme.EQ.'clipping') THEN
84    
85 heimbach 1.6 #ifdef GM_EXCLUDE_CLIPPING
86    
87     STOP 'Need to compile without "#define GM_EXCLUDE_CLIPPING"'
88    
89     #else /* GM_EXCLUDE_CLIPPING */
90    
91 jmc 1.1 C- Original implementation in mitgcmuv
92     C (this turns out to be the same as Cox slope clipping)
93    
94 heimbach 1.3 C-- X-comp
95    
96     #ifdef ALLOW_AUTODIFF_TAMC
97 jmc 1.7 DO j=1-Oly,sNy+Oly
98     DO i=1-Olx+1,sNx+Olx
99     dSigmaDrLtd(i,j) = 0. _d 0
100     ENDDO
101     ENDDO
102 heimbach 1.3 #endif /* ALLOW_AUTODIFF_TAMC */
103    
104 jmc 1.1 C- Cox 1987 "Slope clipping"
105 jmc 1.7 DO j=1-Oly,sNy+Oly
106     DO i=1-Olx+1,sNx+Olx
107     dSigmaDrLtd(i,j) = -(GM_Small_Number+
108     & ABS(SlopeX(i,j))*GM_rMaxSlope)
109 heimbach 1.3 ENDDO
110     ENDDO
111     #ifdef ALLOW_AUTODIFF_TAMC
112 jmc 1.7 CADJ STORE dSigmaDrLtd(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
113 heimbach 1.3 CADJ STORE dSigmaDrW(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
114     #endif
115 jmc 1.7 DO j=1-Oly,sNy+Oly
116     DO i=1-Olx+1,sNx+Olx
117     IF (dSigmaDrW(i,j).GE.dSigmaDrLtd(i,j))
118     & dSigmaDrW(i,j) = dSigmaDrLtd(i,j)
119 heimbach 1.3 ENDDO
120     ENDDO
121     #ifdef ALLOW_AUTODIFF_TAMC
122     CADJ STORE dSigmaDrW(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
123     #endif
124 jmc 1.7 DO j=1-Oly,sNy+Oly
125     DO i=1-Olx+1,sNx+Olx
126 heimbach 1.3 SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)
127     taperX(i,j)=1. _d 0
128     ENDDO
129     ENDDO
130 jmc 1.1
131 heimbach 1.3 C-- Y-comp
132 jmc 1.1
133 heimbach 1.3 #ifdef ALLOW_AUTODIFF_TAMC
134 jmc 1.7 DO j=1-Oly+1,sNy+Oly
135     DO i=1-Olx,sNx+Olx
136     dSigmaDrLtd(i,j) = 0. _d 0
137     ENDDO
138     ENDDO
139 heimbach 1.3 #endif /* ALLOW_AUTODIFF_TAMC */
140 jmc 1.7 DO j=1-Oly+1,sNy+Oly
141     DO i=1-Olx,sNx+Olx
142     dSigmaDrLtd(i,j) = -(GM_Small_Number+
143     & ABS(SlopeY(i,j))*GM_rMaxSlope)
144 heimbach 1.3 ENDDO
145     ENDDO
146     #ifdef ALLOW_AUTODIFF_TAMC
147 jmc 1.7 CADJ STORE dSigmaDrLtd(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
148 heimbach 1.3 CADJ STORE dSigmaDrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
149     #endif
150 jmc 1.7 DO j=1-Oly+1,sNy+Oly
151     DO i=1-Olx,sNx+Olx
152     IF (dSigmaDrS(i,j).GE.dSigmaDrLtd(i,j))
153     & dSigmaDrS(i,j) = dSigmaDrLtd(i,j)
154 heimbach 1.3 ENDDO
155     ENDDO
156     #ifdef ALLOW_AUTODIFF_TAMC
157     CADJ STORE dSigmaDrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
158     #endif
159 jmc 1.7 DO j=1-Oly+1,sNy+Oly
160     DO i=1-Olx,sNx+Olx
161 jmc 1.1 SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)
162     taperY(i,j)=1. _d 0
163     ENDDO
164     ENDDO
165    
166 heimbach 1.6 #endif /* GM_EXCLUDE_CLIPPING */
167    
168 jmc 1.1 ELSE
169    
170 heimbach 1.6 #ifdef GM_EXCLUDE_TAPERING
171    
172     STOP 'Need to compile without "#define GM_EXCLUDE_TAPERING"'
173    
174     #else /* GM_EXCLUDE_TAPERING */
175    
176 heimbach 1.3 #ifdef ALLOW_AUTODIFF_TAMC
177     CADJ STORE slopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
178     CADJ STORE dSigmaDrW(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
179     #endif
180    
181 jmc 1.1 C- Compute the slope, no clipping, but avoid reverse slope in negatively
182     C stratified (Sigma_Z > 0) region :
183 jmc 1.7 DO j=1-Oly,sNy+Oly
184     DO i=1-Olx+1,sNx+Olx
185     IF (dSigmaDrW(i,j).GE.-GM_Small_Number)
186     & dSigmaDrW(i,j) = -GM_Small_Number
187 heimbach 1.3 ENDDO
188     ENDDO
189     #ifdef ALLOW_AUTODIFF_TAMC
190     CADJ STORE dsigmadrW(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
191     #endif
192 jmc 1.7 DO j=1-Oly,sNy+Oly
193     DO i=1-Olx+1,sNx+Olx
194 heimbach 1.3 SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)
195 heimbach 1.6 taperX(i,j)= 1. _d 0
196     ENDDO
197     ENDDO
198     #ifdef ALLOW_AUTODIFF_TAMC
199     CADJ STORE slopex(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
200     #endif
201 jmc 1.7 DO j=1-Oly,sNy+Oly
202     DO i=1-Olx+1,sNx+Olx
203     IF ( ABS(SlopeX(i,j)) .GE. slopeCutoff ) THEN
204     SlopeX(i,j) = SIGN(slopeCutoff,SlopeX(i,j))
205 heimbach 1.6 taperX(i,j) = 0. _d 0
206     ENDIF
207 heimbach 1.3 ENDDO
208     ENDDO
209 jmc 1.1
210 heimbach 1.3 #ifdef ALLOW_AUTODIFF_TAMC
211     CADJ STORE slopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
212     CADJ STORE dSigmaDrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
213     #endif
214 jmc 1.1
215 jmc 1.7 DO j=1-Oly+1,sNy+Oly
216     DO i=1-Olx,sNx+Olx
217     IF (dSigmaDrS(i,j).GE.-GM_Small_Number)
218     & dSigmaDrS(i,j) = -GM_Small_Number
219 heimbach 1.3 ENDDO
220     ENDDO
221     #ifdef ALLOW_AUTODIFF_TAMC
222     CADJ STORE dsigmadrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
223     #endif
224 jmc 1.7 DO j=1-Oly+1,sNy+Oly
225     DO i=1-Olx,sNx+Olx
226 jmc 1.1 SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)
227     taperY(i,j)=1. _d 0
228     ENDDO
229     ENDDO
230 heimbach 1.6 #ifdef ALLOW_AUTODIFF_TAMC
231     CADJ STORE slopey(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
232     #endif
233 jmc 1.7 DO j=1-Oly+1,sNy+Oly
234     DO i=1-Olx,sNx+Olx
235     IF ( ABS(SlopeY(i,j)) .GE. slopeCutoff ) THEN
236     SlopeY(i,j) = SIGN(slopeCutoff,SlopeY(i,j))
237 heimbach 1.6 taperY(i,j) = 0. _d 0
238     ENDIF
239     ENDDO
240     ENDDO
241 jmc 1.1
242     C- Compute the tapering function for the GM+Redi tensor :
243    
244 heimbach 1.3 #ifdef ALLOW_AUTODIFF_TAMC
245     CADJ STORE slopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
246     CADJ STORE slopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
247     #endif
248    
249 jmc 1.1 IF (GM_taper_scheme.EQ.'linear') THEN
250    
251     C- Simplest adiabatic tapering = Smax/Slope (linear)
252 jmc 1.7 DO j=1-Oly,sNy+Oly
253     DO i=1-Olx+1,sNx+Olx
254     Smod = ABS(SlopeX(i,j))
255 heimbach 1.6 IF ( Smod .GT. GM_maxSlope .AND.
256 jmc 1.7 & Smod .LT. slopeCutoff )
257 heimbach 1.6 & taperX(i,j)=GM_maxSlope/(Smod+GM_Small_Number)
258 jmc 1.7 ENDDO
259     ENDDO
260     DO j=1-Oly+1,sNy+Oly
261     DO i=1-Olx,sNx+Olx
262     Smod = ABS(SlopeY(i,j))
263 heimbach 1.6 IF ( Smod .GT. GM_maxSlope .AND.
264 jmc 1.7 & Smod .LT. slopeCutoff )
265 heimbach 1.6 & taperY(i,j)=GM_maxSlope/(Smod+GM_Small_Number)
266 jmc 1.1 ENDDO
267     ENDDO
268    
269     ELSEIF (GM_taper_scheme.EQ.'gkw91') THEN
270    
271     C- Gerdes, Koberle and Willebrand, Clim. Dyn. 1991
272     maxSlopeSqr = GM_maxSlope*GM_maxSlope
273 jmc 1.7 DO j=1-Oly,sNy+Oly
274     DO i=1-Olx+1,sNx+Olx
275     IF ( ABS(SlopeX(i,j)) .GT. GM_maxSlope .AND.
276     & ABS(SlopeX(i,j)) .LT. slopeCutoff )
277 heimbach 1.6 & taperX(i,j)=maxSlopeSqr/
278     & ( SlopeX(i,j)*SlopeX(i,j) + GM_Small_Number )
279 jmc 1.7 ENDDO
280     ENDDO
281     DO j=1-Oly+1,sNy+Oly
282     DO i=1-Olx,sNx+Olx
283     IF ( ABS(SlopeY(i,j)) .GT. GM_maxSlope .AND.
284     & ABS(SlopeY(i,j)) .LT. slopeCutoff )
285 heimbach 1.6 & taperY(i,j)=maxSlopeSqr/
286     & ( SlopeY(i,j)*SlopeY(i,j) + GM_Small_Number )
287 jmc 1.1 ENDDO
288     ENDDO
289    
290     ELSEIF (GM_taper_scheme.EQ.'dm95') THEN
291    
292     C- Danabasoglu and McWilliams, J. Clim. 1995
293 jmc 1.7 DO j=1-Oly,sNy+Oly
294     DO i=1-Olx+1,sNx+Olx
295     Smod = ABS(SlopeX(i,j))
296     taperX(i,j)=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
297     ENDDO
298     ENDDO
299     DO j=1-Oly+1,sNy+Oly
300     DO i=1-Olx,sNx+Olx
301     Smod = ABS(SlopeY(i,j))
302     taperY(i,j)=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
303 jmc 1.1 ENDDO
304     ENDDO
305    
306     ELSEIF (GM_taper_scheme.EQ.'ldd97') THEN
307    
308     C- Large, Danabasoglu and Doney, JPO 1997
309    
310 jmc 1.7 DO j=1-Oly,sNy+Oly
311     DO i=1-Olx+1,sNx+Olx
312     Smod = ABS(SlopeX(i,j))
313     IF ( Smod .LT. slopeCutoff ) THEN
314     f1=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
315     IF (Smod.NE.0.) THEN
316     Rnondim=depthZ/(LrhoW(i,j)*Smod)
317     ELSE
318 jmc 1.1 Rnondim=0.
319 jmc 1.7 ENDIF
320     f2=op5*( 1. _d 0 + SIN( fpi*(Rnondim-op5)))
321 jmc 1.1 taperX(i,j)=f1*f2
322 jmc 1.7 ENDIF
323     ENDDO
324     ENDDO
325 jmc 1.1
326 jmc 1.7 DO j=1-Oly+1,sNy+Oly
327     DO i=1-Olx,sNx+Olx
328     Smod = ABS(SlopeY(i,j))
329     IF ( Smod .LT. slopeCutoff ) THEN
330     f1=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
331     IF (Smod.NE.0.) THEN
332     Rnondim=depthZ/(LrhoS(i,j)*Smod)
333     ELSE
334 jmc 1.1 Rnondim=0.
335 jmc 1.7 ENDIF
336     f2=op5*( 1. _d 0 + SIN( fpi*(Rnondim-op5)))
337 jmc 1.1 taperY(i,j)=f1*f2
338 jmc 1.7 ENDIF
339 jmc 1.1 ENDDO
340     ENDDO
341    
342     ELSEIF (GM_taper_scheme.NE.' ') THEN
343     STOP 'GMREDI_SLOPE_PSI: Bad GM_taper_scheme'
344     ENDIF
345 heimbach 1.6
346     #endif /* GM_EXCLUDE_TAPERING */
347 jmc 1.1
348     ENDIF
349    
350 heimbach 1.3 #endif /* BOLUS_ADVEC */
351 jmc 1.1 #endif /* ALLOW_GMREDI */
352    
353     RETURN
354     END

  ViewVC Help
Powered by ViewVC 1.1.22