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

Last change on this file since 5878 was 5768, checked in by rkazeroni, 5 months ago

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