1 |
! |
2 |
! $Header: /home/cvsroot/LMDZ4/libf/phylmd/orografi.F,v 1.4 2005/12/01 11:27:29 fairhead Exp $ |
3 |
! |
4 |
SUBROUTINE drag_noro (nlon,nlev,dtime,paprs,pplay, |
5 |
e pmea,pstd, psig, pgam, pthe,ppic,pval, |
6 |
e kgwd,kdx,ktest, |
7 |
e t, u, v, |
8 |
s pulow, pvlow, pustr, pvstr, |
9 |
s d_t, d_u, d_v) |
10 |
c |
11 |
use dimens_m |
12 |
use dimphy |
13 |
use YOMCST |
14 |
IMPLICIT none |
15 |
c====================================================================== |
16 |
c Auteur(s): F.Lott (LMD/CNRS) date: 19950201 |
17 |
c Objet: Frottement de la montagne Interface |
18 |
c====================================================================== |
19 |
c Arguments: |
20 |
c dtime---input-R- pas d'integration (s) |
21 |
c paprs---input-R-pression pour chaque inter-couche (en Pa) |
22 |
c pplay---input-R-pression pour le mileu de chaque couche (en Pa) |
23 |
c t-------input-R-temperature (K) |
24 |
c u-------input-R-vitesse horizontale (m/s) |
25 |
c v-------input-R-vitesse horizontale (m/s) |
26 |
c |
27 |
c d_t-----output-R-increment de la temperature |
28 |
c d_u-----output-R-increment de la vitesse u |
29 |
c d_v-----output-R-increment de la vitesse v |
30 |
c====================================================================== |
31 |
c |
32 |
c ARGUMENTS |
33 |
c |
34 |
INTEGER nlon,nlev |
35 |
REAL, intent(in):: dtime |
36 |
REAL, intent(in):: paprs(klon,klev+1) |
37 |
REAL, intent(in):: pplay(klon,klev) |
38 |
REAL pmea(nlon) |
39 |
real, intent(in):: pstd(nlon),psig(nlon) |
40 |
real pgam(nlon),pthe(nlon) |
41 |
REAL ppic(nlon),pval(nlon) |
42 |
REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon) |
43 |
REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev) |
44 |
REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev) |
45 |
c |
46 |
INTEGER i, k, kgwd, kdx(nlon), ktest(nlon) |
47 |
c |
48 |
c Variables locales: |
49 |
c |
50 |
REAL zgeom(klon,klev) |
51 |
REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev) |
52 |
REAL pt(klon,klev), pu(klon,klev), pv(klon,klev) |
53 |
REAL papmf(klon,klev),papmh(klon,klev+1) |
54 |
c |
55 |
c initialiser les variables de sortie (pour securite) |
56 |
c |
57 |
DO i = 1,klon |
58 |
pulow(i) = 0.0 |
59 |
pvlow(i) = 0.0 |
60 |
pustr(i) = 0.0 |
61 |
pvstr(i) = 0.0 |
62 |
ENDDO |
63 |
DO k = 1, klev |
64 |
DO i = 1, klon |
65 |
d_t(i,k) = 0.0 |
66 |
d_u(i,k) = 0.0 |
67 |
d_v(i,k) = 0.0 |
68 |
pdudt(i,k)=0.0 |
69 |
pdvdt(i,k)=0.0 |
70 |
pdtdt(i,k)=0.0 |
71 |
ENDDO |
72 |
ENDDO |
73 |
c |
74 |
c preparer les variables d'entree (attention: l'ordre des niveaux |
75 |
c verticaux augmente du haut vers le bas) |
76 |
c |
77 |
DO k = 1, klev |
78 |
DO i = 1, klon |
79 |
pt(i,k) = t(i,klev-k+1) |
80 |
pu(i,k) = u(i,klev-k+1) |
81 |
pv(i,k) = v(i,klev-k+1) |
82 |
papmf(i,k) = pplay(i,klev-k+1) |
83 |
ENDDO |
84 |
ENDDO |
85 |
DO k = 1, klev+1 |
86 |
DO i = 1, klon |
87 |
papmh(i,k) = paprs(i,klev-k+2) |
88 |
ENDDO |
89 |
ENDDO |
90 |
DO i = 1, klon |
91 |
zgeom(i,klev) = RD * pt(i,klev) |
92 |
. * LOG(papmh(i,klev+1)/papmf(i,klev)) |
93 |
ENDDO |
94 |
DO k = klev-1, 1, -1 |
95 |
DO i = 1, klon |
96 |
zgeom(i,k) = zgeom(i,k+1) + RD * (pt(i,k)+pt(i,k+1))/2.0 |
97 |
. * LOG(papmf(i,k+1)/papmf(i,k)) |
98 |
ENDDO |
99 |
ENDDO |
100 |
c |
101 |
c appeler la routine principale |
102 |
c |
103 |
CALL orodrag(klon,klev,kgwd,kdx,ktest, |
104 |
. dtime, |
105 |
. papmh, papmf, zgeom, |
106 |
. pt, pu, pv, |
107 |
. pmea, pstd, psig, pgam, pthe, ppic,pval, |
108 |
. pulow,pvlow, |
109 |
. pdudt,pdvdt,pdtdt) |
110 |
C |
111 |
DO k = 1, klev |
112 |
DO i = 1, klon |
113 |
d_u(i,klev+1-k) = dtime*pdudt(i,k) |
114 |
d_v(i,klev+1-k) = dtime*pdvdt(i,k) |
115 |
d_t(i,klev+1-k) = dtime*pdtdt(i,k) |
116 |
pustr(i) = pustr(i) |
117 |
cIM BUG . +rg*pdudt(i,k)*(papmh(i,k+1)-papmh(i,k)) |
118 |
. +pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))/RG |
119 |
pvstr(i) = pvstr(i) |
120 |
cIM BUG . +rg*pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k)) |
121 |
. +pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))/RG |
122 |
ENDDO |
123 |
ENDDO |
124 |
c |
125 |
RETURN |
126 |
END |
127 |
SUBROUTINE orodrag( nlon,nlev |
128 |
i , kgwd, kdx, ktest |
129 |
r , ptsphy |
130 |
r , paphm1,papm1,pgeom1,ptm1,pum1,pvm1 |
131 |
r , pmea, pstd, psig, pgamma, ptheta, ppic, pval |
132 |
c outputs |
133 |
r , pulow,pvlow |
134 |
r , pvom,pvol,pte ) |
135 |
|
136 |
use dimens_m |
137 |
use dimphy |
138 |
use YOMCST |
139 |
use yoegwd |
140 |
implicit none |
141 |
|
142 |
c |
143 |
c |
144 |
c**** *gwdrag* - does the gravity wave parametrization. |
145 |
c |
146 |
c purpose. |
147 |
c -------- |
148 |
c |
149 |
c this routine computes the physical tendencies of the |
150 |
c prognostic variables u,v and t due to vertical transports by |
151 |
c subgridscale orographically excited gravity waves |
152 |
c |
153 |
c** interface. |
154 |
c ---------- |
155 |
c called from *callpar*. |
156 |
c |
157 |
c the routine takes its input from the long-term storage: |
158 |
c u,v,t and p at t-1. |
159 |
c |
160 |
c explicit arguments : |
161 |
c -------------------- |
162 |
c ==== inputs === |
163 |
c ==== outputs === |
164 |
c |
165 |
c implicit arguments : none |
166 |
c -------------------- |
167 |
c |
168 |
c implicit logical (l) |
169 |
c |
170 |
c method. |
171 |
c ------- |
172 |
c |
173 |
c externals. |
174 |
c ---------- |
175 |
integer ismin, ismax |
176 |
external ismin, ismax |
177 |
c |
178 |
c reference. |
179 |
c ---------- |
180 |
c |
181 |
c author. |
182 |
c ------- |
183 |
c m.miller + b.ritter e.c.m.w.f. 15/06/86. |
184 |
c |
185 |
c f.lott + m. miller e.c.m.w.f. 22/11/94 |
186 |
c----------------------------------------------------------------------- |
187 |
c |
188 |
c |
189 |
c----------------------------------------------------------------------- |
190 |
c |
191 |
c* 0.1 arguments |
192 |
c --------- |
193 |
c |
194 |
c |
195 |
integer nlon, nlev, klevm1 |
196 |
integer kgwd, jl, ilevp1, jk, ji |
197 |
real zdelp, ztemp, zforc, ztend |
198 |
real rover, zb, zc, zconb, zabsv |
199 |
real zzd1, ratio, zbet, zust,zvst, zdis |
200 |
real pte(nlon,nlev), |
201 |
* pvol(nlon,nlev), |
202 |
* pvom(nlon,nlev), |
203 |
* pulow(klon), |
204 |
* pvlow(klon) |
205 |
real pum1(nlon,nlev), |
206 |
* pvm1(nlon,nlev), |
207 |
* ptm1(nlon,nlev), |
208 |
* pmea(nlon) |
209 |
real, intent(in):: pstd(nlon) |
210 |
real, intent(in):: psig(nlon) |
211 |
real pgamma(nlon),ptheta(nlon),ppic(nlon),pval(nlon), |
212 |
* pgeom1(nlon,nlev), |
213 |
* papm1(nlon,nlev), |
214 |
* paphm1(nlon,nlev+1) |
215 |
c |
216 |
integer kdx(nlon),ktest(nlon) |
217 |
c----------------------------------------------------------------------- |
218 |
c |
219 |
c* 0.2 local arrays |
220 |
c ------------ |
221 |
integer isect(klon), |
222 |
* icrit(klon), |
223 |
* ikcrith(klon), |
224 |
* ikenvh(klon), |
225 |
* iknu(klon), |
226 |
* iknu2(klon), |
227 |
* ikcrit(klon), |
228 |
* ikhlim(klon) |
229 |
c |
230 |
real ztau(klon,klev+1), |
231 |
$ ztauf(klon,klev+1), |
232 |
* zstab(klon,klev+1), |
233 |
* zvph(klon,klev+1), |
234 |
* zrho(klon,klev+1), |
235 |
* zri(klon,klev+1), |
236 |
* zpsi(klon,klev+1), |
237 |
* zzdep(klon,klev) |
238 |
real zdudt(klon), |
239 |
* zdvdt(klon), |
240 |
* zdtdt(klon), |
241 |
* zdedt(klon), |
242 |
* zvidis(klon), |
243 |
* znu(klon), |
244 |
* zd1(klon), |
245 |
* zd2(klon), |
246 |
* zdmod(klon) |
247 |
real ztmst, zrtmst |
248 |
real, intent(in):: ptsphy |
249 |
c |
250 |
c------------------------------------------------------------------ |
251 |
c |
252 |
c* 1. initialization |
253 |
c -------------- |
254 |
c |
255 |
100 continue |
256 |
c |
257 |
c ------------------------------------------------------------------ |
258 |
c |
259 |
c* 1.1 computational constants |
260 |
c ----------------------- |
261 |
c |
262 |
110 continue |
263 |
c |
264 |
c ztmst=twodt |
265 |
c if(nstep.eq.nstart) ztmst=0.5*twodt |
266 |
klevm1=klev-1 |
267 |
ztmst=ptsphy |
268 |
zrtmst=1./ztmst |
269 |
c ------------------------------------------------------------------ |
270 |
c |
271 |
120 continue |
272 |
c |
273 |
c ------------------------------------------------------------------ |
274 |
c |
275 |
c* 1.3 check whether row contains point for printing |
276 |
c --------------------------------------------- |
277 |
c |
278 |
130 continue |
279 |
c |
280 |
c ------------------------------------------------------------------ |
281 |
c |
282 |
c* 2. precompute basic state variables. |
283 |
c* ---------- ----- ----- ---------- |
284 |
c* define low level wind, project winds in plane of |
285 |
c* low level wind, determine sector in which to take |
286 |
c* the variance and set indicator for critical levels. |
287 |
c |
288 |
200 continue |
289 |
c |
290 |
c |
291 |
c |
292 |
call orosetup |
293 |
* ( nlon, ktest |
294 |
* , ikcrit, ikcrith, icrit, ikenvh,iknu,iknu2 |
295 |
* , paphm1, papm1 , pum1 , pvm1 , ptm1 , pgeom1, pstd |
296 |
* , zrho , zri , zstab , ztau , zvph , zpsi, zzdep |
297 |
* , pulow, pvlow |
298 |
* , ptheta,pgamma,pmea,ppic,pval,znu ,zd1, zd2, zdmod ) |
299 |
c |
300 |
c |
301 |
c |
302 |
c*********************************************************** |
303 |
c |
304 |
c |
305 |
c* 3. compute low level stresses using subcritical and |
306 |
c* supercritical forms.computes anisotropy coefficient |
307 |
c* as measure of orographic twodimensionality. |
308 |
c |
309 |
300 continue |
310 |
c |
311 |
call gwstress |
312 |
* ( nlon , nlev |
313 |
* , ktest , icrit, ikenvh, iknu |
314 |
* , zrho , zstab, zvph , pstd, psig, pmea, ppic |
315 |
* , ztau |
316 |
* , pgeom1,zdmod) |
317 |
c |
318 |
c |
319 |
c* 4. compute stress profile. |
320 |
c* ------- ------ -------- |
321 |
c |
322 |
400 continue |
323 |
c |
324 |
c |
325 |
call gwprofil |
326 |
* ( nlon , nlev |
327 |
* , kgwd , kdx , ktest |
328 |
* , ikcrith, icrit |
329 |
* , paphm1, zrho , zstab , zvph |
330 |
* , zri , ztau |
331 |
* , zdmod , psig , pstd) |
332 |
c |
333 |
c |
334 |
c* 5. compute tendencies. |
335 |
c* ------------------- |
336 |
c |
337 |
500 continue |
338 |
c |
339 |
c explicit solution at all levels for the gravity wave |
340 |
c implicit solution for the blocked levels |
341 |
|
342 |
do 510 jl=1,klon |
343 |
zvidis(jl)=0.0 |
344 |
zdudt(jl)=0.0 |
345 |
zdvdt(jl)=0.0 |
346 |
zdtdt(jl)=0.0 |
347 |
510 continue |
348 |
c |
349 |
ilevp1=klev+1 |
350 |
c |
351 |
c |
352 |
do 524 jk=1,klev |
353 |
c |
354 |
c |
355 |
c do 523 jl=1,kgwd |
356 |
c ji=kdx(jl) |
357 |
c Modif vectorisation 02/04/2004 |
358 |
do 523 ji=1,klon |
359 |
if(ktest(ji).eq.1) then |
360 |
|
361 |
zdelp=paphm1(ji,jk+1)-paphm1(ji,jk) |
362 |
ztemp=-rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,ilevp1)*zdelp) |
363 |
zdudt(ji)=(pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji) |
364 |
zdvdt(ji)=(pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji) |
365 |
c |
366 |
c controle des overshoots: |
367 |
c |
368 |
zforc=sqrt(zdudt(ji)**2+zdvdt(ji)**2)+1.E-12 |
369 |
ztend=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/ztmst+1.E-12 |
370 |
rover=0.25 |
371 |
if(zforc.ge.rover*ztend)then |
372 |
zdudt(ji)=rover*ztend/zforc*zdudt(ji) |
373 |
zdvdt(ji)=rover*ztend/zforc*zdvdt(ji) |
374 |
endif |
375 |
c |
376 |
c fin du controle des overshoots |
377 |
c |
378 |
if(jk.ge.ikenvh(ji)) then |
379 |
zb=1.0-0.18*pgamma(ji)-0.04*pgamma(ji)**2 |
380 |
zc=0.48*pgamma(ji)+0.3*pgamma(ji)**2 |
381 |
zconb=2.*ztmst*gkwake*psig(ji)/(4.*pstd(ji)) |
382 |
zabsv=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/2. |
383 |
zzd1=zb*cos(zpsi(ji,jk))**2+zc*sin(zpsi(ji,jk))**2 |
384 |
ratio=(cos(zpsi(ji,jk))**2+pgamma(ji)*sin(zpsi(ji,jk))**2)/ |
385 |
* (pgamma(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2) |
386 |
zbet=max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv |
387 |
c |
388 |
c simplement oppose au vent |
389 |
c |
390 |
zdudt(ji)=-pum1(ji,jk)/ztmst |
391 |
zdvdt(ji)=-pvm1(ji,jk)/ztmst |
392 |
c |
393 |
c projection dans la direction de l'axe principal de l'orographie |
394 |
cmod zdudt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.) |
395 |
cmod * +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.)) |
396 |
cmod * *cos(ptheta(ji)*rpi/180.)/ztmst |
397 |
cmod zdvdt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.) |
398 |
cmod * +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.)) |
399 |
cmod * *sin(ptheta(ji)*rpi/180.)/ztmst |
400 |
zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet)) |
401 |
zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet)) |
402 |
end if |
403 |
pvom(ji,jk)=zdudt(ji) |
404 |
pvol(ji,jk)=zdvdt(ji) |
405 |
zust=pum1(ji,jk)+ztmst*zdudt(ji) |
406 |
zvst=pvm1(ji,jk)+ztmst*zdvdt(ji) |
407 |
zdis=0.5*(pum1(ji,jk)**2+pvm1(ji,jk)**2-zust**2-zvst**2) |
408 |
zdedt(ji)=zdis/ztmst |
409 |
zvidis(ji)=zvidis(ji)+zdis*zdelp |
410 |
zdtdt(ji)=zdedt(ji)/rcpd |
411 |
c pte(ji,jk)=zdtdt(ji) |
412 |
c |
413 |
c ENCORE UN TRUC POUR EVITER LES EXPLOSIONS |
414 |
c |
415 |
pte(ji,jk)=0.0 |
416 |
|
417 |
endif |
418 |
523 continue |
419 |
|
420 |
524 continue |
421 |
c |
422 |
c |
423 |
return |
424 |
end |
425 |
SUBROUTINE orosetup |
426 |
* ( nlon , ktest |
427 |
* , kkcrit, kkcrith, kcrit |
428 |
* , kkenvh, kknu , kknu2 |
429 |
* , paphm1, papm1 , pum1 , pvm1 , ptm1 , pgeom1, pstd |
430 |
* , prho , pri , pstab , ptau , pvph ,ppsi, pzdep |
431 |
* , pulow , pvlow |
432 |
* , ptheta, pgamma, pmea, ppic, pval |
433 |
* , pnu , pd1 , pd2 ,pdmod ) |
434 |
c |
435 |
c**** *gwsetup* |
436 |
c |
437 |
c purpose. |
438 |
c -------- |
439 |
c |
440 |
c** interface. |
441 |
c ---------- |
442 |
c from *orodrag* |
443 |
c |
444 |
c explicit arguments : |
445 |
c -------------------- |
446 |
c ==== inputs === |
447 |
c ==== outputs === |
448 |
c |
449 |
c implicit arguments : none |
450 |
c -------------------- |
451 |
c |
452 |
c method. |
453 |
c ------- |
454 |
c |
455 |
c |
456 |
c externals. |
457 |
c ---------- |
458 |
c |
459 |
c |
460 |
c reference. |
461 |
c ---------- |
462 |
c |
463 |
c see ecmwf research department documentation of the "i.f.s." |
464 |
c |
465 |
c author. |
466 |
c ------- |
467 |
c |
468 |
c modifications. |
469 |
c -------------- |
470 |
c f.lott for the new-gwdrag scheme november 1993 |
471 |
c |
472 |
c----------------------------------------------------------------------- |
473 |
use dimens_m |
474 |
use dimphy |
475 |
use YOMCST |
476 |
use yoegwd |
477 |
implicit none |
478 |
c |
479 |
|
480 |
|
481 |
c----------------------------------------------------------------------- |
482 |
c |
483 |
c* 0.1 arguments |
484 |
c --------- |
485 |
c |
486 |
integer nlon |
487 |
integer jl, jk |
488 |
real zdelp |
489 |
|
490 |
integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon), |
491 |
* ktest(nlon),kkenvh(nlon) |
492 |
|
493 |
c |
494 |
real paphm1(nlon,klev+1),papm1(nlon,klev),pum1(nlon,klev), |
495 |
* pvm1(nlon,klev),ptm1(nlon,klev),pgeom1(nlon,klev), |
496 |
* prho(nlon,klev+1),pri(nlon,klev+1),pstab(nlon,klev+1), |
497 |
* ptau(nlon,klev+1),pvph(nlon,klev+1),ppsi(nlon,klev+1), |
498 |
* pzdep(nlon,klev) |
499 |
real pulow(nlon),pvlow(nlon),ptheta(nlon),pgamma(nlon),pnu(nlon), |
500 |
* pd1(nlon),pd2(nlon),pdmod(nlon) |
501 |
real, intent(in):: pstd(nlon) |
502 |
real pmea(nlon),ppic(nlon),pval(nlon) |
503 |
c |
504 |
c----------------------------------------------------------------------- |
505 |
c |
506 |
c* 0.2 local arrays |
507 |
c ------------ |
508 |
c |
509 |
c |
510 |
integer ilevm1, ilevm2, ilevh |
511 |
real zcons1, zcons2,zcons3, zhgeo |
512 |
real zu, zphi, zvt1,zvt2, zst, zvar, zdwind, zwind |
513 |
real zstabm, zstabp, zrhom, zrhop, alpha |
514 |
real zggeenv, zggeom1,zgvar |
515 |
logical lo |
516 |
logical ll1(klon,klev+1) |
517 |
integer kknu(klon),kknu2(klon),kknub(klon),kknul(klon), |
518 |
* kentp(klon),ncount(klon) |
519 |
c |
520 |
real zhcrit(klon,klev),zvpf(klon,klev), |
521 |
* zdp(klon,klev) |
522 |
real znorm(klon),zb(klon),zc(klon), |
523 |
* zulow(klon),zvlow(klon),znup(klon),znum(klon) |
524 |
c |
525 |
c ------------------------------------------------------------------ |
526 |
c |
527 |
c* 1. initialization |
528 |
c -------------- |
529 |
c |
530 |
c print *,' entree gwsetup' |
531 |
100 continue |
532 |
c |
533 |
c ------------------------------------------------------------------ |
534 |
c |
535 |
c* 1.1 computational constants |
536 |
c ----------------------- |
537 |
c |
538 |
110 continue |
539 |
c |
540 |
ilevm1=klev-1 |
541 |
ilevm2=klev-2 |
542 |
ilevh =klev/3 |
543 |
c |
544 |
zcons1=1./rd |
545 |
cold zcons2=g**2/cpd |
546 |
zcons2=rg**2/rcpd |
547 |
cold zcons3=1.5*api |
548 |
zcons3=1.5*rpi |
549 |
c |
550 |
c |
551 |
c ------------------------------------------------------------------ |
552 |
c |
553 |
c* 2. |
554 |
c -------------- |
555 |
c |
556 |
200 continue |
557 |
c |
558 |
c ------------------------------------------------------------------ |
559 |
c |
560 |
c* 2.1 define low level wind, project winds in plane of |
561 |
c* low level wind, determine sector in which to take |
562 |
c* the variance and set indicator for critical levels. |
563 |
c |
564 |
c |
565 |
c |
566 |
do 2001 jl=1,klon |
567 |
kknu(jl) =klev |
568 |
kknu2(jl) =klev |
569 |
kknub(jl) =klev |
570 |
kknul(jl) =klev |
571 |
pgamma(jl) =max(pgamma(jl),gtsec) |
572 |
ll1(jl,klev+1)=.false. |
573 |
2001 continue |
574 |
c |
575 |
c Ajouter une initialisation (L. Li, le 23fev99): |
576 |
c |
577 |
do jk=klev,ilevh,-1 |
578 |
do jl=1,klon |
579 |
ll1(jl,jk)= .FALSE. |
580 |
ENDDO |
581 |
ENDDO |
582 |
c |
583 |
c* define top of low level flow |
584 |
c ---------------------------- |
585 |
do 2002 jk=klev,ilevh,-1 |
586 |
do 2003 jl=1,klon |
587 |
lo=(paphm1(jl,jk)/paphm1(jl,klev+1)).ge.gsigcr |
588 |
if(lo) then |
589 |
kkcrit(jl)=jk |
590 |
endif |
591 |
zhcrit(jl,jk)=ppic(jl) |
592 |
zhgeo=pgeom1(jl,jk)/rg |
593 |
ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk)) |
594 |
if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then |
595 |
kknu(jl)=jk |
596 |
endif |
597 |
if(.not.ll1(jl,ilevh))kknu(jl)=ilevh |
598 |
2003 continue |
599 |
2002 continue |
600 |
do 2004 jk=klev,ilevh,-1 |
601 |
do 2005 jl=1,klon |
602 |
zhcrit(jl,jk)=ppic(jl)-pval(jl) |
603 |
zhgeo=pgeom1(jl,jk)/rg |
604 |
ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk)) |
605 |
if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then |
606 |
kknu2(jl)=jk |
607 |
endif |
608 |
if(.not.ll1(jl,ilevh))kknu2(jl)=ilevh |
609 |
2005 continue |
610 |
2004 continue |
611 |
do 2006 jk=klev,ilevh,-1 |
612 |
do 2007 jl=1,klon |
613 |
zhcrit(jl,jk)=amax1(ppic(jl)-pmea(jl),pmea(jl)-pval(jl)) |
614 |
zhgeo=pgeom1(jl,jk)/rg |
615 |
ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk)) |
616 |
if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then |
617 |
kknub(jl)=jk |
618 |
endif |
619 |
if(.not.ll1(jl,ilevh))kknub(jl)=ilevh |
620 |
2007 continue |
621 |
2006 continue |
622 |
c |
623 |
do 2010 jl=1,klon |
624 |
kknu(jl)=min(kknu(jl),nktopg) |
625 |
kknu2(jl)=min(kknu2(jl),nktopg) |
626 |
kknub(jl)=min(kknub(jl),nktopg) |
627 |
kknul(jl)=klev |
628 |
2010 continue |
629 |
c |
630 |
|
631 |
210 continue |
632 |
c |
633 |
c |
634 |
cc* initialize various arrays |
635 |
c |
636 |
do 2107 jl=1,klon |
637 |
prho(jl,klev+1) =0.0 |
638 |
pstab(jl,klev+1) =0.0 |
639 |
pstab(jl,1) =0.0 |
640 |
pri(jl,klev+1) =9999.0 |
641 |
ppsi(jl,klev+1) =0.0 |
642 |
pri(jl,1) =0.0 |
643 |
pvph(jl,1) =0.0 |
644 |
pulow(jl) =0.0 |
645 |
pvlow(jl) =0.0 |
646 |
zulow(jl) =0.0 |
647 |
zvlow(jl) =0.0 |
648 |
kkcrith(jl) =klev |
649 |
kkenvh(jl) =klev |
650 |
kentp(jl) =klev |
651 |
kcrit(jl) =1 |
652 |
ncount(jl) =0 |
653 |
ll1(jl,klev+1) =.false. |
654 |
2107 continue |
655 |
c |
656 |
c* define low-level flow |
657 |
c --------------------- |
658 |
c |
659 |
do 223 jk=klev,2,-1 |
660 |
do 222 jl=1,klon |
661 |
if(ktest(jl).eq.1) then |
662 |
zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1) |
663 |
prho(jl,jk)=2.*paphm1(jl,jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1)) |
664 |
pstab(jl,jk)=2.*zcons2/(ptm1(jl,jk)+ptm1(jl,jk-1))* |
665 |
* (1.-rcpd*prho(jl,jk)*(ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk)) |
666 |
pstab(jl,jk)=max(pstab(jl,jk),gssec) |
667 |
endif |
668 |
222 continue |
669 |
223 continue |
670 |
c |
671 |
c******************************************************************** |
672 |
c |
673 |
c* define blocked flow |
674 |
c ------------------- |
675 |
do 2115 jk=klev,ilevh,-1 |
676 |
do 2116 jl=1,klon |
677 |
if(jk.ge.kknub(jl).and.jk.le.kknul(jl)) then |
678 |
pulow(jl)=pulow(jl)+pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk)) |
679 |
pvlow(jl)=pvlow(jl)+pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk)) |
680 |
end if |
681 |
2116 continue |
682 |
2115 continue |
683 |
do 2110 jl=1,klon |
684 |
pulow(jl)=pulow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknub(jl))) |
685 |
pvlow(jl)=pvlow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknub(jl))) |
686 |
znorm(jl)=max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec) |
687 |
pvph(jl,klev+1)=znorm(jl) |
688 |
2110 continue |
689 |
c |
690 |
c******* setup orography axes and define plane of profiles ******* |
691 |
c |
692 |
do 2112 jl=1,klon |
693 |
lo=(pulow(jl).lt.gvsec).and.(pulow(jl).ge.-gvsec) |
694 |
if(lo) then |
695 |
zu=pulow(jl)+2.*gvsec |
696 |
else |
697 |
zu=pulow(jl) |
698 |
endif |
699 |
zphi=atan(pvlow(jl)/zu) |
700 |
ppsi(jl,klev+1)=ptheta(jl)*rpi/180.-zphi |
701 |
zb(jl)=1.-0.18*pgamma(jl)-0.04*pgamma(jl)**2 |
702 |
zc(jl)=0.48*pgamma(jl)+0.3*pgamma(jl)**2 |
703 |
pd1(jl)=zb(jl)-(zb(jl)-zc(jl))*(sin(ppsi(jl,klev+1))**2) |
704 |
pd2(jl)=(zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))*cos(ppsi(jl,klev+1)) |
705 |
pdmod(jl)=sqrt(pd1(jl)**2+pd2(jl)**2) |
706 |
2112 continue |
707 |
c |
708 |
c ************ define flow in plane of lowlevel stress ************* |
709 |
c |
710 |
do 213 jk=1,klev |
711 |
do 212 jl=1,klon |
712 |
if(ktest(jl).eq.1) then |
713 |
zvt1 =pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk) |
714 |
zvt2 =-pvlow(jl)*pum1(jl,jk)+pulow(jl)*pvm1(jl,jk) |
715 |
zvpf(jl,jk)=(zvt1*pd1(jl)+zvt2*pd2(jl))/(znorm(jl)*pdmod(jl)) |
716 |
endif |
717 |
ptau(jl,jk) =0.0 |
718 |
pzdep(jl,jk) =0.0 |
719 |
ppsi(jl,jk) =0.0 |
720 |
ll1(jl,jk) =.false. |
721 |
212 continue |
722 |
213 continue |
723 |
do 215 jk=2,klev |
724 |
do 214 jl=1,klon |
725 |
if(ktest(jl).eq.1) then |
726 |
zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1) |
727 |
pvph(jl,jk)=((paphm1(jl,jk)-papm1(jl,jk-1))*zvpf(jl,jk)+ |
728 |
* (papm1(jl,jk)-paphm1(jl,jk))*zvpf(jl,jk-1)) |
729 |
* /zdp(jl,jk) |
730 |
if(pvph(jl,jk).lt.gvsec) then |
731 |
pvph(jl,jk)=gvsec |
732 |
kcrit(jl)=jk |
733 |
endif |
734 |
endif |
735 |
214 continue |
736 |
215 continue |
737 |
c |
738 |
c |
739 |
c* 2.2 brunt-vaisala frequency and density at half levels. |
740 |
c |
741 |
220 continue |
742 |
c |
743 |
do 2211 jk=ilevh,klev |
744 |
do 221 jl=1,klon |
745 |
if(ktest(jl).eq.1) then |
746 |
if(jk.ge.(kknub(jl)+1).and.jk.le.kknul(jl)) then |
747 |
zst=zcons2/ptm1(jl,jk)*(1.-rcpd*prho(jl,jk)* |
748 |
* (ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk)) |
749 |
pstab(jl,klev+1)=pstab(jl,klev+1)+zst*zdp(jl,jk) |
750 |
pstab(jl,klev+1)=max(pstab(jl,klev+1),gssec) |
751 |
prho(jl,klev+1)=prho(jl,klev+1)+paphm1(jl,jk)*2.*zdp(jl,jk) |
752 |
* *zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1)) |
753 |
endif |
754 |
endif |
755 |
221 continue |
756 |
2211 continue |
757 |
c |
758 |
do 2212 jl=1,klon |
759 |
pstab(jl,klev+1)=pstab(jl,klev+1)/(papm1(jl,kknul(jl)) |
760 |
* -papm1(jl,kknub(jl))) |
761 |
prho(jl,klev+1)=prho(jl,klev+1)/(papm1(jl,kknul(jl)) |
762 |
* -papm1(jl,kknub(jl))) |
763 |
zvar=pstd(jl) |
764 |
2212 continue |
765 |
c |
766 |
c* 2.3 mean flow richardson number. |
767 |
c* and critical height for froude layer |
768 |
c |
769 |
230 continue |
770 |
c |
771 |
do 232 jk=2,klev |
772 |
do 231 jl=1,klon |
773 |
if(ktest(jl).eq.1) then |
774 |
zdwind=max(abs(zvpf(jl,jk)-zvpf(jl,jk-1)),gvsec) |
775 |
pri(jl,jk)=pstab(jl,jk)*(zdp(jl,jk) |
776 |
* /(rg*prho(jl,jk)*zdwind))**2 |
777 |
pri(jl,jk)=max(pri(jl,jk),grcrit) |
778 |
endif |
779 |
231 continue |
780 |
232 continue |
781 |
|
782 |
c |
783 |
c |
784 |
c* define top of 'envelope' layer |
785 |
c ---------------------------- |
786 |
|
787 |
do 233 jl=1,klon |
788 |
pnu (jl)=0.0 |
789 |
znum(jl)=0.0 |
790 |
233 continue |
791 |
|
792 |
do 234 jk=2,klev-1 |
793 |
do 234 jl=1,klon |
794 |
|
795 |
if(ktest(jl).eq.1) then |
796 |
|
797 |
if (jk.ge.kknub(jl)) then |
798 |
|
799 |
znum(jl)=pnu(jl) |
800 |
zwind=(pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/ |
801 |
* max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec) |
802 |
zwind=max(sqrt(zwind**2),gvsec) |
803 |
zdelp=paphm1(jl,jk+1)-paphm1(jl,jk) |
804 |
zstabm=sqrt(max(pstab(jl,jk ),gssec)) |
805 |
zstabp=sqrt(max(pstab(jl,jk+1),gssec)) |
806 |
zrhom=prho(jl,jk ) |
807 |
zrhop=prho(jl,jk+1) |
808 |
pnu(jl) = pnu(jl) + (zdelp/rg)* |
809 |
* ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind |
810 |
if((znum(jl).le.gfrcrit).and.(pnu(jl).gt.gfrcrit) |
811 |
* .and.(kkenvh(jl).eq.klev)) |
812 |
* kkenvh(jl)=jk |
813 |
|
814 |
endif |
815 |
|
816 |
endif |
817 |
|
818 |
234 continue |
819 |
|
820 |
c calculation of a dynamical mixing height for the breaking |
821 |
c of gravity waves: |
822 |
|
823 |
|
824 |
do 235 jl=1,klon |
825 |
znup(jl)=0.0 |
826 |
znum(jl)=0.0 |
827 |
235 continue |
828 |
|
829 |
do 236 jk=klev-1,2,-1 |
830 |
do 236 jl=1,klon |
831 |
|
832 |
if(ktest(jl).eq.1) then |
833 |
|
834 |
znum(jl)=znup(jl) |
835 |
zwind=(pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/ |
836 |
* max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec) |
837 |
zwind=max(sqrt(zwind**2),gvsec) |
838 |
zdelp=paphm1(jl,jk+1)-paphm1(jl,jk) |
839 |
zstabm=sqrt(max(pstab(jl,jk ),gssec)) |
840 |
zstabp=sqrt(max(pstab(jl,jk+1),gssec)) |
841 |
zrhom=prho(jl,jk ) |
842 |
zrhop=prho(jl,jk+1) |
843 |
znup(jl) = znup(jl) + (zdelp/rg)* |
844 |
* ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind |
845 |
if((znum(jl).le.rpi/2.).and.(znup(jl).gt.rpi/2.) |
846 |
* .and.(kkcrith(jl).eq.klev)) |
847 |
* kkcrith(jl)=jk |
848 |
|
849 |
endif |
850 |
|
851 |
236 continue |
852 |
|
853 |
do 237 jl=1,klon |
854 |
kkcrith(jl)=min0(kkcrith(jl),kknu2(jl)) |
855 |
kkcrith(jl)=max0(kkcrith(jl),ilevh*2) |
856 |
237 continue |
857 |
c |
858 |
c directional info for flow blocking ************************* |
859 |
c |
860 |
do 251 jk=ilevh,klev |
861 |
do 252 jl=1,klon |
862 |
if(jk.ge.kkenvh(jl)) then |
863 |
lo=(pum1(jl,jk).lt.gvsec).and.(pum1(jl,jk).ge.-gvsec) |
864 |
if(lo) then |
865 |
zu=pum1(jl,jk)+2.*gvsec |
866 |
else |
867 |
zu=pum1(jl,jk) |
868 |
endif |
869 |
zphi=atan(pvm1(jl,jk)/zu) |
870 |
ppsi(jl,jk)=ptheta(jl)*rpi/180.-zphi |
871 |
end if |
872 |
252 continue |
873 |
251 continue |
874 |
c forms the vertical 'leakiness' ************************** |
875 |
|
876 |
alpha=3. |
877 |
|
878 |
do 254 jk=ilevh,klev |
879 |
do 253 jl=1,klon |
880 |
if(jk.ge.kkenvh(jl)) then |
881 |
zggeenv=amax1(1., |
882 |
* (pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)-1))/2.) |
883 |
zggeom1=amax1(pgeom1(jl,jk),1.) |
884 |
zgvar=amax1(pstd(jl)*rg,1.) |
885 |
cmod pzdep(jl,jk)=sqrt((zggeenv-zggeom1)/(zggeom1+zgvar)) |
886 |
pzdep(jl,jk)=(pgeom1(jl,kkenvh(jl)-1)-pgeom1(jl, jk))/ |
887 |
* (pgeom1(jl,kkenvh(jl)-1)-pgeom1(jl,klev)) |
888 |
end if |
889 |
253 continue |
890 |
254 continue |
891 |
|
892 |
260 continue |
893 |
|
894 |
return |
895 |
end |
896 |
SUBROUTINE gwstress |
897 |
* ( nlon , nlev |
898 |
* , ktest, kcrit, kkenvh |
899 |
* , kknu |
900 |
* , prho , pstab , pvph , pstd, psig |
901 |
* , pmea , ppic , ptau |
902 |
* , pgeom1 , pdmod ) |
903 |
c |
904 |
c**** *gwstress* |
905 |
c |
906 |
c purpose. |
907 |
c -------- |
908 |
c |
909 |
c** interface. |
910 |
c ---------- |
911 |
c call *gwstress* from *gwdrag* |
912 |
c |
913 |
c explicit arguments : |
914 |
c -------------------- |
915 |
c ==== inputs === |
916 |
c ==== outputs === |
917 |
c |
918 |
c implicit arguments : none |
919 |
c -------------------- |
920 |
c |
921 |
c method. |
922 |
c ------- |
923 |
c |
924 |
c |
925 |
c externals. |
926 |
c ---------- |
927 |
c |
928 |
c |
929 |
c reference. |
930 |
c ---------- |
931 |
c |
932 |
c see ecmwf research department documentation of the "i.f.s." |
933 |
c |
934 |
c author. |
935 |
c ------- |
936 |
c |
937 |
c modifications. |
938 |
c -------------- |
939 |
c f. lott put the new gwd on ifs 22/11/93 |
940 |
c |
941 |
c----------------------------------------------------------------------- |
942 |
use dimens_m |
943 |
use dimphy |
944 |
use YOMCST |
945 |
use yoegwd |
946 |
implicit none |
947 |
|
948 |
c----------------------------------------------------------------------- |
949 |
c |
950 |
c* 0.1 arguments |
951 |
c --------- |
952 |
c |
953 |
integer nlon, nlev |
954 |
integer kcrit(nlon), |
955 |
* ktest(nlon),kkenvh(nlon),kknu(nlon) |
956 |
c |
957 |
real prho(nlon,nlev+1),pstab(nlon,nlev+1),ptau(nlon,nlev+1), |
958 |
* pvph(nlon,nlev+1), |
959 |
* pgeom1(nlon,nlev) |
960 |
real, intent(in):: pstd(nlon) |
961 |
c |
962 |
real, intent(in):: psig(nlon) |
963 |
real pmea(nlon),ppic(nlon) |
964 |
real pdmod(nlon) |
965 |
c |
966 |
c----------------------------------------------------------------------- |
967 |
c |
968 |
c* 0.2 local arrays |
969 |
c ------------ |
970 |
integer jl |
971 |
real zblock, zvar, zeff |
972 |
logical lo |
973 |
c |
974 |
c----------------------------------------------------------------------- |
975 |
c |
976 |
c* 0.3 functions |
977 |
c --------- |
978 |
c ------------------------------------------------------------------ |
979 |
c |
980 |
c* 1. initialization |
981 |
c -------------- |
982 |
c |
983 |
100 continue |
984 |
c |
985 |
c* 3.1 gravity wave stress. |
986 |
c |
987 |
300 continue |
988 |
c |
989 |
c |
990 |
do 301 jl=1,klon |
991 |
if(ktest(jl).eq.1) then |
992 |
|
993 |
c effective mountain height above the blocked flow |
994 |
|
995 |
if(kkenvh(jl).eq.klev)then |
996 |
zblock=0.0 |
997 |
else |
998 |
zblock=(pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)+1))/2./rg |
999 |
endif |
1000 |
|
1001 |
zvar=ppic(jl)-pmea(jl) |
1002 |
zeff=amax1(0.,zvar-zblock) |
1003 |
|
1004 |
ptau(jl,klev+1)=prho(jl,klev+1)*gkdrag*psig(jl)*zeff**2 |
1005 |
* /4./pstd(jl)*pvph(jl,klev+1)*pdmod(jl)*sqrt(pstab(jl,klev+1)) |
1006 |
|
1007 |
c too small value of stress or low level flow include critical level |
1008 |
c or low level flow: gravity wave stress nul. |
1009 |
|
1010 |
lo=(ptau(jl,klev+1).lt.gtsec).or.(kcrit(jl).ge.kknu(jl)) |
1011 |
* .or.(pvph(jl,klev+1).lt.gvcrit) |
1012 |
c if(lo) ptau(jl,klev+1)=0.0 |
1013 |
|
1014 |
else |
1015 |
|
1016 |
ptau(jl,klev+1)=0.0 |
1017 |
|
1018 |
endif |
1019 |
|
1020 |
301 continue |
1021 |
c |
1022 |
return |
1023 |
end |
1024 |
SUBROUTINE GWPROFIL |
1025 |
* ( NLON, NLEV |
1026 |
* , kgwd, kdx , ktest |
1027 |
* , KKCRITH, KCRIT |
1028 |
* , PAPHM1, PRHO , PSTAB , PVPH , PRI , PTAU |
1029 |
* , pdmod , psig , pvar) |
1030 |
|
1031 |
C**** *GWPROFIL* |
1032 |
C |
1033 |
C PURPOSE. |
1034 |
C -------- |
1035 |
C |
1036 |
C** INTERFACE. |
1037 |
C ---------- |
1038 |
C FROM *GWDRAG* |
1039 |
C |
1040 |
C EXPLICIT ARGUMENTS : |
1041 |
C -------------------- |
1042 |
C ==== INPUTS === |
1043 |
C ==== OUTPUTS === |
1044 |
C |
1045 |
C IMPLICIT ARGUMENTS : NONE |
1046 |
C -------------------- |
1047 |
C |
1048 |
C METHOD: |
1049 |
C ------- |
1050 |
C THE STRESS PROFILE FOR GRAVITY WAVES IS COMPUTED AS FOLLOWS: |
1051 |
C IT IS CONSTANT (NO GWD) AT THE LEVELS BETWEEN THE GROUND |
1052 |
C AND THE TOP OF THE BLOCKED LAYER (KKENVH). |
1053 |
C IT DECREASES LINEARLY WITH HEIGHTS FROM THE TOP OF THE |
1054 |
C BLOCKED LAYER TO 3*VAROR (kKNU), TO SIMULATES LEE WAVES OR |
1055 |
C NONLINEAR GRAVITY WAVE BREAKING. |
1056 |
C ABOVE IT IS CONSTANT, EXCEPT WHEN THE WAVE ENCOUNTERS A CRITICAL |
1057 |
C LEVEL (KCRIT) OR WHEN IT BREAKS. |
1058 |
C |
1059 |
C |
1060 |
C |
1061 |
C EXTERNALS. |
1062 |
C ---------- |
1063 |
C |
1064 |
C |
1065 |
C REFERENCE. |
1066 |
C ---------- |
1067 |
C |
1068 |
C SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S." |
1069 |
C |
1070 |
C AUTHOR. |
1071 |
C ------- |
1072 |
C |
1073 |
C MODIFICATIONS. |
1074 |
C -------------- |
1075 |
C PASSAGE OF THE NEW GWDRAG TO I.F.S. (F. LOTT, 22/11/93) |
1076 |
C----------------------------------------------------------------------- |
1077 |
use dimens_m |
1078 |
use dimphy |
1079 |
use YOMCST |
1080 |
use yoegwd |
1081 |
implicit none |
1082 |
C |
1083 |
|
1084 |
C |
1085 |
|
1086 |
|
1087 |
C----------------------------------------------------------------------- |
1088 |
C |
1089 |
C* 0.1 ARGUMENTS |
1090 |
C --------- |
1091 |
C |
1092 |
integer nlon,nlev |
1093 |
INTEGER KKCRITH(NLON),KCRIT(NLON) |
1094 |
* ,kdx(nlon) , ktest(nlon) |
1095 |
|
1096 |
C |
1097 |
REAL PAPHM1(NLON,NLEV+1), PSTAB(NLON,NLEV+1), |
1098 |
* PRHO (NLON,NLEV+1), PVPH (NLON,NLEV+1), |
1099 |
* PRI (NLON,NLEV+1), PTAU(NLON,NLEV+1) |
1100 |
|
1101 |
REAL pdmod (NLON) |
1102 |
real, intent(in):: psig(NLON) |
1103 |
real, intent(in):: pvar(NLON) |
1104 |
|
1105 |
C----------------------------------------------------------------------- |
1106 |
C |
1107 |
C* 0.2 LOCAL ARRAYS |
1108 |
C ------------ |
1109 |
C |
1110 |
integer ilevh, ji, kgwd, jl, jk |
1111 |
real zsqr, zalfa, zriw, zdel, zb, zalpha,zdz2n |
1112 |
real zdelp, zdelpt |
1113 |
REAL ZDZ2 (KLON,KLEV) , ZNORM(KLON) , zoro(KLON) |
1114 |
REAL ZTAU (KLON,KLEV+1) |
1115 |
C |
1116 |
C----------------------------------------------------------------------- |
1117 |
C |
1118 |
C* 1. INITIALIZATION |
1119 |
C -------------- |
1120 |
C |
1121 |
c print *,' entree gwprofil' |
1122 |
100 CONTINUE |
1123 |
C |
1124 |
C |
1125 |
C* COMPUTATIONAL CONSTANTS. |
1126 |
C ------------- ---------- |
1127 |
C |
1128 |
ilevh=KLEV/3 |
1129 |
C |
1130 |
c DO 400 ji=1,kgwd |
1131 |
c jl=kdx(ji) |
1132 |
c Modif vectorisation 02/04/2004 |
1133 |
DO 400 jl=1,klon |
1134 |
if (ktest(jl).eq.1) then |
1135 |
Zoro(JL)=Psig(JL)*Pdmod(JL)/4./max(pvar(jl),1.0) |
1136 |
ZTAU(JL,KLEV+1)=PTAU(JL,KLEV+1) |
1137 |
endif |
1138 |
400 CONTINUE |
1139 |
|
1140 |
C |
1141 |
DO 430 JK=KLEV,2,-1 |
1142 |
C |
1143 |
C |
1144 |
C* 4.1 CONSTANT WAVE STRESS UNTIL TOP OF THE |
1145 |
C BLOCKING LAYER. |
1146 |
410 CONTINUE |
1147 |
C |
1148 |
c DO 411 ji=1,kgwd |
1149 |
c jl=kdx(ji) |
1150 |
c Modif vectorisation 02/04/2004 |
1151 |
do 411 jl=1,klon |
1152 |
if (ktest(jl).eq.1) then |
1153 |
IF(JK.GT.KKCRITH(JL)) THEN |
1154 |
PTAU(JL,JK)=ZTAU(JL,KLEV+1) |
1155 |
C ENDIF |
1156 |
C IF(JK.EQ.KKCRITH(JL)) THEN |
1157 |
ELSE |
1158 |
PTAU(JL,JK)=GRAHILO*ZTAU(JL,KLEV+1) |
1159 |
ENDIF |
1160 |
endif |
1161 |
411 CONTINUE |
1162 |
C |
1163 |
C* 4.15 CONSTANT SHEAR STRESS UNTIL THE TOP OF THE |
1164 |
C LOW LEVEL FLOW LAYER. |
1165 |
415 CONTINUE |
1166 |
C |
1167 |
C |
1168 |
C* 4.2 WAVE DISPLACEMENT AT NEXT LEVEL. |
1169 |
C |
1170 |
420 CONTINUE |
1171 |
C |
1172 |
c DO 421 ji=1,kgwd |
1173 |
c jl=kdx(ji) |
1174 |
c Modif vectorisation 02/04/2004 |
1175 |
do 421 jl=1,klon |
1176 |
if(ktest(jl).eq.1) then |
1177 |
IF(JK.LT.KKCRITH(JL)) THEN |
1178 |
ZNORM(JL)=gkdrag*PRHO(JL,JK)*SQRT(PSTAB(JL,JK))*PVPH(JL,JK) |
1179 |
* *zoro(jl) |
1180 |
ZDZ2(JL,JK)=PTAU(JL,JK+1)/max(ZNORM(JL),gssec) |
1181 |
ENDIF |
1182 |
endif |
1183 |
421 CONTINUE |
1184 |
C |
1185 |
C* 4.3 WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT |
1186 |
C* AND STRESS: BREAKING EVALUATION AND CRITICAL |
1187 |
C LEVEL |
1188 |
C |
1189 |
|
1190 |
c DO 431 ji=1,kgwd |
1191 |
c jl=Kdx(ji) |
1192 |
c Modif vectorisation 02/04/2004 |
1193 |
do 431 jl=1,klon |
1194 |
if(ktest(jl).eq.1) then |
1195 |
|
1196 |
IF(JK.LT.KKCRITH(JL)) THEN |
1197 |
IF((PTAU(JL,JK+1).LT.GTSEC).OR.(JK.LE.KCRIT(JL))) THEN |
1198 |
PTAU(JL,JK)=0.0 |
1199 |
ELSE |
1200 |
ZSQR=SQRT(PRI(JL,JK)) |
1201 |
ZALFA=SQRT(PSTAB(JL,JK)*ZDZ2(JL,JK))/PVPH(JL,JK) |
1202 |
ZRIW=PRI(JL,JK)*(1.-ZALFA)/(1+ZALFA*ZSQR)**2 |
1203 |
IF(ZRIW.LT.GRCRIT) THEN |
1204 |
ZDEL=4./ZSQR/GRCRIT+1./GRCRIT**2+4./GRCRIT |
1205 |
ZB=1./GRCRIT+2./ZSQR |
1206 |
ZALPHA=0.5*(-ZB+SQRT(ZDEL)) |
1207 |
ZDZ2N=(PVPH(JL,JK)*ZALPHA)**2/PSTAB(JL,JK) |
1208 |
PTAU(JL,JK)=ZNORM(JL)*ZDZ2N |
1209 |
ELSE |
1210 |
PTAU(JL,JK)=ZNORM(JL)*ZDZ2(JL,JK) |
1211 |
ENDIF |
1212 |
PTAU(JL,JK)=MIN(PTAU(JL,JK),PTAU(JL,JK+1)) |
1213 |
ENDIF |
1214 |
ENDIF |
1215 |
endif |
1216 |
431 CONTINUE |
1217 |
|
1218 |
430 CONTINUE |
1219 |
440 CONTINUE |
1220 |
|
1221 |
C REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL |
1222 |
|
1223 |
c DO 530 ji=1,kgwd |
1224 |
c jl=kdx(ji) |
1225 |
c Modif vectorisation 02/04/2004 |
1226 |
do 530 jl=1,klon |
1227 |
if(ktest(jl).eq.1) then |
1228 |
ZTAU(JL,KKCRITH(JL))=PTAU(JL,KKCRITH(JL)) |
1229 |
ZTAU(JL,NSTRA)=PTAU(JL,NSTRA) |
1230 |
endif |
1231 |
530 CONTINUE |
1232 |
|
1233 |
DO 531 JK=1,KLEV |
1234 |
|
1235 |
c DO 532 ji=1,kgwd |
1236 |
c jl=kdx(ji) |
1237 |
c Modif vectorisation 02/04/2004 |
1238 |
do 532 jl=1,klon |
1239 |
if(ktest(jl).eq.1) then |
1240 |
|
1241 |
|
1242 |
IF(JK.GT.KKCRITH(JL))THEN |
1243 |
|
1244 |
ZDELP=PAPHM1(JL,JK)-PAPHM1(JL,KLEV+1 ) |
1245 |
ZDELPT=PAPHM1(JL,KKCRITH(JL))-PAPHM1(JL,KLEV+1 ) |
1246 |
PTAU(JL,JK)=ZTAU(JL,KLEV+1 ) + |
1247 |
. (ZTAU(JL,KKCRITH(JL))-ZTAU(JL,KLEV+1 ) )* |
1248 |
. ZDELP/ZDELPT |
1249 |
|
1250 |
ENDIF |
1251 |
|
1252 |
endif |
1253 |
532 CONTINUE |
1254 |
|
1255 |
C REORGANISATION IN THE STRATOSPHERE |
1256 |
|
1257 |
c DO 533 ji=1,kgwd |
1258 |
c jl=kdx(ji) |
1259 |
c Modif vectorisation 02/04/2004 |
1260 |
do 533 jl=1,klon |
1261 |
if(ktest(jl).eq.1) then |
1262 |
|
1263 |
|
1264 |
IF(JK.LT.NSTRA)THEN |
1265 |
|
1266 |
ZDELP =PAPHM1(JL,NSTRA) |
1267 |
ZDELPT=PAPHM1(JL,JK) |
1268 |
PTAU(JL,JK)=ZTAU(JL,NSTRA)*ZDELPT/ZDELP |
1269 |
|
1270 |
ENDIF |
1271 |
|
1272 |
endif |
1273 |
533 CONTINUE |
1274 |
|
1275 |
C REORGANISATION IN THE TROPOSPHERE |
1276 |
|
1277 |
c DO 534 ji=1,kgwd |
1278 |
c jl=kdx(ji) |
1279 |
c Modif vectorisation 02/04/2004 |
1280 |
do 534 jl=1,klon |
1281 |
if(ktest(jl).eq.1) then |
1282 |
|
1283 |
|
1284 |
IF(JK.LT.KKCRITH(JL).AND.JK.GT.NSTRA)THEN |
1285 |
|
1286 |
ZDELP=PAPHM1(JL,JK)-PAPHM1(JL,KKCRITH(JL)) |
1287 |
ZDELPT=PAPHM1(JL,NSTRA)-PAPHM1(JL,KKCRITH(JL)) |
1288 |
PTAU(JL,JK)=ZTAU(JL,KKCRITH(JL)) + |
1289 |
* (ZTAU(JL,NSTRA)-ZTAU(JL,KKCRITH(JL)))*ZDELP |
1290 |
. /ZDELPT |
1291 |
|
1292 |
ENDIF |
1293 |
endif |
1294 |
534 CONTINUE |
1295 |
|
1296 |
|
1297 |
531 CONTINUE |
1298 |
|
1299 |
|
1300 |
RETURN |
1301 |
END |
1302 |
SUBROUTINE lift_noro (nlon,nlev,dtime,paprs,pplay, |
1303 |
e plat,pmea,pstd, ppic, |
1304 |
e ktest, |
1305 |
e t, u, v, |
1306 |
s pulow, pvlow, pustr, pvstr, |
1307 |
s d_t, d_u, d_v) |
1308 |
c |
1309 |
use dimens_m |
1310 |
use dimphy |
1311 |
use YOMCST |
1312 |
IMPLICIT none |
1313 |
c====================================================================== |
1314 |
c Auteur(s): F.Lott (LMD/CNRS) date: 19950201 |
1315 |
c Objet: Frottement de la montagne Interface |
1316 |
c====================================================================== |
1317 |
c Arguments: |
1318 |
c dtime---input-R- pas d'integration (s) |
1319 |
c paprs---input-R-pression pour chaque inter-couche (en Pa) |
1320 |
c pplay---input-R-pression pour le mileu de chaque couche (en Pa) |
1321 |
c t-------input-R-temperature (K) |
1322 |
c u-------input-R-vitesse horizontale (m/s) |
1323 |
c v-------input-R-vitesse horizontale (m/s) |
1324 |
c |
1325 |
c d_t-----output-R-increment de la temperature |
1326 |
c d_u-----output-R-increment de la vitesse u |
1327 |
c d_v-----output-R-increment de la vitesse v |
1328 |
c====================================================================== |
1329 |
c |
1330 |
c ARGUMENTS |
1331 |
c |
1332 |
INTEGER nlon,nlev |
1333 |
REAL, intent(in):: dtime |
1334 |
REAL, intent(in):: paprs(klon,klev+1) |
1335 |
REAL, intent(in):: pplay(klon,klev) |
1336 |
REAL, intent(in):: plat(nlon) |
1337 |
real pmea(nlon) |
1338 |
REAL, intent(in):: pstd(nlon) |
1339 |
REAL ppic(nlon) |
1340 |
REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon) |
1341 |
REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev) |
1342 |
REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev) |
1343 |
c |
1344 |
INTEGER i, k, ktest(nlon) |
1345 |
c |
1346 |
c Variables locales: |
1347 |
c |
1348 |
REAL zgeom(klon,klev) |
1349 |
REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev) |
1350 |
REAL pt(klon,klev), pu(klon,klev), pv(klon,klev) |
1351 |
REAL papmf(klon,klev),papmh(klon,klev+1) |
1352 |
c |
1353 |
c initialiser les variables de sortie (pour securite) |
1354 |
c |
1355 |
DO i = 1,klon |
1356 |
pulow(i) = 0.0 |
1357 |
pvlow(i) = 0.0 |
1358 |
pustr(i) = 0.0 |
1359 |
pvstr(i) = 0.0 |
1360 |
ENDDO |
1361 |
DO k = 1, klev |
1362 |
DO i = 1, klon |
1363 |
d_t(i,k) = 0.0 |
1364 |
d_u(i,k) = 0.0 |
1365 |
d_v(i,k) = 0.0 |
1366 |
pdudt(i,k)=0.0 |
1367 |
pdvdt(i,k)=0.0 |
1368 |
pdtdt(i,k)=0.0 |
1369 |
ENDDO |
1370 |
ENDDO |
1371 |
c |
1372 |
c preparer les variables d'entree (attention: l'ordre des niveaux |
1373 |
c verticaux augmente du haut vers le bas) |
1374 |
c |
1375 |
DO k = 1, klev |
1376 |
DO i = 1, klon |
1377 |
pt(i,k) = t(i,klev-k+1) |
1378 |
pu(i,k) = u(i,klev-k+1) |
1379 |
pv(i,k) = v(i,klev-k+1) |
1380 |
papmf(i,k) = pplay(i,klev-k+1) |
1381 |
ENDDO |
1382 |
ENDDO |
1383 |
DO k = 1, klev+1 |
1384 |
DO i = 1, klon |
1385 |
papmh(i,k) = paprs(i,klev-k+2) |
1386 |
ENDDO |
1387 |
ENDDO |
1388 |
DO i = 1, klon |
1389 |
zgeom(i,klev) = RD * pt(i,klev) |
1390 |
. * LOG(papmh(i,klev+1)/papmf(i,klev)) |
1391 |
ENDDO |
1392 |
DO k = klev-1, 1, -1 |
1393 |
DO i = 1, klon |
1394 |
zgeom(i,k) = zgeom(i,k+1) + RD * (pt(i,k)+pt(i,k+1))/2.0 |
1395 |
. * LOG(papmf(i,k+1)/papmf(i,k)) |
1396 |
ENDDO |
1397 |
ENDDO |
1398 |
c |
1399 |
c appeler la routine principale |
1400 |
c |
1401 |
CALL OROLIFT(klon,klev,ktest, |
1402 |
. dtime, |
1403 |
. papmh, zgeom, |
1404 |
. pt, pu, pv, |
1405 |
. plat,pmea, pstd, ppic, |
1406 |
. pulow,pvlow, |
1407 |
. pdudt,pdvdt,pdtdt) |
1408 |
C |
1409 |
DO k = 1, klev |
1410 |
DO i = 1, klon |
1411 |
d_u(i,klev+1-k) = dtime*pdudt(i,k) |
1412 |
d_v(i,klev+1-k) = dtime*pdvdt(i,k) |
1413 |
d_t(i,klev+1-k) = dtime*pdtdt(i,k) |
1414 |
pustr(i) = pustr(i) |
1415 |
cIM BUG . +RG*pdudt(i,k)*(papmh(i,k+1)-papmh(i,k)) |
1416 |
. +pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))/RG |
1417 |
pvstr(i) = pvstr(i) |
1418 |
cIM BUG . +RG*pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k)) |
1419 |
. +pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))/RG |
1420 |
ENDDO |
1421 |
ENDDO |
1422 |
c |
1423 |
RETURN |
1424 |
END |
1425 |
SUBROUTINE OROLIFT( NLON,NLEV |
1426 |
I , KTEST |
1427 |
R , PTSPHY |
1428 |
R , PAPHM1,PGEOM1,PTM1,PUM1,PVM1 |
1429 |
R , PLAT |
1430 |
R , PMEA, PVAROR, ppic |
1431 |
C OUTPUTS |
1432 |
R , PULOW,PVLOW |
1433 |
R , PVOM,PVOL,PTE ) |
1434 |
|
1435 |
C |
1436 |
C**** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT. |
1437 |
C |
1438 |
C PURPOSE. |
1439 |
C -------- |
1440 |
C |
1441 |
C** INTERFACE. |
1442 |
C ---------- |
1443 |
C CALLED FROM *lift_noro |
1444 |
C ---------- |
1445 |
C |
1446 |
C AUTHOR. |
1447 |
C ------- |
1448 |
C F.LOTT LMD 22/11/95 |
1449 |
C |
1450 |
use dimens_m |
1451 |
use dimphy |
1452 |
use YOMCST |
1453 |
use yoegwd |
1454 |
implicit none |
1455 |
C |
1456 |
C |
1457 |
C----------------------------------------------------------------------- |
1458 |
C |
1459 |
C* 0.1 ARGUMENTS |
1460 |
C --------- |
1461 |
C |
1462 |
C |
1463 |
integer nlon, nlev |
1464 |
REAL PTE(NLON,NLEV), |
1465 |
* PVOL(NLON,NLEV), |
1466 |
* PVOM(NLON,NLEV), |
1467 |
* PULOW(NLON), |
1468 |
* PVLOW(NLON) |
1469 |
REAL PUM1(NLON,NLEV), PVM1(NLON,NLEV), PTM1(NLON,NLEV) |
1470 |
real, intent(in):: PLAT(NLON) |
1471 |
real PMEA(NLON) |
1472 |
real, intent(in):: PVAROR(NLON) |
1473 |
real ppic(NLON), PGEOM1(NLON,NLEV), PAPHM1(NLON,NLEV+1) |
1474 |
C |
1475 |
INTEGER KTEST(NLON) |
1476 |
real, intent(in):: ptsphy |
1477 |
C----------------------------------------------------------------------- |
1478 |
C |
1479 |
C* 0.2 LOCAL ARRAYS |
1480 |
C ------------ |
1481 |
logical lifthigh |
1482 |
integer klevm1, jl, ilevh, jk |
1483 |
real zcons1, ztmst, zrtmst,zpi, zhgeo |
1484 |
real zdelp, zslow, zsqua, zscav, zbet |
1485 |
INTEGER |
1486 |
* IKNUB(klon), |
1487 |
* IKNUL(klon) |
1488 |
LOGICAL LL1(KLON,KLEV+1) |
1489 |
C |
1490 |
REAL ZTAU(KLON,KLEV+1), |
1491 |
* ZTAV(KLON,KLEV+1), |
1492 |
* ZRHO(KLON,KLEV+1) |
1493 |
REAL ZDUDT(KLON), |
1494 |
* ZDVDT(KLON) |
1495 |
REAL ZHCRIT(KLON,KLEV) |
1496 |
C----------------------------------------------------------------------- |
1497 |
C |
1498 |
C* 1.1 INITIALIZATIONS |
1499 |
C --------------- |
1500 |
|
1501 |
LIFTHIGH=.FALSE. |
1502 |
|
1503 |
IF(NLON.NE.KLON.OR.NLEV.NE.KLEV)STOP |
1504 |
ZCONS1=1./RD |
1505 |
KLEVM1=KLEV-1 |
1506 |
ZTMST=PTSPHY |
1507 |
ZRTMST=1./ZTMST |
1508 |
ZPI=ACOS(-1.) |
1509 |
C |
1510 |
DO 1001 JL=1,klon |
1511 |
ZRHO(JL,KLEV+1) =0.0 |
1512 |
PULOW(JL) =0.0 |
1513 |
PVLOW(JL) =0.0 |
1514 |
iknub(JL) =klev |
1515 |
iknul(JL) =klev |
1516 |
ilevh=klev/3 |
1517 |
ll1(jl,klev+1)=.false. |
1518 |
DO 1000 JK=1,KLEV |
1519 |
PVOM(JL,JK)=0.0 |
1520 |
PVOL(JL,JK)=0.0 |
1521 |
PTE (JL,JK)=0.0 |
1522 |
1000 CONTINUE |
1523 |
1001 CONTINUE |
1524 |
|
1525 |
C |
1526 |
C* 2.1 DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF |
1527 |
C* LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE |
1528 |
C* THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS. |
1529 |
C |
1530 |
C |
1531 |
C |
1532 |
DO 2006 JK=KLEV,1,-1 |
1533 |
DO 2007 JL=1,klon |
1534 |
IF(KTEST(JL).EQ.1) THEN |
1535 |
ZHCRIT(JL,JK)=amax1(Ppic(JL)-pmea(JL),100.) |
1536 |
ZHGEO=PGEOM1(JL,JK)/RG |
1537 |
ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK)) |
1538 |
IF(ll1(JL,JK).neqv.ll1(JL,JK+1)) THEN |
1539 |
iknub(JL)=JK |
1540 |
ENDIF |
1541 |
ENDIF |
1542 |
2007 CONTINUE |
1543 |
2006 CONTINUE |
1544 |
C |
1545 |
do 2010 jl=1,klon |
1546 |
IF(KTEST(JL).EQ.1) THEN |
1547 |
iknub(jl)=max(iknub(jl),klev/2) |
1548 |
iknul(jl)=max(iknul(jl),2*klev/3) |
1549 |
if(iknub(jl).gt.nktopg) iknub(jl)=nktopg |
1550 |
if(iknub(jl).eq.nktopg) iknul(jl)=klev |
1551 |
if(iknub(jl).eq.iknul(jl)) iknub(jl)=iknul(jl)-1 |
1552 |
ENDIF |
1553 |
2010 continue |
1554 |
|
1555 |
C do 2011 jl=1,klon |
1556 |
C IF(KTEST(JL).EQ.1) THEN |
1557 |
C print *,' iknul= ',iknul(jl),' iknub=',iknub(jl) |
1558 |
C ENDIF |
1559 |
C2011 continue |
1560 |
|
1561 |
C PRINT *,' DANS OROLIFT: 2010' |
1562 |
|
1563 |
DO 223 JK=KLEV,2,-1 |
1564 |
DO 222 JL=1,klon |
1565 |
ZRHO(JL,JK)=2.*PAPHM1(JL,JK)*ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1)) |
1566 |
222 CONTINUE |
1567 |
223 CONTINUE |
1568 |
C PRINT *,' DANS OROLIFT: 223' |
1569 |
|
1570 |
C******************************************************************** |
1571 |
C |
1572 |
C* DEFINE LOW LEVEL FLOW |
1573 |
C ------------------- |
1574 |
DO 2115 JK=klev,1,-1 |
1575 |
DO 2116 JL=1,klon |
1576 |
IF(KTEST(JL).EQ.1) THEN |
1577 |
if(jk.ge.iknub(jl).and.jk.le.iknul(jl)) then |
1578 |
pulow(JL)=pulow(JL)+PUM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK)) |
1579 |
pvlow(JL)=pvlow(JL)+PVM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK)) |
1580 |
zrho(JL,klev+1)=zrho(JL,klev+1) |
1581 |
* +zrho(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK)) |
1582 |
end if |
1583 |
ENDIF |
1584 |
2116 CONTINUE |
1585 |
2115 CONTINUE |
1586 |
DO 2110 JL=1,klon |
1587 |
IF(KTEST(JL).EQ.1) THEN |
1588 |
pulow(JL)=pulow(JL)/(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl))) |
1589 |
pvlow(JL)=pvlow(JL)/(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl))) |
1590 |
zrho(JL,klev+1)=zrho(JL,klev+1) |
1591 |
* /(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl))) |
1592 |
ENDIF |
1593 |
2110 CONTINUE |
1594 |
|
1595 |
|
1596 |
200 CONTINUE |
1597 |
|
1598 |
C*********************************************************** |
1599 |
C |
1600 |
C* 3. COMPUTE MOUNTAIN LIFT |
1601 |
C |
1602 |
300 CONTINUE |
1603 |
C |
1604 |
DO 301 JL=1,klon |
1605 |
IF(KTEST(JL).EQ.1) THEN |
1606 |
ZTAU(JL,KLEV+1)= - GKLIFT*ZRHO(JL,KLEV+1)*2.*ROMEGA* |
1607 |
* 2*PVAROR(JL)* |
1608 |
* SIN(ZPI/180.*PLAT(JL))*PVLOW(JL) |
1609 |
ZTAV(JL,KLEV+1)= GKLIFT*ZRHO(JL,KLEV+1)*2.*ROMEGA* |
1610 |
* 2*PVAROR(JL)* |
1611 |
* SIN(ZPI/180.*PLAT(JL))*PULOW(JL) |
1612 |
ELSE |
1613 |
ZTAU(JL,KLEV+1)=0.0 |
1614 |
ZTAV(JL,KLEV+1)=0.0 |
1615 |
ENDIF |
1616 |
301 CONTINUE |
1617 |
|
1618 |
C |
1619 |
C* 4. COMPUTE LIFT PROFILE |
1620 |
C* -------------------- |
1621 |
C |
1622 |
|
1623 |
400 CONTINUE |
1624 |
|
1625 |
DO 401 JK=1,KLEV |
1626 |
DO 401 JL=1,klon |
1627 |
IF(KTEST(JL).EQ.1) THEN |
1628 |
ZTAU(JL,JK)=ZTAU(JL,KLEV+1)*PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1) |
1629 |
ZTAV(JL,JK)=ZTAV(JL,KLEV+1)*PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1) |
1630 |
ELSE |
1631 |
ZTAU(JL,JK)=0.0 |
1632 |
ZTAV(JL,JK)=0.0 |
1633 |
ENDIF |
1634 |
401 CONTINUE |
1635 |
C |
1636 |
C |
1637 |
C* 5. COMPUTE TENDENCIES. |
1638 |
C* ------------------- |
1639 |
IF(LIFTHIGH)THEN |
1640 |
C |
1641 |
500 CONTINUE |
1642 |
C PRINT *,' DANS OROLIFT: 500' |
1643 |
C |
1644 |
C EXPLICIT SOLUTION AT ALL LEVELS |
1645 |
C |
1646 |
DO 524 JK=1,klev |
1647 |
DO 523 JL=1,KLON |
1648 |
IF(KTEST(JL).EQ.1) THEN |
1649 |
ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK) |
1650 |
ZDUDT(JL)=-RG*(ZTAU(JL,JK+1)-ZTAU(JL,JK))/ZDELP |
1651 |
ZDVDT(JL)=-RG*(ZTAV(JL,JK+1)-ZTAV(JL,JK))/ZDELP |
1652 |
ENDIF |
1653 |
523 CONTINUE |
1654 |
524 CONTINUE |
1655 |
C |
1656 |
C PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY |
1657 |
C |
1658 |
DO 530 JK=1,klev |
1659 |
DO 530 JL=1,KLON |
1660 |
IF(KTEST(JL).EQ.1) THEN |
1661 |
|
1662 |
ZSLOW=SQRT(PULOW(JL)**2+PVLOW(JL)**2) |
1663 |
ZSQUA=AMAX1(SQRT(PUM1(JL,JK)**2+PVM1(JL,JK)**2),GVSEC) |
1664 |
ZSCAV=-ZDUDT(JL)*PVM1(JL,JK)+ZDVDT(JL)*PUM1(JL,JK) |
1665 |
IF(ZSQUA.GT.GVSEC)THEN |
1666 |
PVOM(JL,JK)=-ZSCAV*PVM1(JL,JK)/ZSQUA**2 |
1667 |
PVOL(JL,JK)= ZSCAV*PUM1(JL,JK)/ZSQUA**2 |
1668 |
ELSE |
1669 |
PVOM(JL,JK)=0.0 |
1670 |
PVOL(JL,JK)=0.0 |
1671 |
ENDIF |
1672 |
ZSQUA=SQRT(PUM1(JL,JK)**2+PUM1(JL,JK)**2) |
1673 |
IF(ZSQUA.LT.ZSLOW)THEN |
1674 |
PVOM(JL,JK)=ZSQUA/ZSLOW*PVOM(JL,JK) |
1675 |
PVOL(JL,JK)=ZSQUA/ZSLOW*PVOL(JL,JK) |
1676 |
ENDIF |
1677 |
|
1678 |
ENDIF |
1679 |
530 CONTINUE |
1680 |
C |
1681 |
C 6. LOW LEVEL LIFT, SEMI IMPLICIT: |
1682 |
C ---------------------------------- |
1683 |
|
1684 |
ELSE |
1685 |
|
1686 |
DO 601 JL=1,KLON |
1687 |
IF(KTEST(JL).EQ.1) THEN |
1688 |
DO JK=KLEV,IKNUB(JL),-1 |
1689 |
ZBET=GKLIFT*2.*ROMEGA*SIN(ZPI/180.*PLAT(JL))*ztmst* |
1690 |
* (PGEOM1(JL,IKNUB(JL)-1)-PGEOM1(JL, JK))/ |
1691 |
* (PGEOM1(JL,IKNUB(JL)-1)-PGEOM1(JL,KLEV)) |
1692 |
ZDUDT(JL)=-PUM1(JL,JK)/ztmst/(1+ZBET**2) |
1693 |
ZDVDT(JL)=-PVM1(JL,JK)/ztmst/(1+ZBET**2) |
1694 |
PVOM(JL,JK)= ZBET**2*ZDUDT(JL) - ZBET *ZDVDT(JL) |
1695 |
PVOL(JL,JK)= ZBET *ZDUDT(JL) + ZBET**2*ZDVDT(JL) |
1696 |
ENDDO |
1697 |
ENDIF |
1698 |
601 CONTINUE |
1699 |
|
1700 |
ENDIF |
1701 |
|
1702 |
RETURN |
1703 |
END |
1704 |
SUBROUTINE SUGWD(NLON,NLEV,paprs,pplay) |
1705 |
C |
1706 |
C**** *SUGWD* INITIALIZE COMMON YOEGWD CONTROLLING GRAVITY WAVE DRAG |
1707 |
C |
1708 |
C PURPOSE. |
1709 |
C -------- |
1710 |
C INITIALIZE YOEGWD, THE COMMON THAT CONTROLS THE |
1711 |
C GRAVITY WAVE DRAG PARAMETRIZATION. |
1712 |
C |
1713 |
C** INTERFACE. |
1714 |
C ---------- |
1715 |
C CALL *SUGWD* FROM *SUPHEC* |
1716 |
C ----- ------ |
1717 |
C |
1718 |
C EXPLICIT ARGUMENTS : |
1719 |
C -------------------- |
1720 |
C PSIG : VERTICAL COORDINATE TABLE |
1721 |
C NLEV : NUMBER OF MODEL LEVELS |
1722 |
C |
1723 |
C IMPLICIT ARGUMENTS : |
1724 |
C -------------------- |
1725 |
C COMMON YOEGWD |
1726 |
C |
1727 |
C METHOD. |
1728 |
C ------- |
1729 |
C SEE DOCUMENTATION |
1730 |
C |
1731 |
C EXTERNALS. |
1732 |
C ---------- |
1733 |
C NONE |
1734 |
C |
1735 |
C REFERENCE. |
1736 |
C ---------- |
1737 |
C ECMWF Research Department documentation of the IFS |
1738 |
C |
1739 |
C AUTHOR. |
1740 |
C ------- |
1741 |
C MARTIN MILLER *ECMWF* |
1742 |
C |
1743 |
C MODIFICATIONS. |
1744 |
C -------------- |
1745 |
C ORIGINAL : 90-01-01 |
1746 |
C ------------------------------------------------------------------ |
1747 |
use yoegwd |
1748 |
implicit none |
1749 |
C |
1750 |
C ----------------------------------------------------------------- |
1751 |
C ---------------------------------------------------------------- |
1752 |
C |
1753 |
integer nlon,nlev, jk |
1754 |
REAL, intent(in):: paprs(nlon,nlev+1) |
1755 |
REAL, intent(in):: pplay(nlon,nlev) |
1756 |
real zpr,zstra,zsigt,zpm1r |
1757 |
C |
1758 |
C* 1. SET THE VALUES OF THE PARAMETERS |
1759 |
C -------------------------------- |
1760 |
C |
1761 |
100 CONTINUE |
1762 |
C |
1763 |
PRINT *,' DANS SUGWD NLEV=',NLEV |
1764 |
GHMAX=10000. |
1765 |
C |
1766 |
ZPR=100000. |
1767 |
ZSTRA=0.1 |
1768 |
ZSIGT=0.94 |
1769 |
cold ZPR=80000. |
1770 |
cold ZSIGT=0.85 |
1771 |
C |
1772 |
DO 110 JK=1,NLEV |
1773 |
ZPM1R=pplay(nlon/2,jk)/paprs(nlon/2,1) |
1774 |
IF(ZPM1R.GE.ZSIGT)THEN |
1775 |
nktopg=JK |
1776 |
ENDIF |
1777 |
ZPM1R=pplay(nlon/2,jk)/paprs(nlon/2,1) |
1778 |
IF(ZPM1R.GE.ZSTRA)THEN |
1779 |
NSTRA=JK |
1780 |
ENDIF |
1781 |
110 CONTINUE |
1782 |
c |
1783 |
c inversion car dans orodrag on compte les niveaux a l'envers |
1784 |
nktopg=nlev-nktopg+1 |
1785 |
nstra=nlev-nstra |
1786 |
print *,' DANS SUGWD nktopg=', nktopg |
1787 |
print *,' DANS SUGWD nstra=', nstra |
1788 |
C |
1789 |
GSIGCR=0.80 |
1790 |
C |
1791 |
GKDRAG=0.2 |
1792 |
GRAHILO=1. |
1793 |
GRCRIT=0.01 |
1794 |
GFRCRIT=1.0 |
1795 |
GKWAKE=0.50 |
1796 |
C |
1797 |
GKLIFT=0.50 |
1798 |
GVCRIT =0.0 |
1799 |
C |
1800 |
C |
1801 |
C ---------------------------------------------------------------- |
1802 |
C |
1803 |
C* 2. SET VALUES OF SECURITY PARAMETERS |
1804 |
C --------------------------------- |
1805 |
C |
1806 |
200 CONTINUE |
1807 |
C |
1808 |
GVSEC=0.10 |
1809 |
GSSEC=1.E-12 |
1810 |
C |
1811 |
GTSEC=1.E-07 |
1812 |
C |
1813 |
C ---------------------------------------------------------------- |
1814 |
C |
1815 |
RETURN |
1816 |
END |