source: LMDZ5/trunk/libf/phylmd/orografi.F @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

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