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

Last change on this file since 5392 was 495, checked in by lmdzadmin, 21 years ago

Optimisation de differentes routines, IM, MAF, FH
LF

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