/[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.14 - (show annotations) (download)
Tue Sep 9 22:34:06 2014 UTC (9 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65l, checkpoint65m, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65o, HEAD
Changes since 1.13: +66 -65 lines
Include explicitly AUTODIFF_OPTIONS.h (in case we don't use ECCO_CPPOPTIONS.h)

1 C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_slope_psi.F,v 1.13 2012/09/28 21:12:25 gforget Exp $
2 C $Name: $
3
4 #include "GMREDI_OPTIONS.h"
5 #ifdef ALLOW_AUTODIFF
6 # include "AUTODIFF_OPTIONS.h"
7 #endif
8
9 CStartOfInterface
10 SUBROUTINE GMREDI_SLOPE_PSI(
11 O taperX, taperY,
12 U SlopeX, SlopeY,
13 U dSigmaDrW,dSigmaDrS,
14 I LrhoW, LrhoS, depthZ, K,
15 I bi,bj, myThid )
16 C /==========================================================\
17 C | SUBROUTINE GMREDI_SLOPE_PSI |
18 C | o Calculate slopes for use in GM/Redi tensor |
19 C |==========================================================|
20 C | On entry: |
21 C | dSigmaDrW,S contains the d/dz Sigma |
22 C | SlopeX/Y contains X/Y gradients of sigma |
23 C | depthZ contains the depth (< 0 !) [m] |
24 C | On exit: |
25 C | dSigmaDrW,S contains the effective dSig/dz |
26 C | SlopeX/Y contains X/Y slopes |
27 C | taperFct contains tapering funct. value ; |
28 C | = 1 when using no tapering |
29 C \==========================================================/
30 IMPLICIT NONE
31
32 C == Global variables ==
33 #include "SIZE.h"
34 #include "EEPARAMS.h"
35 #include "GMREDI.h"
36 #include "PARAMS.h"
37
38 #ifdef ALLOW_AUTODIFF_TAMC
39 #include "tamc.h"
40 #include "tamc_keys.h"
41 #endif /* ALLOW_AUTODIFF_TAMC */
42
43 C == Routine arguments ==
44 C
45 _RL taperX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46 _RL taperY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47 _RL SlopeX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48 _RL SlopeY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
49 _RL dSigmaDrW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
50 _RL dSigmaDrS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51 _RL LrhoW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52 _RL LrhoS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
53 _RS depthZ
54 INTEGER K,bi,bj,myThid
55 CEndOfInterface
56
57 #ifdef ALLOW_GMREDI
58 #ifdef GM_BOLUS_ADVEC
59
60 C == Local variables ==
61 _RL dSigmaDrLtd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62 _RL f1,Smod,f2,Rnondim
63 _RL maxSlopeSqr
64 _RL slopeCutoff
65 _RL fpi
66 PARAMETER(fpi=3.141592653589793047592d0)
67 INTEGER i,j
68 #ifdef GMREDI_WITH_STABLE_ADJOINT
69 _RL slopeTmpSpec,slopeMaxSpec
70 #endif
71
72 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
73
74 slopeCutoff = SQRT( GM_slopeSqCutoff )
75
76 #ifdef ALLOW_AUTODIFF_TAMC
77 act1 = bi - myBxLo(myThid)
78 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
79 act2 = bj - myByLo(myThid)
80 max2 = myByHi(myThid) - myByLo(myThid) + 1
81 act3 = myThid - 1
82 max3 = nTx*nTy
83 act4 = ikey_dynamics - 1
84 igmkey = (act1 + 1) + act2*max1
85 & + act3*max1*max2
86 & + act4*max1*max2*max3
87 kkey = (igmkey-1)*Nr + k
88 #endif /* ALLOW_AUTODIFF_TAMC */
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
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 */
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
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 */
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 IF (GM_taper_scheme.NE.'stableGmAdjTap') THEN
214 DO j=1-OLy,sNy+OLy
215 DO i=1-OLx+1,sNx+OLx
216 IF ( ABS(SlopeX(i,j)) .GE. slopeCutoff ) THEN
217 SlopeX(i,j) = SIGN(slopeCutoff,SlopeX(i,j))
218 taperX(i,j) = 0. _d 0
219 ENDIF
220 ENDDO
221 ENDDO
222 ENDIF
223
224 #ifdef ALLOW_AUTODIFF_TAMC
225 CADJ STORE slopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
226 CADJ STORE dSigmaDrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
227 #endif
228
229 DO j=1-OLy+1,sNy+OLy
230 DO i=1-OLx,sNx+OLx
231 IF (dSigmaDrS(i,j).GE.-GM_Small_Number)
232 & dSigmaDrS(i,j) = -GM_Small_Number
233 ENDDO
234 ENDDO
235 #ifdef ALLOW_AUTODIFF_TAMC
236 CADJ STORE dsigmadrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
237 #endif
238 DO j=1-OLy+1,sNy+OLy
239 DO i=1-OLx,sNx+OLx
240 SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)
241 taperY(i,j) = 1. _d 0
242 ENDDO
243 ENDDO
244 #ifdef ALLOW_AUTODIFF_TAMC
245 CADJ STORE slopey(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
246 #endif
247 IF (GM_taper_scheme.NE.'stableGmAdjTap') THEN
248 DO j=1-OLy+1,sNy+OLy
249 DO i=1-OLx,sNx+OLx
250 IF ( ABS(SlopeY(i,j)) .GE. slopeCutoff ) THEN
251 SlopeY(i,j) = SIGN(slopeCutoff,SlopeY(i,j))
252 taperY(i,j) = 0. _d 0
253 ENDIF
254 ENDDO
255 ENDDO
256 ENDIF
257
258 C- Compute the tapering function for the GM+Redi tensor :
259
260 #ifdef ALLOW_AUTODIFF_TAMC
261 CADJ STORE slopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
262 CADJ STORE slopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
263 #endif
264
265 IF (GM_taper_scheme.EQ.'linear') THEN
266
267 C- Simplest adiabatic tapering = Smax/Slope (linear)
268 DO j=1-OLy,sNy+OLy
269 DO i=1-OLx+1,sNx+OLx
270 Smod = ABS(SlopeX(i,j))
271 IF ( Smod .GT. GM_maxSlope .AND.
272 & Smod .LT. slopeCutoff )
273 & taperX(i,j)=GM_maxSlope/(Smod+GM_Small_Number)
274 ENDDO
275 ENDDO
276 DO j=1-OLy+1,sNy+OLy
277 DO i=1-OLx,sNx+OLx
278 Smod = ABS(SlopeY(i,j))
279 IF ( Smod .GT. GM_maxSlope .AND.
280 & Smod .LT. slopeCutoff )
281 & taperY(i,j)=GM_maxSlope/(Smod+GM_Small_Number)
282 ENDDO
283 ENDDO
284
285 ELSEIF (GM_taper_scheme.EQ.'gkw91') THEN
286
287 C- Gerdes, Koberle and Willebrand, Clim. Dyn. 1991
288 maxSlopeSqr = GM_maxSlope*GM_maxSlope
289 DO j=1-OLy,sNy+OLy
290 DO i=1-OLx+1,sNx+OLx
291 IF ( ABS(SlopeX(i,j)) .GT. GM_maxSlope .AND.
292 & ABS(SlopeX(i,j)) .LT. slopeCutoff )
293 & taperX(i,j)=maxSlopeSqr/
294 & ( SlopeX(i,j)*SlopeX(i,j) + GM_Small_Number )
295 ENDDO
296 ENDDO
297 DO j=1-OLy+1,sNy+OLy
298 DO i=1-OLx,sNx+OLx
299 IF ( ABS(SlopeY(i,j)) .GT. GM_maxSlope .AND.
300 & ABS(SlopeY(i,j)) .LT. slopeCutoff )
301 & taperY(i,j)=maxSlopeSqr/
302 & ( SlopeY(i,j)*SlopeY(i,j) + GM_Small_Number )
303 ENDDO
304 ENDDO
305
306 ELSEIF (GM_taper_scheme.EQ.'dm95') THEN
307
308 C- Danabasoglu and McWilliams, J. Clim. 1995
309 DO j=1-OLy,sNy+OLy
310 DO i=1-OLx+1,sNx+OLx
311 Smod = ABS(SlopeX(i,j))
312 taperX(i,j)=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
313 ENDDO
314 ENDDO
315 DO j=1-OLy+1,sNy+OLy
316 DO i=1-OLx,sNx+OLx
317 Smod = ABS(SlopeY(i,j))
318 taperY(i,j)=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
319 ENDDO
320 ENDDO
321
322 ELSEIF (GM_taper_scheme.EQ.'ldd97') THEN
323
324 C- Large, Danabasoglu and Doney, JPO 1997
325
326 DO j=1-OLy,sNy+OLy
327 DO i=1-OLx+1,sNx+OLx
328 Smod = ABS(SlopeX(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/(LrhoW(i,j)*Smod)
333 ELSE
334 Rnondim = 1.
335 ENDIF
336 IF ( Rnondim.GE.1. _d 0 ) THEN
337 f2 = 1. _d 0
338 ELSE
339 f2 = op5*( 1. _d 0 + SIN( fpi*(Rnondim-op5) ))
340 ENDIF
341 taperX(i,j)=f1*f2
342 ENDIF
343 ENDDO
344 ENDDO
345
346 DO j=1-OLy+1,sNy+OLy
347 DO i=1-OLx,sNx+OLx
348 Smod = ABS(SlopeY(i,j))
349 IF ( Smod .LT. slopeCutoff ) THEN
350 f1=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
351 IF (Smod.NE.0.) THEN
352 Rnondim = -depthZ/(LrhoS(i,j)*Smod)
353 ELSE
354 Rnondim = 1.
355 ENDIF
356 IF ( Rnondim.GE.1. _d 0 ) THEN
357 f2 = 1. _d 0
358 ELSE
359 f2 = op5*( 1. _d 0 + SIN( fpi*(Rnondim-op5) ))
360 ENDIF
361 taperY(i,j)=f1*f2
362 ENDIF
363 ENDDO
364 ENDDO
365
366 ELSEIF (GM_taper_scheme.EQ.'stableGmAdjTap') THEN
367
368 #ifndef GMREDI_WITH_STABLE_ADJOINT
369
370 STOP 'Need to compile wth "#define GMREDI_WITH_STABLE_ADJOINT"'
371
372 #else /* GMREDI_WITH_STABLE_ADJOINT */
373
374 c special choice for adjoint/optimization of parameters
375 c (~ strong clipping, reducing non linearity of psi=f(K))
376
377 slopeMaxSpec=1. _d -4
378
379 CADJ STORE slopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
380 CADJ STORE slopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
381
382 DO j=1-OLy,sNy+OLy
383 DO i=1-OLx+1,sNx+OLx
384 slopeTmpSpec=ABS(SlopeX(i,j))
385 IF ( slopeTmpSpec .GT. slopeMaxSpec ) then
386 SlopeX(i,j)=5.*SlopeX(i,j)*slopeMaxSpec/slopeTmpSpec
387 ELSE
388 SlopeX(i,j)=5.*SlopeX(i,j)
389 ENDIF
390 taperX(i,j)=1.
391 ENDDO
392 ENDDO
393 DO j=1-OLy+1,sNy+OLy
394 DO i=1-OLx,sNx+OLx
395 slopeTmpSpec=ABS(SlopeY(i,j))
396 IF ( slopeTmpSpec .GT. slopeMaxSpec ) then
397 SlopeY(i,j)=5.*SlopeY(i,j)*slopeMaxSpec/slopeTmpSpec
398 ELSE
399 SlopeY(i,j)=5.*SlopeY(i,j)
400 ENDIF
401 taperY(i,j)=1.
402 ENDDO
403 ENDDO
404 #endif /* GMREDI_WITH_STABLE_ADJOINT */
405
406 ELSEIF (GM_taper_scheme.NE.' ') THEN
407 STOP 'GMREDI_SLOPE_PSI: Bad GM_taper_scheme'
408 ENDIF
409
410 #endif /* GM_EXCLUDE_TAPERING */
411
412 ENDIF
413
414 #endif /* BOLUS_ADVEC */
415 #endif /* ALLOW_GMREDI */
416
417 RETURN
418 END

  ViewVC Help
Powered by ViewVC 1.1.22