1 |
C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_slope_limit.F,v 1.36 2015/12/29 14:32:40 gforget Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "GMREDI_OPTIONS.h" |
5 |
#ifdef ALLOW_AUTODIFF |
6 |
# include "AUTODIFF_OPTIONS.h" |
7 |
#endif |
8 |
|
9 |
CBOP |
10 |
C !ROUTINE: GMREDI_SLOPE_LIMIT |
11 |
C !INTERFACE: |
12 |
SUBROUTINE GMREDI_SLOPE_LIMIT( |
13 |
O SlopeX, SlopeY, |
14 |
O SlopeSqr, taperFct, |
15 |
U hTransLay, baseSlope, recipLambda, |
16 |
U dSigmaDr, |
17 |
I dSigmaDx, dSigmaDy, |
18 |
I Lrho, hMixLay, depthZ, kLow, |
19 |
I k, bi, bj, myTime, myIter, myThid ) |
20 |
|
21 |
C !DESCRIPTION: \bv |
22 |
C *==========================================================* |
23 |
C | SUBROUTINE GMREDI_SLOPE_LIMIT |
24 |
C | o Calculate slopes for use in GM/Redi tensor |
25 |
C *==========================================================* |
26 |
C | On entry: |
27 |
C | dSigmaDr contains the d/dz Sigma |
28 |
C | dSigmaDx/Dy contains X/Y gradients of sigma |
29 |
C | depthZ contains the depth (< 0 !) [m] |
30 |
C | Lrho |
31 |
C | hMixLay mixed layer depth (> 0) |
32 |
C U hTransLay transition layer depth (> 0) |
33 |
C U baseSlope, recipLambda, |
34 |
C | On exit: |
35 |
C | dSigmaDr contains the effective dSig/dz |
36 |
C | SlopeX/Y contains X/Y slopes |
37 |
C | SlopeSqr contains Sx^2+Sy^2 |
38 |
C | taperFct contains tapering funct. value ; |
39 |
C | = 1 when using no tapering |
40 |
C U hTransLay transition layer depth (> 0) |
41 |
C U baseSlope, recipLambda |
42 |
C *==========================================================* |
43 |
C \ev |
44 |
|
45 |
C !USES: |
46 |
IMPLICIT NONE |
47 |
|
48 |
C == Global variables == |
49 |
#include "SIZE.h" |
50 |
#include "GRID.h" |
51 |
#include "EEPARAMS.h" |
52 |
#include "GMREDI.h" |
53 |
#include "PARAMS.h" |
54 |
#ifdef ALLOW_AUTODIFF_TAMC |
55 |
#include "tamc.h" |
56 |
#include "tamc_keys.h" |
57 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
58 |
|
59 |
C == Routine arguments == |
60 |
C !INPUT PARAMETERS: |
61 |
C dSigmaDx :: zonal gradient of density |
62 |
C dSigmaDy :: meridional gradient of density |
63 |
C Lrho :: |
64 |
C hMixLay :: mixed layer thickness (> 0) [m] |
65 |
C depthZ :: model discretized depth (< 0) [m] |
66 |
C kLow :: bottom level index 2-D array |
67 |
C k :: level index |
68 |
C bi, bj :: tile indices |
69 |
C myTime :: time in simulation |
70 |
C myIter :: iteration number in simulation |
71 |
C myThid :: My Thread Id. number |
72 |
C !OUTPUT PARAMETERS: |
73 |
C SlopeX :: isopycnal slope in X direction |
74 |
C SlopeY :: isopycnal slope in Y direction |
75 |
C SlopeSqr :: square of isopycnal slope |
76 |
C taperFct :: tapering function |
77 |
C hTransLay :: depth of the base of the transition layer (> 0) [m] |
78 |
C baseSlope :: slope at the the base of the transition layer |
79 |
C recipLambda :: Slope vertical gradient at Trans. Layer Base (=recip.Lambda) |
80 |
C dSigmaDr :: vertical gradient of density |
81 |
_RL SlopeX (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
82 |
_RL SlopeY (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
83 |
_RL SlopeSqr (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
84 |
_RL taperFct (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
85 |
_RL hTransLay (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
86 |
_RL baseSlope (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
87 |
_RL recipLambda(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
88 |
_RL dSigmaDr (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
89 |
_RL dSigmaDx (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
90 |
_RL dSigmaDy (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
91 |
_RL Lrho (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
92 |
_RL hMixLay (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
93 |
_RS depthZ(*) |
94 |
INTEGER kLow (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
95 |
INTEGER k, bi,bj |
96 |
_RL myTime |
97 |
INTEGER myIter |
98 |
INTEGER myThid |
99 |
CEOP |
100 |
|
101 |
#ifdef ALLOW_GMREDI |
102 |
|
103 |
C !LOCAL VARIABLES: |
104 |
C == Local variables == |
105 |
#ifdef GMREDI_WITH_STABLE_ADJOINT |
106 |
_RL slopeSqTmp,slopeTmp,slopeMax |
107 |
#endif |
108 |
_RL dSigmMod(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
109 |
_RL SlopeMod(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
110 |
_RL dRdSigmaLtd(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
111 |
_RL tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
112 |
_RL locVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
113 |
_RL f1,Smod,f2,Rnondim |
114 |
_RL maxSlopeSqr |
115 |
_RL fpi |
116 |
PARAMETER( fpi = PI ) |
117 |
INTEGER i,j |
118 |
_RL dTransLay, rLambMin, DoverLamb |
119 |
_RL taperFctLoc, taperFctHat |
120 |
_RL minTransLay |
121 |
|
122 |
C- need to put this one in GM namelist : |
123 |
_RL GM_bigSlope |
124 |
GM_bigSlope = 1. _d +02 |
125 |
|
126 |
#ifdef ALLOW_AUTODIFF_TAMC |
127 |
C TAF thinks for some reason that depthZ is an active variable. |
128 |
C While this does not make the adjoint code wrong, the resulting |
129 |
C code inhibits vectorization in some cases so we tell TAF here |
130 |
C that depthZ is actually a passive grid variable that needs no adjoint |
131 |
CADJ PASSIVE depthz |
132 |
act1 = bi - myBxLo(myThid) |
133 |
max1 = myBxHi(myThid) - myBxLo(myThid) + 1 |
134 |
act2 = bj - myByLo(myThid) |
135 |
max2 = myByHi(myThid) - myByLo(myThid) + 1 |
136 |
act3 = myThid - 1 |
137 |
max3 = nTx*nTy |
138 |
act4 = ikey_dynamics - 1 |
139 |
ikey = (act1 + 1) + act2*max1 |
140 |
& + act3*max1*max2 |
141 |
& + act4*max1*max2*max3 |
142 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
143 |
|
144 |
DO j=1-OLy+1,sNy+OLy-1 |
145 |
DO i=1-OLx+1,sNx+OLx-1 |
146 |
dSigmMod(i,j) = 0. _d 0 |
147 |
tmpFld(i,j) = 0. _d 0 |
148 |
ENDDO |
149 |
ENDDO |
150 |
|
151 |
IF (GM_taper_scheme.EQ.'orig' .OR. |
152 |
& GM_taper_scheme.EQ.'clipping') THEN |
153 |
|
154 |
#ifdef GM_EXCLUDE_CLIPPING |
155 |
|
156 |
STOP 'Need to compile without "#define GM_EXCLUDE_CLIPPING"' |
157 |
|
158 |
#else /* GM_EXCLUDE_CLIPPING */ |
159 |
|
160 |
C- Original implementation in mitgcmuv |
161 |
C (this turns out to be the same as Cox slope clipping) |
162 |
|
163 |
C- Cox 1987 "Slope clipping" |
164 |
DO j=1-OLy+1,sNy+OLy-1 |
165 |
DO i=1-OLx+1,sNx+OLx-1 |
166 |
tmpFld(i,j) = dSigmaDx(i,j)*dSigmaDx(i,j) |
167 |
& + dSigmaDy(i,j)*dSigmaDy(i,j) |
168 |
IF ( tmpFld(i,j) .EQ. 0. ) THEN |
169 |
dSigmMod(i,j) = 0. _d 0 |
170 |
ELSE |
171 |
dSigmMod(i,j) = SQRT( tmpFld(i,j) ) |
172 |
ENDIF |
173 |
ENDDO |
174 |
ENDDO |
175 |
|
176 |
#ifdef ALLOW_AUTODIFF_TAMC |
177 |
cnostore CADJ STORE dSigmMod(:,:) = comlev1_bibj, key=ikey, byte=isbyte |
178 |
cnostore CADJ STORE dSigmaDr(:,:) = comlev1_bibj, key=ikey, byte=isbyte |
179 |
#endif |
180 |
|
181 |
DO j=1-OLy+1,sNy+OLy-1 |
182 |
DO i=1-OLx+1,sNx+OLx-1 |
183 |
IF (dSigmMod(i,j) .NE. 0.) THEN |
184 |
tmpFld(i,j) = -dSigmMod(i,j)*GM_rMaxSlope |
185 |
IF ( dSigmaDr(i,j) .GE. tmpFld(i,j) ) |
186 |
& dSigmaDr(i,j) = tmpFld(i,j) |
187 |
ENDIF |
188 |
ENDDO |
189 |
ENDDO |
190 |
|
191 |
#ifdef ALLOW_AUTODIFF_TAMC |
192 |
cnostore CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte |
193 |
cnostore CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte |
194 |
cnostore CADJ STORE dSigmaDr(:,:) = comlev1_bibj, key=ikey, byte=isbyte |
195 |
#endif |
196 |
|
197 |
DO j=1-OLy+1,sNy+OLy-1 |
198 |
DO i=1-OLx+1,sNx+OLx-1 |
199 |
IF (dSigmMod(i,j) .EQ. 0.) THEN |
200 |
SlopeX(i,j) = 0. _d 0 |
201 |
SlopeY(i,j) = 0. _d 0 |
202 |
ELSE |
203 |
dRdSigmaLtd(i,j) = 1. _d 0/( dSigmaDr(i,j) ) |
204 |
SlopeX(i,j)=-dSigmaDx(i,j)*dRdSigmaLtd(i,j) |
205 |
SlopeY(i,j)=-dSigmaDy(i,j)*dRdSigmaLtd(i,j) |
206 |
ENDIF |
207 |
ENDDO |
208 |
ENDDO |
209 |
|
210 |
#ifdef ALLOW_AUTODIFF_TAMC |
211 |
cnostore CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte |
212 |
cnostore CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte |
213 |
#endif |
214 |
|
215 |
DO j=1-OLy+1,sNy+OLy-1 |
216 |
DO i=1-OLx+1,sNx+OLx-1 |
217 |
SlopeSqr(i,j)=SlopeX(i,j)*SlopeX(i,j) |
218 |
& +SlopeY(i,j)*SlopeY(i,j) |
219 |
taperFct(i,j)=1. _d 0 |
220 |
ENDDO |
221 |
ENDDO |
222 |
|
223 |
#endif /* GM_EXCLUDE_CLIPPING */ |
224 |
|
225 |
ELSEIF (GM_taper_scheme.EQ.'fm07' ) THEN |
226 |
C-- Ferrari & Mc.Williams, 2007: |
227 |
|
228 |
#ifdef GM_EXCLUDE_FM07_TAP |
229 |
|
230 |
STOP 'Need to compile without "#define GM_EXCLUDE_FM07_TAP"' |
231 |
|
232 |
#else /* GM_EXCLUDE_FM07_TAP */ |
233 |
|
234 |
C- a) Calculate separately slope magnitude (SlopeMod) |
235 |
C and slope horiz. direction (Slope_X,Y : normalized) |
236 |
DO j=1-OLy+1,sNy+OLy-1 |
237 |
DO i=1-OLx+1,sNx+OLx-1 |
238 |
IF ( k.GT.kLow(i,j) ) THEN |
239 |
C- Bottom or below: |
240 |
SlopeX (i,j) = 0. _d 0 |
241 |
SlopeY (i,j) = 0. _d 0 |
242 |
SlopeMod(i,j) = 0. _d 0 |
243 |
taperFct(i,j) = 0. _d 0 |
244 |
ELSE |
245 |
C- Above bottom: |
246 |
IF ( dSigmaDr(i,j).GE. -GM_Small_Number ) |
247 |
& dSigmaDr(i,j) = -GM_Small_Number |
248 |
tmpFld(i,j) = dSigmaDx(i,j)*dSigmaDx(i,j) |
249 |
& + dSigmaDy(i,j)*dSigmaDy(i,j) |
250 |
IF ( tmpFld(i,j).GT.0. ) THEN |
251 |
locVar(i,j) = SQRT( tmpFld(i,j) ) |
252 |
SlopeX (i,j) = dSigmaDx(i,j)/locVar(i,j) |
253 |
SlopeY (i,j) = dSigmaDy(i,j)/locVar(i,j) |
254 |
SlopeMod(i,j) = -locVar(i,j)/dSigmaDr(i,j) |
255 |
taperFct(i,j) = 1. _d 0 |
256 |
ELSE |
257 |
SlopeX (i,j) = 0. _d 0 |
258 |
SlopeY (i,j) = 0. _d 0 |
259 |
SlopeMod(i,j) = 0. _d 0 |
260 |
taperFct(i,j) = 0. _d 0 |
261 |
ENDIF |
262 |
ENDIF |
263 |
ENDDO |
264 |
ENDDO |
265 |
|
266 |
C- b) Set Transition Layer Depth: |
267 |
IF ( k.EQ.1 ) THEN |
268 |
minTransLay = GM_facTrL2dz*( depthZ(k) - depthZ(k+1) ) |
269 |
ELSE |
270 |
minTransLay = GM_facTrL2dz*( depthZ(k-1) - depthZ(k) ) |
271 |
ENDIF |
272 |
DO j=1-OLy+1,sNy+OLy-1 |
273 |
DO i=1-OLx+1,sNx+OLx-1 |
274 |
IF ( hTransLay(i,j).LE.0. _d 0 ) THEN |
275 |
C- previously below the transition layer |
276 |
tmpFld(i,j) = Lrho(i,j)*SlopeMod(i,j) |
277 |
C- ensure that transition layer is at least as thick than current level and |
278 |
C not thicker than the larger of the 2 : maxTransLay and facTrL2ML*hMixLay |
279 |
dTransLay = |
280 |
& MIN( MAX( tmpFld(i,j), minTransLay ), |
281 |
& MAX( GM_facTrL2ML*hMixLay(i,j), GM_maxTransLay ) ) |
282 |
IF ( k.GE.kLow(i,j) ) THEN |
283 |
C- bottom & below & 1rst level above the bottom : |
284 |
recipLambda(i,j) = 0. _d 0 |
285 |
baseSlope(i,j) = SlopeMod(i,j) |
286 |
C-- note: do not catch the 1rst level/interface (k=kLow) above the bottom |
287 |
C since no baseSlope has yet been stored (= 0); but might fit |
288 |
C well into transition layer criteria (if locally not stratified) |
289 |
ELSEIF ( dTransLay+hMixLay(i,j)+depthZ(k) .GE. 0. ) THEN |
290 |
C- Found the transition Layer : set depth of base of Transition layer |
291 |
hTransLay(i,j) = -depthZ(k+1) |
292 |
C and compute inverse length scale "1/Lambda" from slope vert. grad |
293 |
IF ( baseSlope(i,j).GT.0. ) THEN |
294 |
recipLambda(i,j) = recipLambda(i,j) |
295 |
& / MIN( baseSlope(i,j), GM_maxSlope ) |
296 |
ELSE |
297 |
recipLambda(i,j) = 0. _d 0 |
298 |
ENDIF |
299 |
C slope^2 & Kwz should remain > 0 : prevent too large negative D/lambda |
300 |
IF ( hMixLay(i,j)+depthZ(k+1).LT.0. ) THEN |
301 |
rLambMin = 1. _d 0 /( hMixLay(i,j)+depthZ(k+1) ) |
302 |
recipLambda(i,j) = MAX( recipLambda(i,j), rLambMin ) |
303 |
ENDIF |
304 |
ELSE |
305 |
C- Still below Trans. layer: store slope & slope vert. grad. |
306 |
recipLambda(i,j) = ( MIN( SlopeMod(i,j), GM_maxSlope ) |
307 |
& - MIN( baseSlope(i,j), GM_maxSlope ) |
308 |
& ) / ( depthZ(k) - depthZ(k+1) ) |
309 |
baseSlope(i,j) = SlopeMod(i,j) |
310 |
ENDIF |
311 |
ENDIF |
312 |
ENDDO |
313 |
ENDDO |
314 |
|
315 |
C- c) Set Slope component according to vertical position |
316 |
C (in Mixed-Layer / in Transition Layer / below Transition Layer) |
317 |
DO j=1-OLy+1,sNy+OLy-1 |
318 |
DO i=1-OLx+1,sNx+OLx-1 |
319 |
IF ( hTransLay(i,j).GT.0. _d 0 ) THEN |
320 |
C- already above base of transition layer: |
321 |
|
322 |
DoverLamb = (hTransLay(i,j)-hMixLay(i,j))*recipLambda(i,j) |
323 |
IF ( -depthZ(k).LE.hMixLay(i,j) ) THEN |
324 |
C- compute tapering function -G(z) in the mixed Layer: |
325 |
taperFctLoc = |
326 |
& ( -depthZ(k)/(hTransLay(i,j)+hMixLay(i,j)) |
327 |
& *( 2. _d 0 + DoverLamb ) |
328 |
& ) |
329 |
C- compute tapering function -G^(z) in the mixed Layer |
330 |
taperFctHat = |
331 |
& ( -depthZ(k)/(hTransLay(i,j)+hMixLay(i,j)) |
332 |
& * 2. _d 0 |
333 |
& *( 1. _d 0 + DoverLamb ) |
334 |
& ) |
335 |
ELSE |
336 |
C- compute tapering function -G(z) in the transition Layer: |
337 |
taperFctLoc = |
338 |
& ( -depthZ(k)/(hTransLay(i,j)+hMixLay(i,j)) |
339 |
& *( 2. _d 0 + DoverLamb ) |
340 |
& ) |
341 |
& - ( (depthZ(k)+hMixLay(i,j))*(depthZ(k)+hMixLay(i,j)) |
342 |
& /( hTransLay(i,j)*hTransLay(i,j) |
343 |
& - hMixLay(i,j)*hMixLay(i,j) ) |
344 |
& *( 1. _d 0 + hTransLay(i,j)*recipLambda(i,j) ) |
345 |
& ) |
346 |
C- compute tapering function -G^(z) in the transition Layer: |
347 |
taperFctHat = |
348 |
& ( -depthZ(k)/(hTransLay(i,j)+hMixLay(i,j)) |
349 |
& * 2. _d 0 |
350 |
& *( 1. _d 0 + DoverLamb ) |
351 |
& ) |
352 |
& - ( (depthZ(k)+hMixLay(i,j))*(depthZ(k)+hMixLay(i,j)) |
353 |
& /( hTransLay(i,j)*hTransLay(i,j) |
354 |
& - hMixLay(i,j)*hMixLay(i,j) ) |
355 |
& *( 1. _d 0 + hTransLay(i,j)*recipLambda(i,j)*2. _d 0 ) |
356 |
& ) |
357 |
ENDIF |
358 |
C- modify the slope (above bottom of transition layer): |
359 |
c Smod = baseSlope(i,j) |
360 |
C- safer to limit the slope (even if it might never exceed GM_maxSlope) |
361 |
Smod = MIN( baseSlope(i,j), GM_maxSlope ) |
362 |
SlopeX(i,j) = SlopeX(i,j)*Smod*taperFctLoc |
363 |
SlopeY(i,j) = SlopeY(i,j)*Smod*taperFctLoc |
364 |
c SlopeSqr(i,j) = Smod*Smod*taperFctHat |
365 |
c SlopeSqr(i,j) = baseSlope(i,j)*Smod*taperFctHat |
366 |
SlopeSqr(i,j) = MIN( baseSlope(i,j), GM_bigSlope ) |
367 |
& *Smod*taperFctHat |
368 |
|
369 |
ELSE |
370 |
C-- Still below base of transition layer: |
371 |
c Smod = SlopeMod(i,j) |
372 |
C- safer to limit the slope: |
373 |
Smod = MIN( SlopeMod(i,j), GM_maxSlope ) |
374 |
SlopeX(i,j) = SlopeX(i,j)*Smod |
375 |
SlopeY(i,j) = SlopeY(i,j)*Smod |
376 |
c SlopeSqr(i,j) = Smod*Smod |
377 |
c SlopeSqr(i,j) = SlopeMod(i,j)*Smod |
378 |
SlopeSqr(i,j) = MIN( SlopeMod(i,j), GM_bigSlope ) |
379 |
& *Smod |
380 |
|
381 |
C-- end if baseSlope > 0 / else => above/below base of Trans. Layer |
382 |
ENDIF |
383 |
|
384 |
ENDDO |
385 |
ENDDO |
386 |
|
387 |
#endif /* GM_EXCLUDE_FM07_TAP */ |
388 |
|
389 |
ELSE IF (GM_taper_scheme.EQ.'ac02') THEN |
390 |
|
391 |
#ifdef GM_EXCLUDE_AC02_TAP |
392 |
|
393 |
STOP 'Need to compile without "#define GM_EXCLUDE_AC02_TAP"' |
394 |
|
395 |
#else /* GM_EXCLUDE_AC02_TAP */ |
396 |
|
397 |
C- New Scheme (A. & C. 2002): relax part of the small slope approximation |
398 |
C compute the true slope (no approximation) |
399 |
C but still neglect Kxy & Kyx (assumed to be zero) |
400 |
|
401 |
maxSlopeSqr = GM_maxSlope*GM_maxSlope |
402 |
DO j=1-OLy+1,sNy+OLy-1 |
403 |
DO i=1-OLx+1,sNx+OLx-1 |
404 |
dRdSigmaLtd(i,j)= |
405 |
& dSigmaDx(i,j)*dSigmaDx(i,j) |
406 |
& + dSigmaDy(i,j)*dSigmaDy(i,j) |
407 |
& + dSigmaDr(i,j)*dSigmaDr(i,j) |
408 |
taperFct(i,j) = 1. _d 0 |
409 |
|
410 |
IF (dRdSigmaLtd(i,j).NE.0.) THEN |
411 |
dRdSigmaLtd(i,j)=1. _d 0 |
412 |
& / ( dRdSigmaLtd(i,j) ) |
413 |
SlopeSqr(i,j)=(dSigmaDx(i,j)*dSigmaDx(i,j) |
414 |
& +dSigmaDy(i,j)*dSigmaDy(i,j))*dRdSigmaLtd(i,j) |
415 |
SlopeX(i,j)=-dSigmaDx(i,j) |
416 |
& *dRdSigmaLtd(i,j)*dSigmaDr(i,j) |
417 |
SlopeY(i,j)=-dSigmaDy(i,j) |
418 |
& *dRdSigmaLtd(i,j)*dSigmaDr(i,j) |
419 |
cph T11(i,j)=(dSigmaDr(i,j)**2)*dRdSigmaLtd(i,j) |
420 |
ENDIF |
421 |
#ifndef ALLOWW_AUTODIFF_TAMC |
422 |
cph-- this part does not adjoint well |
423 |
IF ( SlopeSqr(i,j) .GT. maxSlopeSqr .AND. |
424 |
& SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN |
425 |
taperFct(i,j) = maxSlopeSqr/SlopeSqr(i,j) |
426 |
ELSE IF ( SlopeSqr(i,j) .GT. GM_slopeSqCutoff ) THEN |
427 |
taperFct(i,j) = 0. _d 0 |
428 |
ENDIF |
429 |
#endif |
430 |
ENDDO |
431 |
ENDDO |
432 |
|
433 |
#endif /* GM_EXCLUDE_AC02_TAP */ |
434 |
|
435 |
ELSE |
436 |
|
437 |
#ifdef GM_EXCLUDE_TAPERING |
438 |
|
439 |
STOP 'Need to compile without "#define GM_EXCLUDE_TAPERING"' |
440 |
|
441 |
#else /* GM_EXCLUDE_TAPERING */ |
442 |
|
443 |
C---------------------------------------------------------------------- |
444 |
|
445 |
C- Compute the slope, no clipping, but avoid reverse slope in negatively |
446 |
C stratified (Sigma_Z > 0) region : |
447 |
|
448 |
#ifdef ALLOW_AUTODIFF_TAMC |
449 |
cnostore CADJ STORE dSigmaDr(:,:) = comlev1_bibj, key=ikey, byte=isbyte |
450 |
#endif |
451 |
|
452 |
DO j=1-OLy+1,sNy+OLy-1 |
453 |
DO i=1-OLx+1,sNx+OLx-1 |
454 |
IF ( dSigmaDr(i,j) .NE. 0. ) THEN |
455 |
IF (dSigmaDr(i,j).GE.(-GM_Small_Number)) |
456 |
& dSigmaDr(i,j) = -GM_Small_Number |
457 |
ENDIF |
458 |
ENDDO |
459 |
ENDDO |
460 |
|
461 |
#ifdef ALLOW_AUTODIFF_TAMC |
462 |
cnostore CADJ STORE dSigmaDx(:,:) = comlev1_bibj, key=ikey, byte=isbyte |
463 |
cnostore CADJ STORE dSigmaDy(:,:) = comlev1_bibj, key=ikey, byte=isbyte |
464 |
cnostore CADJ STORE dSigmaDr(:,:) = comlev1_bibj, key=ikey, byte=isbyte |
465 |
#endif |
466 |
|
467 |
DO j=1-OLy+1,sNy+OLy-1 |
468 |
DO i=1-OLx+1,sNx+OLx-1 |
469 |
IF ( dSigmaDr(i,j) .EQ. 0. ) THEN |
470 |
IF ( dSigmaDx(i,j) .NE. 0. ) THEN |
471 |
SlopeX(i,j) = SIGN( GM_bigSlope, dSigmaDx(i,j) ) |
472 |
ELSE |
473 |
SlopeX(i,j) = 0. _d 0 |
474 |
ENDIF |
475 |
IF ( dSigmaDy(i,j) .NE. 0. ) THEN |
476 |
SlopeY(i,j) = SIGN( GM_bigSlope, dSigmaDy(i,j) ) |
477 |
ELSE |
478 |
SlopeY(i,j) = 0. _d 0 |
479 |
ENDIF |
480 |
ELSE |
481 |
dRdSigmaLtd(i,j) = 1. _d 0 / dSigmaDr(i,j) |
482 |
SlopeX(i,j)=-dSigmaDx(i,j)*dRdSigmaLtd(i,j) |
483 |
SlopeY(i,j)=-dSigmaDy(i,j)*dRdSigmaLtd(i,j) |
484 |
c SlopeX(i,j) = -dSigmaDx(i,j)/dSigmaDr(i,j) |
485 |
c SlopeY(i,j) = -dSigmaDy(i,j)/dSigmaDr(i,j) |
486 |
ENDIF |
487 |
ENDDO |
488 |
ENDDO |
489 |
|
490 |
#ifdef ALLOW_AUTODIFF_TAMC |
491 |
cnostore CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte |
492 |
cnostore CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte |
493 |
#endif |
494 |
|
495 |
DO j=1-OLy+1,sNy+OLy-1 |
496 |
DO i=1-OLx+1,sNx+OLx-1 |
497 |
SlopeSqr(i,j) = SlopeX(i,j)*SlopeX(i,j) |
498 |
& +SlopeY(i,j)*SlopeY(i,j) |
499 |
taperFct(i,j) = 1. _d 0 |
500 |
IF ( SlopeSqr(i,j) .GT. GM_slopeSqCutoff ) THEN |
501 |
slopeSqr(i,j) = GM_slopeSqCutoff |
502 |
taperFct(i,j) = 0. _d 0 |
503 |
ENDIF |
504 |
ENDDO |
505 |
ENDDO |
506 |
|
507 |
C- Compute the tapering function for the GM+Redi tensor : |
508 |
|
509 |
IF (GM_taper_scheme.EQ.'linear') THEN |
510 |
|
511 |
C- Simplest adiabatic tapering = Smax/Slope (linear) |
512 |
maxSlopeSqr = GM_maxSlope*GM_maxSlope |
513 |
DO j=1-OLy+1,sNy+OLy-1 |
514 |
DO i=1-OLx+1,sNx+OLx-1 |
515 |
|
516 |
IF ( SlopeSqr(i,j) .EQ. 0. ) THEN |
517 |
taperFct(i,j) = 1. _d 0 |
518 |
ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr .AND. |
519 |
& SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN |
520 |
taperFct(i,j) = SQRT(maxSlopeSqr / SlopeSqr(i,j)) |
521 |
slopeSqr(i,j) = MIN( slopeSqr(i,j),GM_bigSlope*GM_bigSlope ) |
522 |
ENDIF |
523 |
|
524 |
ENDDO |
525 |
ENDDO |
526 |
|
527 |
ELSEIF (GM_taper_scheme.EQ.'gkw91') THEN |
528 |
|
529 |
C- Gerdes, Koberle and Willebrand, Clim. Dyn. 1991 |
530 |
maxSlopeSqr = GM_maxSlope*GM_maxSlope |
531 |
DO j=1-OLy+1,sNy+OLy-1 |
532 |
DO i=1-OLx+1,sNx+OLx-1 |
533 |
|
534 |
IF ( SlopeSqr(i,j) .EQ. 0. ) THEN |
535 |
taperFct(i,j) = 1. _d 0 |
536 |
ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr .AND. |
537 |
& SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN |
538 |
taperFct(i,j) = maxSlopeSqr/SlopeSqr(i,j) |
539 |
ENDIF |
540 |
|
541 |
ENDDO |
542 |
ENDDO |
543 |
|
544 |
ELSEIF (GM_taper_scheme.EQ.'dm95') THEN |
545 |
|
546 |
C- Danabasoglu and McWilliams, J. Clim. 1995 |
547 |
DO j=1-OLy+1,sNy+OLy-1 |
548 |
DO i=1-OLx+1,sNx+OLx-1 |
549 |
|
550 |
IF ( SlopeSqr(i,j) .EQ. 0. ) THEN |
551 |
taperFct(i,j) = 1. _d 0 |
552 |
ELSE IF ( SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN |
553 |
Smod=SQRT(SlopeSqr(i,j)) |
554 |
taperFct(i,j)=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd )) |
555 |
ENDIF |
556 |
ENDDO |
557 |
ENDDO |
558 |
|
559 |
ELSEIF (GM_taper_scheme.EQ.'ldd97') THEN |
560 |
|
561 |
C- Large, Danabasoglu and Doney, JPO 1997 |
562 |
DO j=1-OLy+1,sNy+OLy-1 |
563 |
DO i=1-OLx+1,sNx+OLx-1 |
564 |
|
565 |
IF (SlopeSqr(i,j) .EQ. 0.) THEN |
566 |
taperFct(i,j) = 1. _d 0 |
567 |
ELSEIF ( SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN |
568 |
Smod=SQRT(SlopeSqr(i,j)) |
569 |
f1=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd )) |
570 |
Rnondim= -depthZ(k)/(Lrho(i,j)*Smod) |
571 |
IF ( Rnondim.GE.1. _d 0 ) THEN |
572 |
f2 = 1. _d 0 |
573 |
ELSE |
574 |
f2 = op5*( 1. _d 0 + SIN( fpi*(Rnondim-op5) )) |
575 |
ENDIF |
576 |
taperFct(i,j)=f1*f2 |
577 |
ENDIF |
578 |
|
579 |
ENDDO |
580 |
ENDDO |
581 |
|
582 |
ELSEIF (GM_taper_scheme.EQ.'stableGmAdjTap') THEN |
583 |
|
584 |
#ifndef GMREDI_WITH_STABLE_ADJOINT |
585 |
|
586 |
STOP 'Need to compile wth "#define GMREDI_WITH_STABLE_ADJOINT"' |
587 |
|
588 |
#else /* GMREDI_WITH_STABLE_ADJOINT */ |
589 |
|
590 |
c special choice for adjoint/optimization of parameters |
591 |
c (~ strong clipping, reducing non linearity of kw=f(K)) |
592 |
|
593 |
slopeMax= 2. _d -3 |
594 |
|
595 |
CADJ STORE SlopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte |
596 |
CADJ STORE SlopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte |
597 |
|
598 |
DO j=1-OLy,sNy+OLy |
599 |
DO i=1-OLx+1,sNx+OLx |
600 |
slopeSqTmp=SlopeX(i,j)*SlopeX(i,j) |
601 |
& +SlopeY(i,j)*SlopeY(i,j) |
602 |
|
603 |
IF ( slopeSqTmp .GT. slopeMax**2 ) then |
604 |
slopeTmp=sqrt(slopeSqTmp) |
605 |
SlopeX(i,j)=SlopeX(i,j)*slopeMax/slopeTmp |
606 |
SlopeY(i,j)=SlopeY(i,j)*slopeMax/slopeTmp |
607 |
ENDIF |
608 |
SlopeSqr(i,j) = SlopeX(i,j)*SlopeX(i,j) |
609 |
& +SlopeY(i,j)*SlopeY(i,j) |
610 |
taperFct(i,j) = 1. _d 0 |
611 |
ENDDO |
612 |
ENDDO |
613 |
|
614 |
#endif /* GMREDI_WITH_STABLE_ADJOINT */ |
615 |
|
616 |
ELSEIF (GM_taper_scheme.NE.' ') THEN |
617 |
STOP 'GMREDI_SLOPE_LIMIT: Bad GM_taper_scheme' |
618 |
ENDIF |
619 |
|
620 |
#endif /* GM_EXCLUDE_TAPERING */ |
621 |
|
622 |
ENDIF |
623 |
|
624 |
#endif /* ALLOW_GMREDI */ |
625 |
|
626 |
RETURN |
627 |
END |