source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/orografi.F90 @ 5440

Last change on this file since 5440 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.
  • 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: 44.0 KB
Line 
1
2! $Id: orografi.F90 1992 2014-03-05 13:19:12Z evignon $
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  ! ym#include "dimensions.h"
27  ! ym#include "dimphy.h"
28  include "YOMCST.h"
29
30  ! ARGUMENTS
31
32  INTEGER nlon, nlev
33  REAL dtime
34  REAL paprs(klon, klev+1)
35  REAL pplay(klon, klev)
36  REAL pmea(nlon), pstd(nlon), psig(nlon), pgam(nlon), pthe(nlon)
37  REAL ppic(nlon), pval(nlon)
38  REAL pulow(nlon), pvlow(nlon), pustr(nlon), pvstr(nlon)
39  REAL t(nlon, nlev), u(nlon, nlev), v(nlon, nlev)
40  REAL d_t(nlon, nlev), d_u(nlon, nlev), d_v(nlon, nlev)
41
42  INTEGER i, k, kgwd, kdx(nlon), ktest(nlon)
43
44  ! Variables locales:
45
46  REAL zgeom(klon, klev)
47  REAL pdtdt(klon, klev), pdudt(klon, klev), pdvdt(klon, klev)
48  REAL pt(klon, klev), pu(klon, klev), pv(klon, klev)
49  REAL papmf(klon, klev), papmh(klon, klev+1)
50
51  ! initialiser les variables de sortie (pour securite)
52
53  DO i = 1, klon
54    pulow(i) = 0.0
55    pvlow(i) = 0.0
56    pustr(i) = 0.0
57    pvstr(i) = 0.0
58  END DO
59  DO k = 1, klev
60    DO i = 1, klon
61      d_t(i, k) = 0.0
62      d_u(i, k) = 0.0
63      d_v(i, k) = 0.0
64      pdudt(i, k) = 0.0
65      pdvdt(i, k) = 0.0
66      pdtdt(i, k) = 0.0
67    END DO
68  END DO
69
70  ! preparer les variables d'entree (attention: l'ordre des niveaux
71  ! verticaux augmente du haut vers le bas)
72
73  DO k = 1, klev
74    DO i = 1, klon
75      pt(i, k) = t(i, klev-k+1)
76      pu(i, k) = u(i, klev-k+1)
77      pv(i, k) = v(i, klev-k+1)
78      papmf(i, k) = pplay(i, klev-k+1)
79    END DO
80  END DO
81  DO k = 1, klev + 1
82    DO i = 1, klon
83      papmh(i, k) = paprs(i, klev-k+2)
84    END DO
85  END DO
86  DO i = 1, klon
87    zgeom(i, klev) = rd*pt(i, klev)*log(papmh(i,klev+1)/papmf(i,klev))
88  END DO
89  DO k = klev - 1, 1, -1
90    DO i = 1, klon
91      zgeom(i, k) = zgeom(i, k+1) + rd*(pt(i,k)+pt(i,k+1))/2.0*log(papmf(i,k+ &
92        1)/papmf(i,k))
93    END DO
94  END DO
95
96  ! appeler la routine principale
97
98  CALL orodrag(klon, klev, kgwd, kdx, ktest, dtime, papmh, papmf, zgeom, pt, &
99    pu, pv, pmea, pstd, psig, pgam, pthe, ppic, pval, pulow, pvlow, pdudt, &
100    pdvdt, pdtdt)
101
102  DO k = 1, klev
103    DO i = 1, klon
104      d_u(i, klev+1-k) = dtime*pdudt(i, k)
105      d_v(i, klev+1-k) = dtime*pdvdt(i, k)
106      d_t(i, klev+1-k) = dtime*pdtdt(i, k)
107      pustr(i) = pustr(i) &        ! IM BUG  .
108                                   ! +rg*pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))
109        +pdudt(i, k)*(papmh(i,k+1)-papmh(i,k))/rg
110      pvstr(i) = pvstr(i) &        ! IM BUG  .
111                                   ! +rg*pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))
112        +pdvdt(i, k)*(papmh(i,k+1)-papmh(i,k))/rg
113    END DO
114  END DO
115
116  RETURN
117END SUBROUTINE drag_noro
118SUBROUTINE orodrag(nlon, nlev, kgwd, kdx, ktest, ptsphy, paphm1, papm1, &
119    pgeom1, ptm1, pum1, pvm1, pmea, pstd, psig, pgamma, ptheta, ppic, pval &
120  ! outputs
121    , pulow, pvlow, pvom, pvol, pte)
122
123  USE dimphy
124  IMPLICIT NONE
125
126
127
128  ! **** *gwdrag* - does the gravity wave parametrization.
129
130  ! purpose.
131  ! --------
132
133  ! this routine computes the physical tendencies of the
134  ! prognostic variables u,v  and t due to  vertical transports by
135  ! subgridscale orographically excited gravity waves
136
137  ! **   interface.
138  ! ----------
139  ! called from *callpar*.
140
141  ! the routine takes its input from the long-term storage:
142  ! u,v,t and p at t-1.
143
144  ! explicit arguments :
145  ! --------------------
146  ! ==== inputs ===
147  ! ==== outputs ===
148
149  ! implicit arguments :   none
150  ! --------------------
151
152  ! implicit logical (l)
153
154  ! method.
155  ! -------
156
157  ! externals.
158  ! ----------
159  INTEGER ismin, ismax
160  EXTERNAL ismin, ismax
161
162  ! reference.
163  ! ----------
164
165  ! author.
166  ! -------
167  ! m.miller + b.ritter   e.c.m.w.f.     15/06/86.
168
169  ! f.lott + m. miller    e.c.m.w.f.     22/11/94
170  ! -----------------------------------------------------------------------
171
172
173  ! ym#include "dimensions.h"
174  ! ym#include "dimphy.h"
175  include "YOMCST.h"
176  include "YOEGWD.h"
177  ! -----------------------------------------------------------------------
178
179  ! *       0.1   arguments
180  ! ---------
181
182
183  ! ym      integer nlon, nlev, klevm1
184  INTEGER nlon, nlev
185  INTEGER kgwd, jl, ilevp1, jk, ji
186  REAL zdelp, ztemp, zforc, ztend
187  REAL rover, zb, zc, zconb, zabsv
188  REAL zzd1, ratio, zbet, zust, zvst, zdis
189  REAL pte(nlon, nlev), pvol(nlon, nlev), pvom(nlon, nlev), pulow(klon), &
190    pvlow(klon)
191  REAL pum1(nlon, nlev), pvm1(nlon, nlev), ptm1(nlon, nlev), pmea(nlon), &
192    pstd(nlon), psig(nlon), pgamma(nlon), ptheta(nlon), ppic(nlon), &
193    pval(nlon), pgeom1(nlon, nlev), papm1(nlon, nlev), paphm1(nlon, nlev+1)
194
195  INTEGER kdx(nlon), ktest(nlon)
196  ! -----------------------------------------------------------------------
197
198  ! *       0.2   local arrays
199  ! ------------
200  INTEGER isect(klon), icrit(klon), ikcrith(klon), ikenvh(klon), iknu(klon), &
201    iknu2(klon), ikcrit(klon), ikhlim(klon)
202
203  REAL ztau(klon, klev+1), ztauf(klon, klev+1), zstab(klon, klev+1), &
204    zvph(klon, klev+1), zrho(klon, klev+1), zri(klon, klev+1), &
205    zpsi(klon, klev+1), zzdep(klon, klev)
206  REAL zdudt(klon), zdvdt(klon), zdtdt(klon), zdedt(klon), zvidis(klon), &
207    znu(klon), zd1(klon), zd2(klon), zdmod(klon)
208  REAL ztmst, ptsphy, zrtmst
209
210  ! ------------------------------------------------------------------
211
212  ! *         1.    initialization
213  ! --------------
214
215
216  ! ------------------------------------------------------------------
217
218  ! *         1.1   computational constants
219  ! -----------------------
220
221
222  ! ztmst=twodt
223  ! if(nstep.eq.nstart) ztmst=0.5*twodt
224  ! ym      klevm1=klev-1
225  ztmst = ptsphy
226  zrtmst = 1./ztmst
227
228  ! ------------------------------------------------------------------
229
230  ! *         1.3   check whether row contains point for printing
231  ! ---------------------------------------------
232
233
234  ! ------------------------------------------------------------------
235
236  ! *         2.     precompute basic state variables.
237  ! *                ---------- ----- ----- ----------
238  ! *                define low level wind, project winds in plane of
239  ! *                low level wind, determine sector in which to take
240  ! *                the variance and set indicator for critical levels.
241
242
243
244
245  CALL orosetup(nlon, ktest, ikcrit, ikcrith, icrit, ikenvh, iknu, iknu2, &
246    paphm1, papm1, pum1, pvm1, ptm1, pgeom1, pstd, zrho, zri, zstab, ztau, &
247    zvph, zpsi, zzdep, pulow, pvlow, ptheta, pgamma, pmea, ppic, pval, znu, &
248    zd1, zd2, zdmod)
249
250  ! ***********************************************************
251
252
253  ! *         3.      compute low level stresses using subcritical and
254  ! *                 supercritical forms.computes anisotropy coefficient
255  ! *                 as measure of orographic twodimensionality.
256
257
258  CALL gwstress(nlon, nlev, ktest, icrit, ikenvh, iknu, zrho, zstab, zvph, &
259    pstd, psig, pmea, ppic, ztau, pgeom1, zdmod)
260
261  ! *         4.      compute stress profile.
262  ! *                 ------- ------ --------
263
264
265
266  CALL gwprofil(nlon, nlev, kgwd, kdx, ktest, ikcrith, icrit, paphm1, zrho, &
267    zstab, zvph, zri, ztau, zdmod, psig, pstd)
268
269  ! *         5.      compute tendencies.
270  ! *                 -------------------
271
272
273  ! explicit solution at all levels for the gravity wave
274  ! implicit solution for the blocked levels
275
276  DO jl = kidia, kfdia
277    zvidis(jl) = 0.0
278    zdudt(jl) = 0.0
279    zdvdt(jl) = 0.0
280    zdtdt(jl) = 0.0
281  END DO
282
283  ilevp1 = klev + 1
284
285
286  DO jk = 1, klev
287
288
289    ! do 523 jl=1,kgwd
290    ! ji=kdx(jl)
291    ! Modif vectorisation 02/04/2004
292    DO ji = kidia, kfdia
293      IF (ktest(ji)==1) THEN
294
295        zdelp = paphm1(ji, jk+1) - paphm1(ji, jk)
296        ztemp = -rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,ilevp1)*zdelp)
297        zdudt(ji) = (pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)
298        zdvdt(ji) = (pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)
299
300        ! controle des overshoots:
301
302        zforc = sqrt(zdudt(ji)**2+zdvdt(ji)**2) + 1.E-12
303        ztend = sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/ztmst + 1.E-12
304        rover = 0.25
305        IF (zforc>=rover*ztend) THEN
306          zdudt(ji) = rover*ztend/zforc*zdudt(ji)
307          zdvdt(ji) = rover*ztend/zforc*zdvdt(ji)
308        END IF
309
310        ! fin du controle des overshoots
311
312        IF (jk>=ikenvh(ji)) THEN
313          zb = 1.0 - 0.18*pgamma(ji) - 0.04*pgamma(ji)**2
314          zc = 0.48*pgamma(ji) + 0.3*pgamma(ji)**2
315          zconb = 2.*ztmst*gkwake*psig(ji)/(4.*pstd(ji))
316          zabsv = sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/2.
317          zzd1 = zb*cos(zpsi(ji,jk))**2 + zc*sin(zpsi(ji,jk))**2
318          ratio = (cos(zpsi(ji,jk))**2+pgamma(ji)*sin(zpsi(ji, &
319            jk))**2)/(pgamma(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2)
320          zbet = max(0., 2.-1./ratio)*zconb*zzdep(ji, jk)*zzd1*zabsv
321
322          ! simplement oppose au vent
323
324          zdudt(ji) = -pum1(ji, jk)/ztmst
325          zdvdt(ji) = -pvm1(ji, jk)/ztmst
326
327          ! projection dans la direction de l'axe principal de l'orographie
328          ! mod     zdudt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.)
329          ! mod *              +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.))
330          ! mod *              *cos(ptheta(ji)*rpi/180.)/ztmst
331          ! mod     zdvdt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.)
332          ! mod *              +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.))
333          ! mod *              *sin(ptheta(ji)*rpi/180.)/ztmst
334          zdudt(ji) = zdudt(ji)*(zbet/(1.+zbet))
335          zdvdt(ji) = zdvdt(ji)*(zbet/(1.+zbet))
336        END IF
337        pvom(ji, jk) = zdudt(ji)
338        pvol(ji, jk) = zdvdt(ji)
339        zust = pum1(ji, jk) + ztmst*zdudt(ji)
340        zvst = pvm1(ji, jk) + ztmst*zdvdt(ji)
341        zdis = 0.5*(pum1(ji,jk)**2+pvm1(ji,jk)**2-zust**2-zvst**2)
342        zdedt(ji) = zdis/ztmst
343        zvidis(ji) = zvidis(ji) + zdis*zdelp
344        zdtdt(ji) = zdedt(ji)/rcpd
345        ! pte(ji,jk)=zdtdt(ji)
346
347        ! ENCORE UN TRUC POUR EVITER LES EXPLOSIONS
348
349        pte(ji, jk) = 0.0
350
351      END IF
352    END DO
353
354  END DO
355
356
357  RETURN
358END SUBROUTINE orodrag
359SUBROUTINE orosetup(nlon, ktest, kkcrit, kkcrith, kcrit, kkenvh, kknu, kknu2, &
360    paphm1, papm1, pum1, pvm1, ptm1, pgeom1, pstd, prho, pri, pstab, ptau, &
361    pvph, ppsi, pzdep, pulow, pvlow, ptheta, pgamma, pmea, ppic, pval, pnu, &
362    pd1, pd2, pdmod)
363
364  ! **** *gwsetup*
365
366  ! purpose.
367  ! --------
368
369  ! **   interface.
370  ! ----------
371  ! from *orodrag*
372
373  ! explicit arguments :
374  ! --------------------
375  ! ==== inputs ===
376  ! ==== outputs ===
377
378  ! implicit arguments :   none
379  ! --------------------
380
381  ! method.
382  ! -------
383
384
385  ! externals.
386  ! ----------
387
388
389  ! reference.
390  ! ----------
391
392  ! see ecmwf research department documentation of the "i.f.s."
393
394  ! author.
395  ! -------
396
397  ! modifications.
398  ! --------------
399  ! f.lott  for the new-gwdrag scheme november 1993
400
401  ! -----------------------------------------------------------------------
402  USE dimphy
403  IMPLICIT NONE
404
405
406  ! ym#include "dimensions.h"
407  ! ym#include "dimphy.h"
408  include "YOMCST.h"
409  include "YOEGWD.h"
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  ! -----------------------------------------------------------------------
849  USE dimphy
850  IMPLICIT NONE
851  ! ym#include "dimensions.h"
852  ! ym#include "dimphy.h"
853  include "YOMCST.h"
854  include "YOEGWD.h"
855
856  ! -----------------------------------------------------------------------
857
858  ! *       0.1   arguments
859  ! ---------
860
861  INTEGER nlon, nlev
862  INTEGER kcrit(nlon), ktest(nlon), kkenvh(nlon), kknu(nlon)
863
864  REAL prho(nlon, nlev+1), pstab(nlon, nlev+1), ptau(nlon, nlev+1), &
865    pvph(nlon, nlev+1), pgeom1(nlon, nlev), pstd(nlon)
866
867  REAL psig(nlon)
868  REAL pmea(nlon), ppic(nlon)
869  REAL pdmod(nlon)
870
871  ! -----------------------------------------------------------------------
872
873  ! *       0.2   local arrays
874  ! ------------
875  INTEGER jl
876  REAL zblock, zvar, zeff
877  LOGICAL lo
878
879  ! -----------------------------------------------------------------------
880
881  ! *       0.3   functions
882  ! ---------
883  ! ------------------------------------------------------------------
884
885  ! *         1.    initialization
886  ! --------------
887
888
889  ! *         3.1     gravity wave stress.
890
891
892
893  DO jl = kidia, kfdia
894    IF (ktest(jl)==1) THEN
895
896      ! effective mountain height above the blocked flow
897
898      IF (kkenvh(jl)==klev) THEN
899        zblock = 0.0
900      ELSE
901        zblock = (pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)+1))/2./rg
902      END IF
903
904      zvar = ppic(jl) - pmea(jl)
905      zeff = amax1(0., zvar-zblock)
906
907      ptau(jl, klev+1) = prho(jl, klev+1)*gkdrag*psig(jl)*zeff**2/4./ &
908        pstd(jl)*pvph(jl, klev+1)*pdmod(jl)*sqrt(pstab(jl,klev+1))
909
910      ! too small value of stress or  low level flow include critical level
911      ! or low level flow:  gravity wave stress nul.
912
913      lo = (ptau(jl,klev+1)<gtsec) .OR. (kcrit(jl)>=kknu(jl)) .OR. &
914        (pvph(jl,klev+1)<gvcrit)
915      ! if(lo) ptau(jl,klev+1)=0.0
916
917    ELSE
918
919      ptau(jl, klev+1) = 0.0
920
921    END IF
922
923  END DO
924
925  RETURN
926END SUBROUTINE gwstress
927SUBROUTINE gwprofil(nlon, nlev, kgwd, kdx, ktest, kkcrith, kcrit, paphm1, &
928    prho, pstab, pvph, pri, ptau, pdmod, psig, pvar)
929
930  ! **** *GWPROFIL*
931
932  ! PURPOSE.
933  ! --------
934
935  ! **   INTERFACE.
936  ! ----------
937  ! FROM *GWDRAG*
938
939  ! EXPLICIT ARGUMENTS :
940  ! --------------------
941  ! ==== INPUTS ===
942  ! ==== OUTPUTS ===
943
944  ! IMPLICIT ARGUMENTS :   NONE
945  ! --------------------
946
947  ! METHOD:
948  ! -------
949  ! THE STRESS PROFILE FOR GRAVITY WAVES IS COMPUTED AS FOLLOWS:
950  ! IT IS CONSTANT (NO GWD) AT THE LEVELS BETWEEN THE GROUND
951  ! AND THE TOP OF THE BLOCKED LAYER (KKENVH).
952  ! IT DECREASES LINEARLY WITH HEIGHTS FROM THE TOP OF THE
953  ! BLOCKED LAYER TO 3*VAROR (kKNU), TO SIMULATES LEE WAVES OR
954  ! NONLINEAR GRAVITY WAVE BREAKING.
955  ! ABOVE IT IS CONSTANT, EXCEPT WHEN THE WAVE ENCOUNTERS A CRITICAL
956  ! LEVEL (KCRIT) OR WHEN IT BREAKS.
957
958
959
960  ! EXTERNALS.
961  ! ----------
962
963
964  ! REFERENCE.
965  ! ----------
966
967  ! SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
968
969  ! AUTHOR.
970  ! -------
971
972  ! MODIFICATIONS.
973  ! --------------
974  ! PASSAGE OF THE NEW GWDRAG TO I.F.S. (F. LOTT, 22/11/93)
975  ! -----------------------------------------------------------------------
976  USE dimphy
977  IMPLICIT NONE
978
979
980
981
982  ! ym#include "dimensions.h"
983  ! ym#include "dimphy.h"
984  include "YOMCST.h"
985  include "YOEGWD.h"
986
987  ! -----------------------------------------------------------------------
988
989  ! *       0.1   ARGUMENTS
990  ! ---------
991
992  INTEGER nlon, nlev
993  INTEGER kkcrith(nlon), kcrit(nlon), kdx(nlon), ktest(nlon)
994
995
996  REAL paphm1(nlon, nlev+1), pstab(nlon, nlev+1), prho(nlon, nlev+1), &
997    pvph(nlon, nlev+1), pri(nlon, nlev+1), ptau(nlon, nlev+1)
998
999  REAL pdmod(nlon), psig(nlon), pvar(nlon)
1000
1001  ! -----------------------------------------------------------------------
1002
1003  ! *       0.2   LOCAL ARRAYS
1004  ! ------------
1005
1006  INTEGER ilevh, ji, kgwd, jl, jk
1007  REAL zsqr, zalfa, zriw, zdel, zb, zalpha, zdz2n
1008  REAL zdelp, zdelpt
1009  REAL zdz2(klon, klev), znorm(klon), zoro(klon)
1010  REAL ztau(klon, klev+1)
1011
1012  ! -----------------------------------------------------------------------
1013
1014  ! *         1.    INITIALIZATION
1015  ! --------------
1016
1017  ! print *,' entree gwprofil'
1018
1019
1020  ! *    COMPUTATIONAL CONSTANTS.
1021  ! ------------- ----------
1022
1023  ilevh = klev/3
1024
1025  ! DO 400 ji=1,kgwd
1026  ! jl=kdx(ji)
1027  ! Modif vectorisation 02/04/2004
1028  DO jl = kidia, kfdia
1029    IF (ktest(jl)==1) THEN
1030      zoro(jl) = psig(jl)*pdmod(jl)/4./max(pvar(jl), 1.0)
1031      ztau(jl, klev+1) = ptau(jl, klev+1)
1032    END IF
1033  END DO
1034
1035
1036  DO jk = klev, 2, -1
1037
1038    ! *         4.1    CONSTANT WAVE STRESS UNTIL TOP OF THE
1039    ! BLOCKING LAYER.
1040
1041    ! DO 411 ji=1,kgwd
1042    ! jl=kdx(ji)
1043    ! Modif vectorisation 02/04/2004
1044    DO jl = kidia, kfdia
1045      IF (ktest(jl)==1) THEN
1046        IF (jk>kkcrith(jl)) THEN
1047          ptau(jl, jk) = ztau(jl, klev+1)
1048          ! ENDIF
1049          ! IF(JK.EQ.KKCRITH(JL)) THEN
1050        ELSE
1051          ptau(jl, jk) = grahilo*ztau(jl, klev+1)
1052        END IF
1053      END IF
1054    END DO
1055
1056    ! *         4.15   CONSTANT SHEAR STRESS UNTIL THE TOP OF THE
1057    ! LOW LEVEL FLOW LAYER.
1058
1059
1060    ! *         4.2    WAVE DISPLACEMENT AT NEXT LEVEL.
1061
1062
1063    ! DO 421 ji=1,kgwd
1064    ! jl=kdx(ji)
1065    ! Modif vectorisation 02/04/2004
1066    DO jl = kidia, kfdia
1067      IF (ktest(jl)==1) THEN
1068        IF (jk<kkcrith(jl)) THEN
1069          znorm(jl) = gkdrag*prho(jl, jk)*sqrt(pstab(jl,jk))*pvph(jl, jk)* &
1070            zoro(jl)
1071          zdz2(jl, jk) = ptau(jl, jk+1)/max(znorm(jl), gssec)
1072        END IF
1073      END IF
1074    END DO
1075
1076    ! *         4.3    WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT
1077    ! *                AND STRESS:  BREAKING EVALUATION AND CRITICAL
1078    ! LEVEL
1079
1080
1081    ! DO 431 ji=1,kgwd
1082    ! jl=Kdx(ji)
1083    ! Modif vectorisation 02/04/2004
1084    DO jl = kidia, kfdia
1085      IF (ktest(jl)==1) THEN
1086
1087        IF (jk<kkcrith(jl)) THEN
1088          IF ((ptau(jl,jk+1)<gtsec) .OR. (jk<=kcrit(jl))) THEN
1089            ptau(jl, jk) = 0.0
1090          ELSE
1091            zsqr = sqrt(pri(jl,jk))
1092            zalfa = sqrt(pstab(jl,jk)*zdz2(jl,jk))/pvph(jl, jk)
1093            zriw = pri(jl, jk)*(1.-zalfa)/(1+zalfa*zsqr)**2
1094            IF (zriw<grcrit) THEN
1095              zdel = 4./zsqr/grcrit + 1./grcrit**2 + 4./grcrit
1096              zb = 1./grcrit + 2./zsqr
1097              zalpha = 0.5*(-zb+sqrt(zdel))
1098              zdz2n = (pvph(jl,jk)*zalpha)**2/pstab(jl, jk)
1099              ptau(jl, jk) = znorm(jl)*zdz2n
1100            ELSE
1101              ptau(jl, jk) = znorm(jl)*zdz2(jl, jk)
1102            END IF
1103            ptau(jl, jk) = min(ptau(jl,jk), ptau(jl,jk+1))
1104          END IF
1105        END IF
1106      END IF
1107    END DO
1108
1109  END DO
1110
1111  ! REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL
1112
1113  ! DO 530 ji=1,kgwd
1114  ! jl=kdx(ji)
1115  ! Modif vectorisation 02/04/2004
1116  DO jl = kidia, kfdia
1117    IF (ktest(jl)==1) THEN
1118      ztau(jl, kkcrith(jl)) = ptau(jl, kkcrith(jl))
1119      ztau(jl, nstra) = ptau(jl, nstra)
1120    END IF
1121  END DO
1122
1123  DO jk = 1, klev
1124
1125    ! DO 532 ji=1,kgwd
1126    ! jl=kdx(ji)
1127    ! Modif vectorisation 02/04/2004
1128    DO jl = kidia, kfdia
1129      IF (ktest(jl)==1) THEN
1130
1131
1132        IF (jk>kkcrith(jl)) THEN
1133
1134          zdelp = paphm1(jl, jk) - paphm1(jl, klev+1)
1135          zdelpt = paphm1(jl, kkcrith(jl)) - paphm1(jl, klev+1)
1136          ptau(jl, jk) = ztau(jl, klev+1) + (ztau(jl,kkcrith(jl))-ztau(jl, &
1137            klev+1))*zdelp/zdelpt
1138
1139        END IF
1140
1141      END IF
1142    END DO
1143
1144    ! REORGANISATION IN THE STRATOSPHERE
1145
1146    ! DO 533 ji=1,kgwd
1147    ! jl=kdx(ji)
1148    ! Modif vectorisation 02/04/2004
1149    DO jl = kidia, kfdia
1150      IF (ktest(jl)==1) THEN
1151
1152
1153        IF (jk<nstra) THEN
1154
1155          zdelp = paphm1(jl, nstra)
1156          zdelpt = paphm1(jl, jk)
1157          ptau(jl, jk) = ztau(jl, nstra)*zdelpt/zdelp
1158
1159        END IF
1160
1161      END IF
1162    END DO
1163
1164    ! REORGANISATION IN THE TROPOSPHERE
1165
1166    ! DO 534 ji=1,kgwd
1167    ! jl=kdx(ji)
1168    ! Modif vectorisation 02/04/2004
1169    DO jl = kidia, kfdia
1170      IF (ktest(jl)==1) THEN
1171
1172
1173        IF (jk<kkcrith(jl) .AND. jk>nstra) THEN
1174
1175          zdelp = paphm1(jl, jk) - paphm1(jl, kkcrith(jl))
1176          zdelpt = paphm1(jl, nstra) - paphm1(jl, kkcrith(jl))
1177          ptau(jl, jk) = ztau(jl, kkcrith(jl)) + (ztau(jl,nstra)-ztau(jl, &
1178            kkcrith(jl)))*zdelp/zdelpt
1179
1180        END IF
1181      END IF
1182    END DO
1183
1184
1185  END DO
1186
1187
1188  RETURN
1189END SUBROUTINE gwprofil
1190SUBROUTINE lift_noro(nlon, nlev, dtime, paprs, pplay, plat, pmea, pstd, ppic, &
1191    ktest, t, u, v, pulow, pvlow, pustr, pvstr, d_t, d_u, d_v)
1192
1193  USE dimphy
1194  IMPLICIT NONE
1195  ! ======================================================================
1196  ! Auteur(s): F.Lott (LMD/CNRS) date: 19950201
1197  ! Objet: Frottement de la montagne Interface
1198  ! ======================================================================
1199  ! Arguments:
1200  ! dtime---input-R- pas d'integration (s)
1201  ! paprs---input-R-pression pour chaque inter-couche (en Pa)
1202  ! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
1203  ! t-------input-R-temperature (K)
1204  ! u-------input-R-vitesse horizontale (m/s)
1205  ! v-------input-R-vitesse horizontale (m/s)
1206
1207  ! d_t-----output-R-increment de la temperature
1208  ! d_u-----output-R-increment de la vitesse u
1209  ! d_v-----output-R-increment de la vitesse v
1210  ! ======================================================================
1211  ! ym#include "dimensions.h"
1212  ! ym#include "dimphy.h"
1213  include "YOMCST.h"
1214
1215  ! ARGUMENTS
1216
1217  INTEGER nlon, nlev
1218  REAL dtime
1219  REAL paprs(klon, klev+1)
1220  REAL pplay(klon, klev)
1221  REAL plat(nlon), pmea(nlon)
1222  REAL pstd(nlon)
1223  REAL ppic(nlon)
1224  REAL pulow(nlon), pvlow(nlon), pustr(nlon), pvstr(nlon)
1225  REAL t(nlon, nlev), u(nlon, nlev), v(nlon, nlev)
1226  REAL d_t(nlon, nlev), d_u(nlon, nlev), d_v(nlon, nlev)
1227
1228  INTEGER i, k, ktest(nlon)
1229
1230  ! Variables locales:
1231
1232  REAL zgeom(klon, klev)
1233  REAL pdtdt(klon, klev), pdudt(klon, klev), pdvdt(klon, klev)
1234  REAL pt(klon, klev), pu(klon, klev), pv(klon, klev)
1235  REAL papmf(klon, klev), papmh(klon, klev+1)
1236
1237  ! initialiser les variables de sortie (pour securite)
1238
1239  DO i = 1, klon
1240    pulow(i) = 0.0
1241    pvlow(i) = 0.0
1242    pustr(i) = 0.0
1243    pvstr(i) = 0.0
1244  END DO
1245  DO k = 1, klev
1246    DO i = 1, klon
1247      d_t(i, k) = 0.0
1248      d_u(i, k) = 0.0
1249      d_v(i, k) = 0.0
1250      pdudt(i, k) = 0.0
1251      pdvdt(i, k) = 0.0
1252      pdtdt(i, k) = 0.0
1253    END DO
1254  END DO
1255
1256  ! preparer les variables d'entree (attention: l'ordre des niveaux
1257  ! verticaux augmente du haut vers le bas)
1258
1259  DO k = 1, klev
1260    DO i = 1, klon
1261      pt(i, k) = t(i, klev-k+1)
1262      pu(i, k) = u(i, klev-k+1)
1263      pv(i, k) = v(i, klev-k+1)
1264      papmf(i, k) = pplay(i, klev-k+1)
1265    END DO
1266  END DO
1267  DO k = 1, klev + 1
1268    DO i = 1, klon
1269      papmh(i, k) = paprs(i, klev-k+2)
1270    END DO
1271  END DO
1272  DO i = 1, klon
1273    zgeom(i, klev) = rd*pt(i, klev)*log(papmh(i,klev+1)/papmf(i,klev))
1274  END DO
1275  DO k = klev - 1, 1, -1
1276    DO i = 1, klon
1277      zgeom(i, k) = zgeom(i, k+1) + rd*(pt(i,k)+pt(i,k+1))/2.0*log(papmf(i,k+ &
1278        1)/papmf(i,k))
1279    END DO
1280  END DO
1281
1282  ! appeler la routine principale
1283
1284  CALL orolift(klon, klev, ktest, dtime, papmh, zgeom, pt, pu, pv, plat, &
1285    pmea, pstd, ppic, pulow, pvlow, pdudt, pdvdt, pdtdt)
1286
1287  DO k = 1, klev
1288    DO i = 1, klon
1289      d_u(i, klev+1-k) = dtime*pdudt(i, k)
1290      d_v(i, klev+1-k) = dtime*pdvdt(i, k)
1291      d_t(i, klev+1-k) = dtime*pdtdt(i, k)
1292      pustr(i) = pustr(i) &        ! IM BUG .
1293                                   ! +RG*pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))
1294        +pdudt(i, k)*(papmh(i,k+1)-papmh(i,k))/rg
1295      pvstr(i) = pvstr(i) &        ! IM BUG .
1296                                   ! +RG*pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))
1297        +pdvdt(i, k)*(papmh(i,k+1)-papmh(i,k))/rg
1298    END DO
1299  END DO
1300
1301  RETURN
1302END SUBROUTINE lift_noro
1303SUBROUTINE orolift(nlon, nlev, ktest, ptsphy, paphm1, pgeom1, ptm1, pum1, &
1304    pvm1, plat, pmea, pvaror, ppic & ! OUTPUTS
1305    , pulow, pvlow, pvom, pvol, pte)
1306
1307
1308  ! **** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT.
1309
1310  ! PURPOSE.
1311  ! --------
1312
1313  ! **   INTERFACE.
1314  ! ----------
1315  ! CALLED FROM *lift_noro
1316  ! ----------
1317
1318  ! AUTHOR.
1319  ! -------
1320  ! F.LOTT  LMD 22/11/95
1321
1322  USE dimphy
1323  IMPLICIT NONE
1324
1325
1326  ! ym#include "dimensions.h"
1327  ! ym#include "dimphy.h"
1328  include "YOMCST.h"
1329  include "YOEGWD.h"
1330  ! -----------------------------------------------------------------------
1331
1332  ! *       0.1   ARGUMENTS
1333  ! ---------
1334
1335
1336  INTEGER nlon, nlev
1337  REAL pte(nlon, nlev), pvol(nlon, nlev), pvom(nlon, nlev), pulow(nlon), &
1338    pvlow(nlon)
1339  REAL pum1(nlon, nlev), pvm1(nlon, nlev), ptm1(nlon, nlev), plat(nlon), &
1340    pmea(nlon), pvaror(nlon), ppic(nlon), pgeom1(nlon, nlev), &
1341    paphm1(nlon, nlev+1)
1342
1343  INTEGER ktest(nlon)
1344  REAL ptsphy
1345  ! -----------------------------------------------------------------------
1346
1347  ! *       0.2   LOCAL ARRAYS
1348  ! ------------
1349  LOGICAL lifthigh
1350  ! ym      integer klevm1, jl, ilevh, jk
1351  INTEGER jl, ilevh, jk
1352  REAL zcons1, ztmst, zrtmst, zpi, zhgeo
1353  REAL zdelp, zslow, zsqua, zscav, zbet
1354  INTEGER iknub(klon), iknul(klon)
1355  LOGICAL ll1(klon, klev+1)
1356
1357  REAL ztau(klon, klev+1), ztav(klon, klev+1), zrho(klon, klev+1)
1358  REAL zdudt(klon), zdvdt(klon)
1359  REAL zhcrit(klon, klev)
1360  CHARACTER (LEN=20) :: modname = 'orografi'
1361  CHARACTER (LEN=80) :: abort_message
1362  ! -----------------------------------------------------------------------
1363
1364  ! *         1.1  INITIALIZATIONS
1365  ! ---------------
1366
1367  lifthigh = .FALSE.
1368
1369  IF (nlon/=klon .OR. nlev/=klev) THEN
1370    abort_message = 'pb dimension'
1371    CALL abort_gcm(modname, abort_message, 1)
1372  END IF
1373  zcons1 = 1./rd
1374  ! ym      KLEVM1=KLEV-1
1375  ztmst = ptsphy
1376  zrtmst = 1./ztmst
1377  zpi = acos(-1.)
1378
1379  DO jl = kidia, kfdia
1380    zrho(jl, klev+1) = 0.0
1381    pulow(jl) = 0.0
1382    pvlow(jl) = 0.0
1383    iknub(jl) = klev
1384    iknul(jl) = klev
1385    ilevh = klev/3
1386    ll1(jl, klev+1) = .FALSE.
1387    DO jk = 1, klev
1388      pvom(jl, jk) = 0.0
1389      pvol(jl, jk) = 0.0
1390      pte(jl, jk) = 0.0
1391    END DO
1392  END DO
1393
1394
1395  ! *         2.1     DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
1396  ! *                 LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
1397  ! *                 THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
1398
1399
1400
1401  DO jk = klev, 1, -1
1402    DO jl = kidia, kfdia
1403      IF (ktest(jl)==1) THEN
1404        zhcrit(jl, jk) = amax1(ppic(jl)-pmea(jl), 100.)
1405        zhgeo = pgeom1(jl, jk)/rg
1406        ll1(jl, jk) = (zhgeo>zhcrit(jl,jk))
1407        IF (ll1(jl,jk) .NEQV. ll1(jl,jk+1)) THEN
1408          iknub(jl) = jk
1409        END IF
1410      END IF
1411    END DO
1412  END DO
1413
1414  DO jl = kidia, kfdia
1415    IF (ktest(jl)==1) THEN
1416      iknub(jl) = max(iknub(jl), klev/2)
1417      iknul(jl) = max(iknul(jl), 2*klev/3)
1418      IF (iknub(jl)>nktopg) iknub(jl) = nktopg
1419      IF (iknub(jl)==nktopg) iknul(jl) = klev
1420      IF (iknub(jl)==iknul(jl)) iknub(jl) = iknul(jl) - 1
1421    END IF
1422  END DO
1423
1424  ! do 2011 jl=kidia,kfdia
1425  ! IF(KTEST(JL).EQ.1) THEN
1426  ! print *,' iknul= ',iknul(jl),'  iknub=',iknub(jl)
1427  ! ENDIF
1428  ! 2011 continue
1429
1430  ! PRINT *,'  DANS OROLIFT: 2010'
1431
1432  DO jk = klev, 2, -1
1433    DO jl = kidia, kfdia
1434      zrho(jl, jk) = 2.*paphm1(jl, jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
1435    END DO
1436  END DO
1437  ! PRINT *,'  DANS OROLIFT: 223'
1438
1439  ! ********************************************************************
1440
1441  ! *     DEFINE LOW LEVEL FLOW
1442  ! -------------------
1443  DO jk = klev, 1, -1
1444    DO jl = kidia, kfdia
1445      IF (ktest(jl)==1) THEN
1446        IF (jk>=iknub(jl) .AND. jk<=iknul(jl)) THEN
1447          pulow(jl) = pulow(jl) + pum1(jl, jk)*(paphm1(jl,jk+1)-paphm1(jl,jk) &
1448            )
1449          pvlow(jl) = pvlow(jl) + pvm1(jl, jk)*(paphm1(jl,jk+1)-paphm1(jl,jk) &
1450            )
1451          zrho(jl, klev+1) = zrho(jl, klev+1) + zrho(jl, jk)*(paphm1(jl,jk+1) &
1452            -paphm1(jl,jk))
1453        END IF
1454      END IF
1455    END DO
1456  END DO
1457  DO jl = kidia, kfdia
1458    IF (ktest(jl)==1) THEN
1459      pulow(jl) = pulow(jl)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
1460      pvlow(jl) = pvlow(jl)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
1461      zrho(jl, klev+1) = zrho(jl, klev+1)/(paphm1(jl,iknul(jl)+1)-paphm1(jl, &
1462        iknub(jl)))
1463    END IF
1464  END DO
1465
1466  ! ***********************************************************
1467
1468  ! *         3.      COMPUTE MOUNTAIN LIFT
1469
1470
1471  DO jl = kidia, kfdia
1472    IF (ktest(jl)==1) THEN
1473      ztau(jl, klev+1) = -gklift*zrho(jl, klev+1)*2.*romega* & ! *
1474                                                               ! (2*PVAROR(JL)+PMEA(JL))*
1475        2*pvaror(jl)*sin(zpi/180.*plat(jl))*pvlow(jl)
1476      ztav(jl, klev+1) = gklift*zrho(jl, klev+1)*2.*romega* & ! *
1477                                                              ! (2*PVAROR(JL)+PMEA(JL))*
1478        2*pvaror(jl)*sin(zpi/180.*plat(jl))*pulow(jl)
1479    ELSE
1480      ztau(jl, klev+1) = 0.0
1481      ztav(jl, klev+1) = 0.0
1482    END IF
1483  END DO
1484
1485  ! *         4.      COMPUTE LIFT PROFILE
1486  ! *                 --------------------
1487
1488
1489
1490  DO jk = 1, klev
1491    DO jl = kidia, kfdia
1492      IF (ktest(jl)==1) THEN
1493        ztau(jl, jk) = ztau(jl, klev+1)*paphm1(jl, jk)/paphm1(jl, klev+1)
1494        ztav(jl, jk) = ztav(jl, klev+1)*paphm1(jl, jk)/paphm1(jl, klev+1)
1495      ELSE
1496        ztau(jl, jk) = 0.0
1497        ztav(jl, jk) = 0.0
1498      END IF
1499    END DO
1500  END DO
1501
1502
1503  ! *         5.      COMPUTE TENDENCIES.
1504  ! *                 -------------------
1505  IF (lifthigh) THEN
1506    ! PRINT *,'  DANS OROLIFT: 500'
1507
1508    ! EXPLICIT SOLUTION AT ALL LEVELS
1509
1510    DO jk = 1, klev
1511      DO jl = kidia, kfdia
1512        IF (ktest(jl)==1) THEN
1513          zdelp = paphm1(jl, jk+1) - paphm1(jl, jk)
1514          zdudt(jl) = -rg*(ztau(jl,jk+1)-ztau(jl,jk))/zdelp
1515          zdvdt(jl) = -rg*(ztav(jl,jk+1)-ztav(jl,jk))/zdelp
1516        END IF
1517      END DO
1518    END DO
1519
1520    ! PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY
1521
1522    DO jk = 1, klev
1523      DO jl = kidia, kfdia
1524        IF (ktest(jl)==1) THEN
1525
1526          zslow = sqrt(pulow(jl)**2+pvlow(jl)**2)
1527          zsqua = amax1(sqrt(pum1(jl,jk)**2+pvm1(jl,jk)**2), gvsec)
1528          zscav = -zdudt(jl)*pvm1(jl, jk) + zdvdt(jl)*pum1(jl, jk)
1529          IF (zsqua>gvsec) THEN
1530            pvom(jl, jk) = -zscav*pvm1(jl, jk)/zsqua**2
1531            pvol(jl, jk) = zscav*pum1(jl, jk)/zsqua**2
1532          ELSE
1533            pvom(jl, jk) = 0.0
1534            pvol(jl, jk) = 0.0
1535          END IF
1536          zsqua = sqrt(pum1(jl,jk)**2+pum1(jl,jk)**2)
1537          IF (zsqua<zslow) THEN
1538            pvom(jl, jk) = zsqua/zslow*pvom(jl, jk)
1539            pvol(jl, jk) = zsqua/zslow*pvol(jl, jk)
1540          END IF
1541
1542        END IF
1543      END DO
1544    END DO
1545
1546    ! 6.  LOW LEVEL LIFT, SEMI IMPLICIT:
1547    ! ----------------------------------
1548
1549  ELSE
1550
1551    DO jl = kidia, kfdia
1552      IF (ktest(jl)==1) THEN
1553        DO jk = klev, iknub(jl), -1
1554          zbet = gklift*2.*romega*sin(zpi/180.*plat(jl))*ztmst* &
1555            (pgeom1(jl,iknub(jl)-1)-pgeom1(jl,jk))/ &
1556            (pgeom1(jl,iknub(jl)-1)-pgeom1(jl,klev))
1557          zdudt(jl) = -pum1(jl, jk)/ztmst/(1+zbet**2)
1558          zdvdt(jl) = -pvm1(jl, jk)/ztmst/(1+zbet**2)
1559          pvom(jl, jk) = zbet**2*zdudt(jl) - zbet*zdvdt(jl)
1560          pvol(jl, jk) = zbet*zdudt(jl) + zbet**2*zdvdt(jl)
1561        END DO
1562      END IF
1563    END DO
1564
1565  END IF
1566
1567  RETURN
1568END SUBROUTINE orolift
1569
1570
1571SUBROUTINE sugwd(nlon, nlev, paprs, pplay)
1572  USE dimphy
1573  USE mod_phys_lmdz_para
1574  USE mod_grid_phy_lmdz
1575  ! USE parallel
1576
1577  ! **** *SUGWD* INITIALIZE COMMON YOEGWD CONTROLLING GRAVITY WAVE DRAG
1578
1579  ! PURPOSE.
1580  ! --------
1581  ! INITIALIZE YOEGWD, THE COMMON THAT CONTROLS THE
1582  ! GRAVITY WAVE DRAG PARAMETRIZATION.
1583
1584  ! **   INTERFACE.
1585  ! ----------
1586  ! CALL *SUGWD* FROM *SUPHEC*
1587  ! -----        ------
1588
1589  ! EXPLICIT ARGUMENTS :
1590  ! --------------------
1591  ! PSIG        : VERTICAL COORDINATE TABLE
1592  ! NLEV        : NUMBER OF MODEL LEVELS
1593
1594  ! IMPLICIT ARGUMENTS :
1595  ! --------------------
1596  ! COMMON YOEGWD
1597
1598  ! METHOD.
1599  ! -------
1600  ! SEE DOCUMENTATION
1601
1602  ! EXTERNALS.
1603  ! ----------
1604  ! NONE
1605
1606  ! REFERENCE.
1607  ! ----------
1608  ! ECMWF Research Department documentation of the IFS
1609
1610  ! AUTHOR.
1611  ! -------
1612  ! MARTIN MILLER             *ECMWF*
1613
1614  ! MODIFICATIONS.
1615  ! --------------
1616  ! ORIGINAL : 90-01-01
1617  ! ------------------------------------------------------------------
1618  IMPLICIT NONE
1619
1620  ! -----------------------------------------------------------------
1621  include "YOEGWD.h"
1622  ! ----------------------------------------------------------------
1623
1624  INTEGER nlon, nlev, jk
1625  REAL paprs(nlon, nlev+1)
1626  REAL pplay(nlon, nlev)
1627  REAL zpr, zstra, zsigt, zpm1r
1628  REAL :: pplay_glo(klon_glo, nlev)
1629  REAL :: paprs_glo(klon_glo, nlev+1)
1630
1631  ! *       1.    SET THE VALUES OF THE PARAMETERS
1632  ! --------------------------------
1633
1634
1635  PRINT *, ' DANS SUGWD NLEV=', nlev
1636  ghmax = 10000.
1637
1638  zpr = 100000.
1639  zstra = 0.1
1640  zsigt = 0.94
1641  ! old  ZPR=80000.
1642  ! old  ZSIGT=0.85
1643
1644
1645  CALL gather(pplay, pplay_glo)
1646  CALL bcast(pplay_glo)
1647  CALL gather(paprs, paprs_glo)
1648  CALL bcast(paprs_glo)
1649
1650
1651  DO jk = 1, nlev
1652    zpm1r = pplay_glo((klon_glo/2)+1, jk)/paprs_glo((klon_glo/2)+1, 1)
1653    IF (zpm1r>=zsigt) THEN
1654      nktopg = jk
1655    END IF
1656    zpm1r = pplay_glo((klon_glo/2)+1, jk)/paprs_glo((klon_glo/2)+1, 1)
1657    IF (zpm1r>=zstra) THEN
1658      nstra = jk
1659    END IF
1660  END DO
1661
1662
1663
1664  ! inversion car dans orodrag on compte les niveaux a l'envers
1665  nktopg = nlev - nktopg + 1
1666  nstra = nlev - nstra
1667  PRINT *, ' DANS SUGWD nktopg=', nktopg
1668  PRINT *, ' DANS SUGWD nstra=', nstra
1669
1670  gsigcr = 0.80
1671
1672  gkdrag = 0.2
1673  grahilo = 1.
1674  grcrit = 0.01
1675  gfrcrit = 1.0
1676  gkwake = 0.50
1677
1678  gklift = 0.50
1679  gvcrit = 0.0
1680
1681  ! ----------------------------------------------------------------
1682
1683  ! *       2.    SET VALUES OF SECURITY PARAMETERS
1684  ! ---------------------------------
1685
1686
1687  gvsec = 0.10
1688  gssec = 1.E-12
1689
1690  gtsec = 1.E-07
1691
1692  ! ----------------------------------------------------------------
1693
1694  RETURN
1695END SUBROUTINE sugwd
Note: See TracBrowser for help on using the repository browser.