source: LMDZ6/trunk/libf/phylmd/orografi.f90 @ 5300

Last change on this file since 5300 was 5285, checked in by abarral, 5 weeks ago

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