source: LMDZ6/branches/Amaury_dev/libf/phylmd/orografi.F90 @ 5144

Last change on this file since 5144 was 5144, checked in by abarral, 3 months ago

Put YOMCST.h into modules

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