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