source: LMDZ4/branches/LMDZ4_V2_patch/libf/phylmd/orografi.F @ 757

Last change on this file since 757 was 644, checked in by Laurent Fairhead, 20 years ago

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