1 |
|
44659008 |
SUBROUTINE drag_noro_strato(partdrag, nlon, nlev, dtime, paprs, pplay, pmea, pstd, & |
2 |
|
576 |
psig, pgam, pthe, ppic, pval, kgwd, kdx, ktest, t, u, v, pulow, pvlow, & |
3 |
|
576 |
pustr, pvstr, d_t, d_u, d_v) |
4 |
|
|
|
5 |
|
|
USE dimphy |
6 |
|
|
IMPLICIT NONE |
7 |
|
|
! ====================================================================== |
8 |
|
|
! Auteur(s): F.Lott (LMD/CNRS) date: 19950201 |
9 |
|
|
! Object: Mountain drag interface. Made necessary because: |
10 |
|
|
! 1. in the LMD-GCM Layers are from bottom to top, |
11 |
|
|
! contrary to most European GCM. |
12 |
|
|
! 2. the altitude above ground of each model layers |
13 |
|
|
! needs to be known (variable zgeom) |
14 |
|
|
! ====================================================================== |
15 |
|
|
! Explicit Arguments: |
16 |
|
|
! ================== |
17 |
|
|
! partdrag-input-I-control which part of the drag we consider (total part or GW part) |
18 |
|
|
! nlon----input-I-Total number of horizontal points that get into physics |
19 |
|
|
! nlev----input-I-Number of vertical levels |
20 |
|
|
! dtime---input-R-Time-step (s) |
21 |
|
|
! paprs---input-R-Pressure in semi layers (Pa) |
22 |
|
|
! pplay---input-R-Pressure model-layers (Pa) |
23 |
|
|
! t-------input-R-temperature (K) |
24 |
|
|
! u-------input-R-Horizontal wind (m/s) |
25 |
|
|
! v-------input-R-Meridional wind (m/s) |
26 |
|
|
! pmea----input-R-Mean Orography (m) |
27 |
|
|
! pstd----input-R-SSO standard deviation (m) |
28 |
|
|
! psig----input-R-SSO slope |
29 |
|
|
! pgam----input-R-SSO Anisotropy |
30 |
|
|
! pthe----input-R-SSO Angle |
31 |
|
|
! ppic----input-R-SSO Peacks elevation (m) |
32 |
|
|
! pval----input-R-SSO Valleys elevation (m) |
33 |
|
|
|
34 |
|
|
! kgwd- -input-I: Total nb of points where the orography schemes are active |
35 |
|
|
! ktest--input-I: Flags to indicate active points |
36 |
|
|
! kdx----input-I: Locate the physical location of an active point. |
37 |
|
|
|
38 |
|
|
! pulow, pvlow -output-R: Low-level wind |
39 |
|
|
! pustr, pvstr -output-R: Surface stress due to SSO drag (Pa) |
40 |
|
|
|
41 |
|
|
! d_t-----output-R: T increment |
42 |
|
|
! d_u-----output-R: U increment |
43 |
|
|
! d_v-----output-R: V increment |
44 |
|
|
|
45 |
|
|
! Implicit Arguments: |
46 |
|
|
! =================== |
47 |
|
|
|
48 |
|
|
! iim--common-I: Number of longitude intervals |
49 |
|
|
! jjm--common-I: Number of latitude intervals |
50 |
|
|
! klon-common-I: Number of points seen by the physics |
51 |
|
|
! (iim+1)*(jjm+1) for instance |
52 |
|
|
! klev-common-I: Number of vertical layers |
53 |
|
|
! ====================================================================== |
54 |
|
|
! Local Variables: |
55 |
|
|
! ================ |
56 |
|
|
|
57 |
|
|
! zgeom-----R: Altitude of layer above ground |
58 |
|
|
! pt, pu, pv --R: t u v from top to bottom |
59 |
|
|
! pdtdt, pdudt, pdvdt --R: t u v tendencies (from top to bottom) |
60 |
|
|
! papmf: pressure at model layer (from top to bottom) |
61 |
|
|
! papmh: pressure at model 1/2 layer (from top to bottom) |
62 |
|
|
|
63 |
|
|
! ====================================================================== |
64 |
|
|
include "YOMCST.h" |
65 |
|
|
include "YOEGWD.h" |
66 |
|
|
|
67 |
|
|
! ARGUMENTS |
68 |
|
|
|
69 |
|
|
INTEGER partdrag,nlon, nlev |
70 |
|
|
REAL dtime |
71 |
|
|
REAL paprs(nlon, nlev+1) |
72 |
|
|
REAL pplay(nlon, nlev) |
73 |
|
|
REAL pmea(nlon), pstd(nlon), psig(nlon), pgam(nlon), pthe(nlon) |
74 |
|
|
REAL ppic(nlon), pval(nlon) |
75 |
|
|
REAL pulow(nlon), pvlow(nlon), pustr(nlon), pvstr(nlon) |
76 |
|
|
REAL t(nlon, nlev), u(nlon, nlev), v(nlon, nlev) |
77 |
|
|
REAL d_t(nlon, nlev), d_u(nlon, nlev), d_v(nlon, nlev) |
78 |
|
|
|
79 |
|
|
INTEGER i, k, kgwd, kdx(nlon), ktest(nlon) |
80 |
|
|
|
81 |
|
|
! LOCAL VARIABLES: |
82 |
|
|
|
83 |
|
1152 |
REAL zgeom(klon, klev) |
84 |
|
1152 |
REAL pdtdt(klon, klev), pdudt(klon, klev), pdvdt(klon, klev) |
85 |
|
1152 |
REAL pt(klon, klev), pu(klon, klev), pv(klon, klev) |
86 |
|
576 |
REAL papmf(klon, klev), papmh(klon, klev+1) |
87 |
|
|
CHARACTER (LEN=20) :: modname = 'orografi_strato' |
88 |
|
|
CHARACTER (LEN=80) :: abort_message |
89 |
|
|
|
90 |
|
|
! INITIALIZE OUTPUT VARIABLES |
91 |
|
|
|
92 |
✓✓ |
573120 |
DO i = 1, klon |
93 |
|
572544 |
pulow(i) = 0.0 |
94 |
|
572544 |
pvlow(i) = 0.0 |
95 |
|
572544 |
pustr(i) = 0.0 |
96 |
|
573120 |
pvstr(i) = 0.0 |
97 |
|
|
END DO |
98 |
✓✓ |
23040 |
DO k = 1, klev |
99 |
✓✓ |
22352256 |
DO i = 1, klon |
100 |
|
22329216 |
d_t(i, k) = 0.0 |
101 |
|
22329216 |
d_u(i, k) = 0.0 |
102 |
|
22329216 |
d_v(i, k) = 0.0 |
103 |
|
22329216 |
pdudt(i, k) = 0.0 |
104 |
|
22329216 |
pdvdt(i, k) = 0.0 |
105 |
|
22351680 |
pdtdt(i, k) = 0.0 |
106 |
|
|
END DO |
107 |
|
|
END DO |
108 |
|
|
|
109 |
|
|
! PREPARE INPUT VARIABLES FOR ORODRAG (i.e., ORDERED FROM TOP TO BOTTOM) |
110 |
|
|
! CALCULATE LAYERS HEIGHT ABOVE GROUND) |
111 |
|
|
|
112 |
✓✓ |
23040 |
DO k = 1, klev |
113 |
✓✓ |
22352256 |
DO i = 1, klon |
114 |
|
22329216 |
pt(i, k) = t(i, klev-k+1) |
115 |
|
22329216 |
pu(i, k) = u(i, klev-k+1) |
116 |
|
22329216 |
pv(i, k) = v(i, klev-k+1) |
117 |
|
22351680 |
papmf(i, k) = pplay(i, klev-k+1) |
118 |
|
|
END DO |
119 |
|
|
END DO |
120 |
✓✓ |
23616 |
DO k = 1, klev + 1 |
121 |
✓✓ |
22925376 |
DO i = 1, klon |
122 |
|
22924800 |
papmh(i, k) = paprs(i, klev-k+2) |
123 |
|
|
END DO |
124 |
|
|
END DO |
125 |
✓✓ |
573120 |
DO i = 1, klon |
126 |
|
573120 |
zgeom(i, klev) = rd*pt(i, klev)*log(papmh(i,klev+1)/papmf(i,klev)) |
127 |
|
|
END DO |
128 |
✓✓ |
22464 |
DO k = klev - 1, 1, -1 |
129 |
✓✓ |
21779136 |
DO i = 1, klon |
130 |
|
|
zgeom(i, k) = zgeom(i, k+1) + rd*(pt(i,k)+pt(i,k+1))/2.0*log(papmf(i,k+ & |
131 |
|
21778560 |
1)/papmf(i,k)) |
132 |
|
|
END DO |
133 |
|
|
END DO |
134 |
|
|
|
135 |
|
|
! CALL SSO DRAG ROUTINES |
136 |
|
|
|
137 |
|
|
CALL orodrag_strato(partdrag,klon, klev, kgwd, kdx, ktest, dtime, papmh, papmf, & |
138 |
|
|
zgeom, pt, pu, pv, pmea, pstd, psig, pgam, pthe, ppic, pval, pulow, & |
139 |
|
576 |
pvlow, pdudt, pdvdt, pdtdt) |
140 |
|
|
|
141 |
|
|
! COMPUTE INCREMENTS AND STRESS FROM TENDENCIES |
142 |
|
|
|
143 |
✓✓ |
23040 |
DO k = 1, klev |
144 |
✓✓ |
22352256 |
DO i = 1, klon |
145 |
|
22329216 |
d_u(i, klev+1-k) = dtime*pdudt(i, k) |
146 |
|
22329216 |
d_v(i, klev+1-k) = dtime*pdvdt(i, k) |
147 |
|
22329216 |
d_t(i, klev+1-k) = dtime*pdtdt(i, k) |
148 |
|
22329216 |
pustr(i) = pustr(i) + pdudt(i, k)*(papmh(i,k+1)-papmh(i,k))/rg |
149 |
|
22351680 |
pvstr(i) = pvstr(i) + pdvdt(i, k)*(papmh(i,k+1)-papmh(i,k))/rg |
150 |
|
|
END DO |
151 |
|
|
END DO |
152 |
|
|
|
153 |
|
576 |
RETURN |
154 |
|
|
END SUBROUTINE drag_noro_strato |
155 |
|
|
|
156 |
|
20847168 |
SUBROUTINE orodrag_strato(partdrag,nlon, nlev, kgwd, kdx, ktest, ptsphy, paphm1, & |
157 |
|
|
papm1, pgeom1, ptm1, pum1, pvm1, pmea, pstd, psig, pgam, pthe, ppic, pval & |
158 |
|
|
! outputs |
159 |
|
|
, pulow, pvlow, pvom, pvol, pte) |
160 |
|
|
|
161 |
|
35716417 |
USE dimphy |
162 |
|
|
IMPLICIT NONE |
163 |
|
|
|
164 |
|
|
|
165 |
|
|
! **** *orodrag* - does the SSO drag parametrization. |
166 |
|
|
|
167 |
|
|
! purpose. |
168 |
|
|
! -------- |
169 |
|
|
|
170 |
|
|
! this routine computes the physical tendencies of the |
171 |
|
|
! prognostic variables u,v and t due to vertical transports by |
172 |
|
|
! subgridscale orographically excited gravity waves, and to |
173 |
|
|
! low level blocked flow drag. |
174 |
|
|
|
175 |
|
|
! ** interface. |
176 |
|
|
! ---------- |
177 |
|
|
! called from *drag_noro*. |
178 |
|
|
|
179 |
|
|
! the routine takes its input from the long-term storage: |
180 |
|
|
! u,v,t and p at t-1. |
181 |
|
|
|
182 |
|
|
! explicit arguments : |
183 |
|
|
! -------------------- |
184 |
|
|
! ==== inputs === |
185 |
|
|
! partdrag-input-I-control which part of the drag we consider (total part or GW part) |
186 |
|
|
! nlon----input-I-Total number of horizontal points that get into physics |
187 |
|
|
! nlev----input-I-Number of vertical levels |
188 |
|
|
|
189 |
|
|
! kgwd- -input-I: Total nb of points where the orography schemes are active |
190 |
|
|
! ktest--input-I: Flags to indicate active points |
191 |
|
|
! kdx----input-I: Locate the physical location of an active point. |
192 |
|
|
! ptsphy--input-R-Time-step (s) |
193 |
|
|
! paphm1--input-R: pressure at model 1/2 layer |
194 |
|
|
! papm1---input-R: pressure at model layer |
195 |
|
|
! pgeom1--input-R: Altitude of layer above ground |
196 |
|
|
! ptm1, pum1, pvm1--R-: t, u and v |
197 |
|
|
! pmea----input-R-Mean Orography (m) |
198 |
|
|
! pstd----input-R-SSO standard deviation (m) |
199 |
|
|
! psig----input-R-SSO slope |
200 |
|
|
! pgam----input-R-SSO Anisotropy |
201 |
|
|
! pthe----input-R-SSO Angle |
202 |
|
|
! ppic----input-R-SSO Peacks elevation (m) |
203 |
|
|
! pval----input-R-SSO Valleys elevation (m) |
204 |
|
|
|
205 |
|
|
INTEGER nlon, nlev, kgwd |
206 |
|
|
REAL ptsphy |
207 |
|
|
|
208 |
|
|
! ==== outputs === |
209 |
|
|
! pulow, pvlow -output-R: Low-level wind |
210 |
|
|
|
211 |
|
|
! pte -----output-R: T tendency |
212 |
|
|
! pvom-----output-R: U tendency |
213 |
|
|
! pvol-----output-R: V tendency |
214 |
|
|
|
215 |
|
|
|
216 |
|
|
! Implicit Arguments: |
217 |
|
|
! =================== |
218 |
|
|
|
219 |
|
|
! klon-common-I: Number of points seen by the physics |
220 |
|
|
! klev-common-I: Number of vertical layers |
221 |
|
|
|
222 |
|
|
! method. |
223 |
|
|
! ------- |
224 |
|
|
|
225 |
|
|
! externals. |
226 |
|
|
! ---------- |
227 |
|
|
INTEGER ismin, ismax |
228 |
|
|
EXTERNAL ismin, ismax |
229 |
|
|
|
230 |
|
|
! reference. |
231 |
|
|
! ---------- |
232 |
|
|
|
233 |
|
|
! author. |
234 |
|
|
! ------- |
235 |
|
|
! m.miller + b.ritter e.c.m.w.f. 15/06/86. |
236 |
|
|
|
237 |
|
|
! f.lott + m. miller e.c.m.w.f. 22/11/94 |
238 |
|
|
! ----------------------------------------------------------------------- |
239 |
|
|
|
240 |
|
|
|
241 |
|
|
include "YOMCST.h" |
242 |
|
|
include "YOEGWD.h" |
243 |
|
|
|
244 |
|
|
! ----------------------------------------------------------------------- |
245 |
|
|
|
246 |
|
|
! * 0.1 arguments |
247 |
|
|
! --------- |
248 |
|
|
|
249 |
|
|
INTEGER partdrag |
250 |
|
|
REAL pte(nlon, nlev), pvol(nlon, nlev), pvom(nlon, nlev), pulow(nlon), & |
251 |
|
|
pvlow(nlon) |
252 |
|
|
REAL pum1(nlon, nlev), pvm1(nlon, nlev), ptm1(nlon, nlev), pmea(nlon), & |
253 |
|
|
pstd(nlon), psig(nlon), pgam(nlon), pthe(nlon), ppic(nlon), pval(nlon), & |
254 |
|
|
pgeom1(nlon, nlev), papm1(nlon, nlev), paphm1(nlon, nlev+1) |
255 |
|
|
|
256 |
|
|
INTEGER kdx(nlon), ktest(nlon) |
257 |
|
|
! ----------------------------------------------------------------------- |
258 |
|
|
|
259 |
|
|
! * 0.2 local arrays |
260 |
|
|
! ------------ |
261 |
|
1152 |
INTEGER isect(klon), icrit(klon), ikcrith(klon), ikenvh(klon), iknu(klon), & |
262 |
|
1152 |
iknu2(klon), ikcrit(klon), ikhlim(klon) |
263 |
|
|
|
264 |
|
1152 |
REAL ztau(klon, klev+1), zstab(klon, klev+1), zvph(klon, klev+1), & |
265 |
|
1152 |
zrho(klon, klev+1), zri(klon, klev+1), zpsi(klon, klev+1), & |
266 |
|
1152 |
zzdep(klon, klev) |
267 |
|
1152 |
REAL zdudt(klon), zdvdt(klon), zdtdt(klon), zdedt(klon), zvidis(klon), & |
268 |
|
1152 |
ztfr(klon), znu(klon), zd1(klon), zd2(klon), zdmod(klon) |
269 |
|
|
|
270 |
|
|
|
271 |
|
|
! local quantities: |
272 |
|
|
|
273 |
|
|
INTEGER jl, jk, ji |
274 |
|
|
REAL ztmst, zdelp, ztemp, zforc, ztend, rover, facpart |
275 |
|
|
REAL zb, zc, zconb, zabsv, zzd1, ratio, zbet, zust, zvst, zdis |
276 |
|
|
|
277 |
|
|
! ------------------------------------------------------------------ |
278 |
|
|
|
279 |
|
|
! * 1. initialization |
280 |
|
|
! -------------- |
281 |
|
|
|
282 |
|
|
! print *,' in orodrag' |
283 |
|
|
|
284 |
|
|
! ------------------------------------------------------------------ |
285 |
|
|
|
286 |
|
|
! * 1.1 computational constants |
287 |
|
|
! ----------------------- |
288 |
|
|
|
289 |
|
|
|
290 |
|
|
! ztmst=twodt |
291 |
|
|
! if(nstep.eq.nstart) ztmst=0.5*twodt |
292 |
|
576 |
ztmst = ptsphy |
293 |
|
|
|
294 |
|
|
! ------------------------------------------------------------------ |
295 |
|
|
|
296 |
|
|
! * 1.3 check whether row contains point for printing |
297 |
|
|
! --------------------------------------------- |
298 |
|
|
|
299 |
|
|
|
300 |
|
|
! ------------------------------------------------------------------ |
301 |
|
|
|
302 |
|
|
! * 2. precompute basic state variables. |
303 |
|
|
! * ---------- ----- ----- ---------- |
304 |
|
|
! * define low level wind, project winds in plane of |
305 |
|
|
! * low level wind, determine sector in which to take |
306 |
|
|
! * the variance and set indicator for critical levels. |
307 |
|
|
|
308 |
|
|
|
309 |
|
|
|
310 |
|
|
|
311 |
|
|
|
312 |
|
|
CALL orosetup_strato(nlon, nlev, ktest, ikcrit, ikcrith, icrit, isect, & |
313 |
|
|
ikhlim, ikenvh, iknu, iknu2, paphm1, papm1, pum1, pvm1, ptm1, pgeom1, & |
314 |
|
|
pstd, zrho, zri, zstab, ztau, zvph, zpsi, zzdep, pulow, pvlow, pthe, & |
315 |
|
576 |
pgam, pmea, ppic, pval, znu, zd1, zd2, zdmod) |
316 |
|
|
|
317 |
|
|
! *********************************************************** |
318 |
|
|
|
319 |
|
|
|
320 |
|
|
! * 3. compute low level stresses using subcritical and |
321 |
|
|
! * supercritical forms.computes anisotropy coefficient |
322 |
|
|
! * as measure of orographic twodimensionality. |
323 |
|
|
|
324 |
|
|
|
325 |
|
|
CALL gwstress_strato(nlon, nlev, ikcrit, isect, ikhlim, ktest, ikcrith, & |
326 |
|
|
icrit, ikenvh, iknu, zrho, zstab, zvph, pstd, psig, pmea, ppic, pval, & |
327 |
|
576 |
ztfr, ztau, pgeom1, pgam, zd1, zd2, zdmod, znu) |
328 |
|
|
|
329 |
|
|
! * 4. compute stress profile including |
330 |
|
|
! trapped waves, wave breaking, |
331 |
|
|
! linear decay in stratosphere. |
332 |
|
|
|
333 |
|
|
|
334 |
|
|
|
335 |
|
|
|
336 |
|
|
CALL gwprofil_strato(nlon, nlev, kgwd, kdx, ktest, ikcrit, ikcrith, icrit, & |
337 |
|
|
ikenvh, iknu, iknu2, paphm1, zrho, zstab, ztfr, zvph, zri, ztau & |
338 |
|
576 |
, zdmod, znu, psig, pgam, pstd, ppic, pval) |
339 |
|
|
|
340 |
|
|
! * 5. Compute tendencies from waves stress profile. |
341 |
|
|
! Compute low level blocked flow drag. |
342 |
|
|
! * -------------------------------------------- |
343 |
|
|
|
344 |
|
|
|
345 |
|
|
|
346 |
|
|
|
347 |
|
|
! explicit solution at all levels for the gravity wave |
348 |
|
|
! implicit solution for the blocked levels |
349 |
|
|
|
350 |
✓✓ |
573120 |
DO jl = kidia, kfdia |
351 |
|
572544 |
zvidis(jl) = 0.0 |
352 |
|
572544 |
zdudt(jl) = 0.0 |
353 |
|
572544 |
zdvdt(jl) = 0.0 |
354 |
|
573120 |
zdtdt(jl) = 0.0 |
355 |
|
|
END DO |
356 |
|
|
|
357 |
|
|
|
358 |
✓✓ |
23040 |
DO jk = 1, klev |
359 |
|
|
|
360 |
|
|
|
361 |
|
|
! WAVE STRESS |
362 |
|
|
! ------------- |
363 |
|
|
|
364 |
|
|
|
365 |
✓✓ |
22352256 |
DO ji = kidia, kfdia |
366 |
|
|
|
367 |
✓✓ |
22351680 |
IF (ktest(ji)==1) THEN |
368 |
|
|
|
369 |
|
10423296 |
zdelp = paphm1(ji, jk+1) - paphm1(ji, jk) |
370 |
|
10423296 |
ztemp = -rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,klev+1)*zdelp) |
371 |
|
|
|
372 |
|
10423296 |
zdudt(ji) = (pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji) |
373 |
|
10423296 |
zdvdt(ji) = (pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji) |
374 |
|
|
|
375 |
|
|
! Control Overshoots |
376 |
|
|
|
377 |
|
|
|
378 |
✓✗ |
10423296 |
IF (jk>=nstra) THEN |
379 |
|
|
rover = 0.10 |
380 |
✓✓ |
10423296 |
IF (abs(zdudt(ji))>rover*abs(pum1(ji,jk))/ztmst) zdudt(ji) = rover* & |
381 |
|
11152 |
abs(pum1(ji,jk))/ztmst*zdudt(ji)/(abs(zdudt(ji))+1.E-10) |
382 |
✓✓ |
10423296 |
IF (abs(zdvdt(ji))>rover*abs(pvm1(ji,jk))/ztmst) zdvdt(ji) = rover* & |
383 |
|
23454 |
abs(pvm1(ji,jk))/ztmst*zdvdt(ji)/(abs(zdvdt(ji))+1.E-10) |
384 |
|
|
END IF |
385 |
|
|
|
386 |
|
|
rover = 0.25 |
387 |
|
10423296 |
zforc = sqrt(zdudt(ji)**2+zdvdt(ji)**2) |
388 |
|
10423296 |
ztend = sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/ztmst |
389 |
|
|
|
390 |
✗✓ |
10423296 |
IF (zforc>=rover*ztend) THEN |
391 |
|
|
zdudt(ji) = rover*ztend/zforc*zdudt(ji) |
392 |
|
|
zdvdt(ji) = rover*ztend/zforc*zdvdt(ji) |
393 |
|
|
END IF |
394 |
|
|
|
395 |
|
|
! BLOCKED FLOW DRAG: |
396 |
|
|
! ----------------- |
397 |
|
|
|
398 |
✓✓ |
10423296 |
IF (partdrag .GE. 2) THEN |
399 |
|
|
facpart=0. |
400 |
|
|
ELSE |
401 |
|
5211648 |
facpart=gkwake |
402 |
|
|
ENDIF |
403 |
|
|
|
404 |
|
|
|
405 |
✓✓ |
10423296 |
IF (jk>ikenvh(ji)) THEN |
406 |
|
1599026 |
zb = 1.0 - 0.18*pgam(ji) - 0.04*pgam(ji)**2 |
407 |
|
1599026 |
zc = 0.48*pgam(ji) + 0.3*pgam(ji)**2 |
408 |
|
1599026 |
zconb = 2.*ztmst*facpart*psig(ji)/(4.*pstd(ji)) |
409 |
|
1599026 |
zabsv = sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/2. |
410 |
|
1599026 |
zzd1 = zb*cos(zpsi(ji,jk))**2 + zc*sin(zpsi(ji,jk))**2 |
411 |
|
|
ratio = (cos(zpsi(ji,jk))**2+pgam(ji)*sin(zpsi(ji, & |
412 |
|
1599026 |
jk))**2)/(pgam(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2) |
413 |
|
1599026 |
zbet = max(0., 2.-1./ratio)*zconb*zzdep(ji, jk)*zzd1*zabsv |
414 |
|
|
|
415 |
|
|
! OPPOSED TO THE WIND |
416 |
|
|
|
417 |
|
1599026 |
zdudt(ji) = -pum1(ji, jk)/ztmst |
418 |
|
1599026 |
zdvdt(ji) = -pvm1(ji, jk)/ztmst |
419 |
|
|
|
420 |
|
|
! PERPENDICULAR TO THE SSO MAIN AXIS: |
421 |
|
|
|
422 |
|
|
! mod zdudt(ji)=-(pum1(ji,jk)*cos(pthe(ji)*rpi/180.) |
423 |
|
|
! mod * +pvm1(ji,jk)*sin(pthe(ji)*rpi/180.)) |
424 |
|
|
! mod * *cos(pthe(ji)*rpi/180.)/ztmst |
425 |
|
|
! mod zdvdt(ji)=-(pum1(ji,jk)*cos(pthe(ji)*rpi/180.) |
426 |
|
|
! mod * +pvm1(ji,jk)*sin(pthe(ji)*rpi/180.)) |
427 |
|
|
! mod * *sin(pthe(ji)*rpi/180.)/ztmst |
428 |
|
|
|
429 |
|
1599026 |
zdudt(ji) = zdudt(ji)*(zbet/(1.+zbet)) |
430 |
|
1599026 |
zdvdt(ji) = zdvdt(ji)*(zbet/(1.+zbet)) |
431 |
|
|
END IF |
432 |
|
10423296 |
pvom(ji, jk) = zdudt(ji) |
433 |
|
10423296 |
pvol(ji, jk) = zdvdt(ji) |
434 |
|
10423296 |
zust = pum1(ji, jk) + ztmst*zdudt(ji) |
435 |
|
10423296 |
zvst = pvm1(ji, jk) + ztmst*zdvdt(ji) |
436 |
|
10423296 |
zdis = 0.5*(pum1(ji,jk)**2+pvm1(ji,jk)**2-zust**2-zvst**2) |
437 |
|
10423296 |
zdedt(ji) = zdis/ztmst |
438 |
|
10423296 |
zvidis(ji) = zvidis(ji) + zdis*zdelp |
439 |
|
10423296 |
zdtdt(ji) = zdedt(ji)/rcpd |
440 |
|
|
|
441 |
|
|
! NO TENDENCIES ON TEMPERATURE ..... |
442 |
|
|
|
443 |
|
|
! Instead of, pte(ji,jk)=zdtdt(ji), due to mechanical dissipation |
444 |
|
|
|
445 |
|
10423296 |
pte(ji, jk) = 0.0 |
446 |
|
|
|
447 |
|
|
END IF |
448 |
|
|
|
449 |
|
|
END DO |
450 |
|
|
END DO |
451 |
|
|
|
452 |
|
576 |
RETURN |
453 |
|
|
END SUBROUTINE orodrag_strato |
454 |
|
576 |
SUBROUTINE orosetup_strato(nlon, nlev, ktest, kkcrit, kkcrith, kcrit, ksect, & |
455 |
|
576 |
kkhlim, kkenvh, kknu, kknu2, paphm1, papm1, pum1, pvm1, ptm1, pgeom1, & |
456 |
|
|
pstd, prho, pri, pstab, ptau, pvph, ppsi, pzdep, pulow, pvlow, ptheta, & |
457 |
|
576 |
pgam, pmea, ppic, pval, pnu, pd1, pd2, pdmod) |
458 |
|
|
|
459 |
|
|
! **** *gwsetup* |
460 |
|
|
|
461 |
|
|
! purpose. |
462 |
|
|
! -------- |
463 |
|
|
! SET-UP THE ESSENTIAL PARAMETERS OF THE SSO DRAG SCHEME: |
464 |
|
|
! DEPTH OF LOW WBLOCKED LAYER, LOW-LEVEL FLOW, BACKGROUND |
465 |
|
|
! STRATIFICATION..... |
466 |
|
|
|
467 |
|
|
! ** interface. |
468 |
|
|
! ---------- |
469 |
|
|
! from *orodrag* |
470 |
|
|
|
471 |
|
|
! explicit arguments : |
472 |
|
|
! -------------------- |
473 |
|
|
! ==== inputs === |
474 |
|
|
|
475 |
|
|
! nlon----input-I-Total number of horizontal points that get into physics |
476 |
|
|
! nlev----input-I-Number of vertical levels |
477 |
|
|
! ktest--input-I: Flags to indicate active points |
478 |
|
|
|
479 |
|
|
! ptsphy--input-R-Time-step (s) |
480 |
|
|
! paphm1--input-R: pressure at model 1/2 layer |
481 |
|
|
! papm1---input-R: pressure at model layer |
482 |
|
|
! pgeom1--input-R: Altitude of layer above ground |
483 |
|
|
! ptm1, pum1, pvm1--R-: t, u and v |
484 |
|
|
! pmea----input-R-Mean Orography (m) |
485 |
|
|
! pstd----input-R-SSO standard deviation (m) |
486 |
|
|
! psig----input-R-SSO slope |
487 |
|
|
! pgam----input-R-SSO Anisotropy |
488 |
|
|
! pthe----input-R-SSO Angle |
489 |
|
|
! ppic----input-R-SSO Peacks elevation (m) |
490 |
|
|
! pval----input-R-SSO Valleys elevation (m) |
491 |
|
|
|
492 |
|
|
! ==== outputs === |
493 |
|
|
! pulow, pvlow -output-R: Low-level wind |
494 |
|
|
! kkcrit----I-: Security value for top of low level flow |
495 |
|
|
! kcrit-----I-: Critical level |
496 |
|
|
! ksect-----I-: Not used |
497 |
|
|
! kkhlim----I-: Not used |
498 |
|
|
! kkenvh----I-: Top of blocked flow layer |
499 |
|
|
! kknu------I-: Layer that sees mountain peacks |
500 |
|
|
! kknu2-----I-: Layer that sees mountain peacks above mountain mean |
501 |
|
|
! kknub-----I-: Layer that sees mountain mean above valleys |
502 |
|
|
! prho------R-: Density at 1/2 layers |
503 |
|
|
! pri-------R-: Background Richardson Number, Wind shear measured along GW |
504 |
|
|
! stress |
505 |
|
|
! pstab-----R-: Brunt-Vaisala freq. at 1/2 layers |
506 |
|
|
! pvph------R-: Wind in plan of GW stress, Half levels. |
507 |
|
|
! ppsi------R-: Angle between low level wind and SS0 main axis. |
508 |
|
|
! pd1-------R-| Compared the ratio of the stress |
509 |
|
|
! pd2-------R-| that is along the wind to that Normal to it. |
510 |
|
|
! pdi define the plane of low level stress |
511 |
|
|
! compared to the low level wind. |
512 |
|
|
! see p. 108 Lott & Miller (1997). |
513 |
|
|
! pdmod-----R-: Norme of pdi |
514 |
|
|
|
515 |
|
|
! === local arrays === |
516 |
|
|
|
517 |
|
|
! zvpf------R-: Wind projected in the plan of the low-level stress. |
518 |
|
|
|
519 |
|
|
! ==== outputs === |
520 |
|
|
|
521 |
|
|
! implicit arguments : none |
522 |
|
|
! -------------------- |
523 |
|
|
|
524 |
|
|
! method. |
525 |
|
|
! ------- |
526 |
|
|
|
527 |
|
|
|
528 |
|
|
! externals. |
529 |
|
|
! ---------- |
530 |
|
|
|
531 |
|
|
|
532 |
|
|
! reference. |
533 |
|
|
! ---------- |
534 |
|
|
|
535 |
|
|
! see ecmwf research department documentation of the "i.f.s." |
536 |
|
|
|
537 |
|
|
! author. |
538 |
|
|
! ------- |
539 |
|
|
|
540 |
|
|
! modifications. |
541 |
|
|
! -------------- |
542 |
|
|
! f.lott for the new-gwdrag scheme november 1993 |
543 |
|
|
|
544 |
|
|
! ----------------------------------------------------------------------- |
545 |
|
|
USE dimphy |
546 |
|
|
IMPLICIT NONE |
547 |
|
|
|
548 |
|
|
|
549 |
|
|
include "YOMCST.h" |
550 |
|
|
include "YOEGWD.h" |
551 |
|
|
|
552 |
|
|
! ----------------------------------------------------------------------- |
553 |
|
|
|
554 |
|
|
! * 0.1 arguments |
555 |
|
|
! --------- |
556 |
|
|
|
557 |
|
|
INTEGER nlon, nlev |
558 |
|
|
INTEGER kkcrit(nlon), kkcrith(nlon), kcrit(nlon), ksect(nlon), & |
559 |
|
|
kkhlim(nlon), ktest(nlon), kkenvh(nlon) |
560 |
|
|
|
561 |
|
|
|
562 |
|
|
REAL paphm1(nlon, klev+1), papm1(nlon, klev), pum1(nlon, klev), & |
563 |
|
|
pvm1(nlon, klev), ptm1(nlon, klev), pgeom1(nlon, klev), & |
564 |
|
|
prho(nlon, klev+1), pri(nlon, klev+1), pstab(nlon, klev+1), & |
565 |
|
|
ptau(nlon, klev+1), pvph(nlon, klev+1), ppsi(nlon, klev+1), & |
566 |
|
|
pzdep(nlon, klev) |
567 |
|
|
REAL pulow(nlon), pvlow(nlon), ptheta(nlon), pgam(nlon), pnu(nlon), & |
568 |
|
|
pd1(nlon), pd2(nlon), pdmod(nlon) |
569 |
|
|
REAL pstd(nlon), pmea(nlon), ppic(nlon), pval(nlon) |
570 |
|
|
|
571 |
|
|
! ----------------------------------------------------------------------- |
572 |
|
|
|
573 |
|
|
! * 0.2 local arrays |
574 |
|
|
! ------------ |
575 |
|
|
|
576 |
|
|
|
577 |
|
|
INTEGER ilevh, jl, jk |
578 |
|
|
REAL zcons1, zcons2, zhgeo, zu, zphi |
579 |
|
|
REAL zvt1, zvt2, zdwind, zwind, zdelp |
580 |
|
|
REAL zstabm, zstabp, zrhom, zrhop |
581 |
|
|
LOGICAL lo |
582 |
|
1152 |
LOGICAL ll1(klon, klev+1) |
583 |
|
1152 |
INTEGER kknu(klon), kknu2(klon), kknub(klon), kknul(klon), kentp(klon), & |
584 |
|
1152 |
ncount(klon) |
585 |
|
|
|
586 |
|
1152 |
REAL zhcrit(klon, klev), zvpf(klon, klev), zdp(klon, klev) |
587 |
|
1152 |
REAL znorm(klon), zb(klon), zc(klon), zulow(klon), zvlow(klon), znup(klon), & |
588 |
|
1152 |
znum(klon) |
589 |
|
|
|
590 |
|
|
! ------------------------------------------------------------------ |
591 |
|
|
|
592 |
|
|
! * 1. initialization |
593 |
|
|
! -------------- |
594 |
|
|
|
595 |
|
|
! PRINT *,' in orosetup' |
596 |
|
|
|
597 |
|
|
! ------------------------------------------------------------------ |
598 |
|
|
|
599 |
|
|
! * 1.1 computational constants |
600 |
|
|
! ----------------------- |
601 |
|
|
|
602 |
|
|
|
603 |
|
576 |
ilevh = klev/3 |
604 |
|
|
|
605 |
|
576 |
zcons1 = 1./rd |
606 |
|
576 |
zcons2 = rg**2/rcpd |
607 |
|
|
|
608 |
|
|
! ------------------------------------------------------------------ |
609 |
|
|
|
610 |
|
|
! * 2. |
611 |
|
|
! -------------- |
612 |
|
|
|
613 |
|
|
|
614 |
|
|
! ------------------------------------------------------------------ |
615 |
|
|
|
616 |
|
|
! * 2.1 define low level wind, project winds in plane of |
617 |
|
|
! * low level wind, determine sector in which to take |
618 |
|
|
! * the variance and set indicator for critical levels. |
619 |
|
|
|
620 |
|
|
|
621 |
|
|
|
622 |
✓✓ |
573120 |
DO jl = kidia, kfdia |
623 |
|
572544 |
kknu(jl) = klev |
624 |
|
572544 |
kknu2(jl) = klev |
625 |
|
572544 |
kknub(jl) = klev |
626 |
|
572544 |
kknul(jl) = klev |
627 |
|
572544 |
pgam(jl) = max(pgam(jl), gtsec) |
628 |
|
576 |
ll1(jl, klev+1) = .FALSE. |
629 |
|
|
END DO |
630 |
|
|
|
631 |
|
|
! Ajouter une initialisation (L. Li, le 23fev99): |
632 |
|
|
|
633 |
✓✓ |
16128 |
DO jk = klev, ilevh, -1 |
634 |
✓✓ |
15474816 |
DO jl = kidia, kfdia |
635 |
|
15474240 |
ll1(jl, jk) = .FALSE. |
636 |
|
|
END DO |
637 |
|
|
END DO |
638 |
|
|
|
639 |
|
|
! * define top of low level flow |
640 |
|
|
! ---------------------------- |
641 |
✓✓ |
16128 |
DO jk = klev, ilevh, -1 |
642 |
✓✓ |
15474816 |
DO jl = kidia, kfdia |
643 |
✓✓ |
15474240 |
IF (ktest(jl)==1) THEN |
644 |
|
7216128 |
lo = (paphm1(jl,jk)/paphm1(jl,klev+1)) >= gsigcr |
645 |
|
7216128 |
IF (lo) THEN |
646 |
|
1906012 |
kkcrit(jl) = jk |
647 |
|
|
END IF |
648 |
|
7216128 |
zhcrit(jl, jk) = ppic(jl) - pval(jl) |
649 |
|
7216128 |
zhgeo = pgeom1(jl, jk)/rg |
650 |
|
7216128 |
ll1(jl, jk) = (zhgeo>zhcrit(jl,jk)) |
651 |
✓✓ |
7216128 |
IF (ll1(jl,jk) .NEQV. ll1(jl,jk+1)) THEN |
652 |
|
267264 |
kknu(jl) = jk |
653 |
|
|
END IF |
654 |
✓✓ |
7216128 |
IF (.NOT. ll1(jl,ilevh)) kknu(jl) = ilevh |
655 |
|
|
END IF |
656 |
|
|
END DO |
657 |
|
|
END DO |
658 |
✓✓ |
16128 |
DO jk = klev, ilevh, -1 |
659 |
✓✓ |
15474816 |
DO jl = kidia, kfdia |
660 |
✓✓ |
15474240 |
IF (ktest(jl)==1) THEN |
661 |
|
7216128 |
zhcrit(jl, jk) = ppic(jl) - pmea(jl) |
662 |
|
7216128 |
zhgeo = pgeom1(jl, jk)/rg |
663 |
|
7216128 |
ll1(jl, jk) = (zhgeo>zhcrit(jl,jk)) |
664 |
✓✓ |
7216128 |
IF (ll1(jl,jk) .NEQV. ll1(jl,jk+1)) THEN |
665 |
|
267264 |
kknu2(jl) = jk |
666 |
|
|
END IF |
667 |
✗✓ |
7216128 |
IF (.NOT. ll1(jl,ilevh)) kknu2(jl) = ilevh |
668 |
|
|
END IF |
669 |
|
|
END DO |
670 |
|
|
END DO |
671 |
✓✓ |
16128 |
DO jk = klev, ilevh, -1 |
672 |
✓✓ |
15474816 |
DO jl = kidia, kfdia |
673 |
✓✓ |
15474240 |
IF (ktest(jl)==1) THEN |
674 |
|
7216128 |
zhcrit(jl, jk) = amin1(ppic(jl)-pmea(jl), pmea(jl)-pval(jl)) |
675 |
|
7216128 |
zhgeo = pgeom1(jl, jk)/rg |
676 |
|
7216128 |
ll1(jl, jk) = (zhgeo>zhcrit(jl,jk)) |
677 |
✓✓ |
7216128 |
IF (ll1(jl,jk) .NEQV. ll1(jl,jk+1)) THEN |
678 |
|
267264 |
kknub(jl) = jk |
679 |
|
|
END IF |
680 |
✗✓ |
7216128 |
IF (.NOT. ll1(jl,ilevh)) kknub(jl) = ilevh |
681 |
|
|
END IF |
682 |
|
|
END DO |
683 |
|
|
END DO |
684 |
|
|
|
685 |
✓✓ |
573120 |
DO jl = kidia, kfdia |
686 |
✓✓ |
573120 |
IF (ktest(jl)==1) THEN |
687 |
|
267264 |
kknu(jl) = min(kknu(jl), nktopg) |
688 |
|
267264 |
kknu2(jl) = min(kknu2(jl), nktopg) |
689 |
|
267264 |
kknub(jl) = min(kknub(jl), nktopg) |
690 |
|
267264 |
kknul(jl) = klev |
691 |
|
|
END IF |
692 |
|
|
END DO |
693 |
|
|
|
694 |
|
|
! c* initialize various arrays |
695 |
|
|
|
696 |
✓✓ |
573120 |
DO jl = kidia, kfdia |
697 |
|
572544 |
prho(jl, klev+1) = 0.0 |
698 |
|
|
! ym correction en attendant mieux |
699 |
|
572544 |
prho(jl, 1) = 0.0 |
700 |
|
572544 |
pstab(jl, klev+1) = 0.0 |
701 |
|
572544 |
pstab(jl, 1) = 0.0 |
702 |
|
572544 |
pri(jl, klev+1) = 9999.0 |
703 |
|
572544 |
ppsi(jl, klev+1) = 0.0 |
704 |
|
572544 |
pri(jl, 1) = 0.0 |
705 |
|
572544 |
pvph(jl, 1) = 0.0 |
706 |
|
572544 |
pvph(jl, klev+1) = 0.0 |
707 |
|
|
! ym correction en attendant mieux |
708 |
|
|
! ym pvph(jl,klev) =0.0 |
709 |
|
572544 |
pulow(jl) = 0.0 |
710 |
|
572544 |
pvlow(jl) = 0.0 |
711 |
|
572544 |
zulow(jl) = 0.0 |
712 |
|
572544 |
zvlow(jl) = 0.0 |
713 |
|
572544 |
kkcrith(jl) = klev |
714 |
|
572544 |
kkenvh(jl) = klev |
715 |
|
572544 |
kentp(jl) = klev |
716 |
|
572544 |
kcrit(jl) = 1 |
717 |
|
572544 |
ncount(jl) = 0 |
718 |
|
573120 |
ll1(jl, klev+1) = .FALSE. |
719 |
|
|
END DO |
720 |
|
|
|
721 |
|
|
! * define flow density and stratification (rho and N2) |
722 |
|
|
! at semi layers. |
723 |
|
|
! ------------------------------------------------------- |
724 |
|
|
|
725 |
✓✓ |
22464 |
DO jk = klev, 2, -1 |
726 |
✓✓ |
21779136 |
DO jl = kidia, kfdia |
727 |
✓✓ |
21778560 |
IF (ktest(jl)==1) THEN |
728 |
|
10156032 |
zdp(jl, jk) = papm1(jl, jk) - papm1(jl, jk-1) |
729 |
|
10156032 |
prho(jl, jk) = 2.*paphm1(jl, jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1)) |
730 |
|
|
pstab(jl, jk) = 2.*zcons2/(ptm1(jl,jk)+ptm1(jl,jk-1))* & |
731 |
|
10156032 |
(1.-rcpd*prho(jl,jk)*(ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk)) |
732 |
|
10156032 |
pstab(jl, jk) = max(pstab(jl,jk), gssec) |
733 |
|
|
END IF |
734 |
|
|
END DO |
735 |
|
|
END DO |
736 |
|
|
|
737 |
|
|
! ******************************************************************** |
738 |
|
|
|
739 |
|
|
! * define Low level flow (between ground and peacks-valleys) |
740 |
|
|
! --------------------------------------------------------- |
741 |
✓✓ |
16128 |
DO jk = klev, ilevh, -1 |
742 |
✓✓ |
15474816 |
DO jl = kidia, kfdia |
743 |
✓✓ |
15474240 |
IF (ktest(jl)==1) THEN |
744 |
✓✓✓✗
|
7216128 |
IF (jk>=kknu2(jl) .AND. jk<=kknul(jl)) THEN |
745 |
|
|
pulow(jl) = pulow(jl) + pum1(jl, jk)*(paphm1(jl,jk+1)-paphm1(jl,jk) & |
746 |
|
1936914 |
) |
747 |
|
|
pvlow(jl) = pvlow(jl) + pvm1(jl, jk)*(paphm1(jl,jk+1)-paphm1(jl,jk) & |
748 |
|
1936914 |
) |
749 |
|
|
pstab(jl, klev+1) = pstab(jl, klev+1) + pstab(jl, jk)*(paphm1(jl,jk & |
750 |
|
1936914 |
+1)-paphm1(jl,jk)) |
751 |
|
|
prho(jl, klev+1) = prho(jl, klev+1) + prho(jl, jk)*(paphm1(jl,jk+1) & |
752 |
|
1936914 |
-paphm1(jl,jk)) |
753 |
|
|
END IF |
754 |
|
|
END IF |
755 |
|
|
END DO |
756 |
|
|
END DO |
757 |
✓✓ |
573120 |
DO jl = kidia, kfdia |
758 |
✓✓ |
573120 |
IF (ktest(jl)==1) THEN |
759 |
|
267264 |
pulow(jl) = pulow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl))) |
760 |
|
267264 |
pvlow(jl) = pvlow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl))) |
761 |
|
267264 |
znorm(jl) = max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec) |
762 |
|
267264 |
pvph(jl, klev+1) = znorm(jl) |
763 |
|
|
pstab(jl, klev+1) = pstab(jl, klev+1)/(paphm1(jl,kknul(jl)+1)-paphm1(jl & |
764 |
|
267264 |
,kknu2(jl))) |
765 |
|
|
prho(jl, klev+1) = prho(jl, klev+1)/(paphm1(jl,kknul(jl)+1)-paphm1(jl, & |
766 |
|
267264 |
kknu2(jl))) |
767 |
|
|
END IF |
768 |
|
|
END DO |
769 |
|
|
|
770 |
|
|
|
771 |
|
|
! ******* setup orography orientation relative to the low level |
772 |
|
|
! wind and define parameters of the Anisotropic wave stress. |
773 |
|
|
|
774 |
✓✓ |
573120 |
DO jl = kidia, kfdia |
775 |
✓✓ |
573120 |
IF (ktest(jl)==1) THEN |
776 |
✓✓ |
122938 |
lo = (pulow(jl)<gvsec) .AND. (pulow(jl)>=-gvsec) |
777 |
|
|
IF (lo) THEN |
778 |
|
4252 |
zu = pulow(jl) + 2.*gvsec |
779 |
|
|
ELSE |
780 |
|
|
zu = pulow(jl) |
781 |
|
|
END IF |
782 |
|
267264 |
zphi = atan(pvlow(jl)/zu) |
783 |
|
267264 |
ppsi(jl, klev+1) = ptheta(jl)*rpi/180. - zphi |
784 |
|
267264 |
zb(jl) = 1. - 0.18*pgam(jl) - 0.04*pgam(jl)**2 |
785 |
|
267264 |
zc(jl) = 0.48*pgam(jl) + 0.3*pgam(jl)**2 |
786 |
|
267264 |
pd1(jl) = zb(jl) - (zb(jl)-zc(jl))*(sin(ppsi(jl,klev+1))**2) |
787 |
|
267264 |
pd2(jl) = (zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))*cos(ppsi(jl,klev+1)) |
788 |
|
267264 |
pdmod(jl) = sqrt(pd1(jl)**2+pd2(jl)**2) |
789 |
|
|
END IF |
790 |
|
|
END DO |
791 |
|
|
|
792 |
|
|
! ************ projet flow in plane of lowlevel stress ************* |
793 |
|
|
! ************ Find critical levels... ************* |
794 |
|
|
|
795 |
✓✓ |
23040 |
DO jk = 1, klev |
796 |
✓✓ |
22352256 |
DO jl = kidia, kfdia |
797 |
✓✓ |
22329216 |
IF (ktest(jl)==1) THEN |
798 |
|
10423296 |
zvt1 = pulow(jl)*pum1(jl, jk) + pvlow(jl)*pvm1(jl, jk) |
799 |
|
10423296 |
zvt2 = -pvlow(jl)*pum1(jl, jk) + pulow(jl)*pvm1(jl, jk) |
800 |
|
10423296 |
zvpf(jl, jk) = (zvt1*pd1(jl)+zvt2*pd2(jl))/(znorm(jl)*pdmod(jl)) |
801 |
|
|
END IF |
802 |
|
22329216 |
ptau(jl, jk) = 0.0 |
803 |
|
22329216 |
pzdep(jl, jk) = 0.0 |
804 |
|
22329216 |
ppsi(jl, jk) = 0.0 |
805 |
|
22351680 |
ll1(jl, jk) = .FALSE. |
806 |
|
|
END DO |
807 |
|
|
END DO |
808 |
✓✓ |
22464 |
DO jk = 2, klev |
809 |
✓✓ |
21779136 |
DO jl = kidia, kfdia |
810 |
✓✓ |
21778560 |
IF (ktest(jl)==1) THEN |
811 |
|
10156032 |
zdp(jl, jk) = papm1(jl, jk) - papm1(jl, jk-1) |
812 |
|
|
pvph(jl, jk) = ((paphm1(jl,jk)-papm1(jl,jk-1))*zvpf(jl,jk)+(papm1(jl, & |
813 |
|
10156032 |
jk)-paphm1(jl,jk))*zvpf(jl,jk-1))/zdp(jl, jk) |
814 |
|
10156032 |
IF (pvph(jl,jk)<gvsec) THEN |
815 |
|
3030215 |
pvph(jl, jk) = gvsec |
816 |
|
3030215 |
kcrit(jl) = jk |
817 |
|
|
END IF |
818 |
|
|
END IF |
819 |
|
|
END DO |
820 |
|
|
END DO |
821 |
|
|
|
822 |
|
|
! * 2.3 mean flow richardson number. |
823 |
|
|
|
824 |
|
|
|
825 |
✓✓ |
22464 |
DO jk = 2, klev |
826 |
✓✓ |
21779136 |
DO jl = kidia, kfdia |
827 |
✓✓ |
21778560 |
IF (ktest(jl)==1) THEN |
828 |
|
10156032 |
zdwind = max(abs(zvpf(jl,jk)-zvpf(jl,jk-1)), gvsec) |
829 |
|
10156032 |
pri(jl, jk) = pstab(jl, jk)*(zdp(jl,jk)/(rg*prho(jl,jk)*zdwind))**2 |
830 |
|
10156032 |
pri(jl, jk) = max(pri(jl,jk), grcrit) |
831 |
|
|
END IF |
832 |
|
|
END DO |
833 |
|
|
END DO |
834 |
|
|
|
835 |
|
|
|
836 |
|
|
|
837 |
|
|
! * define top of 'envelope' layer |
838 |
|
|
! ---------------------------- |
839 |
|
|
|
840 |
✓✓ |
573120 |
DO jl = kidia, kfdia |
841 |
|
572544 |
pnu(jl) = 0.0 |
842 |
|
573120 |
znum(jl) = 0.0 |
843 |
|
|
END DO |
844 |
|
|
|
845 |
✓✓ |
21888 |
DO jk = 2, klev - 1 |
846 |
✓✓ |
21206016 |
DO jl = kidia, kfdia |
847 |
|
|
|
848 |
✓✓ |
21205440 |
IF (ktest(jl)==1) THEN |
849 |
|
|
|
850 |
✓✓ |
9888768 |
IF (jk>=kknu2(jl)) THEN |
851 |
|
|
|
852 |
|
1669650 |
znum(jl) = pnu(jl) |
853 |
|
|
zwind = (pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/ & |
854 |
|
1669650 |
max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec) |
855 |
|
1669650 |
zwind = max(sqrt(zwind**2), gvsec) |
856 |
|
1669650 |
zdelp = paphm1(jl, jk+1) - paphm1(jl, jk) |
857 |
|
1669650 |
zstabm = sqrt(max(pstab(jl,jk),gssec)) |
858 |
|
1669650 |
zstabp = sqrt(max(pstab(jl,jk+1),gssec)) |
859 |
|
1669650 |
zrhom = prho(jl, jk) |
860 |
|
1669650 |
zrhop = prho(jl, jk+1) |
861 |
|
|
pnu(jl) = pnu(jl) + (zdelp/rg)*((zstabp/zrhop+zstabm/zrhom)/2.)/ & |
862 |
|
1669650 |
zwind |
863 |
✓✓✓✗
|
336931 |
IF ((znum(jl)<=gfrcrit) .AND. (pnu(jl)>gfrcrit) .AND. (kkenvh( & |
864 |
|
266307 |
jl)==klev)) kkenvh(jl) = jk |
865 |
|
|
|
866 |
|
|
END IF |
867 |
|
|
|
868 |
|
|
END IF |
869 |
|
|
|
870 |
|
|
END DO |
871 |
|
|
END DO |
872 |
|
|
|
873 |
|
|
! calculation of a dynamical mixing height for when the waves |
874 |
|
|
! BREAK AT LOW LEVEL: The drag will be repartited over |
875 |
|
|
! a depths that depends on waves vertical wavelength, |
876 |
|
|
! not just between two adjacent model layers. |
877 |
|
|
! of gravity waves: |
878 |
|
|
|
879 |
✓✓ |
573120 |
DO jl = kidia, kfdia |
880 |
|
572544 |
znup(jl) = 0.0 |
881 |
|
573120 |
znum(jl) = 0.0 |
882 |
|
|
END DO |
883 |
|
|
|
884 |
✓✓ |
21888 |
DO jk = klev - 1, 2, -1 |
885 |
✓✓ |
21206016 |
DO jl = kidia, kfdia |
886 |
|
|
|
887 |
✓✓ |
21205440 |
IF (ktest(jl)==1) THEN |
888 |
|
|
|
889 |
|
9888768 |
znum(jl) = znup(jl) |
890 |
|
|
zwind = (pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/ & |
891 |
|
9888768 |
max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec) |
892 |
|
9888768 |
zwind = max(sqrt(zwind**2), gvsec) |
893 |
|
9888768 |
zdelp = paphm1(jl, jk+1) - paphm1(jl, jk) |
894 |
|
9888768 |
zstabm = sqrt(max(pstab(jl,jk),gssec)) |
895 |
|
9888768 |
zstabp = sqrt(max(pstab(jl,jk+1),gssec)) |
896 |
|
9888768 |
zrhom = prho(jl, jk) |
897 |
|
9888768 |
zrhop = prho(jl, jk+1) |
898 |
|
|
znup(jl) = znup(jl) + (zdelp/rg)*((zstabp/zrhop+zstabm/zrhom)/2.)/ & |
899 |
|
9888768 |
zwind |
900 |
✓✓✓✗
|
698087 |
IF ((znum(jl)<=rpi/4.) .AND. (znup(jl)>rpi/4.) .AND. (kkcrith( & |
901 |
|
267264 |
jl)==klev)) kkcrith(jl) = jk |
902 |
|
|
|
903 |
|
|
END IF |
904 |
|
|
|
905 |
|
|
END DO |
906 |
|
|
END DO |
907 |
|
|
|
908 |
✓✓ |
573120 |
DO jl = kidia, kfdia |
909 |
✓✓ |
573120 |
IF (ktest(jl)==1) THEN |
910 |
|
267264 |
kkcrith(jl) = max0(kkcrith(jl), ilevh*2) |
911 |
|
267264 |
kkcrith(jl) = max0(kkcrith(jl), kknu(jl)) |
912 |
✓✓ |
267264 |
IF (kcrit(jl)>=kkcrith(jl)) kcrit(jl) = 1 |
913 |
|
|
END IF |
914 |
|
|
END DO |
915 |
|
|
|
916 |
|
|
! directional info for flow blocking ************************* |
917 |
|
|
|
918 |
✓✓ |
23040 |
DO jk = 1, klev |
919 |
✓✓ |
22352256 |
DO jl = kidia, kfdia |
920 |
✓✓ |
22351680 |
IF (ktest(jl)==1) THEN |
921 |
✓✓ |
2754132 |
lo = (pum1(jl,jk)<gvsec) .AND. (pum1(jl,jk)>=-gvsec) |
922 |
|
|
IF (lo) THEN |
923 |
|
62361 |
zu = pum1(jl, jk) + 2.*gvsec |
924 |
|
|
ELSE |
925 |
|
|
zu = pum1(jl, jk) |
926 |
|
|
END IF |
927 |
|
10423296 |
zphi = atan(pvm1(jl,jk)/zu) |
928 |
|
10423296 |
ppsi(jl, jk) = ptheta(jl)*rpi/180. - zphi |
929 |
|
|
END IF |
930 |
|
|
END DO |
931 |
|
|
END DO |
932 |
|
|
|
933 |
|
|
! forms the vertical 'leakiness' ************************** |
934 |
|
|
|
935 |
✓✓ |
16128 |
DO jk = ilevh, klev |
936 |
✓✓ |
15474816 |
DO jl = kidia, kfdia |
937 |
✓✓ |
15474240 |
IF (ktest(jl)==1) THEN |
938 |
|
7216128 |
pzdep(jl, jk) = 0 |
939 |
✓✓✓✓
|
7216128 |
IF (jk>=kkenvh(jl) .AND. kkenvh(jl)/=klev) THEN |
940 |
|
|
pzdep(jl, jk) = (pgeom1(jl,kkenvh(jl))-pgeom1(jl,jk))/ & |
941 |
|
1865333 |
(pgeom1(jl,kkenvh(jl))-pgeom1(jl,klev)) |
942 |
|
|
END IF |
943 |
|
|
END IF |
944 |
|
|
END DO |
945 |
|
|
END DO |
946 |
|
|
|
947 |
|
576 |
RETURN |
948 |
|
|
END SUBROUTINE orosetup_strato |
949 |
|
576 |
SUBROUTINE gwstress_strato(nlon, nlev, kkcrit, ksect, kkhlim, ktest, kkcrith, & |
950 |
|
576 |
kcrit, kkenvh, kknu, prho, pstab, pvph, pstd, psig, pmea, ppic, pval, & |
951 |
|
|
ptfr, ptau, pgeom1, pgamma, pd1, pd2, pdmod, pnu) |
952 |
|
|
|
953 |
|
|
! **** *gwstress* |
954 |
|
|
|
955 |
|
|
! purpose. |
956 |
|
|
! -------- |
957 |
|
|
! Compute the surface stress due to Gravity Waves, according |
958 |
|
|
! to the Phillips (1979) theory of 3-D flow above |
959 |
|
|
! anisotropic elliptic ridges. |
960 |
|
|
|
961 |
|
|
! The stress is reduced two account for cut-off flow over |
962 |
|
|
! hill. The flow only see that part of the ridge located |
963 |
|
|
! above the blocked layer (see zeff). |
964 |
|
|
|
965 |
|
|
! ** interface. |
966 |
|
|
! ---------- |
967 |
|
|
! call *gwstress* from *gwdrag* |
968 |
|
|
|
969 |
|
|
! explicit arguments : |
970 |
|
|
! -------------------- |
971 |
|
|
! ==== inputs === |
972 |
|
|
! ==== outputs === |
973 |
|
|
|
974 |
|
|
! implicit arguments : none |
975 |
|
|
! -------------------- |
976 |
|
|
|
977 |
|
|
! method. |
978 |
|
|
! ------- |
979 |
|
|
|
980 |
|
|
|
981 |
|
|
! externals. |
982 |
|
|
! ---------- |
983 |
|
|
|
984 |
|
|
|
985 |
|
|
! reference. |
986 |
|
|
! ---------- |
987 |
|
|
|
988 |
|
|
! LOTT and MILLER (1997) & LOTT (1999) |
989 |
|
|
|
990 |
|
|
! author. |
991 |
|
|
! ------- |
992 |
|
|
|
993 |
|
|
! modifications. |
994 |
|
|
! -------------- |
995 |
|
|
! f. lott put the new gwd on ifs 22/11/93 |
996 |
|
|
|
997 |
|
|
! ----------------------------------------------------------------------- |
998 |
✓✓✓✓
✓✓✓✓
✓✓✓✓
|
63917920 |
USE dimphy |
999 |
|
|
IMPLICIT NONE |
1000 |
|
|
|
1001 |
|
|
include "YOMCST.h" |
1002 |
|
|
include "YOEGWD.h" |
1003 |
|
|
|
1004 |
|
|
! ----------------------------------------------------------------------- |
1005 |
|
|
|
1006 |
|
|
! * 0.1 arguments |
1007 |
|
|
! --------- |
1008 |
|
|
|
1009 |
|
|
INTEGER nlon, nlev |
1010 |
|
|
INTEGER kkcrit(nlon), kkcrith(nlon), kcrit(nlon), ksect(nlon), & |
1011 |
|
|
kkhlim(nlon), ktest(nlon), kkenvh(nlon), kknu(nlon) |
1012 |
|
|
|
1013 |
|
|
REAL prho(nlon, nlev+1), pstab(nlon, nlev+1), ptau(nlon, nlev+1), & |
1014 |
|
|
pvph(nlon, nlev+1), ptfr(nlon), pgeom1(nlon, nlev), pstd(nlon) |
1015 |
|
|
|
1016 |
|
|
REAL pd1(nlon), pd2(nlon), pnu(nlon), psig(nlon), pgamma(nlon) |
1017 |
|
|
REAL pmea(nlon), ppic(nlon), pval(nlon) |
1018 |
|
|
REAL pdmod(nlon) |
1019 |
|
|
|
1020 |
|
|
! ----------------------------------------------------------------------- |
1021 |
|
|
|
1022 |
|
|
! * 0.2 local arrays |
1023 |
|
|
! ------------ |
1024 |
|
|
! zeff--real: effective height seen by the flow when there is blocking |
1025 |
|
|
|
1026 |
|
|
INTEGER jl |
1027 |
|
|
REAL zeff |
1028 |
|
|
|
1029 |
|
|
! ----------------------------------------------------------------------- |
1030 |
|
|
|
1031 |
|
|
! * 0.3 functions |
1032 |
|
|
! --------- |
1033 |
|
|
! ------------------------------------------------------------------ |
1034 |
|
|
|
1035 |
|
|
! * 1. initialization |
1036 |
|
|
! -------------- |
1037 |
|
|
|
1038 |
|
|
! PRINT *,' in gwstress' |
1039 |
|
|
|
1040 |
|
|
! * 3.1 gravity wave stress. |
1041 |
|
|
|
1042 |
|
|
|
1043 |
|
|
|
1044 |
✓✓ |
573120 |
DO jl = kidia, kfdia |
1045 |
✓✓ |
573120 |
IF (ktest(jl)==1) THEN |
1046 |
|
|
|
1047 |
|
|
! effective mountain height above the blocked flow |
1048 |
|
|
|
1049 |
|
267264 |
zeff = ppic(jl) - pval(jl) |
1050 |
✓✓ |
267264 |
IF (kkenvh(jl)<klev) THEN |
1051 |
|
266307 |
zeff = amin1(gfrcrit*pvph(jl,klev+1)/sqrt(pstab(jl,klev+1)), zeff) |
1052 |
|
|
END IF |
1053 |
|
|
|
1054 |
|
|
|
1055 |
|
|
ptau(jl, klev+1) = gkdrag*prho(jl, klev+1)*psig(jl)*pdmod(jl)/4./ & |
1056 |
|
267264 |
pstd(jl)*pvph(jl, klev+1)*sqrt(pstab(jl,klev+1))*zeff**2 |
1057 |
|
|
|
1058 |
|
|
|
1059 |
|
|
! too small value of stress or low level flow include critical level |
1060 |
|
|
! or low level flow: gravity wave stress nul. |
1061 |
|
|
|
1062 |
|
|
! lo=(ptau(jl,klev+1).lt.gtsec).or.(kcrit(jl).ge.kknu(jl)) |
1063 |
|
|
! * .or.(pvph(jl,klev+1).lt.gvcrit) |
1064 |
|
|
! if(lo) ptau(jl,klev+1)=0.0 |
1065 |
|
|
|
1066 |
|
|
! print *,jl,ptau(jl,klev+1) |
1067 |
|
|
|
1068 |
|
|
ELSE |
1069 |
|
|
|
1070 |
|
305280 |
ptau(jl, klev+1) = 0.0 |
1071 |
|
|
|
1072 |
|
|
END IF |
1073 |
|
|
|
1074 |
|
|
END DO |
1075 |
|
|
|
1076 |
|
|
! write(21)(ptau(jl,klev+1),jl=kidia,kfdia) |
1077 |
|
|
|
1078 |
|
576 |
RETURN |
1079 |
|
|
END SUBROUTINE gwstress_strato |
1080 |
|
|
|
1081 |
|
576 |
SUBROUTINE gwprofil_strato(nlon, nlev, kgwd, kdx, ktest, kkcrit, kkcrith, & |
1082 |
|
576 |
kcrit, kkenvh, kknu, kknu2, paphm1, prho, pstab, ptfr, pvph, pri, ptau, & |
1083 |
|
|
pdmod, pnu, psig, pgamma, pstd, ppic, pval) |
1084 |
|
|
|
1085 |
|
|
! **** *gwprofil* |
1086 |
|
|
|
1087 |
|
|
! purpose. |
1088 |
|
|
! -------- |
1089 |
|
|
|
1090 |
|
|
! ** interface. |
1091 |
|
|
! ---------- |
1092 |
|
|
! from *gwdrag* |
1093 |
|
|
|
1094 |
|
|
! explicit arguments : |
1095 |
|
|
! -------------------- |
1096 |
|
|
! ==== inputs === |
1097 |
|
|
|
1098 |
|
|
! ==== outputs === |
1099 |
|
|
|
1100 |
|
|
! implicit arguments : none |
1101 |
|
|
! -------------------- |
1102 |
|
|
|
1103 |
|
|
! method: |
1104 |
|
|
! ------- |
1105 |
|
|
! the stress profile for gravity waves is computed as follows: |
1106 |
|
|
! it decreases linearly with heights from the ground |
1107 |
|
|
! to the low-level indicated by kkcrith, |
1108 |
|
|
! to simulates lee waves or |
1109 |
|
|
! low-level gravity wave breaking. |
1110 |
|
|
! above it is constant, except when the waves encounter a critical |
1111 |
|
|
! level (kcrit) or when they break. |
1112 |
|
|
! The stress is also uniformly distributed above the level |
1113 |
|
|
! nstra. |
1114 |
|
|
|
1115 |
|
267264 |
USE dimphy |
1116 |
|
|
IMPLICIT NONE |
1117 |
|
|
|
1118 |
|
|
include "YOMCST.h" |
1119 |
|
|
include "YOEGWD.h" |
1120 |
|
|
|
1121 |
|
|
! ----------------------------------------------------------------------- |
1122 |
|
|
|
1123 |
|
|
! * 0.1 ARGUMENTS |
1124 |
|
|
! --------- |
1125 |
|
|
|
1126 |
|
|
INTEGER nlon, nlev, kgwd |
1127 |
|
|
INTEGER kkcrit(nlon), kkcrith(nlon), kcrit(nlon), kdx(nlon), ktest(nlon), & |
1128 |
|
|
kkenvh(nlon), kknu(nlon), kknu2(nlon) |
1129 |
|
|
|
1130 |
|
|
REAL paphm1(nlon, nlev+1), pstab(nlon, nlev+1), prho(nlon, nlev+1), & |
1131 |
|
|
pvph(nlon, nlev+1), pri(nlon, nlev+1), ptfr(nlon), ptau(nlon, nlev+1) |
1132 |
|
|
|
1133 |
|
|
REAL pdmod(nlon), pnu(nlon), psig(nlon), pgamma(nlon), pstd(nlon), & |
1134 |
|
|
ppic(nlon), pval(nlon) |
1135 |
|
|
|
1136 |
|
|
! ----------------------------------------------------------------------- |
1137 |
|
|
|
1138 |
|
|
! * 0.2 local arrays |
1139 |
|
|
! ------------ |
1140 |
|
|
|
1141 |
|
|
INTEGER jl, jk |
1142 |
|
|
REAL zsqr, zalfa, zriw, zdel, zb, zalpha, zdz2n, zdelp, zdelpt |
1143 |
|
|
|
1144 |
|
1152 |
REAL zdz2(klon, klev), znorm(klon), zoro(klon) |
1145 |
|
1152 |
REAL ztau(klon, klev+1) |
1146 |
|
|
|
1147 |
|
|
! ----------------------------------------------------------------------- |
1148 |
|
|
|
1149 |
|
|
! * 1. INITIALIZATION |
1150 |
|
|
! -------------- |
1151 |
|
|
|
1152 |
|
|
! print *,' entree gwprofil' |
1153 |
|
|
|
1154 |
|
|
|
1155 |
|
|
! * COMPUTATIONAL CONSTANTS. |
1156 |
|
|
! ------------- ---------- |
1157 |
|
|
|
1158 |
✓✓ |
573120 |
DO jl = kidia, kfdia |
1159 |
✓✓ |
573120 |
IF (ktest(jl)==1) THEN |
1160 |
|
267264 |
zoro(jl) = psig(jl)*pdmod(jl)/4./pstd(jl) |
1161 |
|
267264 |
ztau(jl, klev+1) = ptau(jl, klev+1) |
1162 |
|
|
! print *,jl,ptau(jl,klev+1) |
1163 |
|
267264 |
ztau(jl, kkcrith(jl)) = grahilo*ptau(jl, klev+1) |
1164 |
|
|
END IF |
1165 |
|
|
END DO |
1166 |
|
|
|
1167 |
|
|
|
1168 |
✓✓ |
23616 |
DO jk = klev + 1, 1, -1 |
1169 |
|
|
! * 4.1 constant shear stress until top of the |
1170 |
|
|
! low-level breaking/trapped layer |
1171 |
|
|
|
1172 |
✓✓ |
22925376 |
DO jl = kidia, kfdia |
1173 |
✓✓ |
22924800 |
IF (ktest(jl)==1) THEN |
1174 |
✓✓ |
10690560 |
IF (jk>kkcrith(jl)) THEN |
1175 |
|
965351 |
zdelp = paphm1(jl, jk) - paphm1(jl, klev+1) |
1176 |
|
965351 |
zdelpt = paphm1(jl, kkcrith(jl)) - paphm1(jl, klev+1) |
1177 |
|
|
ptau(jl, jk) = ztau(jl, klev+1) + zdelp/zdelpt*(ztau(jl,kkcrith(jl) & |
1178 |
|
965351 |
)-ztau(jl,klev+1)) |
1179 |
|
|
ELSE |
1180 |
|
9725209 |
ptau(jl, jk) = ztau(jl, kkcrith(jl)) |
1181 |
|
|
END IF |
1182 |
|
|
END IF |
1183 |
|
|
END DO |
1184 |
|
|
|
1185 |
|
|
! * 4.15 constant shear stress until the top of the |
1186 |
|
|
! low level flow layer. |
1187 |
|
|
|
1188 |
|
|
|
1189 |
|
|
! * 4.2 wave displacement at next level. |
1190 |
|
|
|
1191 |
|
|
|
1192 |
|
|
END DO |
1193 |
|
|
|
1194 |
|
|
|
1195 |
|
|
! * 4.4 wave richardson number, new wave displacement |
1196 |
|
|
! * and stress: breaking evaluation and critical |
1197 |
|
|
! level |
1198 |
|
|
|
1199 |
|
|
|
1200 |
✓✓ |
23040 |
DO jk = klev, 1, -1 |
1201 |
|
|
|
1202 |
✓✓ |
22351680 |
DO jl = kidia, kfdia |
1203 |
✓✓ |
22351680 |
IF (ktest(jl)==1) THEN |
1204 |
|
10423296 |
znorm(jl) = prho(jl, jk)*sqrt(pstab(jl,jk))*pvph(jl, jk) |
1205 |
|
10423296 |
zdz2(jl, jk) = ptau(jl, jk)/amax1(znorm(jl), gssec)/zoro(jl) |
1206 |
|
|
END IF |
1207 |
|
|
END DO |
1208 |
|
|
|
1209 |
✓✓ |
22352256 |
DO jl = kidia, kfdia |
1210 |
✓✓ |
22351680 |
IF (ktest(jl)==1) THEN |
1211 |
✓✓ |
10423296 |
IF (jk<kkcrith(jl)) THEN |
1212 |
✓✓✓✓
|
9457945 |
IF ((ptau(jl,jk+1)<gtsec) .OR. (jk<=kcrit(jl))) THEN |
1213 |
|
4580593 |
ptau(jl, jk) = 0.0 |
1214 |
|
|
ELSE |
1215 |
|
4877352 |
zsqr = sqrt(pri(jl,jk)) |
1216 |
|
4877352 |
zalfa = sqrt(pstab(jl,jk)*zdz2(jl,jk))/pvph(jl, jk) |
1217 |
|
4877352 |
zriw = pri(jl, jk)*(1.-zalfa)/(1+zalfa*zsqr)**2 |
1218 |
✓✓ |
4877352 |
IF (zriw<grcrit) THEN |
1219 |
|
|
! print *,' breaking!!!',ptau(jl,jk) |
1220 |
|
801697 |
zdel = 4./zsqr/grcrit + 1./grcrit**2 + 4./grcrit |
1221 |
|
801697 |
zb = 1./grcrit + 2./zsqr |
1222 |
|
801697 |
zalpha = 0.5*(-zb+sqrt(zdel)) |
1223 |
|
801697 |
zdz2n = (pvph(jl,jk)*zalpha)**2/pstab(jl, jk) |
1224 |
|
801697 |
ptau(jl, jk) = znorm(jl)*zdz2n*zoro(jl) |
1225 |
|
|
END IF |
1226 |
|
|
|
1227 |
|
4877352 |
ptau(jl, jk) = amin1(ptau(jl,jk), ptau(jl,jk+1)) |
1228 |
|
|
|
1229 |
|
|
END IF |
1230 |
|
|
END IF |
1231 |
|
|
END IF |
1232 |
|
|
END DO |
1233 |
|
|
END DO |
1234 |
|
|
|
1235 |
|
|
! REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL |
1236 |
|
|
|
1237 |
✓✓ |
573120 |
DO jl = kidia, kfdia |
1238 |
✓✓ |
573120 |
IF (ktest(jl)==1) THEN |
1239 |
|
267264 |
ztau(jl, kkcrith(jl)-1) = ptau(jl, kkcrith(jl)-1) |
1240 |
|
267264 |
ztau(jl, nstra) = ptau(jl, nstra) |
1241 |
|
|
END IF |
1242 |
|
|
END DO |
1243 |
|
|
|
1244 |
✓✓ |
23040 |
DO jk = 1, klev |
1245 |
|
|
|
1246 |
✓✓ |
22351680 |
DO jl = kidia, kfdia |
1247 |
✓✓ |
22351680 |
IF (ktest(jl)==1) THEN |
1248 |
|
|
|
1249 |
✓✓ |
10423296 |
IF (jk>kkcrith(jl)-1) THEN |
1250 |
|
|
|
1251 |
|
965351 |
zdelp = paphm1(jl, jk) - paphm1(jl, klev+1) |
1252 |
|
965351 |
zdelpt = paphm1(jl, kkcrith(jl)-1) - paphm1(jl, klev+1) |
1253 |
|
|
ptau(jl, jk) = ztau(jl, klev+1) + (ztau(jl,kkcrith(jl)-1)-ztau(jl, & |
1254 |
|
965351 |
klev+1))*zdelp/zdelpt |
1255 |
|
|
|
1256 |
|
|
END IF |
1257 |
|
|
END IF |
1258 |
|
|
|
1259 |
|
|
END DO |
1260 |
|
|
|
1261 |
|
|
! REORGANISATION AT THE MODEL TOP.... |
1262 |
|
|
|
1263 |
✓✓ |
22352256 |
DO jl = kidia, kfdia |
1264 |
✓✓ |
22351680 |
IF (ktest(jl)==1) THEN |
1265 |
|
|
|
1266 |
✗✓ |
10423296 |
IF (jk<nstra) THEN |
1267 |
|
|
|
1268 |
|
|
zdelp = paphm1(jl, nstra) |
1269 |
|
|
zdelpt = paphm1(jl, jk) |
1270 |
|
|
ptau(jl, jk) = ztau(jl, nstra)*zdelpt/zdelp |
1271 |
|
|
! ptau(jl,jk)=ztau(jl,nstra) |
1272 |
|
|
|
1273 |
|
|
END IF |
1274 |
|
|
|
1275 |
|
|
END IF |
1276 |
|
|
|
1277 |
|
|
END DO |
1278 |
|
|
|
1279 |
|
|
|
1280 |
|
|
END DO |
1281 |
|
|
|
1282 |
|
|
|
1283 |
|
|
123 FORMAT (I4, 1X, 20(F6.3,1X)) |
1284 |
|
|
|
1285 |
|
|
|
1286 |
|
576 |
RETURN |
1287 |
|
|
END SUBROUTINE gwprofil_strato |
1288 |
|
22329504 |
SUBROUTINE lift_noro_strato(nlon, nlev, dtime, paprs, pplay, plat, pmea, & |
1289 |
|
288 |
pstd, psig, pgam, pthe, ppic, pval, kgwd, kdx, ktest, t, u, v, pulow, & |
1290 |
|
288 |
pvlow, pustr, pvstr, d_t, d_u, d_v) |
1291 |
|
|
|
1292 |
|
|
USE dimphy |
1293 |
|
|
IMPLICIT NONE |
1294 |
|
|
! ====================================================================== |
1295 |
|
|
! Auteur(s): F.Lott (LMD/CNRS) date: 19950201 |
1296 |
|
|
! Object: Mountain lift interface (enhanced vortex stretching). |
1297 |
|
|
! Made necessary because: |
1298 |
|
|
! 1. in the LMD-GCM Layers are from bottom to top, |
1299 |
|
|
! contrary to most European GCM. |
1300 |
|
|
! 2. the altitude above ground of each model layers |
1301 |
|
|
! needs to be known (variable zgeom) |
1302 |
|
|
! ====================================================================== |
1303 |
|
|
! Explicit Arguments: |
1304 |
|
|
! ================== |
1305 |
|
|
! nlon----input-I-Total number of horizontal points that get into physics |
1306 |
|
|
! nlev----input-I-Number of vertical levels |
1307 |
|
|
! dtime---input-R-Time-step (s) |
1308 |
|
|
! paprs---input-R-Pressure in semi layers (Pa) |
1309 |
|
|
! pplay---input-R-Pressure model-layers (Pa) |
1310 |
|
|
! t-------input-R-temperature (K) |
1311 |
|
|
! u-------input-R-Horizontal wind (m/s) |
1312 |
|
|
! v-------input-R-Meridional wind (m/s) |
1313 |
|
|
! pmea----input-R-Mean Orography (m) |
1314 |
|
|
! pstd----input-R-SSO standard deviation (m) |
1315 |
|
|
! psig----input-R-SSO slope |
1316 |
|
|
! pgam----input-R-SSO Anisotropy |
1317 |
|
|
! pthe----input-R-SSO Angle |
1318 |
|
|
! ppic----input-R-SSO Peacks elevation (m) |
1319 |
|
|
! pval----input-R-SSO Valleys elevation (m) |
1320 |
|
|
|
1321 |
|
|
! kgwd- -input-I: Total nb of points where the orography schemes are active |
1322 |
|
|
! ktest--input-I: Flags to indicate active points |
1323 |
|
|
! kdx----input-I: Locate the physical location of an active point. |
1324 |
|
|
|
1325 |
|
|
! pulow, pvlow -output-R: Low-level wind |
1326 |
|
|
! pustr, pvstr -output-R: Surface stress due to SSO drag (Pa) |
1327 |
|
|
|
1328 |
|
|
! d_t-----output-R: T increment |
1329 |
|
|
! d_u-----output-R: U increment |
1330 |
|
|
! d_v-----output-R: V increment |
1331 |
|
|
|
1332 |
|
|
! Implicit Arguments: |
1333 |
|
|
! =================== |
1334 |
|
|
|
1335 |
|
|
! iim--common-I: Number of longitude intervals |
1336 |
|
|
! jjm--common-I: Number of latitude intervals |
1337 |
|
|
! klon-common-I: Number of points seen by the physics |
1338 |
|
|
! (iim+1)*(jjm+1) for instance |
1339 |
|
|
! klev-common-I: Number of vertical layers |
1340 |
|
|
! ====================================================================== |
1341 |
|
|
! Local Variables: |
1342 |
|
|
! ================ |
1343 |
|
|
|
1344 |
|
|
! zgeom-----R: Altitude of layer above ground |
1345 |
|
|
! pt, pu, pv --R: t u v from top to bottom |
1346 |
|
|
! pdtdt, pdudt, pdvdt --R: t u v tendencies (from top to bottom) |
1347 |
|
|
! papmf: pressure at model layer (from top to bottom) |
1348 |
|
|
! papmh: pressure at model 1/2 layer (from top to bottom) |
1349 |
|
|
|
1350 |
|
|
! ====================================================================== |
1351 |
|
|
|
1352 |
|
|
include "YOMCST.h" |
1353 |
|
|
include "YOEGWD.h" |
1354 |
|
|
|
1355 |
|
|
! ARGUMENTS |
1356 |
|
|
|
1357 |
|
|
INTEGER nlon, nlev |
1358 |
|
|
REAL dtime |
1359 |
|
|
REAL paprs(klon, klev+1) |
1360 |
|
|
REAL pplay(klon, klev) |
1361 |
|
|
REAL plat(nlon), pmea(nlon) |
1362 |
|
|
REAL pstd(nlon), psig(nlon), pgam(nlon), pthe(nlon) |
1363 |
|
|
REAL ppic(nlon), pval(nlon) |
1364 |
|
|
REAL pulow(nlon), pvlow(nlon), pustr(nlon), pvstr(nlon) |
1365 |
|
|
REAL t(nlon, nlev), u(nlon, nlev), v(nlon, nlev) |
1366 |
|
|
REAL d_t(nlon, nlev), d_u(nlon, nlev), d_v(nlon, nlev) |
1367 |
|
|
|
1368 |
|
|
INTEGER i, k, kgwd, kdx(nlon), ktest(nlon) |
1369 |
|
|
|
1370 |
|
|
! Variables locales: |
1371 |
|
|
|
1372 |
|
576 |
REAL zgeom(klon, klev) |
1373 |
|
576 |
REAL pdtdt(klon, klev), pdudt(klon, klev), pdvdt(klon, klev) |
1374 |
|
576 |
REAL pt(klon, klev), pu(klon, klev), pv(klon, klev) |
1375 |
|
288 |
REAL papmf(klon, klev), papmh(klon, klev+1) |
1376 |
|
|
|
1377 |
|
|
! initialiser les variables de sortie (pour securite) |
1378 |
|
|
|
1379 |
|
|
|
1380 |
|
|
! print *,'in lift_noro' |
1381 |
✓✓ |
286560 |
DO i = 1, klon |
1382 |
|
286272 |
pulow(i) = 0.0 |
1383 |
|
286272 |
pvlow(i) = 0.0 |
1384 |
|
286272 |
pustr(i) = 0.0 |
1385 |
|
286560 |
pvstr(i) = 0.0 |
1386 |
|
|
END DO |
1387 |
✓✓ |
11520 |
DO k = 1, klev |
1388 |
✓✓ |
11176128 |
DO i = 1, klon |
1389 |
|
11164608 |
d_t(i, k) = 0.0 |
1390 |
|
11164608 |
d_u(i, k) = 0.0 |
1391 |
|
11164608 |
d_v(i, k) = 0.0 |
1392 |
|
11164608 |
pdudt(i, k) = 0.0 |
1393 |
|
11164608 |
pdvdt(i, k) = 0.0 |
1394 |
|
11175840 |
pdtdt(i, k) = 0.0 |
1395 |
|
|
END DO |
1396 |
|
|
END DO |
1397 |
|
|
|
1398 |
|
|
! preparer les variables d'entree (attention: l'ordre des niveaux |
1399 |
|
|
! verticaux augmente du haut vers le bas) |
1400 |
|
|
|
1401 |
✓✓ |
11520 |
DO k = 1, klev |
1402 |
✓✓ |
11176128 |
DO i = 1, klon |
1403 |
|
11164608 |
pt(i, k) = t(i, klev-k+1) |
1404 |
|
11164608 |
pu(i, k) = u(i, klev-k+1) |
1405 |
|
11164608 |
pv(i, k) = v(i, klev-k+1) |
1406 |
|
11175840 |
papmf(i, k) = pplay(i, klev-k+1) |
1407 |
|
|
END DO |
1408 |
|
|
END DO |
1409 |
✓✓ |
11808 |
DO k = 1, klev + 1 |
1410 |
✓✓ |
11462688 |
DO i = 1, klon |
1411 |
|
11462400 |
papmh(i, k) = paprs(i, klev-k+2) |
1412 |
|
|
END DO |
1413 |
|
|
END DO |
1414 |
✓✓ |
286560 |
DO i = 1, klon |
1415 |
|
286560 |
zgeom(i, klev) = rd*pt(i, klev)*log(papmh(i,klev+1)/papmf(i,klev)) |
1416 |
|
|
END DO |
1417 |
✓✓ |
11232 |
DO k = klev - 1, 1, -1 |
1418 |
✓✓ |
10889568 |
DO i = 1, klon |
1419 |
|
|
zgeom(i, k) = zgeom(i, k+1) + rd*(pt(i,k)+pt(i,k+1))/2.0*log(papmf(i,k+ & |
1420 |
|
10889280 |
1)/papmf(i,k)) |
1421 |
|
|
END DO |
1422 |
|
|
END DO |
1423 |
|
|
|
1424 |
|
|
! appeler la routine principale |
1425 |
|
|
|
1426 |
|
|
|
1427 |
|
|
CALL orolift_strato(klon, klev, kgwd, kdx, ktest, dtime, papmh, papmf, & |
1428 |
|
|
zgeom, pt, pu, pv, plat, pmea, pstd, psig, pgam, pthe, ppic, pval, pulow, & |
1429 |
|
288 |
pvlow, pdudt, pdvdt, pdtdt) |
1430 |
|
|
|
1431 |
✓✓ |
11520 |
DO k = 1, klev |
1432 |
✓✓ |
11176128 |
DO i = 1, klon |
1433 |
|
11164608 |
d_u(i, klev+1-k) = dtime*pdudt(i, k) |
1434 |
|
11164608 |
d_v(i, klev+1-k) = dtime*pdvdt(i, k) |
1435 |
|
11164608 |
d_t(i, klev+1-k) = dtime*pdtdt(i, k) |
1436 |
|
11164608 |
pustr(i) = pustr(i) + pdudt(i, k)*(papmh(i,k+1)-papmh(i,k))/rg |
1437 |
|
11175840 |
pvstr(i) = pvstr(i) + pdvdt(i, k)*(papmh(i,k+1)-papmh(i,k))/rg |
1438 |
|
|
END DO |
1439 |
|
|
END DO |
1440 |
|
|
|
1441 |
|
|
! print *,' out lift_noro' |
1442 |
|
|
|
1443 |
|
288 |
RETURN |
1444 |
|
|
END SUBROUTINE lift_noro_strato |
1445 |
|
6574890 |
SUBROUTINE orolift_strato(nlon, nlev, kgwd, kdx, ktest, ptsphy, paphm1, & |
1446 |
|
|
papm1, pgeom1, ptm1, pum1, pvm1, plat, pmea, pstd, psig, pgam, pthe, & |
1447 |
|
|
ppic, pval & ! OUTPUTS |
1448 |
|
288 |
, pulow, pvlow, pvom, pvol, pte) |
1449 |
|
|
|
1450 |
|
|
|
1451 |
|
|
! **** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT. |
1452 |
|
|
|
1453 |
|
|
! PURPOSE. |
1454 |
|
|
! -------- |
1455 |
|
|
! this routine computes the physical tendencies of the |
1456 |
|
|
! prognostic variables u,v when enhanced vortex stretching |
1457 |
|
|
! is needed. |
1458 |
|
|
|
1459 |
|
|
! ** INTERFACE. |
1460 |
|
|
! ---------- |
1461 |
|
|
! CALLED FROM *lift_noro |
1462 |
|
|
! explicit arguments : |
1463 |
|
|
! -------------------- |
1464 |
|
|
! ==== inputs === |
1465 |
|
|
! nlon----input-I-Total number of horizontal points that get into physics |
1466 |
|
|
! nlev----input-I-Number of vertical levels |
1467 |
|
|
|
1468 |
|
|
! kgwd- -input-I: Total nb of points where the orography schemes are active |
1469 |
|
|
! ktest--input-I: Flags to indicate active points |
1470 |
|
|
! kdx----input-I: Locate the physical location of an active point. |
1471 |
|
|
! ptsphy--input-R-Time-step (s) |
1472 |
|
|
! paphm1--input-R: pressure at model 1/2 layer |
1473 |
|
|
! papm1---input-R: pressure at model layer |
1474 |
|
|
! pgeom1--input-R: Altitude of layer above ground |
1475 |
|
|
! ptm1, pum1, pvm1--R-: t, u and v |
1476 |
|
|
! pmea----input-R-Mean Orography (m) |
1477 |
|
|
! pstd----input-R-SSO standard deviation (m) |
1478 |
|
|
! psig----input-R-SSO slope |
1479 |
|
|
! pgam----input-R-SSO Anisotropy |
1480 |
|
|
! pthe----input-R-SSO Angle |
1481 |
|
|
! ppic----input-R-SSO Peacks elevation (m) |
1482 |
|
|
! pval----input-R-SSO Valleys elevation (m) |
1483 |
|
|
! plat----input-R-Latitude (degree) |
1484 |
|
|
|
1485 |
|
|
! ==== outputs === |
1486 |
|
|
! pulow, pvlow -output-R: Low-level wind |
1487 |
|
|
|
1488 |
|
|
! pte -----output-R: T tendency |
1489 |
|
|
! pvom-----output-R: U tendency |
1490 |
|
|
! pvol-----output-R: V tendency |
1491 |
|
|
|
1492 |
|
|
|
1493 |
|
|
! Implicit Arguments: |
1494 |
|
|
! =================== |
1495 |
|
|
|
1496 |
|
|
! klon-common-I: Number of points seen by the physics |
1497 |
|
|
! klev-common-I: Number of vertical layers |
1498 |
|
|
|
1499 |
|
|
|
1500 |
|
|
! ---------- |
1501 |
|
|
|
1502 |
|
|
! AUTHOR. |
1503 |
|
|
! ------- |
1504 |
|
|
! F.LOTT LMD 22/11/95 |
1505 |
|
|
|
1506 |
|
|
USE dimphy |
1507 |
|
|
IMPLICIT NONE |
1508 |
|
|
|
1509 |
|
|
|
1510 |
|
|
include "YOMCST.h" |
1511 |
|
|
include "YOEGWD.h" |
1512 |
|
|
! ----------------------------------------------------------------------- |
1513 |
|
|
|
1514 |
|
|
! * 0.1 ARGUMENTS |
1515 |
|
|
! --------- |
1516 |
|
|
|
1517 |
|
|
|
1518 |
|
|
INTEGER nlon, nlev, kgwd |
1519 |
|
|
REAL ptsphy |
1520 |
|
|
REAL pte(nlon, nlev), pvol(nlon, nlev), pvom(nlon, nlev), pulow(nlon), & |
1521 |
|
|
pvlow(nlon) |
1522 |
|
|
REAL pum1(nlon, nlev), pvm1(nlon, nlev), ptm1(nlon, nlev), plat(nlon), & |
1523 |
|
|
pmea(nlon), pstd(nlon), psig(nlon), pgam(nlon), pthe(nlon), ppic(nlon), & |
1524 |
|
|
pval(nlon), pgeom1(nlon, nlev), papm1(nlon, nlev), paphm1(nlon, nlev+1) |
1525 |
|
|
|
1526 |
|
|
INTEGER kdx(nlon), ktest(nlon) |
1527 |
|
|
! ----------------------------------------------------------------------- |
1528 |
|
|
|
1529 |
|
|
! * 0.2 local arrays |
1530 |
|
|
|
1531 |
|
|
INTEGER jl, ilevh, jk |
1532 |
|
|
REAL zhgeo, zdelp, zslow, zsqua, zscav, zbet |
1533 |
|
|
! ------------ |
1534 |
|
576 |
INTEGER iknub(klon), iknul(klon) |
1535 |
|
576 |
LOGICAL ll1(klon, klev+1) |
1536 |
|
|
|
1537 |
|
576 |
REAL ztau(klon, klev+1), ztav(klon, klev+1), zrho(klon, klev+1) |
1538 |
|
576 |
REAL zdudt(klon), zdvdt(klon) |
1539 |
|
576 |
REAL zhcrit(klon, klev) |
1540 |
|
|
|
1541 |
|
|
LOGICAL lifthigh |
1542 |
|
|
REAL zcons1, ztmst |
1543 |
|
|
CHARACTER (LEN=20) :: modname = 'orolift_strato' |
1544 |
|
|
CHARACTER (LEN=80) :: abort_message |
1545 |
|
|
|
1546 |
|
|
|
1547 |
|
|
! ----------------------------------------------------------------------- |
1548 |
|
|
|
1549 |
|
|
! * 1.1 initialisations |
1550 |
|
|
! --------------- |
1551 |
|
|
|
1552 |
|
|
lifthigh = .FALSE. |
1553 |
|
|
|
1554 |
✓✗✗✓
|
288 |
IF (nlon/=klon .OR. nlev/=klev) THEN |
1555 |
|
|
abort_message = 'pb dimension' |
1556 |
|
|
CALL abort_physic(modname, abort_message, 1) |
1557 |
|
|
END IF |
1558 |
|
288 |
zcons1 = 1./rd |
1559 |
|
288 |
ztmst = ptsphy |
1560 |
|
|
|
1561 |
✓✓ |
286560 |
DO jl = kidia, kfdia |
1562 |
|
286272 |
zrho(jl, klev+1) = 0.0 |
1563 |
|
286272 |
pulow(jl) = 0.0 |
1564 |
|
286272 |
pvlow(jl) = 0.0 |
1565 |
|
286272 |
iknub(jl) = klev |
1566 |
|
286272 |
iknul(jl) = klev |
1567 |
|
|
ilevh = klev/3 |
1568 |
|
286272 |
ll1(jl, klev+1) = .FALSE. |
1569 |
✓✓ |
11451168 |
DO jk = 1, klev |
1570 |
|
11164608 |
pvom(jl, jk) = 0.0 |
1571 |
|
11164608 |
pvol(jl, jk) = 0.0 |
1572 |
|
11450880 |
pte(jl, jk) = 0.0 |
1573 |
|
|
END DO |
1574 |
|
|
END DO |
1575 |
|
|
|
1576 |
|
|
|
1577 |
|
|
! * 2.1 DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF |
1578 |
|
|
! * LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE |
1579 |
|
|
! * THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS. |
1580 |
|
|
|
1581 |
|
|
|
1582 |
|
|
|
1583 |
✓✓ |
11520 |
DO jk = klev, 1, -1 |
1584 |
✓✓ |
11176128 |
DO jl = kidia, kfdia |
1585 |
✓✓ |
11175840 |
IF (ktest(jl)==1) THEN |
1586 |
|
5211648 |
zhcrit(jl, jk) = amax1(ppic(jl)-pval(jl), 100.) |
1587 |
|
5211648 |
zhgeo = pgeom1(jl, jk)/rg |
1588 |
|
5211648 |
ll1(jl, jk) = (zhgeo>zhcrit(jl,jk)) |
1589 |
✓✓ |
5211648 |
IF (ll1(jl,jk) .NEQV. ll1(jl,jk+1)) THEN |
1590 |
|
133632 |
iknub(jl) = jk |
1591 |
|
|
END IF |
1592 |
|
|
END IF |
1593 |
|
|
END DO |
1594 |
|
|
END DO |
1595 |
|
|
|
1596 |
|
|
|
1597 |
✓✓ |
286560 |
DO jl = kidia, kfdia |
1598 |
✓✓ |
286560 |
IF (ktest(jl)==1) THEN |
1599 |
|
133632 |
iknub(jl) = max(iknub(jl), klev/2) |
1600 |
|
133632 |
iknul(jl) = max(iknul(jl), 2*klev/3) |
1601 |
✗✓ |
133632 |
IF (iknub(jl)>nktopg) iknub(jl) = nktopg |
1602 |
✓✓ |
133632 |
IF (iknub(jl)==nktopg) iknul(jl) = klev |
1603 |
✗✓ |
133632 |
IF (iknub(jl)==iknul(jl)) iknub(jl) = iknul(jl) - 1 |
1604 |
|
|
END IF |
1605 |
|
|
END DO |
1606 |
|
|
|
1607 |
✓✓ |
11232 |
DO jk = klev, 2, -1 |
1608 |
✓✓ |
10889568 |
DO jl = kidia, kfdia |
1609 |
|
10889280 |
zrho(jl, jk) = 2.*paphm1(jl, jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1)) |
1610 |
|
|
END DO |
1611 |
|
|
END DO |
1612 |
|
|
! print *,' dans orolift: 223' |
1613 |
|
|
|
1614 |
|
|
! ******************************************************************** |
1615 |
|
|
|
1616 |
|
|
! * define low level flow |
1617 |
|
|
! ------------------- |
1618 |
✓✓ |
11520 |
DO jk = klev, 1, -1 |
1619 |
✓✓ |
11176128 |
DO jl = kidia, kfdia |
1620 |
✓✓ |
11175840 |
IF (ktest(jl)==1) THEN |
1621 |
✓✓✓✗
|
5211648 |
IF (jk>=iknub(jl) .AND. jk<=iknul(jl)) THEN |
1622 |
|
|
pulow(jl) = pulow(jl) + pum1(jl, jk)*(paphm1(jl,jk+1)-paphm1(jl,jk) & |
1623 |
|
1095402 |
) |
1624 |
|
|
pvlow(jl) = pvlow(jl) + pvm1(jl, jk)*(paphm1(jl,jk+1)-paphm1(jl,jk) & |
1625 |
|
1095402 |
) |
1626 |
|
|
zrho(jl, klev+1) = zrho(jl, klev+1) + zrho(jl, jk)*(paphm1(jl,jk+1) & |
1627 |
|
1095402 |
-paphm1(jl,jk)) |
1628 |
|
|
END IF |
1629 |
|
|
END IF |
1630 |
|
|
END DO |
1631 |
|
|
END DO |
1632 |
✓✓ |
286560 |
DO jl = kidia, kfdia |
1633 |
✓✓ |
286560 |
IF (ktest(jl)==1) THEN |
1634 |
|
133632 |
pulow(jl) = pulow(jl)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl))) |
1635 |
|
133632 |
pvlow(jl) = pvlow(jl)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl))) |
1636 |
|
|
zrho(jl, klev+1) = zrho(jl, klev+1)/(paphm1(jl,iknul(jl)+1)-paphm1(jl, & |
1637 |
|
133632 |
iknub(jl))) |
1638 |
|
|
END IF |
1639 |
|
|
END DO |
1640 |
|
|
|
1641 |
|
|
! *********************************************************** |
1642 |
|
|
|
1643 |
|
|
! * 3. COMPUTE MOUNTAIN LIFT |
1644 |
|
|
|
1645 |
|
|
|
1646 |
✓✓ |
286560 |
DO jl = kidia, kfdia |
1647 |
✓✓ |
286560 |
IF (ktest(jl)==1) THEN |
1648 |
|
|
ztau(jl, klev+1) = -gklift*zrho(jl, klev+1)*2.*romega* & ! * |
1649 |
|
|
! (2*pstd(jl)+pmea(jl))* |
1650 |
|
133632 |
2*pstd(jl)*sin(rpi/180.*plat(jl))*pvlow(jl) |
1651 |
|
|
ztav(jl, klev+1) = gklift*zrho(jl, klev+1)*2.*romega* & ! * |
1652 |
|
|
! (2*pstd(jl)+pmea(jl))* |
1653 |
|
133632 |
2*pstd(jl)*sin(rpi/180.*plat(jl))*pulow(jl) |
1654 |
|
|
ELSE |
1655 |
|
152640 |
ztau(jl, klev+1) = 0.0 |
1656 |
|
152640 |
ztav(jl, klev+1) = 0.0 |
1657 |
|
|
END IF |
1658 |
|
|
END DO |
1659 |
|
|
|
1660 |
|
|
! * 4. COMPUTE LIFT PROFILE |
1661 |
|
|
! * -------------------- |
1662 |
|
|
|
1663 |
|
|
|
1664 |
|
|
|
1665 |
✓✓ |
11520 |
DO jk = 1, klev |
1666 |
✓✓ |
11176128 |
DO jl = kidia, kfdia |
1667 |
✓✓ |
11175840 |
IF (ktest(jl)==1) THEN |
1668 |
|
5211648 |
ztau(jl, jk) = ztau(jl, klev+1)*paphm1(jl, jk)/paphm1(jl, klev+1) |
1669 |
|
5211648 |
ztav(jl, jk) = ztav(jl, klev+1)*paphm1(jl, jk)/paphm1(jl, klev+1) |
1670 |
|
|
ELSE |
1671 |
|
5952960 |
ztau(jl, jk) = 0.0 |
1672 |
|
5952960 |
ztav(jl, jk) = 0.0 |
1673 |
|
|
END IF |
1674 |
|
|
END DO |
1675 |
|
|
END DO |
1676 |
|
|
|
1677 |
|
|
|
1678 |
|
|
! * 5. COMPUTE TENDENCIES. |
1679 |
|
|
! * ------------------- |
1680 |
|
|
IF (lifthigh) THEN |
1681 |
|
|
! EXPLICIT SOLUTION AT ALL LEVELS |
1682 |
|
|
|
1683 |
|
|
DO jk = 1, klev |
1684 |
|
|
DO jl = kidia, kfdia |
1685 |
|
|
IF (ktest(jl)==1) THEN |
1686 |
|
|
zdelp = paphm1(jl, jk+1) - paphm1(jl, jk) |
1687 |
|
|
zdudt(jl) = -rg*(ztau(jl,jk+1)-ztau(jl,jk))/zdelp |
1688 |
|
|
zdvdt(jl) = -rg*(ztav(jl,jk+1)-ztav(jl,jk))/zdelp |
1689 |
|
|
END IF |
1690 |
|
|
END DO |
1691 |
|
|
END DO |
1692 |
|
|
|
1693 |
|
|
! PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY |
1694 |
|
|
|
1695 |
|
|
DO jk = 1, klev |
1696 |
|
|
DO jl = kidia, kfdia |
1697 |
|
|
IF (ktest(jl)==1) THEN |
1698 |
|
|
|
1699 |
|
|
zslow = sqrt(pulow(jl)**2+pvlow(jl)**2) |
1700 |
|
|
zsqua = amax1(sqrt(pum1(jl,jk)**2+pvm1(jl,jk)**2), gvsec) |
1701 |
|
|
zscav = -zdudt(jl)*pvm1(jl, jk) + zdvdt(jl)*pum1(jl, jk) |
1702 |
|
|
IF (zsqua>gvsec) THEN |
1703 |
|
|
pvom(jl, jk) = -zscav*pvm1(jl, jk)/zsqua**2 |
1704 |
|
|
pvol(jl, jk) = zscav*pum1(jl, jk)/zsqua**2 |
1705 |
|
|
ELSE |
1706 |
|
|
pvom(jl, jk) = 0.0 |
1707 |
|
|
pvol(jl, jk) = 0.0 |
1708 |
|
|
END IF |
1709 |
|
|
zsqua = sqrt(pum1(jl,jk)**2+pum1(jl,jk)**2) |
1710 |
|
|
IF (zsqua<zslow) THEN |
1711 |
|
|
pvom(jl, jk) = zsqua/zslow*pvom(jl, jk) |
1712 |
|
|
pvol(jl, jk) = zsqua/zslow*pvol(jl, jk) |
1713 |
|
|
END IF |
1714 |
|
|
|
1715 |
|
|
END IF |
1716 |
|
|
END DO |
1717 |
|
|
END DO |
1718 |
|
|
|
1719 |
|
|
! 6. LOW LEVEL LIFT, SEMI IMPLICIT: |
1720 |
|
|
! ---------------------------------- |
1721 |
|
|
|
1722 |
|
|
ELSE |
1723 |
|
|
|
1724 |
✓✓ |
286560 |
DO jl = kidia, kfdia |
1725 |
✓✓ |
286560 |
IF (ktest(jl)==1) THEN |
1726 |
✓✓ |
1229034 |
DO jk = klev, iknub(jl), -1 |
1727 |
|
|
zbet = gklift*2.*romega*sin(rpi/180.*plat(jl))*ztmst* & |
1728 |
|
|
(pgeom1(jl,iknub(jl)-1)-pgeom1(jl,jk))/ & |
1729 |
|
1095402 |
(pgeom1(jl,iknub(jl)-1)-pgeom1(jl,klev)) |
1730 |
|
1095402 |
zdudt(jl) = -pum1(jl, jk)/ztmst/(1+zbet**2) |
1731 |
|
1095402 |
zdvdt(jl) = -pvm1(jl, jk)/ztmst/(1+zbet**2) |
1732 |
|
1095402 |
pvom(jl, jk) = zbet**2*zdudt(jl) - zbet*zdvdt(jl) |
1733 |
|
1229034 |
pvol(jl, jk) = zbet*zdudt(jl) + zbet**2*zdvdt(jl) |
1734 |
|
|
END DO |
1735 |
|
|
END IF |
1736 |
|
|
END DO |
1737 |
|
|
|
1738 |
|
|
END IF |
1739 |
|
|
|
1740 |
|
|
! print *,' out orolift' |
1741 |
|
|
|
1742 |
|
288 |
RETURN |
1743 |
|
|
END SUBROUTINE orolift_strato |
1744 |
|
7 |
SUBROUTINE sugwd_strato(nlon, nlev, paprs, pplay) |
1745 |
|
|
|
1746 |
|
|
|
1747 |
|
|
! **** *SUGWD* INITIALIZE COMMON YOEGWD CONTROLLING GRAVITY WAVE DRAG |
1748 |
|
|
|
1749 |
|
|
! PURPOSE. |
1750 |
|
|
! -------- |
1751 |
|
|
! INITIALIZE YOEGWD, THE COMMON THAT CONTROLS THE |
1752 |
|
|
! GRAVITY WAVE DRAG PARAMETRIZATION. |
1753 |
|
|
! VERY IMPORTANT: |
1754 |
|
|
! ______________ |
1755 |
|
|
! THIS ROUTINE SET_UP THE "TUNABLE PARAMETERS" OF THE |
1756 |
|
|
! VARIOUS SSO SCHEMES |
1757 |
|
|
|
1758 |
|
|
! ** INTERFACE. |
1759 |
|
|
! ---------- |
1760 |
|
|
! CALL *SUGWD* FROM *SUPHEC* |
1761 |
|
|
! ----- ------ |
1762 |
|
|
|
1763 |
|
|
! EXPLICIT ARGUMENTS : |
1764 |
|
|
! -------------------- |
1765 |
|
|
! PAPRS,PPLAY : Pressure at semi and full model levels |
1766 |
|
|
! NLEV : number of model levels |
1767 |
|
|
! NLON : number of points treated in the physics |
1768 |
|
|
|
1769 |
|
|
! IMPLICIT ARGUMENTS : |
1770 |
|
|
! -------------------- |
1771 |
|
|
! COMMON YOEGWD |
1772 |
|
|
! -GFRCRIT-R: Critical Non-dimensional mountain Height |
1773 |
|
|
! (HNC in (1), LOTT 1999) |
1774 |
|
|
! -GKWAKE--R: Bluff-body drag coefficient for low level wake |
1775 |
|
|
! (Cd in (2), LOTT 1999) |
1776 |
|
|
! -GRCRIT--R: Critical Richardson Number |
1777 |
|
|
! (Ric, End of first column p791 of LOTT 1999) |
1778 |
|
|
! -GKDRAG--R: Gravity wave drag coefficient |
1779 |
|
|
! (G in (3), LOTT 1999) |
1780 |
|
|
! -GKLIFT--R: Mountain Lift coefficient |
1781 |
|
|
! (Cl in (4), LOTT 1999) |
1782 |
|
|
! -GHMAX---R: Not used |
1783 |
|
|
! -GRAHILO-R: Set-up the trapped waves fraction |
1784 |
|
|
! (Beta , End of first column, LOTT 1999) |
1785 |
|
|
|
1786 |
|
|
! -GSIGCR--R: Security value for blocked flow depth |
1787 |
|
|
! -NKTOPG--I: Security value for blocked flow level |
1788 |
|
|
! -nstra----I: An estimate to qualify the upper levels of |
1789 |
|
|
! the model where one wants to impose strees |
1790 |
|
|
! profiles |
1791 |
|
|
! -GSSECC--R: Security min value for low-level B-V frequency |
1792 |
|
|
! -GTSEC---R: Security min value for anisotropy and GW stress. |
1793 |
|
|
! -GVSEC---R: Security min value for ulow |
1794 |
|
|
|
1795 |
|
|
|
1796 |
|
|
! METHOD. |
1797 |
|
|
! ------- |
1798 |
|
|
! SEE DOCUMENTATION |
1799 |
|
|
|
1800 |
|
|
! EXTERNALS. |
1801 |
|
|
! ---------- |
1802 |
|
|
! NONE |
1803 |
|
|
|
1804 |
|
|
! REFERENCE. |
1805 |
|
|
! ---------- |
1806 |
|
|
! Lott, 1999: Alleviation of stationary biases in a GCM through... |
1807 |
|
|
! Monthly Weather Review, 127, pp 788-801. |
1808 |
|
|
|
1809 |
|
|
! AUTHOR. |
1810 |
|
|
! ------- |
1811 |
|
|
! FRANCOIS LOTT *LMD* |
1812 |
|
|
|
1813 |
|
|
! MODIFICATIONS. |
1814 |
|
|
! -------------- |
1815 |
|
|
! ORIGINAL : 90-01-01 (MARTIN MILLER, ECMWF) |
1816 |
|
|
! LAST: 99-07-09 (FRANCOIS LOTT,LMD) |
1817 |
|
|
! ------------------------------------------------------------------ |
1818 |
|
|
USE dimphy |
1819 |
|
|
USE mod_phys_lmdz_para |
1820 |
|
|
USE mod_grid_phy_lmdz |
1821 |
|
|
USE geometry_mod |
1822 |
|
|
IMPLICIT NONE |
1823 |
|
|
|
1824 |
|
|
! ----------------------------------------------------------------- |
1825 |
|
|
include "YOEGWD.h" |
1826 |
|
|
! ---------------------------------------------------------------- |
1827 |
|
|
|
1828 |
|
|
! ARGUMENTS |
1829 |
|
|
INTEGER nlon, nlev |
1830 |
|
|
REAL paprs(nlon, nlev+1) |
1831 |
|
|
REAL pplay(nlon, nlev) |
1832 |
|
|
|
1833 |
|
|
INTEGER jk |
1834 |
|
|
REAL zpr, ztop, zsigt, zpm1r |
1835 |
|
|
INTEGER :: cell,ij,nstra_tmp,nktopg_tmp |
1836 |
|
|
REAL :: current_dist, dist_min,dist_min_glo |
1837 |
|
|
|
1838 |
|
|
! * 1. SET THE VALUES OF THE PARAMETERS |
1839 |
|
|
! -------------------------------- |
1840 |
|
|
|
1841 |
|
|
|
1842 |
|
1 |
PRINT *, ' DANS SUGWD NLEV=', nlev |
1843 |
|
1 |
ghmax = 10000. |
1844 |
|
|
|
1845 |
|
|
zpr = 100000. |
1846 |
|
|
ZTOP=0.00005 |
1847 |
|
|
zsigt = 0.94 |
1848 |
|
|
! old ZPR=80000. |
1849 |
|
|
! old ZSIGT=0.85 |
1850 |
|
|
|
1851 |
|
|
|
1852 |
|
|
!ym Take the point at equator close to (0,0) coordinates. |
1853 |
|
1 |
dist_min=360 |
1854 |
|
1 |
dist_min_glo=360. |
1855 |
|
1 |
cell=-1 |
1856 |
✓✓ |
995 |
DO ij=1,klon |
1857 |
|
994 |
current_dist=sqrt(longitude_deg(ij)**2+latitude_deg(ij)**2) |
1858 |
|
994 |
current_dist=current_dist*(1+(1e-10*ind_cell_glo(ij))/klon_glo) ! For point unicity |
1859 |
✓✓ |
995 |
IF (dist_min>current_dist) THEN |
1860 |
|
39 |
dist_min=current_dist |
1861 |
|
39 |
cell=ij |
1862 |
|
|
ENDIF |
1863 |
|
|
ENDDO |
1864 |
|
|
|
1865 |
|
|
!PRINT *, 'SUGWD distmin cell=', dist_min,cell |
1866 |
|
1 |
CALL reduce_min(dist_min,dist_min_glo) |
1867 |
|
1 |
CALL bcast(dist_min_glo) |
1868 |
✗✓ |
1 |
IF (dist_min/=dist_min_glo) cell=-1 |
1869 |
|
|
!ym in future find the point at equator close to (0,0) coordinates. |
1870 |
|
1 |
PRINT *, 'SUGWD distmin dist_min_glo cell=', dist_min,dist_min_glo,cell |
1871 |
|
|
|
1872 |
|
1 |
nktopg_tmp=nktopg |
1873 |
|
1 |
nstra_tmp=nstra |
1874 |
|
|
|
1875 |
✓✗ |
1 |
IF (cell/=-1) THEN |
1876 |
|
|
|
1877 |
|
|
!print*,'SUGWD shape ',shape(pplay),cell+1 |
1878 |
|
|
|
1879 |
✓✓ |
40 |
DO jk = 1, nlev |
1880 |
|
|
!zpm1r = pplay(cell+1, jk)/paprs(cell+1, 1) |
1881 |
|
39 |
zpm1r = pplay(cell, jk)/paprs(cell, 1) |
1882 |
✓✓ |
39 |
IF (zpm1r>=zsigt) THEN |
1883 |
|
4 |
nktopg_tmp = jk |
1884 |
|
|
END IF |
1885 |
✓✓ |
40 |
IF (zpm1r>=ztop) THEN |
1886 |
|
38 |
nstra_tmp = jk |
1887 |
|
|
END IF |
1888 |
|
|
END DO |
1889 |
|
|
ELSE |
1890 |
|
|
nktopg_tmp=0 |
1891 |
|
|
nstra_tmp=0 |
1892 |
|
|
ENDIF |
1893 |
|
|
|
1894 |
|
1 |
CALL reduce_sum(nktopg_tmp,nktopg) |
1895 |
|
1 |
CALL bcast(nktopg) |
1896 |
|
1 |
CALL reduce_sum(nstra_tmp,nstra) |
1897 |
|
1 |
CALL bcast(nstra) |
1898 |
|
|
|
1899 |
|
|
! inversion car dans orodrag on compte les niveaux a l'envers |
1900 |
|
1 |
nktopg = nlev - nktopg + 1 |
1901 |
|
1 |
nstra = nlev - nstra |
1902 |
|
1 |
PRINT *, ' DANS SUGWD nktopg=', nktopg |
1903 |
|
1 |
PRINT *, ' DANS SUGWD nstra=', nstra |
1904 |
✗✓ |
1 |
if (nstra == 0) call abort_physic("sugwd_strato", "no level in stratosphere", 1) |
1905 |
|
|
|
1906 |
|
|
! Valeurs lues dans les .def, ou attribues dans conf_phys |
1907 |
|
|
!gkdrag = 0.2 |
1908 |
|
|
!grahilo = 0.1 |
1909 |
|
|
!grcrit = 1.00 |
1910 |
|
|
!gfrcrit = 0.70 |
1911 |
|
|
!gkwake = 0.40 |
1912 |
|
|
!gklift = 0.25 |
1913 |
|
|
|
1914 |
|
1 |
gsigcr = 0.80 ! Top of low level flow |
1915 |
|
1 |
gvcrit = 0.1 |
1916 |
|
|
|
1917 |
|
1 |
WRITE (UNIT=6, FMT='('' *** SSO essential constants ***'')') |
1918 |
|
1 |
WRITE (UNIT=6, FMT='('' *** SPECIFIED IN SUGWD ***'')') |
1919 |
|
1 |
WRITE (UNIT=6, FMT='('' Gravity wave ct '',E13.7,'' '')') gkdrag |
1920 |
|
1 |
WRITE (UNIT=6, FMT='('' Trapped/total wave dag '',E13.7,'' '')') grahilo |
1921 |
|
1 |
WRITE (UNIT=6, FMT='('' Critical Richardson = '',E13.7,'' '')') grcrit |
1922 |
|
1 |
WRITE (UNIT=6, FMT='('' Critical Froude'',e13.7)') gfrcrit |
1923 |
|
1 |
WRITE (UNIT=6, FMT='('' Low level Wake bluff cte'',e13.7)') gkwake |
1924 |
|
1 |
WRITE (UNIT=6, FMT='('' Low level lift cte'',e13.7)') gklift |
1925 |
|
|
|
1926 |
|
|
! ---------------------------------------------------------------- |
1927 |
|
|
|
1928 |
|
|
! * 2. SET VALUES OF SECURITY PARAMETERS |
1929 |
|
|
! --------------------------------- |
1930 |
|
|
|
1931 |
|
|
|
1932 |
|
1 |
gvsec = 0.10 |
1933 |
|
1 |
gssec = 0.0001 |
1934 |
|
|
|
1935 |
|
1 |
gtsec = 0.00001 |
1936 |
|
|
|
1937 |
|
1 |
RETURN |
1938 |
|
|
END SUBROUTINE sugwd_strato |