/[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.19 - (hide annotations) (download)
Tue Jan 21 19:34:13 2003 UTC (21 years, 4 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint50c_post, checkpoint48e_post, checkpoint50c_pre, checkpoint48i_post, checkpoint51, checkpoint50, checkpoint50d_post, checkpoint50b_pre, checkpoint48b_post, checkpoint51d_post, checkpoint48c_pre, checkpoint48d_pre, checkpoint48d_post, checkpoint48f_post, checkpoint48h_post, checkpoint51b_pre, checkpoint48a_post, checkpoint50f_post, checkpoint50a_post, checkpoint50f_pre, checkpoint47j_post, branchpoint-genmake2, checkpoint48c_post, checkpoint51b_post, checkpoint51c_post, checkpoint50g_post, checkpoint50h_post, checkpoint50e_pre, checkpoint50i_post, checkpoint50e_post, checkpoint50d_pre, checkpoint51e_post, checkpoint48, checkpoint49, checkpoint51f_pre, checkpoint48g_post, checkpoint50b_post, checkpoint51a_post
Branch point for: branch-genmake2
Changes since 1.18: +9 -7 lines
Yet more changes:
o adgmredi_calc_tensor
  avoiding all recomputation of gmredi_slope_limit
o adgmredi_x/y/rtransport
  added flag for excessive storing to avoid recomp. of
  u/v/rtans, dTdx/y/z
  -> this is not really necessary and very memory-consuming
o adgmredi_slope_psi:
  consistency with gmredi_slope_limit in treatment of GM_slopeSqCutoff
o gmredi_slope_limit
  re-activated full calculation of taperfct for case 'ac02'

1 heimbach 1.19 C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_slope_limit.F,v 1.18 2003/01/13 19:02:45 jmc Exp $
2 jmc 1.17 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 heimbach 1.12 _RL Small_Taper
59     PARAMETER(Small_Taper=1.D+03)
60 heimbach 1.15
61 heimbach 1.10 _RL gradSmod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
62     _RL dSigmaDrLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
63     _RL dRdSigmaLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
64 heimbach 1.15 _RL tmpfld(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
65 heimbach 1.10 _RL f1,Smod,f2,Rnondim,Cspd,Lrho
66 jmc 1.8 _RL maxSlopeSqr
67 adcroft 1.1 _RL fpi
68 heimbach 1.6 PARAMETER(fpi=3.141592653589793047592d0)
69 adcroft 1.1 INTEGER i,j
70    
71 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
72 heimbach 1.12 act1 = bi - myBxLo(myThid)
73     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
74     act2 = bj - myByLo(myThid)
75     max2 = myByHi(myThid) - myByLo(myThid) + 1
76     act3 = myThid - 1
77     max3 = nTx*nTy
78     act4 = ikey_dynamics - 1
79     ikey = (act1 + 1) + act2*max1
80     & + act3*max1*max2
81     & + act4*max1*max2*max3
82 heimbach 1.10 #endif /* ALLOW_AUTODIFF_TAMC */
83    
84 heimbach 1.12 DO j=1-Oly+1,sNy+Oly-1
85     DO i=1-Olx+1,sNx+Olx-1
86     gradSmod(i,j) = 0. _d 0
87     dSigmaDrLtd(i,j) = 0. _d 0
88 heimbach 1.15 tmpfld(i,j) = 0. _d 0
89 heimbach 1.12 ENDDO
90     ENDDO
91    
92 jmc 1.8 IF (GM_taper_scheme.EQ.'orig' .OR.
93     & GM_taper_scheme.EQ.'clipping') THEN
94 adcroft 1.1
95 jmc 1.18 #ifdef GM_EXCLUDE_CLIPPING
96    
97     STOP 'Need to compile without "#define GM_EXCLUDE_CLIPPING"'
98    
99     #else /* GM_EXCLUDE_CLIPPING */
100 heimbach 1.13
101 adcroft 1.1 C- Original implementation in mitgcmuv
102     C (this turns out to be the same as Cox slope clipping)
103 jmc 1.8
104     C- Cox 1987 "Slope clipping"
105 adcroft 1.1 DO j=1-Oly+1,sNy+Oly-1
106     DO i=1-Olx+1,sNx+Olx-1
107 jmc 1.17 tmpfld(i,j) = dSigmaDx(i,j)*dSigmaDx(i,j) +
108     & dSigmaDy(i,j)*dSigmaDy(i,j)
109 heimbach 1.15 IF ( tmpfld(i,j) .EQ. 0. ) THEN
110 heimbach 1.12 gradSmod(i,j) = 0. _d 0
111     ELSE
112 heimbach 1.15 gradSmod(i,j) = sqrt( tmpfld(i,j) )
113 heimbach 1.12 ENDIF
114 heimbach 1.10 ENDDO
115     ENDDO
116 adcroft 1.1
117 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
118 heimbach 1.13 cnostore CADJ STORE gradSmod(:,:) = comlev1_bibj, key=ikey, byte=isbyte
119     cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
120 heimbach 1.10 #endif
121    
122     DO j=1-Oly+1,sNy+Oly-1
123     DO i=1-Olx+1,sNx+Olx-1
124 heimbach 1.12 IF (gradSmod(i,j) .NE. 0.) THEN
125     dSigmaDrLtd(i,j) = -gradSmod(i,j)*GM_rMaxSlope
126 heimbach 1.15 IF ( dSigmaDrReal(i,j) .GE. dSigmaDrLtd(i,j) )
127     & dSigmaDrReal(i,j) = dSigmaDrLtd(i,j)
128 heimbach 1.12 ENDIF
129 heimbach 1.10 ENDDO
130     ENDDO
131    
132     #ifdef ALLOW_AUTODIFF_TAMC
133 heimbach 1.13 cnostore CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
134     cnostore CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
135     cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
136 heimbach 1.10 #endif
137    
138     DO j=1-Oly+1,sNy+Oly-1
139     DO i=1-Olx+1,sNx+Olx-1
140 heimbach 1.12 IF (gradSmod(i,j) .EQ. 0.) THEN
141     SlopeX(i,j) = 0. _d 0
142     SlopeY(i,j) = 0. _d 0
143     ELSE
144 heimbach 1.13 dRdSigmaLtd(i,j) = 1./( dSigmaDrReal(i,j) )
145 jmc 1.17 SlopeX(i,j)=-dSigmaDx(i,j)*dRdSigmaLtd(i,j)
146     SlopeY(i,j)=-dSigmaDy(i,j)*dRdSigmaLtd(i,j)
147 heimbach 1.12 ENDIF
148 heimbach 1.10 ENDDO
149     ENDDO
150 adcroft 1.1
151 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
152 heimbach 1.13 cnostore CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
153     cnostore CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
154 heimbach 1.10 #endif
155 jmc 1.8
156 heimbach 1.10 DO j=1-Oly+1,sNy+Oly-1
157     DO i=1-Olx+1,sNx+Olx-1
158 jmc 1.8 SlopeSqr(i,j)=SlopeX(i,j)*SlopeX(i,j)
159     & +SlopeY(i,j)*SlopeY(i,j)
160     taperFct(i,j)=1. _d 0
161 adcroft 1.1 ENDDO
162     ENDDO
163    
164 jmc 1.18 #endif /* GM_EXCLUDE_CLIPPING */
165 heimbach 1.13
166 jmc 1.18 ELSE IF (GM_taper_scheme.EQ.'ac02') THEN
167 heimbach 1.13
168 jmc 1.18 #ifdef GM_EXCLUDE_AC02_TAP
169 heimbach 1.13
170 jmc 1.18 STOP 'Need to compile without "#define GM_EXCLUDE_AC02_TAP"'
171    
172     #else /* GM_EXCLUDE_AC02_TAP */
173 heimbach 1.13
174 jmc 1.18 C- New Scheme (A. & C. 2002): relax part of the small slope approximation
175     C compute the true slope (no approximation)
176     C but still neglect Kxy & Kyx (assumed to be zero)
177 heimbach 1.13
178 heimbach 1.15 maxSlopeSqr = GM_maxSlope*GM_maxSlope
179 heimbach 1.13 DO j=1-Oly+1,sNy+Oly-1
180     DO i=1-Olx+1,sNx+Olx-1
181     dRdSigmaLtd(i,j)=
182 heimbach 1.15 & dSigmaDx(i,j)*dSigmaDx(i,j)
183 heimbach 1.13 & + dSigmaDy(i,j)*dSigmaDy(i,j)
184 heimbach 1.15 & + dSigmaDrReal(i,j)*dSigmaDrReal(i,j)
185     taperFct(i,j) = 1. _d 0
186     c
187 heimbach 1.13 IF (dRdSigmaLtd(i,j).NE.0.) then
188     dRdSigmaLtd(i,j)=1. _d 0
189     & / ( dRdSigmaLtd(i,j) )
190     SlopeSqr(i,j)=(dSigmaDx(i,j)*dSigmaDx(i,j)
191     & +dSigmaDy(i,j)*dSigmaDy(i,j))*dRdSigmaLtd(i,j)
192     SlopeX(i,j)=-dSigmaDx(i,j)
193     & *dRdSigmaLtd(i,j)*dSigmaDrReal(i,j)
194     SlopeY(i,j)=-dSigmaDy(i,j)
195     & *dRdSigmaLtd(i,j)*dSigmaDrReal(i,j)
196     cph T11(i,j)=(dSigmaDrReal(i,j)**2)*dRdSigmaLtd(i,j)
197     ENDIF
198 heimbach 1.19 #ifndef ALLOWW_AUTODIFF_TAMC
199 heimbach 1.15 cph-- this part doesn't adjoint well
200 heimbach 1.19 IF ( SlopeSqr(i,j) .GT. maxSlopeSqr .AND.
201     & SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN
202     taperFct(i,j) = maxSlopeSqr/SlopeSqr(i,j)
203     ELSE IF ( SlopeSqr(i,j) .GT. GM_slopeSqCutoff ) THEN
204     taperFct(i,j) = 0. _d 0
205     ENDIF
206     #endif
207 heimbach 1.13 ENDDO
208     ENDDO
209    
210 jmc 1.18 #endif /* GM_EXCLUDE_AC02_TAP */
211 heimbach 1.13
212 jmc 1.18 ELSE
213    
214     #ifdef GM_EXCLUDE_TAPERING
215 heimbach 1.13
216 jmc 1.18 STOP 'Need to compile without "#define GM_EXCLUDE_TAPERING"'
217 adcroft 1.1
218 jmc 1.18 #else /* GM_EXCLUDE_TAPERING */
219 heimbach 1.13
220 heimbach 1.12 C----------------------------------------------------------------------
221    
222 heimbach 1.10 C- Compute the slope, no clipping, but avoid reverse slope in negatively
223 jmc 1.8 C stratified (Sigma_Z > 0) region :
224 heimbach 1.10
225     #ifdef ALLOW_AUTODIFF_TAMC
226 heimbach 1.13 cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
227 heimbach 1.10 #endif
228    
229     DO j=1-Oly+1,sNy+Oly-1
230     DO i=1-Olx+1,sNx+Olx-1
231 heimbach 1.12 IF ( dSigmaDrReal(i,j) .NE. 0. ) THEN
232 jmc 1.18 IF (dSigmaDrReal(i,j).GE.(-GM_Small_Number))
233     & dSigmaDrReal(i,j) = -GM_Small_Number
234 heimbach 1.12 ENDIF
235 heimbach 1.10 ENDDO
236     ENDDO
237    
238     #ifdef ALLOW_AUTODIFF_TAMC
239 heimbach 1.13 cnostore CADJ STORE dSigmaDx(:,:) = comlev1_bibj, key=ikey, byte=isbyte
240     cnostore CADJ STORE dSigmaDy(:,:) = comlev1_bibj, key=ikey, byte=isbyte
241     cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
242 heimbach 1.10 #endif
243    
244 adcroft 1.1 DO j=1-Oly+1,sNy+Oly-1
245     DO i=1-Olx+1,sNx+Olx-1
246 heimbach 1.12 IF ( dSigmaDrReal(i,j) .EQ. 0. ) THEN
247 heimbach 1.14 IF ( dSigmaDx(i,j) .NE. 0. ) THEN
248     SlopeX(i,j) = SIGN(Small_taper,dSigmaDx(i,j))
249     ELSE
250     SlopeX(i,j) = 0. _d 0
251     ENDIF
252     IF ( dSigmaDy(i,j) .NE. 0. ) THEN
253     SlopeY(i,j) = SIGN(Small_taper,dSigmaDy(i,j))
254     ELSE
255     SlopeY(i,j) = 0. _d 0
256     ENDIF
257 heimbach 1.12 ELSE
258 jmc 1.17 dRdSigmaLtd(i,j) = 1. _d 0 / dSigmaDrReal(i,j)
259     SlopeX(i,j)=-dSigmaDx(i,j)*dRdSigmaLtd(i,j)
260     SlopeY(i,j)=-dSigmaDy(i,j)*dRdSigmaLtd(i,j)
261     c SlopeX(i,j) = -dSigmaDx(i,j)/dSigmaDrReal(i,j)
262     c SlopeY(i,j) = -dSigmaDy(i,j)/dSigmaDrReal(i,j)
263 heimbach 1.12 ENDIF
264 heimbach 1.10 ENDDO
265     ENDDO
266 adcroft 1.1
267 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
268 heimbach 1.13 cnostore CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
269     cnostore CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
270 heimbach 1.10 #endif
271 jmc 1.8
272 heimbach 1.10 DO j=1-Oly+1,sNy+Oly-1
273     DO i=1-Olx+1,sNx+Olx-1
274 heimbach 1.12 SlopeSqr(i,j) = SlopeX(i,j)*SlopeX(i,j)
275     & +SlopeY(i,j)*SlopeY(i,j)
276     taperFct(i,j) = 1. _d 0
277 jmc 1.18 IF ( SlopeSqr(i,j) .GT. GM_slopeSqCutoff ) THEN
278     slopeSqr(i,j) = GM_slopeSqCutoff
279 heimbach 1.15 taperFct(i,j) = 0. _d 0
280 jmc 1.17 ENDIF
281 jmc 1.8 ENDDO
282     ENDDO
283    
284     C- Compute the tapering function for the GM+Redi tensor :
285    
286     IF (GM_taper_scheme.EQ.'linear') THEN
287    
288     C- Simplest adiabatic tapering = Smax/Slope (linear)
289     maxSlopeSqr = GM_maxSlope*GM_maxSlope
290     DO j=1-Oly+1,sNy+Oly-1
291     DO i=1-Olx+1,sNx+Olx-1
292 adcroft 1.1
293 heimbach 1.12 IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
294     taperFct(i,j) = 1. _d 0
295 heimbach 1.15 ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr .AND.
296 jmc 1.18 & SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN
297 heimbach 1.12 taperFct(i,j) = sqrt(maxSlopeSqr / SlopeSqr(i,j))
298     ENDIF
299 adcroft 1.1
300     ENDDO
301     ENDDO
302    
303 jmc 1.8 ELSEIF (GM_taper_scheme.EQ.'gkw91') THEN
304 adcroft 1.1
305     C- Gerdes, Koberle and Willebrand, Clim. Dyn. 1991
306 jmc 1.8 maxSlopeSqr = GM_maxSlope*GM_maxSlope
307 adcroft 1.1 DO j=1-Oly+1,sNy+Oly-1
308     DO i=1-Olx+1,sNx+Olx-1
309    
310 heimbach 1.12 IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
311     taperFct(i,j) = 1. _d 0
312 heimbach 1.15 ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr .AND.
313 jmc 1.18 & SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN
314 heimbach 1.12 taperFct(i,j) = maxSlopeSqr/SlopeSqr(i,j)
315     ENDIF
316 adcroft 1.1
317     ENDDO
318     ENDDO
319    
320 jmc 1.8 ELSEIF (GM_taper_scheme.EQ.'dm95') THEN
321 adcroft 1.1
322     C- Danabasoglu and McWilliams, J. Clim. 1995
323     DO j=1-Oly+1,sNy+Oly-1
324     DO i=1-Olx+1,sNx+Olx-1
325    
326 heimbach 1.12 IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
327     taperFct(i,j) = 1. _d 0
328 jmc 1.18 ELSE IF ( SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN
329 heimbach 1.12 Smod=sqrt(SlopeSqr(i,j))
330 heimbach 1.16 taperFct(i,j)=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))
331 heimbach 1.12 ENDIF
332 adcroft 1.1 ENDDO
333     ENDDO
334    
335 jmc 1.8 ELSEIF (GM_taper_scheme.EQ.'ldd97') THEN
336 adcroft 1.1
337     C- Large, Danabasoglu and Doney, JPO 1997
338     DO j=1-Oly+1,sNy+Oly-1
339     DO i=1-Olx+1,sNx+Olx-1
340    
341 heimbach 1.12 IF (SlopeSqr(i,j) .EQ. 0.) THEN
342     taperFct(i,j) = 1. _d 0
343 jmc 1.18 ELSE IF ( SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN
344 heimbach 1.12 Smod=sqrt(SlopeSqr(i,j))
345 heimbach 1.16 f1=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))
346 jmc 1.17 Cspd=2. _d 0
347 heimbach 1.12 Lrho=100. _d 3
348 jmc 1.17 if (fCori(i,j,bi,bj).NE.0.) Lrho=Cspd/abs(fCori(i,j,bi,bj))
349 heimbach 1.12 Lrho=min(Lrho , 100. _d 3)
350     Lrho=max(Lrho , 15. _d 3)
351     Rnondim=depthZ/(Lrho*Smod)
352 heimbach 1.16 f2=op5*( 1. _d 0 + sin( fpi*(Rnondim-op5)))
353 heimbach 1.12 taperFct(i,j)=f1*f2
354     ENDIF
355 adcroft 1.1
356     ENDDO
357     ENDDO
358    
359 jmc 1.9 ELSEIF (GM_taper_scheme.NE.' ') THEN
360     STOP 'GMREDI_SLOPE_LIMIT: Bad GM_taper_scheme'
361 jmc 1.8 ENDIF
362 heimbach 1.13
363 jmc 1.18 #endif /* GM_EXCLUDE_TAPERING */
364 adcroft 1.1
365 jmc 1.8 ENDIF
366 adcroft 1.1
367    
368     #endif /* ALLOW_GMREDI */
369    
370     RETURN
371     END

  ViewVC Help
Powered by ViewVC 1.1.22