source: LMDZ6/trunk/libf/phylmd/orografi.F90 @ 5006

Last change on this file since 5006 was 2357, checked in by lguez, 9 years ago

The default value for gwd_rando_ruwmax is changed. This is ok because
the Lott GWD rando parameterization has not been much used yet.

New parameters GWD_FRONT_RUWMAX, GWD_FRONT_SAT, that can be chosen at
run-time, for frontal gravity waves.

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