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