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

Last change on this file since 5447 was 5309, checked in by abarral, 7 weeks ago

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