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

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

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


Revision 1.12 - (show annotations) (download)
Tue Dec 8 21:42:22 2009 UTC (14 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62c, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62w, checkpoint62z, checkpoint62y, checkpoint62x, checkpoint63g, checkpoint62, checkpoint63, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint62b
Changes since 1.11: +9 -5 lines
avoid un-used variables

1 C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_slope_psi.F,v 1.11 2008/05/30 21:10:25 gforget Exp $
2 C $Name: $
3
4 #include "GMREDI_OPTIONS.h"
5
6 CStartOfInterface
7 SUBROUTINE GMREDI_SLOPE_PSI(
8 O taperX, taperY,
9 U SlopeX, SlopeY,
10 U dSigmaDrW,dSigmaDrS,
11 I LrhoW, LrhoS, depthZ, K,
12 I bi,bj, myThid )
13 C /==========================================================\
14 C | SUBROUTINE GMREDI_SLOPE_PSI |
15 C | o Calculate slopes for use in GM/Redi tensor |
16 C |==========================================================|
17 C | On entry: |
18 C | dSigmaDrW,S contains the d/dz Sigma |
19 C | SlopeX/Y contains X/Y gradients of sigma |
20 C | depthZ contains the depth (< 0 !) [m] |
21 C | On exit: |
22 C | dSigmaDrW,S contains the effective dSig/dz |
23 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 #ifdef ALLOW_AUTODIFF_TAMC
36 #include "tamc.h"
37 #include "tamc_keys.h"
38 #endif /* ALLOW_AUTODIFF_TAMC */
39
40 C == Routine arguments ==
41 C
42 _RL taperX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
43 _RL taperY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
44 _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 _RL LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49 _RL LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50 _RS depthZ
51 INTEGER K,bi,bj,myThid
52 CEndOfInterface
53
54 #ifdef ALLOW_GMREDI
55 #ifdef GM_BOLUS_ADVEC
56
57 C == Local variables ==
58 _RL dSigmaDrLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59 _RL f1,Smod,f2,Rnondim
60 _RL maxSlopeSqr
61 _RL slopeCutoff
62 _RL fpi
63 PARAMETER(fpi=3.141592653589793047592d0)
64 INTEGER i,j
65 #ifdef GMREDI_WITH_STABLE_ADJOINT
66 _RL slopeTmpSpec,slopeMaxSpec
67 #endif
68
69 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
70
71 slopeCutoff = SQRT( GM_slopeSqCutoff )
72
73 #ifdef ALLOW_AUTODIFF_TAMC
74 act1 = bi - myBxLo(myThid)
75 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
76 act2 = bj - myByLo(myThid)
77 max2 = myByHi(myThid) - myByLo(myThid) + 1
78 act3 = myThid - 1
79 max3 = nTx*nTy
80 act4 = ikey_dynamics - 1
81 igmkey = (act1 + 1) + act2*max1
82 & + act3*max1*max2
83 & + act4*max1*max2*max3
84 kkey = (igmkey-1)*Nr + k
85 #endif /* ALLOW_AUTODIFF_TAMC */
86
87 #ifndef GMREDI_WITH_STABLE_ADJOINT
88 c common case:
89
90 IF (GM_taper_scheme.EQ.'orig' .OR.
91 & GM_taper_scheme.EQ.'clipping') THEN
92
93 #ifdef GM_EXCLUDE_CLIPPING
94
95 STOP 'Need to compile without "#define GM_EXCLUDE_CLIPPING"'
96
97 #else /* GM_EXCLUDE_CLIPPING */
98
99 C- Original implementation in mitgcmuv
100 C (this turns out to be the same as Cox slope clipping)
101
102 C-- X-comp
103
104 #ifdef ALLOW_AUTODIFF_TAMC
105 DO j=1-Oly,sNy+Oly
106 DO i=1-Olx+1,sNx+Olx
107 dSigmaDrLtd(i,j) = 0. _d 0
108 ENDDO
109 ENDDO
110 #endif /* ALLOW_AUTODIFF_TAMC */
111
112 C- Cox 1987 "Slope clipping"
113 DO j=1-Oly,sNy+Oly
114 DO i=1-Olx+1,sNx+Olx
115 dSigmaDrLtd(i,j) = -(GM_Small_Number+
116 & ABS(SlopeX(i,j))*GM_rMaxSlope)
117 ENDDO
118 ENDDO
119 #ifdef ALLOW_AUTODIFF_TAMC
120 CADJ STORE dSigmaDrLtd(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
121 CADJ STORE dSigmaDrW(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
122 #endif
123 DO j=1-Oly,sNy+Oly
124 DO i=1-Olx+1,sNx+Olx
125 IF (dSigmaDrW(i,j).GE.dSigmaDrLtd(i,j))
126 & dSigmaDrW(i,j) = dSigmaDrLtd(i,j)
127 ENDDO
128 ENDDO
129 #ifdef ALLOW_AUTODIFF_TAMC
130 CADJ STORE dSigmaDrW(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
131 #endif
132 DO j=1-Oly,sNy+Oly
133 DO i=1-Olx+1,sNx+Olx
134 SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)
135 taperX(i,j) = 1. _d 0
136 ENDDO
137 ENDDO
138
139 C-- Y-comp
140
141 #ifdef ALLOW_AUTODIFF_TAMC
142 DO j=1-Oly+1,sNy+Oly
143 DO i=1-Olx,sNx+Olx
144 dSigmaDrLtd(i,j) = 0. _d 0
145 ENDDO
146 ENDDO
147 #endif /* ALLOW_AUTODIFF_TAMC */
148 DO j=1-Oly+1,sNy+Oly
149 DO i=1-Olx,sNx+Olx
150 dSigmaDrLtd(i,j) = -(GM_Small_Number+
151 & ABS(SlopeY(i,j))*GM_rMaxSlope)
152 ENDDO
153 ENDDO
154 #ifdef ALLOW_AUTODIFF_TAMC
155 CADJ STORE dSigmaDrLtd(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
156 CADJ STORE dSigmaDrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
157 #endif
158 DO j=1-Oly+1,sNy+Oly
159 DO i=1-Olx,sNx+Olx
160 IF (dSigmaDrS(i,j).GE.dSigmaDrLtd(i,j))
161 & dSigmaDrS(i,j) = dSigmaDrLtd(i,j)
162 ENDDO
163 ENDDO
164 #ifdef ALLOW_AUTODIFF_TAMC
165 CADJ STORE dSigmaDrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
166 #endif
167 DO j=1-Oly+1,sNy+Oly
168 DO i=1-Olx,sNx+Olx
169 SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)
170 taperY(i,j) = 1. _d 0
171 ENDDO
172 ENDDO
173
174 #endif /* GM_EXCLUDE_CLIPPING */
175
176 ELSEIF (GM_taper_scheme.EQ.'fm07') THEN
177
178 STOP 'GMREDI_SLOPE_PSI: AdvForm not yet implemented for fm07'
179
180 ELSE
181
182 #ifdef GM_EXCLUDE_TAPERING
183
184 STOP 'Need to compile without "#define GM_EXCLUDE_TAPERING"'
185
186 #else /* GM_EXCLUDE_TAPERING */
187
188 #ifdef ALLOW_AUTODIFF_TAMC
189 CADJ STORE slopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
190 CADJ STORE dSigmaDrW(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
191 #endif
192
193 C- Compute the slope, no clipping, but avoid reverse slope in negatively
194 C stratified (Sigma_Z > 0) region :
195 DO j=1-Oly,sNy+Oly
196 DO i=1-Olx+1,sNx+Olx
197 IF (dSigmaDrW(i,j).GE.-GM_Small_Number)
198 & dSigmaDrW(i,j) = -GM_Small_Number
199 ENDDO
200 ENDDO
201 #ifdef ALLOW_AUTODIFF_TAMC
202 CADJ STORE dsigmadrW(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
203 #endif
204 DO j=1-Oly,sNy+Oly
205 DO i=1-Olx+1,sNx+Olx
206 SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)
207 taperX(i,j) = 1. _d 0
208 ENDDO
209 ENDDO
210 #ifdef ALLOW_AUTODIFF_TAMC
211 CADJ STORE slopex(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
212 #endif
213 DO j=1-Oly,sNy+Oly
214 DO i=1-Olx+1,sNx+Olx
215 IF ( ABS(SlopeX(i,j)) .GE. slopeCutoff ) THEN
216 SlopeX(i,j) = SIGN(slopeCutoff,SlopeX(i,j))
217 taperX(i,j) = 0. _d 0
218 ENDIF
219 ENDDO
220 ENDDO
221
222 #ifdef ALLOW_AUTODIFF_TAMC
223 CADJ STORE slopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
224 CADJ STORE dSigmaDrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
225 #endif
226
227 DO j=1-Oly+1,sNy+Oly
228 DO i=1-Olx,sNx+Olx
229 IF (dSigmaDrS(i,j).GE.-GM_Small_Number)
230 & dSigmaDrS(i,j) = -GM_Small_Number
231 ENDDO
232 ENDDO
233 #ifdef ALLOW_AUTODIFF_TAMC
234 CADJ STORE dsigmadrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
235 #endif
236 DO j=1-Oly+1,sNy+Oly
237 DO i=1-Olx,sNx+Olx
238 SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)
239 taperY(i,j) = 1. _d 0
240 ENDDO
241 ENDDO
242 #ifdef ALLOW_AUTODIFF_TAMC
243 CADJ STORE slopey(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
244 #endif
245 DO j=1-Oly+1,sNy+Oly
246 DO i=1-Olx,sNx+Olx
247 IF ( ABS(SlopeY(i,j)) .GE. slopeCutoff ) THEN
248 SlopeY(i,j) = SIGN(slopeCutoff,SlopeY(i,j))
249 taperY(i,j) = 0. _d 0
250 ENDIF
251 ENDDO
252 ENDDO
253
254 C- Compute the tapering function for the GM+Redi tensor :
255
256 #ifdef ALLOW_AUTODIFF_TAMC
257 CADJ STORE slopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
258 CADJ STORE slopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
259 #endif
260
261 IF (GM_taper_scheme.EQ.'linear') THEN
262
263 C- Simplest adiabatic tapering = Smax/Slope (linear)
264 DO j=1-Oly,sNy+Oly
265 DO i=1-Olx+1,sNx+Olx
266 Smod = ABS(SlopeX(i,j))
267 IF ( Smod .GT. GM_maxSlope .AND.
268 & Smod .LT. slopeCutoff )
269 & taperX(i,j)=GM_maxSlope/(Smod+GM_Small_Number)
270 ENDDO
271 ENDDO
272 DO j=1-Oly+1,sNy+Oly
273 DO i=1-Olx,sNx+Olx
274 Smod = ABS(SlopeY(i,j))
275 IF ( Smod .GT. GM_maxSlope .AND.
276 & Smod .LT. slopeCutoff )
277 & taperY(i,j)=GM_maxSlope/(Smod+GM_Small_Number)
278 ENDDO
279 ENDDO
280
281 ELSEIF (GM_taper_scheme.EQ.'gkw91') THEN
282
283 C- Gerdes, Koberle and Willebrand, Clim. Dyn. 1991
284 maxSlopeSqr = GM_maxSlope*GM_maxSlope
285 DO j=1-Oly,sNy+Oly
286 DO i=1-Olx+1,sNx+Olx
287 IF ( ABS(SlopeX(i,j)) .GT. GM_maxSlope .AND.
288 & ABS(SlopeX(i,j)) .LT. slopeCutoff )
289 & taperX(i,j)=maxSlopeSqr/
290 & ( SlopeX(i,j)*SlopeX(i,j) + GM_Small_Number )
291 ENDDO
292 ENDDO
293 DO j=1-Oly+1,sNy+Oly
294 DO i=1-Olx,sNx+Olx
295 IF ( ABS(SlopeY(i,j)) .GT. GM_maxSlope .AND.
296 & ABS(SlopeY(i,j)) .LT. slopeCutoff )
297 & taperY(i,j)=maxSlopeSqr/
298 & ( SlopeY(i,j)*SlopeY(i,j) + GM_Small_Number )
299 ENDDO
300 ENDDO
301
302 ELSEIF (GM_taper_scheme.EQ.'dm95') THEN
303
304 C- Danabasoglu and McWilliams, J. Clim. 1995
305 DO j=1-Oly,sNy+Oly
306 DO i=1-Olx+1,sNx+Olx
307 Smod = ABS(SlopeX(i,j))
308 taperX(i,j)=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
309 ENDDO
310 ENDDO
311 DO j=1-Oly+1,sNy+Oly
312 DO i=1-Olx,sNx+Olx
313 Smod = ABS(SlopeY(i,j))
314 taperY(i,j)=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
315 ENDDO
316 ENDDO
317
318 ELSEIF (GM_taper_scheme.EQ.'ldd97') THEN
319
320 C- Large, Danabasoglu and Doney, JPO 1997
321
322 DO j=1-Oly,sNy+Oly
323 DO i=1-Olx+1,sNx+Olx
324 Smod = ABS(SlopeX(i,j))
325 IF ( Smod .LT. slopeCutoff ) THEN
326 f1=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
327 IF (Smod.NE.0.) THEN
328 Rnondim = -depthZ/(LrhoW(i,j)*Smod)
329 ELSE
330 Rnondim = 1.
331 ENDIF
332 IF ( Rnondim.GE.1. _d 0 ) THEN
333 f2 = 1. _d 0
334 ELSE
335 f2 = op5*( 1. _d 0 + SIN( fpi*(Rnondim-op5) ))
336 ENDIF
337 taperX(i,j)=f1*f2
338 ENDIF
339 ENDDO
340 ENDDO
341
342 DO j=1-Oly+1,sNy+Oly
343 DO i=1-Olx,sNx+Olx
344 Smod = ABS(SlopeY(i,j))
345 IF ( Smod .LT. slopeCutoff ) THEN
346 f1=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
347 IF (Smod.NE.0.) THEN
348 Rnondim = -depthZ/(LrhoS(i,j)*Smod)
349 ELSE
350 Rnondim = 1.
351 ENDIF
352 IF ( Rnondim.GE.1. _d 0 ) THEN
353 f2 = 1. _d 0
354 ELSE
355 f2 = op5*( 1. _d 0 + SIN( fpi*(Rnondim-op5) ))
356 ENDIF
357 taperY(i,j)=f1*f2
358 ENDIF
359 ENDDO
360 ENDDO
361
362 ELSEIF (GM_taper_scheme.NE.' ') THEN
363 STOP 'GMREDI_SLOPE_PSI: Bad GM_taper_scheme'
364 ENDIF
365
366 #endif /* GM_EXCLUDE_TAPERING */
367
368 ENDIF
369
370
371 #else /* GMREDI_WITH_STABLE_ADJOINT */
372 c special choice for adjoint/optimization of parameters
373 c (~ strong clipping, reducing non linearity of psi=f(K))
374
375
376 #ifdef ALLOW_AUTODIFF_TAMC
377 CADJ STORE slopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
378 CADJ STORE dSigmaDrW(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
379 #endif
380 DO j=1-Oly,sNy+Oly
381 DO i=1-Olx+1,sNx+Olx
382 IF (dSigmaDrW(i,j).GE.-GM_Small_Number)
383 & dSigmaDrW(i,j) = -GM_Small_Number
384 ENDDO
385 ENDDO
386 #ifdef ALLOW_AUTODIFF_TAMC
387 CADJ STORE dsigmadrW(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
388 #endif
389 DO j=1-Oly,sNy+Oly
390 DO i=1-Olx+1,sNx+Olx
391 SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)
392 taperX(i,j) = 1. _d 0
393 ENDDO
394 ENDDO
395
396 #ifdef ALLOW_AUTODIFF_TAMC
397 CADJ STORE slopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
398 CADJ STORE dSigmaDrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
399 #endif
400
401 DO j=1-Oly+1,sNy+Oly
402 DO i=1-Olx,sNx+Olx
403 IF (dSigmaDrS(i,j).GE.-GM_Small_Number)
404 & dSigmaDrS(i,j) = -GM_Small_Number
405 ENDDO
406 ENDDO
407 #ifdef ALLOW_AUTODIFF_TAMC
408 CADJ STORE dsigmadrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
409 #endif
410 DO j=1-Oly+1,sNy+Oly
411 DO i=1-Olx,sNx+Olx
412 SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)
413 taperY(i,j) = 1. _d 0
414 ENDDO
415 ENDDO
416
417 slopeMaxSpec=1. _d -4
418
419 CADJ STORE slopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
420 CADJ STORE slopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
421
422 DO j=1-Oly,sNy+Oly
423 DO i=1-Olx+1,sNx+Olx
424 slopeTmpSpec=ABS(SlopeX(i,j))
425 IF ( slopeTmpSpec .GT. slopeMaxSpec ) then
426 SlopeX(i,j)=5.*SlopeX(i,j)*slopeMaxSpec/slopeTmpSpec
427 ELSE
428 SlopeX(i,j)=5.*SlopeX(i,j)
429 ENDIF
430 taperX(i,j)=1.
431 ENDDO
432 ENDDO
433 DO j=1-Oly+1,sNy+Oly
434 DO i=1-Olx,sNx+Olx
435 slopeTmpSpec=ABS(SlopeY(i,j))
436 IF ( slopeTmpSpec .GT. slopeMaxSpec ) then
437 SlopeY(i,j)=5.*SlopeY(i,j)*slopeMaxSpec/slopeTmpSpec
438 ELSE
439 SlopeY(i,j)=5.*SlopeY(i,j)
440 ENDIF
441 taperY(i,j)=1.
442 ENDDO
443 ENDDO
444
445 #endif /* GMREDI_WITH_STABLE_ADJOINT */
446
447 #endif /* BOLUS_ADVEC */
448 #endif /* ALLOW_GMREDI */
449
450 RETURN
451 END

  ViewVC Help
Powered by ViewVC 1.1.22