source: LMDZ4/branches/LMDZ4_AR5/libf/phylmd/orografi.F @ 3536

Last change on this file since 3536 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 48.3 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)
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)
1499C-----------------------------------------------------------------------
1500C
1501C*         1.1  INITIALIZATIONS
1502C               ---------------
1503
1504      LIFTHIGH=.FALSE.
1505
1506      IF(NLON.NE.KLON.OR.NLEV.NE.KLEV)STOP
1507      ZCONS1=1./RD
1508cym      KLEVM1=KLEV-1
1509      ZTMST=PTSPHY
1510      ZRTMST=1./ZTMST
1511      ZPI=ACOS(-1.)
1512C
1513      DO 1001 JL=kidia,kfdia
1514      ZRHO(JL,KLEV+1)  =0.0
1515      PULOW(JL)        =0.0
1516      PVLOW(JL)        =0.0
1517      iknub(JL)   =klev
1518      iknul(JL)   =klev
1519      ilevh=klev/3
1520      ll1(jl,klev+1)=.false.
1521      DO 1000 JK=1,KLEV
1522      PVOM(JL,JK)=0.0
1523      PVOL(JL,JK)=0.0
1524      PTE (JL,JK)=0.0
1525 1000 CONTINUE
1526 1001 CONTINUE
1527
1528C
1529C*         2.1     DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
1530C*                 LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
1531C*                 THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
1532C
1533C
1534C
1535      DO 2006 JK=KLEV,1,-1
1536      DO 2007 JL=kidia,kfdia
1537      IF(KTEST(JL).EQ.1) THEN
1538      ZHCRIT(JL,JK)=amax1(Ppic(JL)-pmea(JL),100.)
1539      ZHGEO=PGEOM1(JL,JK)/RG
1540      ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
1541      IF(ll1(JL,JK).neqv.ll1(JL,JK+1)) THEN
1542        iknub(JL)=JK
1543      ENDIF
1544      ENDIF
1545 2007 CONTINUE
1546 2006 CONTINUE
1547C
1548      do 2010 jl=kidia,kfdia
1549      IF(KTEST(JL).EQ.1) THEN
1550      iknub(jl)=max(iknub(jl),klev/2)
1551      iknul(jl)=max(iknul(jl),2*klev/3)
1552      if(iknub(jl).gt.nktopg) iknub(jl)=nktopg
1553      if(iknub(jl).eq.nktopg) iknul(jl)=klev
1554      if(iknub(jl).eq.iknul(jl)) iknub(jl)=iknul(jl)-1
1555      ENDIF
1556 2010 continue
1557
1558C     do 2011 jl=kidia,kfdia
1559C     IF(KTEST(JL).EQ.1) THEN
1560C       print *,' iknul= ',iknul(jl),'  iknub=',iknub(jl)
1561C     ENDIF
1562C2011 continue
1563
1564C     PRINT *,'  DANS OROLIFT: 2010'
1565
1566      DO 223 JK=KLEV,2,-1
1567      DO 222 JL=kidia,kfdia
1568        ZRHO(JL,JK)=2.*PAPHM1(JL,JK)*ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1))
1569  222 CONTINUE
1570  223 CONTINUE
1571C     PRINT *,'  DANS OROLIFT: 223'
1572
1573C********************************************************************
1574C
1575C*     DEFINE LOW LEVEL FLOW
1576C      -------------------
1577      DO 2115 JK=klev,1,-1
1578      DO 2116 JL=kidia,kfdia
1579      IF(KTEST(JL).EQ.1) THEN
1580      if(jk.ge.iknub(jl).and.jk.le.iknul(jl)) then
1581        pulow(JL)=pulow(JL)+PUM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
1582        pvlow(JL)=pvlow(JL)+PVM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
1583        zrho(JL,klev+1)=zrho(JL,klev+1)
1584     *                 +zrho(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
1585      end if
1586      ENDIF
1587 2116 CONTINUE
1588 2115 CONTINUE
1589      DO 2110 JL=kidia,kfdia
1590      IF(KTEST(JL).EQ.1) THEN
1591      pulow(JL)=pulow(JL)/(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl)))
1592      pvlow(JL)=pvlow(JL)/(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl)))
1593      zrho(JL,klev+1)=zrho(JL,klev+1)
1594     *               /(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl)))
1595      ENDIF
1596 2110 CONTINUE
1597
1598
1599200   CONTINUE
1600
1601C***********************************************************
1602C
1603C*         3.      COMPUTE MOUNTAIN LIFT
1604C
1605  300 CONTINUE
1606C
1607      DO 301 JL=kidia,kfdia
1608      IF(KTEST(JL).EQ.1) THEN
1609       ZTAU(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))*PVLOW(JL)
1613       ZTAV(JL,KLEV+1)=   GKLIFT*ZRHO(JL,KLEV+1)*2.*ROMEGA*
1614C    *                 (2*PVAROR(JL)+PMEA(JL))*
1615     *                 2*PVAROR(JL)*
1616     *                 SIN(ZPI/180.*PLAT(JL))*PULOW(JL)
1617      ELSE
1618       ZTAU(JL,KLEV+1)=0.0
1619       ZTAV(JL,KLEV+1)=0.0
1620      ENDIF
1621301   CONTINUE
1622
1623C
1624C*         4.      COMPUTE LIFT PROFILE         
1625C*                 --------------------   
1626C
1627
1628  400 CONTINUE
1629
1630      DO 401 JK=1,KLEV
1631      DO 401 JL=kidia,kfdia
1632      IF(KTEST(JL).EQ.1) THEN
1633      ZTAU(JL,JK)=ZTAU(JL,KLEV+1)*PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1)
1634      ZTAV(JL,JK)=ZTAV(JL,KLEV+1)*PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1)
1635      ELSE
1636      ZTAU(JL,JK)=0.0
1637      ZTAV(JL,JK)=0.0
1638      ENDIF
1639401   CONTINUE
1640C
1641C
1642C*         5.      COMPUTE TENDENCIES.
1643C*                 -------------------
1644      IF(LIFTHIGH)THEN
1645C
1646  500 CONTINUE
1647C     PRINT *,'  DANS OROLIFT: 500'
1648C
1649C  EXPLICIT SOLUTION AT ALL LEVELS
1650C
1651      DO 524 JK=1,klev
1652      DO 523 JL=KIDIA,KFDIA
1653      IF(KTEST(JL).EQ.1) THEN
1654      ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK)
1655      ZDUDT(JL)=-RG*(ZTAU(JL,JK+1)-ZTAU(JL,JK))/ZDELP
1656      ZDVDT(JL)=-RG*(ZTAV(JL,JK+1)-ZTAV(JL,JK))/ZDELP
1657      ENDIF 
1658  523 CONTINUE
1659  524 CONTINUE
1660C
1661C  PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY
1662C
1663      DO 530 JK=1,klev
1664      DO 530 JL=KIDIA,KFDIA
1665      IF(KTEST(JL).EQ.1) THEN
1666
1667        ZSLOW=SQRT(PULOW(JL)**2+PVLOW(JL)**2)
1668        ZSQUA=AMAX1(SQRT(PUM1(JL,JK)**2+PVM1(JL,JK)**2),GVSEC)
1669        ZSCAV=-ZDUDT(JL)*PVM1(JL,JK)+ZDVDT(JL)*PUM1(JL,JK)
1670        IF(ZSQUA.GT.GVSEC)THEN
1671          PVOM(JL,JK)=-ZSCAV*PVM1(JL,JK)/ZSQUA**2
1672          PVOL(JL,JK)= ZSCAV*PUM1(JL,JK)/ZSQUA**2
1673        ELSE
1674          PVOM(JL,JK)=0.0
1675          PVOL(JL,JK)=0.0     
1676        ENDIF 
1677        ZSQUA=SQRT(PUM1(JL,JK)**2+PUM1(JL,JK)**2)               
1678        IF(ZSQUA.LT.ZSLOW)THEN
1679          PVOM(JL,JK)=ZSQUA/ZSLOW*PVOM(JL,JK)
1680          PVOL(JL,JK)=ZSQUA/ZSLOW*PVOL(JL,JK)
1681        ENDIF
1682
1683      ENDIF 
1684530   CONTINUE
1685C
1686C  6.  LOW LEVEL LIFT, SEMI IMPLICIT:
1687C  ----------------------------------
1688
1689      ELSE
1690
1691        DO 601 JL=KIDIA,KFDIA
1692        IF(KTEST(JL).EQ.1) THEN
1693          DO JK=KLEV,IKNUB(JL),-1
1694          ZBET=GKLIFT*2.*ROMEGA*SIN(ZPI/180.*PLAT(JL))*ztmst*
1695     *        (PGEOM1(JL,IKNUB(JL)-1)-PGEOM1(JL,  JK))/
1696     *        (PGEOM1(JL,IKNUB(JL)-1)-PGEOM1(JL,KLEV))
1697          ZDUDT(JL)=-PUM1(JL,JK)/ztmst/(1+ZBET**2)
1698          ZDVDT(JL)=-PVM1(JL,JK)/ztmst/(1+ZBET**2)
1699          PVOM(JL,JK)= ZBET**2*ZDUDT(JL) - ZBET   *ZDVDT(JL)
1700          PVOL(JL,JK)= ZBET   *ZDUDT(JL) + ZBET**2*ZDVDT(JL)   
1701          ENDDO
1702        ENDIF
1703 601    CONTINUE
1704
1705      ENDIF
1706
1707      RETURN
1708      END
1709
1710
1711      SUBROUTINE SUGWD(NLON,NLEV,paprs,pplay)
1712      USE dimphy
1713      USE mod_phys_lmdz_para
1714      USE mod_grid_phy_lmdz
1715c      USE parallel
1716C
1717C**** *SUGWD* INITIALIZE COMMON YOEGWD CONTROLLING GRAVITY WAVE DRAG
1718C
1719C     PURPOSE.
1720C     --------
1721C           INITIALIZE YOEGWD, THE COMMON THAT CONTROLS THE
1722C           GRAVITY WAVE DRAG PARAMETRIZATION.
1723C
1724C**   INTERFACE.
1725C     ----------
1726C        CALL *SUGWD* FROM *SUPHEC*
1727C              -----        ------
1728C
1729C        EXPLICIT ARGUMENTS :
1730C        --------------------
1731C        PSIG        : VERTICAL COORDINATE TABLE
1732C        NLEV        : NUMBER OF MODEL LEVELS
1733C
1734C        IMPLICIT ARGUMENTS :
1735C        --------------------
1736C        COMMON YOEGWD
1737C
1738C     METHOD.
1739C     -------
1740C        SEE DOCUMENTATION
1741C
1742C     EXTERNALS.
1743C     ----------
1744C        NONE
1745C
1746C     REFERENCE.
1747C     ----------
1748C        ECMWF Research Department documentation of the IFS
1749C
1750C     AUTHOR.
1751C     -------
1752C        MARTIN MILLER             *ECMWF*
1753C
1754C     MODIFICATIONS.
1755C     --------------
1756C        ORIGINAL : 90-01-01
1757C     ------------------------------------------------------------------
1758      implicit none
1759C
1760C     -----------------------------------------------------------------
1761#include "YOEGWD.h"
1762C      ----------------------------------------------------------------
1763C
1764      integer nlon,nlev, jk
1765      REAL paprs(nlon,nlev+1)
1766      REAL pplay(nlon,nlev)
1767      real zpr,zstra,zsigt,zpm1r
1768      REAL :: pplay_glo(klon_glo,nlev)
1769      REAL :: paprs_glo(klon_glo,nlev+1)
1770
1771C
1772C*       1.    SET THE VALUES OF THE PARAMETERS
1773C              --------------------------------
1774C
1775 100  CONTINUE
1776C
1777      PRINT *,' DANS SUGWD NLEV=',NLEV
1778      GHMAX=10000.
1779C
1780      ZPR=100000.
1781      ZSTRA=0.1
1782      ZSIGT=0.94
1783cold  ZPR=80000.
1784cold  ZSIGT=0.85
1785C
1786     
1787      CALL gather(pplay,pplay_glo)
1788      CALL bcast(pplay_glo)
1789      CALL gather(paprs,paprs_glo)
1790      CALL bcast(paprs_glo)
1791     
1792           
1793      DO 110 JK=1,NLEV
1794      ZPM1R=pplay_glo((klon_glo/2)+1,jk)/paprs_glo((klon_glo/2)+1,1)
1795      IF(ZPM1R.GE.ZSIGT)THEN
1796         nktopg=JK
1797      ENDIF
1798      ZPM1R=pplay_glo((klon_glo/2)+1,jk)/paprs_glo((klon_glo/2)+1,1)
1799      IF(ZPM1R.GE.ZSTRA)THEN
1800         NSTRA=JK
1801      ENDIF
1802  110 CONTINUE
1803
1804
1805c
1806c  inversion car dans orodrag on compte les niveaux a l'envers
1807      nktopg=nlev-nktopg+1
1808      nstra=nlev-nstra
1809      print *,' DANS SUGWD nktopg=', nktopg
1810      print *,' DANS SUGWD nstra=', nstra
1811C
1812      GSIGCR=0.80
1813C
1814      GKDRAG=0.2
1815      GRAHILO=1.   
1816      GRCRIT=0.01
1817      GFRCRIT=1.0
1818      GKWAKE=0.50
1819C
1820      GKLIFT=0.50 
1821      GVCRIT =0.0
1822C
1823C
1824C      ----------------------------------------------------------------
1825C
1826C*       2.    SET VALUES OF SECURITY PARAMETERS
1827C              ---------------------------------
1828C
1829 200  CONTINUE
1830C
1831      GVSEC=0.10
1832      GSSEC=1.E-12
1833C
1834      GTSEC=1.E-07
1835C
1836C      ----------------------------------------------------------------
1837C
1838      RETURN
1839      END
Note: See TracBrowser for help on using the repository browser.