/[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.11 - (show annotations) (download)
Fri May 30 21:10:25 2008 UTC (15 years, 11 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint60, checkpoint61, checkpoint61f, checkpoint61n, checkpoint61q, checkpoint61e, checkpoint61g, checkpoint61d, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61z, checkpoint61x, checkpoint61y
Changes since 1.10: +82 -1 lines
GMREDI_WITH_STABLE_ADJOINT CPP option: special setup of gmredi
for which the adjoint is stable enough for parameter optimization.

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

  ViewVC Help
Powered by ViewVC 1.1.22