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

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

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


Revision 1.7.2.3 - (show annotations) (download)
Mon Apr 8 20:27:13 2002 UTC (22 years, 2 months ago) by heimbach
Branch: release1
CVS Tags: release1_p13_pre, release1_p13, release1_p8, release1_p9, release1_p1, release1_p2, release1_p3, release1_p4, release1_p5, release1_p6, release1_p7, release1_p11, release1_p12, release1_p10, release1_p16, release1_p17, release1_p14, release1_p15, release1_p12_pre
Branch point for: release1_50yr
Changes since 1.7.2.2: +89 -61 lines
These changes constitute a patch to release1.
They were made on a development branch called "release1_final"
and are on the main trunk between chkpt44d_post and
checkpoint44h_post along with other changes.

This code is equivalent to chkpt44d_post with the following patches:
  - AD-related changes for GMRedi
  - fixes i KPP (delZ -> drF)
  - hook to OBCS songe layer code in external_forcing
  - errorMessageUnit non-zero in eeboot.F
  - modified test cost function and carbon verif.

1 C $Header$
2 C $Name$
3
4 #include "GMREDI_OPTIONS.h"
5
6 CStartOfInterface
7 SUBROUTINE GMREDI_SLOPE_LIMIT(
8 I dSigmaDrReal,
9 I depthZ,
10 U SlopeX, SlopeY,
11 O SlopeSqr, taperFct,
12 I bi,bj, myThid )
13 C /==========================================================\
14 C | SUBROUTINE GMREDI_SLOPE_LIMIT |
15 C | o Calculate slopes for use in GM/Redi tensor |
16 C |==========================================================|
17 C | On entry: |
18 C | dSigmaDrReal conatins the d/dz Sigma |
19 C | SlopeX/Y contains X/Y gradients of sigma |
20 C | depthZ conatins the height (m) of level |
21 C | On exit: |
22 C | dSigmaDrReal conatins the effective dSig/dz |
23 C | SlopeX/Y contains X/Y slopes |
24 C | SlopeSqr contains Sx^2+Sy^2 |
25 C | taperFct contains tapering funct. value ; |
26 C | = 1 when using no tapering |
27 C \==========================================================/
28 IMPLICIT NONE
29
30 C == Global variables ==
31 #include "SIZE.h"
32 #include "GRID.h"
33 #include "EEPARAMS.h"
34 #include "GMREDI.h"
35 #include "PARAMS.h"
36 #ifdef ALLOW_AUTODIFF_TAMC
37 #include "tamc.h"
38 #include "tamc_keys.h"
39 #endif /* ALLOW_AUTODIFF_TAMC */
40
41 C == Routine arguments ==
42 C
43 _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
44 _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
45 _RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
46 _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47 _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48 _RL depthZ
49 INTEGER bi,bj,myThid
50 CEndOfInterface
51
52 #ifdef ALLOW_GMREDI
53
54 C == Local variables ==
55 _RL Small_Number
56 _RL Small_Taper
57 PARAMETER(Small_Number=1.D-12)
58 PARAMETER(Small_Taper=1.D+03)
59 _RL gradSmod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60 _RL dSigmaDrLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
61 _RL dRdSigmaLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
62 _RL f1,Smod,f2,Rnondim,Cspd,Lrho
63 _RL maxSlopeSqr
64 _RL fpi
65 PARAMETER(fpi=3.141592653589793047592d0)
66 INTEGER i,j
67
68 #ifdef ALLOW_AUTODIFF_TAMC
69 act1 = bi - myBxLo(myThid)
70 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
71 act2 = bj - myByLo(myThid)
72 max2 = myByHi(myThid) - myByLo(myThid) + 1
73 act3 = myThid - 1
74 max3 = nTx*nTy
75 act4 = ikey_dynamics - 1
76 ikey = (act1 + 1) + act2*max1
77 & + act3*max1*max2
78 & + act4*max1*max2*max3
79 #endif /* ALLOW_AUTODIFF_TAMC */
80
81 DO j=1-Oly+1,sNy+Oly-1
82 DO i=1-Olx+1,sNx+Olx-1
83 gradSmod(i,j) = 0. _d 0
84 dSigmaDrLtd(i,j) = 0. _d 0
85 ENDDO
86 ENDDO
87
88 IF (GM_taper_scheme.EQ.'orig' .OR.
89 & GM_taper_scheme.EQ.'clipping') THEN
90
91 C- Original implementation in mitgcmuv
92 C (this turns out to be the same as Cox slope clipping)
93
94 C- Cox 1987 "Slope clipping"
95 DO j=1-Oly+1,sNy+Oly-1
96 DO i=1-Olx+1,sNx+Olx-1
97 IF ( SlopeX(i,j)*SlopeX(i,j) + SlopeY(i,j)*SlopeY(i,j)
98 & .EQ. 0. ) THEN
99 gradSmod(i,j) = 0. _d 0
100 ELSE
101 gradSmod(i,j) = sqrt(SlopeX(i,j)*SlopeX(i,j)
102 & +SlopeY(i,j)*SlopeY(i,j))
103 ENDIF
104 ENDDO
105 ENDDO
106
107 #ifdef ALLOW_AUTODIFF_TAMC
108 CADJ STORE gradSmod(:,:) = comlev1_bibj, key=ikey, byte=isbyte
109 CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
110 #endif
111
112 DO j=1-Oly+1,sNy+Oly-1
113 DO i=1-Olx+1,sNx+Olx-1
114 IF (gradSmod(i,j) .NE. 0.) THEN
115 dSigmaDrLtd(i,j) = -gradSmod(i,j)*GM_rMaxSlope
116 IF (dSigmaDrReal(i,j) .GE. dSigmaDrLtd(i,j))
117 & dSigmaDrReal(i,j) = dSigmaDrLtd(i,j)
118 ENDIF
119 ENDDO
120 ENDDO
121
122 #ifdef ALLOW_AUTODIFF_TAMC
123 CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
124 CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
125 CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
126 #endif
127
128 DO j=1-Oly+1,sNy+Oly-1
129 DO i=1-Olx+1,sNx+Olx-1
130 IF (gradSmod(i,j) .EQ. 0.) THEN
131 SlopeX(i,j) = 0. _d 0
132 SlopeY(i,j) = 0. _d 0
133 ELSE
134 dRdSigmaLtd(i,j) = 1./dSigmaDrReal(i,j)
135 SlopeX(i,j)=-SlopeX(i,j)*dRdSigmaLtd(i,j)
136 SlopeY(i,j)=-SlopeY(i,j)*dRdSigmaLtd(i,j)
137 ENDIF
138 ENDDO
139 ENDDO
140
141 #ifdef ALLOW_AUTODIFF_TAMC
142 CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
143 CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
144 #endif
145
146 DO j=1-Oly+1,sNy+Oly-1
147 DO i=1-Olx+1,sNx+Olx-1
148 SlopeSqr(i,j)=SlopeX(i,j)*SlopeX(i,j)
149 & +SlopeY(i,j)*SlopeY(i,j)
150 taperFct(i,j)=1. _d 0
151 ENDDO
152 ENDDO
153
154 ELSE
155
156 C----------------------------------------------------------------------
157
158 C- Compute the slope, no clipping, but avoid reverse slope in negatively
159 C stratified (Sigma_Z > 0) region :
160
161 #ifdef ALLOW_AUTODIFF_TAMC
162 CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
163 #endif
164
165 DO j=1-Oly+1,sNy+Oly-1
166 DO i=1-Olx+1,sNx+Olx-1
167 IF ( dSigmaDrReal(i,j) .NE. 0. ) THEN
168 IF (dSigmaDrReal(i,j).GE.(-Small_Number))
169 & dSigmaDrReal(i,j) = -Small_Number
170 ENDIF
171 ENDDO
172 ENDDO
173
174 #ifdef ALLOW_AUTODIFF_TAMC
175 CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
176 CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
177 CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
178 #endif
179
180 DO j=1-Oly+1,sNy+Oly-1
181 DO i=1-Olx+1,sNx+Olx-1
182 IF ( dSigmaDrReal(i,j) .EQ. 0. ) THEN
183 IF ( SlopeX(i,j) .NE. 0. )
184 & SlopeX(i,j) = SIGN(Small_taper,SlopeX(i,j))
185 IF ( SlopeY(i,j) .NE. 0. )
186 & SlopeY(i,j) = SIGN(Small_taper,SlopeY(i,j))
187 ELSE
188 dRdSigmaLtd(i,j) = 1./dSigmaDrReal(i,j)
189 SlopeX(i,j) = -SlopeX(i,j)*dRdSigmaLtd(i,j)
190 SlopeY(i,j) = -SlopeY(i,j)*dRdSigmaLtd(i,j)
191 ENDIF
192 ENDDO
193 ENDDO
194
195 #ifdef ALLOW_AUTODIFF_TAMC
196 CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
197 CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
198 #endif
199
200 DO j=1-Oly+1,sNy+Oly-1
201 DO i=1-Olx+1,sNx+Olx-1
202 SlopeSqr(i,j) = SlopeX(i,j)*SlopeX(i,j)
203 & +SlopeY(i,j)*SlopeY(i,j)
204 taperFct(i,j) = 1. _d 0
205 ENDDO
206 ENDDO
207
208 C- Compute the tapering function for the GM+Redi tensor :
209
210 IF (GM_taper_scheme.EQ.'linear') THEN
211
212 C- Simplest adiabatic tapering = Smax/Slope (linear)
213 maxSlopeSqr = GM_maxSlope*GM_maxSlope
214 DO j=1-Oly+1,sNy+Oly-1
215 DO i=1-Olx+1,sNx+Olx-1
216
217 IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
218 taperFct(i,j) = 1. _d 0
219 ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr ) THEN
220 taperFct(i,j) = sqrt(maxSlopeSqr / SlopeSqr(i,j))
221 ENDIF
222
223 ENDDO
224 ENDDO
225
226 ELSEIF (GM_taper_scheme.EQ.'gkw91') THEN
227
228 C- Gerdes, Koberle and Willebrand, Clim. Dyn. 1991
229 maxSlopeSqr = GM_maxSlope*GM_maxSlope
230 DO j=1-Oly+1,sNy+Oly-1
231 DO i=1-Olx+1,sNx+Olx-1
232
233 IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
234 taperFct(i,j) = 1. _d 0
235 ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr ) THEN
236 taperFct(i,j) = maxSlopeSqr/SlopeSqr(i,j)
237 ENDIF
238
239 ENDDO
240 ENDDO
241
242 ELSEIF (GM_taper_scheme.EQ.'dm95') THEN
243
244 C- Danabasoglu and McWilliams, J. Clim. 1995
245 DO j=1-Oly+1,sNy+Oly-1
246 DO i=1-Olx+1,sNx+Olx-1
247
248 IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
249 taperFct(i,j) = 1. _d 0
250 ELSE
251 Smod=sqrt(SlopeSqr(i,j))
252 taperFct(i,j)=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))
253 ENDIF
254 ENDDO
255 ENDDO
256
257 ELSEIF (GM_taper_scheme.EQ.'ldd97') THEN
258
259 C- Large, Danabasoglu and Doney, JPO 1997
260 DO j=1-Oly+1,sNy+Oly-1
261 DO i=1-Olx+1,sNx+Olx-1
262
263 IF (SlopeSqr(i,j) .EQ. 0.) THEN
264 taperFct(i,j) = 1. _d 0
265 ELSE
266 Smod=sqrt(SlopeSqr(i,j))
267 f1=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))
268 Cspd=2.
269 Lrho=100. _d 3
270 if (FCori(i,j,bi,bj).NE.0.) Lrho=Cspd/abs(Fcori(i,j,bi,bj))
271 Lrho=min(Lrho , 100. _d 3)
272 Lrho=max(Lrho , 15. _d 3)
273 Rnondim=depthZ/(Lrho*Smod)
274 f2=0.5*(1.+sin( fpi*(Rnondim-0.5)))
275 taperFct(i,j)=f1*f2
276 ENDIF
277
278 ENDDO
279 ENDDO
280
281 ELSEIF (GM_taper_scheme.NE.' ') THEN
282 STOP 'GMREDI_SLOPE_LIMIT: Bad GM_taper_scheme'
283 ENDIF
284
285 ENDIF
286
287
288 #endif /* ALLOW_GMREDI */
289
290 RETURN
291 END

  ViewVC Help
Powered by ViewVC 1.1.22