/[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.7 - (hide annotations) (download)
Sun Nov 21 15:57:17 2004 UTC (19 years, 6 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint57m_post, checkpoint57g_pre, checkpoint57s_post, checkpoint57b_post, checkpoint57g_post, checkpoint56b_post, checkpoint57r_post, checkpoint57d_post, checkpoint57i_post, checkpoint57, checkpoint57n_post, checkpoint57l_post, checkpoint57t_post, checkpoint57v_post, checkpoint57f_post, checkpoint57a_post, checkpoint57h_pre, checkpoint57h_post, checkpoint57c_post, checkpoint57c_pre, checkpoint57e_post, checkpoint57p_post, checkpint57u_post, checkpoint57q_post, eckpoint57e_pre, checkpoint56a_post, checkpoint57h_done, checkpoint57j_post, checkpoint57f_pre, checkpoint56c_post, checkpoint57a_pre, checkpoint57o_post, checkpoint57k_post, checkpoint57w_post
Changes since 1.6: +121 -122 lines
fix ldd97 slope limit ; extend Psi-Bolus valid domain ; clean-up time-ave

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

  ViewVC Help
Powered by ViewVC 1.1.22