/[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.14 - (hide annotations) (download)
Fri Nov 15 02:57:47 2002 UTC (21 years, 6 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint47e_post, checkpoint47c_post, checkpoint47d_pre, checkpoint47a_post, checkpoint47d_post, branch-exfmods-tag, checkpoint47b_post, checkpoint47f_post, checkpoint47
Branch point for: branch-exfmods-curt
Changes since 1.13: +16 -6 lines
o * "clean" adjoint code (in terms of extensive recomputations)
    can now be obtained for all GMREDI options (i.e. for
    - GM_VISBECK_VARIABLE_K
    - GM_NON_UNITY_DIAGONAL
    - GM_EXTRA_DIAGONAL
    - GM_BOLUS_ADVEC )
  * However, wrong gradient check problem remains unsolved.
  * New CPP options have been introduced for different
    tapering schemes

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

  ViewVC Help
Powered by ViewVC 1.1.22