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

Last change on this file since 5116 was 5116, checked in by abarral, 4 months ago

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