source: LMDZ.3.3/trunk/libf/phylmd/orografi.F @ 343

Last change on this file since 343 was 2, checked in by lmdz, 25 years ago

Initial revision

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