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

Annotation of /MITgcm/pkg/gmredi/gmredi_slope_limit.F

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


Revision 1.17 - (hide annotations) (download)
Sun Jan 12 21:35:27 2003 UTC (21 years, 5 months ago) by jmc
Branch: MAIN
Changes since 1.16: +29 -19 lines
fix few bugs and restore parameter value (e.g., Small_Number=1.D-12)
and scheme (e.g., Large_SlopeSqr=1.D+48) of checkpoint47f_post

1 jmc 1.17 C $Header: $
2     C $Name: $
3 adcroft 1.1
4     #include "GMREDI_OPTIONS.h"
5    
6     CStartOfInterface
7     SUBROUTINE GMREDI_SLOPE_LIMIT(
8 jmc 1.8 I dSigmaDrReal,
9 heimbach 1.13 I depthZ,K,
10 adcroft 1.1 U SlopeX, SlopeY,
11 heimbach 1.13 U dSigmaDx, dSigmaDy,
12 jmc 1.8 O SlopeSqr, taperFct,
13 adcroft 1.1 I bi,bj, myThid )
14     C /==========================================================\
15     C | SUBROUTINE GMREDI_SLOPE_LIMIT |
16     C | o Calculate slopes for use in GM/Redi tensor |
17     C |==========================================================|
18     C | On entry: |
19 jmc 1.8 C | dSigmaDrReal conatins the d/dz Sigma |
20 adcroft 1.1 C | SlopeX/Y contains X/Y gradients of sigma |
21     C | depthZ conatins the height (m) of level |
22     C | On exit: |
23 jmc 1.8 C | dSigmaDrReal conatins the effective dSig/dz |
24 adcroft 1.1 C | SlopeX/Y contains X/Y slopes |
25 jmc 1.8 C | SlopeSqr contains Sx^2+Sy^2 |
26     C | taperFct contains tapering funct. value ; |
27     C | = 1 when using no tapering |
28 adcroft 1.1 C \==========================================================/
29     IMPLICIT NONE
30    
31     C == Global variables ==
32     #include "SIZE.h"
33 heimbach 1.12 #include "GRID.h"
34 adcroft 1.1 #include "EEPARAMS.h"
35     #include "GMREDI.h"
36     #include "PARAMS.h"
37 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
38     #include "tamc.h"
39     #include "tamc_keys.h"
40     #endif /* ALLOW_AUTODIFF_TAMC */
41 adcroft 1.1
42     C == Routine arguments ==
43     C
44     _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
45     _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
46 heimbach 1.13 _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47     _RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48 adcroft 1.1 _RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49 jmc 1.8 _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50     _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51 adcroft 1.1 _RL depthZ
52 heimbach 1.13 INTEGER bi,bj,K,myThid
53 adcroft 1.1 CEndOfInterface
54    
55     #ifdef ALLOW_GMREDI
56    
57     C == Local variables ==
58     _RL Small_Number
59 heimbach 1.12 _RL Small_Taper
60 heimbach 1.15 _RL Large_SlopeSqr
61 jmc 1.17 PARAMETER(Small_Number=1.D-12)
62 heimbach 1.12 PARAMETER(Small_Taper=1.D+03)
63 jmc 1.17 PARAMETER(Large_SlopeSqr=1.D+48)
64 heimbach 1.15
65 heimbach 1.10 _RL gradSmod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
66     _RL dSigmaDrLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
67     _RL dRdSigmaLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
68 heimbach 1.15 _RL tmpfld(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
69 heimbach 1.10 _RL f1,Smod,f2,Rnondim,Cspd,Lrho
70 jmc 1.8 _RL maxSlopeSqr
71 adcroft 1.1 _RL fpi
72 heimbach 1.6 PARAMETER(fpi=3.141592653589793047592d0)
73 adcroft 1.1 INTEGER i,j
74 jmc 1.17 c Small_Number=GM_Small_Number
75 adcroft 1.1
76 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
77 heimbach 1.12 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     ikey = (act1 + 1) + act2*max1
85     & + act3*max1*max2
86     & + act4*max1*max2*max3
87 heimbach 1.10 #endif /* ALLOW_AUTODIFF_TAMC */
88    
89 heimbach 1.12 DO j=1-Oly+1,sNy+Oly-1
90     DO i=1-Olx+1,sNx+Olx-1
91     gradSmod(i,j) = 0. _d 0
92     dSigmaDrLtd(i,j) = 0. _d 0
93 heimbach 1.15 tmpfld(i,j) = 0. _d 0
94 heimbach 1.12 ENDDO
95     ENDDO
96    
97 jmc 1.8 IF (GM_taper_scheme.EQ.'orig' .OR.
98     & GM_taper_scheme.EQ.'clipping') THEN
99 adcroft 1.1
100 jmc 1.17 #ifdef GM_TAPER_ORIG_CLIPPING
101 heimbach 1.13
102 adcroft 1.1 C- Original implementation in mitgcmuv
103     C (this turns out to be the same as Cox slope clipping)
104 jmc 1.8
105     C- Cox 1987 "Slope clipping"
106 adcroft 1.1 DO j=1-Oly+1,sNy+Oly-1
107     DO i=1-Olx+1,sNx+Olx-1
108 jmc 1.17 tmpfld(i,j) = dSigmaDx(i,j)*dSigmaDx(i,j) +
109     & dSigmaDy(i,j)*dSigmaDy(i,j)
110 heimbach 1.15 IF ( tmpfld(i,j) .EQ. 0. ) THEN
111 heimbach 1.12 gradSmod(i,j) = 0. _d 0
112     ELSE
113 heimbach 1.15 gradSmod(i,j) = sqrt( tmpfld(i,j) )
114 heimbach 1.12 ENDIF
115 heimbach 1.10 ENDDO
116     ENDDO
117 adcroft 1.1
118 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
119 heimbach 1.13 cnostore CADJ STORE gradSmod(:,:) = comlev1_bibj, key=ikey, byte=isbyte
120     cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
121 heimbach 1.10 #endif
122    
123     DO j=1-Oly+1,sNy+Oly-1
124     DO i=1-Olx+1,sNx+Olx-1
125 heimbach 1.12 IF (gradSmod(i,j) .NE. 0.) THEN
126     dSigmaDrLtd(i,j) = -gradSmod(i,j)*GM_rMaxSlope
127 heimbach 1.15 IF ( dSigmaDrReal(i,j) .GE. dSigmaDrLtd(i,j) )
128     & dSigmaDrReal(i,j) = dSigmaDrLtd(i,j)
129 heimbach 1.12 ENDIF
130 heimbach 1.10 ENDDO
131     ENDDO
132    
133     #ifdef ALLOW_AUTODIFF_TAMC
134 heimbach 1.13 cnostore CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
135     cnostore CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
136     cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
137 heimbach 1.10 #endif
138    
139     DO j=1-Oly+1,sNy+Oly-1
140     DO i=1-Olx+1,sNx+Olx-1
141 heimbach 1.12 IF (gradSmod(i,j) .EQ. 0.) THEN
142     SlopeX(i,j) = 0. _d 0
143     SlopeY(i,j) = 0. _d 0
144     ELSE
145 heimbach 1.13 dRdSigmaLtd(i,j) = 1./( dSigmaDrReal(i,j) )
146 jmc 1.17 SlopeX(i,j)=-dSigmaDx(i,j)*dRdSigmaLtd(i,j)
147     SlopeY(i,j)=-dSigmaDy(i,j)*dRdSigmaLtd(i,j)
148 heimbach 1.12 ENDIF
149 heimbach 1.10 ENDDO
150     ENDDO
151 adcroft 1.1
152 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
153 heimbach 1.13 cnostore CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
154     cnostore CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
155 heimbach 1.10 #endif
156 jmc 1.8
157 heimbach 1.10 DO j=1-Oly+1,sNy+Oly-1
158     DO i=1-Olx+1,sNx+Olx-1
159 jmc 1.8 SlopeSqr(i,j)=SlopeX(i,j)*SlopeX(i,j)
160     & +SlopeY(i,j)*SlopeY(i,j)
161     taperFct(i,j)=1. _d 0
162 adcroft 1.1 ENDDO
163     ENDDO
164    
165 jmc 1.17 #else /* GM_TAPER_ORIG_CLIPPING */
166 heimbach 1.13
167 jmc 1.17 STOP 'Need to compile with "#define GM_TAPER_ORIG_CLIPPING"'
168 heimbach 1.13
169     #endif /* GM_TAPER_ORIG_CLIPPING */
170    
171     ELSE IF (GM_taper_scheme.EQ.'ac02') THEN
172    
173 jmc 1.17 #ifdef GM_TAPER_AC02
174 heimbach 1.13
175 heimbach 1.15 maxSlopeSqr = GM_maxSlope*GM_maxSlope
176 heimbach 1.13 DO j=1-Oly+1,sNy+Oly-1
177     DO i=1-Olx+1,sNx+Olx-1
178     dRdSigmaLtd(i,j)=
179 heimbach 1.15 & dSigmaDx(i,j)*dSigmaDx(i,j)
180 heimbach 1.13 & + dSigmaDy(i,j)*dSigmaDy(i,j)
181 heimbach 1.15 & + dSigmaDrReal(i,j)*dSigmaDrReal(i,j)
182     taperFct(i,j) = 1. _d 0
183     c
184 heimbach 1.13 IF (dRdSigmaLtd(i,j).NE.0.) then
185     dRdSigmaLtd(i,j)=1. _d 0
186     & / ( dRdSigmaLtd(i,j) )
187     SlopeSqr(i,j)=(dSigmaDx(i,j)*dSigmaDx(i,j)
188     & +dSigmaDy(i,j)*dSigmaDy(i,j))*dRdSigmaLtd(i,j)
189     SlopeX(i,j)=-dSigmaDx(i,j)
190     & *dRdSigmaLtd(i,j)*dSigmaDrReal(i,j)
191     SlopeY(i,j)=-dSigmaDy(i,j)
192     & *dRdSigmaLtd(i,j)*dSigmaDrReal(i,j)
193     cph T11(i,j)=(dSigmaDrReal(i,j)**2)*dRdSigmaLtd(i,j)
194     ENDIF
195 heimbach 1.15 cph-- this part doesn't adjoint well
196     cph IF ( SlopeSqr(i,j) .GT. maxSlopeSqr .AND.
197     cph & SlopeSqr(i,j) .LT. Large_SlopeSqr ) THEN
198     cph taperFct(i,j) = maxSlopeSqr/SlopeSqr(i,j)
199     cph ELSE IF ( SlopeSqr(i,j) .GT. Large_SlopeSqr ) THEN
200     cph taperFct(i,j) = 0. _d 0
201 heimbach 1.13 cph ENDIF
202     ENDDO
203     ENDDO
204    
205 jmc 1.17 #else /* GM_TAPER_AC02 */
206 heimbach 1.13
207 jmc 1.17 STOP 'Need to compile with "#define GM_TAPER_AC02"'
208 heimbach 1.13
209     #endif /* GM_TAPER_AC02 */
210    
211 jmc 1.8 ELSE
212 adcroft 1.1
213 jmc 1.17 #ifdef GM_TAPER_REST
214 heimbach 1.13
215 heimbach 1.12 C----------------------------------------------------------------------
216    
217 heimbach 1.10 C- Compute the slope, no clipping, but avoid reverse slope in negatively
218 jmc 1.8 C stratified (Sigma_Z > 0) region :
219 heimbach 1.10
220     #ifdef ALLOW_AUTODIFF_TAMC
221 heimbach 1.13 cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
222 heimbach 1.10 #endif
223    
224     DO j=1-Oly+1,sNy+Oly-1
225     DO i=1-Olx+1,sNx+Olx-1
226 heimbach 1.12 IF ( dSigmaDrReal(i,j) .NE. 0. ) THEN
227     IF (dSigmaDrReal(i,j).GE.(-Small_Number))
228     & dSigmaDrReal(i,j) = -Small_Number
229     ENDIF
230 heimbach 1.10 ENDDO
231     ENDDO
232    
233     #ifdef ALLOW_AUTODIFF_TAMC
234 heimbach 1.13 cnostore CADJ STORE dSigmaDx(:,:) = comlev1_bibj, key=ikey, byte=isbyte
235     cnostore CADJ STORE dSigmaDy(:,:) = comlev1_bibj, key=ikey, byte=isbyte
236     cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
237 heimbach 1.10 #endif
238    
239 adcroft 1.1 DO j=1-Oly+1,sNy+Oly-1
240     DO i=1-Olx+1,sNx+Olx-1
241 heimbach 1.12 IF ( dSigmaDrReal(i,j) .EQ. 0. ) THEN
242 heimbach 1.14 IF ( dSigmaDx(i,j) .NE. 0. ) THEN
243     SlopeX(i,j) = SIGN(Small_taper,dSigmaDx(i,j))
244     ELSE
245     SlopeX(i,j) = 0. _d 0
246     ENDIF
247     IF ( dSigmaDy(i,j) .NE. 0. ) THEN
248     SlopeY(i,j) = SIGN(Small_taper,dSigmaDy(i,j))
249     ELSE
250     SlopeY(i,j) = 0. _d 0
251     ENDIF
252 heimbach 1.12 ELSE
253 jmc 1.17 dRdSigmaLtd(i,j) = 1. _d 0 / dSigmaDrReal(i,j)
254     SlopeX(i,j)=-dSigmaDx(i,j)*dRdSigmaLtd(i,j)
255     SlopeY(i,j)=-dSigmaDy(i,j)*dRdSigmaLtd(i,j)
256     c SlopeX(i,j) = -dSigmaDx(i,j)/dSigmaDrReal(i,j)
257     c SlopeY(i,j) = -dSigmaDy(i,j)/dSigmaDrReal(i,j)
258 heimbach 1.12 ENDIF
259 heimbach 1.10 ENDDO
260     ENDDO
261 adcroft 1.1
262 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
263 heimbach 1.13 cnostore CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
264     cnostore CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
265 heimbach 1.10 #endif
266 jmc 1.8
267 heimbach 1.10 DO j=1-Oly+1,sNy+Oly-1
268     DO i=1-Olx+1,sNx+Olx-1
269 heimbach 1.12 SlopeSqr(i,j) = SlopeX(i,j)*SlopeX(i,j)
270     & +SlopeY(i,j)*SlopeY(i,j)
271     taperFct(i,j) = 1. _d 0
272 heimbach 1.15 IF ( SlopeSqr(i,j) .GT. Large_SlopeSqr ) THEN
273     slopeSqr(i,j) = Large_SlopeSqr
274     taperFct(i,j) = 0. _d 0
275 jmc 1.17 ENDIF
276 jmc 1.8 ENDDO
277     ENDDO
278    
279     C- Compute the tapering function for the GM+Redi tensor :
280    
281     IF (GM_taper_scheme.EQ.'linear') THEN
282    
283     C- Simplest adiabatic tapering = Smax/Slope (linear)
284     maxSlopeSqr = GM_maxSlope*GM_maxSlope
285     DO j=1-Oly+1,sNy+Oly-1
286     DO i=1-Olx+1,sNx+Olx-1
287 adcroft 1.1
288 heimbach 1.12 IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
289     taperFct(i,j) = 1. _d 0
290 heimbach 1.15 ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr .AND.
291     & SlopeSqr(i,j) .LT. Large_SlopeSqr ) THEN
292 heimbach 1.12 taperFct(i,j) = sqrt(maxSlopeSqr / SlopeSqr(i,j))
293     ENDIF
294 adcroft 1.1
295     ENDDO
296     ENDDO
297    
298 jmc 1.8 ELSEIF (GM_taper_scheme.EQ.'gkw91') THEN
299 adcroft 1.1
300     C- Gerdes, Koberle and Willebrand, Clim. Dyn. 1991
301 jmc 1.8 maxSlopeSqr = GM_maxSlope*GM_maxSlope
302 adcroft 1.1 DO j=1-Oly+1,sNy+Oly-1
303     DO i=1-Olx+1,sNx+Olx-1
304    
305 heimbach 1.12 IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
306     taperFct(i,j) = 1. _d 0
307 heimbach 1.15 ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr .AND.
308     & SlopeSqr(i,j) .LT. Large_SlopeSqr ) THEN
309 heimbach 1.12 taperFct(i,j) = maxSlopeSqr/SlopeSqr(i,j)
310     ENDIF
311 adcroft 1.1
312     ENDDO
313     ENDDO
314    
315 jmc 1.8 ELSEIF (GM_taper_scheme.EQ.'dm95') THEN
316 adcroft 1.1
317     C- Danabasoglu and McWilliams, J. Clim. 1995
318     DO j=1-Oly+1,sNy+Oly-1
319     DO i=1-Olx+1,sNx+Olx-1
320    
321 heimbach 1.12 IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
322     taperFct(i,j) = 1. _d 0
323 heimbach 1.15 ELSE IF ( SlopeSqr(i,j) .LT. Large_SlopeSqr ) THEN
324 heimbach 1.12 Smod=sqrt(SlopeSqr(i,j))
325 heimbach 1.16 taperFct(i,j)=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))
326 heimbach 1.12 ENDIF
327 adcroft 1.1 ENDDO
328     ENDDO
329    
330 jmc 1.8 ELSEIF (GM_taper_scheme.EQ.'ldd97') THEN
331 adcroft 1.1
332     C- Large, Danabasoglu and Doney, JPO 1997
333     DO j=1-Oly+1,sNy+Oly-1
334     DO i=1-Olx+1,sNx+Olx-1
335    
336 heimbach 1.12 IF (SlopeSqr(i,j) .EQ. 0.) THEN
337     taperFct(i,j) = 1. _d 0
338 heimbach 1.15 ELSE IF ( SlopeSqr(i,j) .LT. Large_SlopeSqr ) THEN
339 heimbach 1.12 Smod=sqrt(SlopeSqr(i,j))
340 heimbach 1.16 f1=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))
341 jmc 1.17 Cspd=2. _d 0
342 heimbach 1.12 Lrho=100. _d 3
343 jmc 1.17 if (fCori(i,j,bi,bj).NE.0.) Lrho=Cspd/abs(fCori(i,j,bi,bj))
344 heimbach 1.12 Lrho=min(Lrho , 100. _d 3)
345     Lrho=max(Lrho , 15. _d 3)
346     Rnondim=depthZ/(Lrho*Smod)
347 heimbach 1.16 f2=op5*( 1. _d 0 + sin( fpi*(Rnondim-op5)))
348 heimbach 1.12 taperFct(i,j)=f1*f2
349     ENDIF
350 adcroft 1.1
351     ENDDO
352     ENDDO
353    
354 jmc 1.9 ELSEIF (GM_taper_scheme.NE.' ') THEN
355     STOP 'GMREDI_SLOPE_LIMIT: Bad GM_taper_scheme'
356 jmc 1.8 ENDIF
357 heimbach 1.13
358 jmc 1.17 #else /* GM_TAPER_REST */
359    
360     STOP 'Need to compile with "#define GM_TAPER_REST"'
361    
362     #endif /* GM_TAPER_REST */
363 adcroft 1.1
364 jmc 1.8 ENDIF
365 adcroft 1.1
366    
367     #endif /* ALLOW_GMREDI */
368    
369     RETURN
370     END

  ViewVC Help
Powered by ViewVC 1.1.22