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

Last change on this file since 161 was 160, checked in by lmdzadmin, 24 years ago

zb et zc declares deux fois
LF

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