source: LMDZ4/branches/LMDZ4_par_0/libf/phylmd/orografi.F @ 3872

Last change on this file since 3872 was 634, checked in by Laurent Fairhead, 19 years ago

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