/[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.10 - (hide annotations) (download)
Tue Jun 26 15:34:15 2007 UTC (16 years, 11 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59e, checkpoint59g, checkpoint59f, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint59j
Changes since 1.9: +5 -1 lines
Modify routines to restore adjoint.

1 heimbach 1.10 C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_slope_psi.F,v 1.9 2005/12/08 21:40:16 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.9 C | depthZ contains the depth (< 0 !) [m] |
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.9 _RS 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 jmc 1.9 taperX(i,j) = 1. _d 0
128 heimbach 1.3 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 jmc 1.9 taperY(i,j) = 1. _d 0
163 jmc 1.1 ENDDO
164     ENDDO
165    
166 heimbach 1.6 #endif /* GM_EXCLUDE_CLIPPING */
167    
168 heimbach 1.10 ELSEIF (GM_taper_scheme.EQ.'fm07') THEN
169    
170     STOP 'GMREDI_SLOPE_PSI: AdvForm not yet implemented for fm07'
171    
172 jmc 1.1 ELSE
173    
174 heimbach 1.6 #ifdef GM_EXCLUDE_TAPERING
175    
176     STOP 'Need to compile without "#define GM_EXCLUDE_TAPERING"'
177    
178     #else /* GM_EXCLUDE_TAPERING */
179    
180 heimbach 1.3 #ifdef ALLOW_AUTODIFF_TAMC
181     CADJ STORE slopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
182     CADJ STORE dSigmaDrW(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
183     #endif
184    
185 jmc 1.1 C- Compute the slope, no clipping, but avoid reverse slope in negatively
186     C stratified (Sigma_Z > 0) region :
187 jmc 1.7 DO j=1-Oly,sNy+Oly
188     DO i=1-Olx+1,sNx+Olx
189     IF (dSigmaDrW(i,j).GE.-GM_Small_Number)
190     & dSigmaDrW(i,j) = -GM_Small_Number
191 heimbach 1.3 ENDDO
192     ENDDO
193     #ifdef ALLOW_AUTODIFF_TAMC
194     CADJ STORE dsigmadrW(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
195     #endif
196 jmc 1.7 DO j=1-Oly,sNy+Oly
197     DO i=1-Olx+1,sNx+Olx
198 heimbach 1.3 SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)
199 jmc 1.9 taperX(i,j) = 1. _d 0
200 heimbach 1.6 ENDDO
201     ENDDO
202     #ifdef ALLOW_AUTODIFF_TAMC
203     CADJ STORE slopex(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
204     #endif
205 jmc 1.7 DO j=1-Oly,sNy+Oly
206     DO i=1-Olx+1,sNx+Olx
207     IF ( ABS(SlopeX(i,j)) .GE. slopeCutoff ) THEN
208     SlopeX(i,j) = SIGN(slopeCutoff,SlopeX(i,j))
209 heimbach 1.6 taperX(i,j) = 0. _d 0
210     ENDIF
211 heimbach 1.3 ENDDO
212     ENDDO
213 jmc 1.1
214 heimbach 1.3 #ifdef ALLOW_AUTODIFF_TAMC
215     CADJ STORE slopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
216     CADJ STORE dSigmaDrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
217     #endif
218 jmc 1.1
219 jmc 1.7 DO j=1-Oly+1,sNy+Oly
220     DO i=1-Olx,sNx+Olx
221     IF (dSigmaDrS(i,j).GE.-GM_Small_Number)
222     & dSigmaDrS(i,j) = -GM_Small_Number
223 heimbach 1.3 ENDDO
224     ENDDO
225     #ifdef ALLOW_AUTODIFF_TAMC
226     CADJ STORE dsigmadrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
227     #endif
228 jmc 1.7 DO j=1-Oly+1,sNy+Oly
229     DO i=1-Olx,sNx+Olx
230 jmc 1.1 SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)
231 jmc 1.9 taperY(i,j) = 1. _d 0
232 jmc 1.1 ENDDO
233     ENDDO
234 heimbach 1.6 #ifdef ALLOW_AUTODIFF_TAMC
235     CADJ STORE slopey(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
236     #endif
237 jmc 1.7 DO j=1-Oly+1,sNy+Oly
238     DO i=1-Olx,sNx+Olx
239     IF ( ABS(SlopeY(i,j)) .GE. slopeCutoff ) THEN
240     SlopeY(i,j) = SIGN(slopeCutoff,SlopeY(i,j))
241 heimbach 1.6 taperY(i,j) = 0. _d 0
242     ENDIF
243     ENDDO
244     ENDDO
245 jmc 1.1
246     C- Compute the tapering function for the GM+Redi tensor :
247    
248 heimbach 1.3 #ifdef ALLOW_AUTODIFF_TAMC
249     CADJ STORE slopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
250     CADJ STORE slopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
251     #endif
252    
253 jmc 1.1 IF (GM_taper_scheme.EQ.'linear') THEN
254    
255     C- Simplest adiabatic tapering = Smax/Slope (linear)
256 jmc 1.7 DO j=1-Oly,sNy+Oly
257     DO i=1-Olx+1,sNx+Olx
258     Smod = ABS(SlopeX(i,j))
259 heimbach 1.6 IF ( Smod .GT. GM_maxSlope .AND.
260 jmc 1.7 & Smod .LT. slopeCutoff )
261 heimbach 1.6 & taperX(i,j)=GM_maxSlope/(Smod+GM_Small_Number)
262 jmc 1.7 ENDDO
263     ENDDO
264     DO j=1-Oly+1,sNy+Oly
265     DO i=1-Olx,sNx+Olx
266     Smod = ABS(SlopeY(i,j))
267 heimbach 1.6 IF ( Smod .GT. GM_maxSlope .AND.
268 jmc 1.7 & Smod .LT. slopeCutoff )
269 heimbach 1.6 & taperY(i,j)=GM_maxSlope/(Smod+GM_Small_Number)
270 jmc 1.1 ENDDO
271     ENDDO
272    
273     ELSEIF (GM_taper_scheme.EQ.'gkw91') THEN
274    
275     C- Gerdes, Koberle and Willebrand, Clim. Dyn. 1991
276     maxSlopeSqr = GM_maxSlope*GM_maxSlope
277 jmc 1.7 DO j=1-Oly,sNy+Oly
278     DO i=1-Olx+1,sNx+Olx
279     IF ( ABS(SlopeX(i,j)) .GT. GM_maxSlope .AND.
280     & ABS(SlopeX(i,j)) .LT. slopeCutoff )
281 heimbach 1.6 & taperX(i,j)=maxSlopeSqr/
282     & ( SlopeX(i,j)*SlopeX(i,j) + GM_Small_Number )
283 jmc 1.7 ENDDO
284     ENDDO
285     DO j=1-Oly+1,sNy+Oly
286     DO i=1-Olx,sNx+Olx
287     IF ( ABS(SlopeY(i,j)) .GT. GM_maxSlope .AND.
288     & ABS(SlopeY(i,j)) .LT. slopeCutoff )
289 heimbach 1.6 & taperY(i,j)=maxSlopeSqr/
290     & ( SlopeY(i,j)*SlopeY(i,j) + GM_Small_Number )
291 jmc 1.1 ENDDO
292     ENDDO
293    
294     ELSEIF (GM_taper_scheme.EQ.'dm95') THEN
295    
296     C- Danabasoglu and McWilliams, J. Clim. 1995
297 jmc 1.7 DO j=1-Oly,sNy+Oly
298     DO i=1-Olx+1,sNx+Olx
299     Smod = ABS(SlopeX(i,j))
300     taperX(i,j)=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
301     ENDDO
302     ENDDO
303     DO j=1-Oly+1,sNy+Oly
304     DO i=1-Olx,sNx+Olx
305     Smod = ABS(SlopeY(i,j))
306     taperY(i,j)=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
307 jmc 1.1 ENDDO
308     ENDDO
309    
310     ELSEIF (GM_taper_scheme.EQ.'ldd97') THEN
311    
312     C- Large, Danabasoglu and Doney, JPO 1997
313    
314 jmc 1.7 DO j=1-Oly,sNy+Oly
315     DO i=1-Olx+1,sNx+Olx
316     Smod = ABS(SlopeX(i,j))
317     IF ( Smod .LT. slopeCutoff ) THEN
318 jmc 1.9 f1=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
319     IF (Smod.NE.0.) THEN
320     Rnondim = -depthZ/(LrhoW(i,j)*Smod)
321     ELSE
322     Rnondim = 1.
323     ENDIF
324     IF ( Rnondim.GE.1. _d 0 ) THEN
325     f2 = 1. _d 0
326     ELSE
327     f2 = op5*( 1. _d 0 + SIN( fpi*(Rnondim-op5) ))
328     ENDIF
329     taperX(i,j)=f1*f2
330 jmc 1.7 ENDIF
331     ENDDO
332     ENDDO
333 jmc 1.1
334 jmc 1.7 DO j=1-Oly+1,sNy+Oly
335     DO i=1-Olx,sNx+Olx
336     Smod = ABS(SlopeY(i,j))
337     IF ( Smod .LT. slopeCutoff ) THEN
338 jmc 1.9 f1=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
339     IF (Smod.NE.0.) THEN
340     Rnondim = -depthZ/(LrhoS(i,j)*Smod)
341     ELSE
342     Rnondim = 1.
343     ENDIF
344     IF ( Rnondim.GE.1. _d 0 ) THEN
345     f2 = 1. _d 0
346     ELSE
347     f2 = op5*( 1. _d 0 + SIN( fpi*(Rnondim-op5) ))
348     ENDIF
349     taperY(i,j)=f1*f2
350 jmc 1.7 ENDIF
351 jmc 1.1 ENDDO
352     ENDDO
353    
354     ELSEIF (GM_taper_scheme.NE.' ') THEN
355     STOP 'GMREDI_SLOPE_PSI: Bad GM_taper_scheme'
356     ENDIF
357 heimbach 1.6
358     #endif /* GM_EXCLUDE_TAPERING */
359 jmc 1.1
360     ENDIF
361    
362 heimbach 1.3 #endif /* BOLUS_ADVEC */
363 jmc 1.1 #endif /* ALLOW_GMREDI */
364    
365     RETURN
366     END

  ViewVC Help
Powered by ViewVC 1.1.22