source: LMDZ.3.3/branches/rel-LF/libf/phylmd/orografi.F @ 338

Last change on this file since 338 was 177, checked in by lmdzadmin, 24 years ago

Lots of stuff, plus particulierement:

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