!
! $Id: orografi.F 1403 2010-07-01 09:02:53Z fairhead $
!
MODULE orografi

  CONTAINS

      SUBROUTINE drag_noro (nlon,nlev,dtime,paprs,pplay,                             &
       &                   pmea,pstd, psig, pgam, pthe,ppic,pval,                    &
       &                   kgwd,kdx,ktest,                                           &
       &                   t, u, v,                                                  &
       &                   pulow, pvlow, pustr, pvstr,                               &
       &                   d_t, d_u, d_v)
!c
      USE dimphy
      IMPLICIT none
!c======================================================================
!c Auteur(s): F.Lott (LMD/CNRS) date: 19950201
!c Objet: Frottement de la montagne Interface
!c======================================================================
!c Arguments:
!c dtime---input-R- pas d'integration (s)
!c paprs---input-R-pression pour chaque inter-couche (en Pa)
!c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
!c t-------input-R-temperature (K)
!c u-------input-R-vitesse horizontale (m/s)
!c v-------input-R-vitesse horizontale (m/s)
!c
!c d_t-----output-R-increment de la temperature             
!c d_u-----output-R-increment de la vitesse u
!c d_v-----output-R-increment de la vitesse v
!c======================================================================
!cym#include "dimensions.h"
!cym#include "dimphy.h"
#include "YOMCST.h"
!c
!c ARGUMENTS
!c
      INTEGER nlon,nlev
      REAL dtime
      REAL paprs(klon,klev+1)
      REAL pplay(klon,klev)
      REAL pmea(nlon),pstd(nlon),psig(nlon),pgam(nlon),pthe(nlon)
      REAL ppic(nlon),pval(nlon)
      REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon)
      REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev)
      REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev)
!c
      INTEGER i, k, kgwd, kdx(nlon), ktest(nlon)
!c
!c Variables locales:
!c
      REAL zgeom(klon,klev)
      REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev)
      REAL pt(klon,klev), pu(klon,klev), pv(klon,klev)
      REAL papmf(klon,klev),papmh(klon,klev+1)
!c
!c initialiser les variables de sortie (pour securite)
!c
      DO i = 1,klon
         pulow(i) = 0.0
         pvlow(i) = 0.0
         pustr(i) = 0.0
         pvstr(i) = 0.0
      ENDDO
      DO k = 1, klev
      DO i = 1, klon
         d_t(i,k) = 0.0
         d_u(i,k) = 0.0
         d_v(i,k) = 0.0
         pdudt(i,k)=0.0
         pdvdt(i,k)=0.0
         pdtdt(i,k)=0.0
      ENDDO
      ENDDO
!c
!c preparer les variables d'entree (attention: l'ordre des niveaux 
!c verticaux augmente du haut vers le bas)
!c
      DO k = 1, klev
      DO i = 1, klon
         pt(i,k) = t(i,klev-k+1) 
         pu(i,k) = u(i,klev-k+1)
         pv(i,k) = v(i,klev-k+1)
         papmf(i,k) = pplay(i,klev-k+1)
      ENDDO
      ENDDO
      DO k = 1, klev+1
      DO i = 1, klon
         papmh(i,k) = paprs(i,klev-k+2)
      ENDDO
      ENDDO
      DO i = 1, klon
         zgeom(i,klev) = RD * pt(i,klev)                                             &
       &                  * LOG(papmh(i,klev+1)/papmf(i,klev))
      ENDDO
      DO k = klev-1, 1, -1
      DO i = 1, klon
         zgeom(i,k) = zgeom(i,k+1) + RD * (pt(i,k)+pt(i,k+1))/2.0                    &
       &               * LOG(papmf(i,k+1)/papmf(i,k))
      ENDDO
      ENDDO
!c
!c appeler la routine principale
!c
      CALL orodrag(klon,klev,kgwd,kdx,ktest,                                         &
       &            dtime,                                                           &
       &            papmh, papmf, zgeom,                                             &
       &            pt, pu, pv,                                                      &
       &            pmea, pstd, psig, pgam, pthe, ppic,pval,                         &
       &            pulow,pvlow,                                                     &
       &            pdudt,pdvdt,pdtdt)
!C
      DO k = 1, klev
      DO i = 1, klon
         d_u(i,klev+1-k) = dtime*pdudt(i,k)
         d_v(i,klev+1-k) = dtime*pdvdt(i,k)
         d_t(i,klev+1-k) = dtime*pdtdt(i,k)
         pustr(i)        = pustr(i)                                                  &
!IM BUG  .                +rg*pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))                   &
       &                    +pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))/RG
         pvstr(i)        = pvstr(i)                                                  &
!IM BUG  .                +rg*pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))                   &
       &                    +pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))/RG
      ENDDO
      ENDDO
!c
      RETURN
      END SUBROUTINE drag_noro

      SUBROUTINE orodrag( nlon,nlev                                                  &
       &                 , kgwd, kdx, ktest                                          &
       &                 , ptsphy                                                    &
       &                 , paphm1,papm1,pgeom1,ptm1,pum1,pvm1                        &
       &                 , pmea, pstd, psig, pgamma, ptheta, ppic, pval              &
!c outputs                                                                           &
       &                 , pulow,pvlow                                               &
       &                 , pvom,pvol,pte )

      USE dimphy
      implicit none

!c
!c
!c**** *gwdrag* - does the gravity wave parametrization.
!c
!c     purpose.
!c     --------
!c
!c          this routine computes the physical tendencies of the
!c     prognostic variables u,v  and t due to  vertical transports by
!c     subgridscale orographically excited gravity waves
!c
!c**   interface.
!c     ----------
!c          called from *callpar*.
!c
!c          the routine takes its input from the long-term storage:
!c          u,v,t and p at t-1.
!c
!c        explicit arguments :
!c        --------------------
!c     ==== inputs ===
!c     ==== outputs ===
!c
!c        implicit arguments :   none
!c        --------------------
!c
!c      implicit logical (l)
!c
!c     method.
!c     -------
!c
!c     externals.
!c     ----------
      integer ismin, ismax
      external ismin, ismax
!c
!c     reference.
!c     ----------
!c
!c     author.
!c     -------
!c     m.miller + b.ritter   e.c.m.w.f.     15/06/86.
!c
!c     f.lott + m. miller    e.c.m.w.f.     22/11/94
!c-----------------------------------------------------------------------
!c
!c
!cym#include "dimensions.h"
!cym#include "dimphy.h"
#include "YOMCST.h"
#include "YOEGWD.h"
!c-----------------------------------------------------------------------
!c
!c*       0.1   arguments
!c              ---------
!c
!c
!cym      integer nlon, nlev, klevm1
      integer nlon, nlev
      integer kgwd, jl, ilevp1, jk, ji
      real zdelp, ztemp, zforc, ztend
      real rover, zb, zc, zconb, zabsv
      real zzd1, ratio, zbet, zust,zvst, zdis
      real  pte(nlon,nlev),                                                          &
       &      pvol(nlon,nlev),                                                       &
       &      pvom(nlon,nlev),                                                       &
       &      pulow(klon),                                                           &
       &      pvlow(klon)
      real  pum1(nlon,nlev),                                                         &
       &      pvm1(nlon,nlev),                                                       &
       &      ptm1(nlon,nlev),                                                       &
       &      pmea(nlon),pstd(nlon),psig(nlon),                                      &
       &      pgamma(nlon),ptheta(nlon),ppic(nlon),pval(nlon),                       &
       &      pgeom1(nlon,nlev),                                                     &
       &      papm1(nlon,nlev),                                                      &
       &      paphm1(nlon,nlev+1)
!c
      integer  kdx(nlon),ktest(nlon)
!c-----------------------------------------------------------------------
!c
!c*       0.2   local arrays
!c              ------------
      integer  isect(klon),                                                          &
       &         icrit(klon),                                                        &
       &         ikcrith(klon),                                                      &
       &         ikenvh(klon),                                                       &
       &         iknu(klon),                                                         &
       &         iknu2(klon),                                                        &
       &         ikcrit(klon),                                                       &
       &         ikhlim(klon)
!c
      real   ztau(klon,klev+1),                                                      &
       &       ztauf(klon,klev+1),                                                   &
       &       zstab(klon,klev+1),                                                   &
       &       zvph(klon,klev+1),                                                    &
       &       zrho(klon,klev+1),                                                    &
       &       zri(klon,klev+1),                                                     &
       &       zpsi(klon,klev+1),                                                    &
       &       zzdep(klon,klev)
      real   zdudt(klon),                                                            &
       &       zdvdt(klon),                                                          &
       &       zdtdt(klon),                                                          &
       &       zdedt(klon),                                                          &
       &       zvidis(klon),                                                         &
       &       znu(klon),                                                            &
       &       zd1(klon),                                                            &
       &       zd2(klon),                                                            &
       &       zdmod(klon)
      real ztmst, ptsphy, zrtmst 
!c
!c------------------------------------------------------------------
!c
!c*         1.    initialization
!c                --------------
!c
 100  continue
!c
!c     ------------------------------------------------------------------
!c
!c*         1.1   computational constants
!c                -----------------------
!c
 110  continue
!c
!c     ztmst=twodt
!c     if(nstep.eq.nstart) ztmst=0.5*twodt
!cym      klevm1=klev-1
      ztmst=ptsphy
      zrtmst=1./ztmst
!c     ------------------------------------------------------------------
!c
 120  continue
!c
!c     ------------------------------------------------------------------
!c
!c*         1.3   check whether row contains point for printing
!c                ---------------------------------------------
!c
 130  continue
!c
!c     ------------------------------------------------------------------
!c
!c*         2.     precompute basic state variables.
!c*                ---------- ----- ----- ----------
!c*                define low level wind, project winds in plane of
!c*                low level wind, determine sector in which to take
!c*                the variance and set indicator for critical levels.
!c
  200 continue
!c
!c
!c
      call orosetup                                                                  &
       &     ( nlon, ktest                                                           &
       &     , ikcrit, ikcrith, icrit,  ikenvh,iknu,iknu2                            &
       &     , paphm1, papm1 , pum1   , pvm1 , ptm1 , pgeom1, pstd                   &
       &     , zrho  , zri   , zstab  , ztau , zvph , zpsi, zzdep                    &
       &     , pulow, pvlow                                                          &
       &     , ptheta,pgamma,pmea,ppic,pval,znu  ,zd1,  zd2,  zdmod )
!c
!c
!c
!c***********************************************************
!c
!c
!c*         3.      compute low level stresses using subcritical and
!c*                 supercritical forms.computes anisotropy coefficient
!c*                 as measure of orographic twodimensionality.
!c
  300 continue
!c
      call gwstress                                                                  &
       &    ( nlon  , nlev                                                           &
       &    , ktest , icrit, ikenvh, iknu                                            &
       &    , zrho  , zstab, zvph  , pstd,  psig, pmea, ppic                         &
       &    , ztau                                                                   &
       &    , pgeom1,zdmod)
!c
!c
!c*         4.      compute stress profile.
!c*                 ------- ------ --------
!c
  400 continue
!c
!c
      call gwprofil                                                                  &
       &       (  nlon , nlev                                                        &
       &       , kgwd   , kdx , ktest                                                &
       &       , ikcrith, icrit                                                      &
       &       , paphm1, zrho   , zstab ,  zvph                                      &
       &       , zri   , ztau                                                        &
       &       , zdmod , psig  , pstd)
!c
!c
!c*         5.      compute tendencies.
!c*                 -------------------
!c
  500 continue
!c
!c  explicit solution at all levels for the gravity wave
!c  implicit solution for the blocked levels

      do 510 jl=kidia,kfdia
      zvidis(jl)=0.0
      zdudt(jl)=0.0
      zdvdt(jl)=0.0
      zdtdt(jl)=0.0
  510 continue
!c
      ilevp1=klev+1
!c
!c
      do 524 jk=1,klev
!c
!c
!c     do 523 jl=1,kgwd
!c     ji=kdx(jl)
!c  Modif vectorisation 02/04/2004
      do 523 ji=kidia,kfdia
      if(ktest(ji).eq.1) then

      zdelp=paphm1(ji,jk+1)-paphm1(ji,jk)
      ztemp=-rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,ilevp1)*zdelp)
      zdudt(ji)=(pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)
      zdvdt(ji)=(pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)
!c
!c controle des overshoots:
!c
      zforc=sqrt(zdudt(ji)**2+zdvdt(ji)**2)+1.E-12
      ztend=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/ztmst+1.E-12
      rover=0.25
      if(zforc.ge.rover*ztend)then
        zdudt(ji)=rover*ztend/zforc*zdudt(ji)
        zdvdt(ji)=rover*ztend/zforc*zdvdt(ji)
      endif
!c
!c fin du controle des overshoots
!c
      if(jk.ge.ikenvh(ji)) then
         zb=1.0-0.18*pgamma(ji)-0.04*pgamma(ji)**2
         zc=0.48*pgamma(ji)+0.3*pgamma(ji)**2
         zconb=2.*ztmst*gkwake*psig(ji)/(4.*pstd(ji))
         zabsv=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/2.
         zzd1=zb*cos(zpsi(ji,jk))**2+zc*sin(zpsi(ji,jk))**2
	     ratio=(cos(zpsi(ji,jk))**2+pgamma(ji)*sin(zpsi(ji,jk))**2)/                    &
       &   (pgamma(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2)
         zbet=max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv
!c
!c simplement oppose au vent
!c
         zdudt(ji)=-pum1(ji,jk)/ztmst
         zdvdt(ji)=-pvm1(ji,jk)/ztmst
!c
!c  projection dans la direction de l'axe principal de l'orographie
!cmod     zdudt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.)
!cmod *              +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.))
!cmod *              *cos(ptheta(ji)*rpi/180.)/ztmst
!cmod     zdvdt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.)
!cmod *              +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.))
!cmod *              *sin(ptheta(ji)*rpi/180.)/ztmst
         zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet))
         zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet))
      end if
      pvom(ji,jk)=zdudt(ji)
      pvol(ji,jk)=zdvdt(ji)
      zust=pum1(ji,jk)+ztmst*zdudt(ji)
      zvst=pvm1(ji,jk)+ztmst*zdvdt(ji)
      zdis=0.5*(pum1(ji,jk)**2+pvm1(ji,jk)**2-zust**2-zvst**2)
      zdedt(ji)=zdis/ztmst
      zvidis(ji)=zvidis(ji)+zdis*zdelp
      zdtdt(ji)=zdedt(ji)/rcpd
!c     pte(ji,jk)=zdtdt(ji)
!c
!c  ENCORE UN TRUC POUR EVITER LES EXPLOSIONS
!c
      pte(ji,jk)=0.0

      endif
  523 continue

  524 continue
!c
!c
      return
      end SUBROUTINE orodrag

      SUBROUTINE orosetup                                                            &
       &         ( nlon   , ktest                                                    &
       &         , kkcrit, kkcrith, kcrit                                            &
       &         , kkenvh, kknu  , kknu2                                             &
       &         , paphm1, papm1 , pum1   , pvm1 , ptm1  , pgeom1, pstd              &
       &         , prho  , pri   , pstab  , ptau , pvph  ,ppsi, pzdep                &
       &         , pulow , pvlow                                                     &
       &         , ptheta, pgamma, pmea, ppic, pval                                  &
       &         , pnu  ,  pd1  ,  pd2  ,pdmod  )
!c
!c**** *gwsetup*
!c
!c     purpose.
!c     --------
!c
!c**   interface.
!c     ----------
!c          from *orodrag*
!c
!c        explicit arguments :
!c        --------------------
!c     ==== inputs ===
!c     ==== outputs ===
!c
!c        implicit arguments :   none
!c        --------------------
!c
!c     method.
!c     -------
!c
!c
!c     externals.
!c     ----------
!c
!c
!c     reference.
!c     ----------
!c
!c        see ecmwf research department documentation of the "i.f.s."
!c
!c     author.
!c     -------
!c
!c     modifications.
!c     --------------
!c     f.lott  for the new-gwdrag scheme november 1993
!c
!c-----------------------------------------------------------------------
      USE dimphy
      implicit none
!c

!cym#include "dimensions.h"
!cym#include "dimphy.h"
#include "YOMCST.h"
#include "YOEGWD.h"

!c-----------------------------------------------------------------------
!c
!c*       0.1   arguments
!c              ---------
!c
      integer nlon
      integer jl, jk
      real zdelp

      integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon),                                &
       &        ktest(nlon),kkenvh(nlon)

!c
      real paphm1(nlon,klev+1),papm1(nlon,klev),pum1(nlon,klev),                     &
       &     pvm1(nlon,klev),ptm1(nlon,klev),pgeom1(nlon,klev),                      &
       &     prho(nlon,klev+1),pri(nlon,klev+1),pstab(nlon,klev+1),                  &
       &     ptau(nlon,klev+1),pvph(nlon,klev+1),ppsi(nlon,klev+1),                  &
       &     pzdep(nlon,klev)
       real pulow(nlon),pvlow(nlon),ptheta(nlon),pgamma(nlon),pnu(nlon),             &
       &     pd1(nlon),pd2(nlon),pdmod(nlon)
      real pstd(nlon),pmea(nlon),ppic(nlon),pval(nlon)
!c
!c-----------------------------------------------------------------------
!c
!c*       0.2   local arrays
!c              ------------
!c
!c
      integer ilevm1, ilevm2, ilevh
      real zcons1, zcons2,zcons3, zhgeo
      real zu, zphi, zvt1,zvt2, zst, zvar, zdwind, zwind
      real zstabm, zstabp, zrhom,  zrhop, alpha
      real zggeenv, zggeom1,zgvar 
      logical lo 
      logical ll1(klon,klev+1)
      integer kknu(klon),kknu2(klon),kknub(klon),kknul(klon),                        &
       &        kentp(klon),ncount(klon)  
!c
      real zhcrit(klon,klev),zvpf(klon,klev),                                        &
       &     zdp(klon,klev)
      real znorm(klon),zb(klon),zc(klon),                                            &
       &      zulow(klon),zvlow(klon),znup(klon),znum(klon)
!c
!c     ------------------------------------------------------------------
!c
!c*         1.    initialization
!c                --------------
!c
!c     print *,' entree gwsetup'
 100  continue
!c
!c     ------------------------------------------------------------------
!c
!c*         1.1   computational constants
!c                -----------------------
!c
 110  continue
!c
      ilevm1=klev-1
      ilevm2=klev-2
      ilevh =klev/3
!c
      zcons1=1./rd
!cold  zcons2=g**2/cpd
      zcons2=rg**2/rcpd
!cold  zcons3=1.5*api
      zcons3=1.5*rpi
!c
!c
!c     ------------------------------------------------------------------
!c
!c*         2.
!c                --------------
!c
 200  continue
!c
!c     ------------------------------------------------------------------
!c
!c*         2.1     define low level wind, project winds in plane of
!c*                 low level wind, determine sector in which to take
!c*                 the variance and set indicator for critical levels.
!c
!c
!c
      do 2001 jl=kidia,kfdia
      kknu(jl)    =klev
      kknu2(jl)   =klev
      kknub(jl)   =klev
      kknul(jl)   =klev
      pgamma(jl) =max(pgamma(jl),gtsec)
      ll1(jl,klev+1)=.false.
 2001 continue
!c
!c Ajouter une initialisation (L. Li, le 23fev99):
!c
      do jk=klev,ilevh,-1
      do jl=kidia,kfdia
      ll1(jl,jk)= .FALSE.
      ENDDO
      ENDDO
!c
!c*      define top of low level flow
!c       ----------------------------
      do 2002 jk=klev,ilevh,-1
      do 2003 jl=kidia,kfdia
      lo=(paphm1(jl,jk)/paphm1(jl,klev+1)).ge.gsigcr
      if(lo) then
        kkcrit(jl)=jk
      endif
      zhcrit(jl,jk)=ppic(jl)
      zhgeo=pgeom1(jl,jk)/rg
      ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
      if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
        kknu(jl)=jk
      endif
      if(.not.ll1(jl,ilevh))kknu(jl)=ilevh
 2003 continue
 2002 continue
      do 2004 jk=klev,ilevh,-1
      do 2005 jl=kidia,kfdia
      zhcrit(jl,jk)=ppic(jl)-pval(jl)
      zhgeo=pgeom1(jl,jk)/rg
      ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
      if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
        kknu2(jl)=jk
      endif
      if(.not.ll1(jl,ilevh))kknu2(jl)=ilevh
 2005 continue
 2004 continue
      do 2006 jk=klev,ilevh,-1
      do 2007 jl=kidia,kfdia
      zhcrit(jl,jk)=amax1(ppic(jl)-pmea(jl),pmea(jl)-pval(jl))
      zhgeo=pgeom1(jl,jk)/rg
      ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
      if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
        kknub(jl)=jk
      endif
      if(.not.ll1(jl,ilevh))kknub(jl)=ilevh
 2007 continue
 2006 continue
!c
      do 2010 jl=kidia,kfdia  
      kknu(jl)=min(kknu(jl),nktopg)
      kknu2(jl)=min(kknu2(jl),nktopg)
      kknub(jl)=min(kknub(jl),nktopg)
      kknul(jl)=klev
 2010 continue      
!c

 210  continue
!c
!c
!cc*     initialize various arrays
!c
      do 2107 jl=kidia,kfdia
      prho(jl,klev+1)  =0.0
      pstab(jl,klev+1) =0.0
      pstab(jl,1)      =0.0
      pri(jl,klev+1)   =9999.0
      ppsi(jl,klev+1)  =0.0
      pri(jl,1)        =0.0
      pvph(jl,1)       =0.0
      pulow(jl)        =0.0
      pvlow(jl)        =0.0
      zulow(jl)        =0.0
      zvlow(jl)        =0.0
      kkcrith(jl)      =klev
      kkenvh(jl)       =klev
      kentp(jl)        =klev
      kcrit(jl)        =1
      ncount(jl)       =0
      ll1(jl,klev+1)   =.false.
 2107 continue
!c
!c*     define low-level flow
!c      ---------------------
!c
      do 223 jk=klev,2,-1
      do 222 jl=kidia,kfdia
      if(ktest(jl).eq.1) then
        zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
        prho(jl,jk)=2.*paphm1(jl,jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
        pstab(jl,jk)=2.*zcons2/(ptm1(jl,jk)+ptm1(jl,jk-1))*                          &
       &  (1.-rcpd*prho(jl,jk)*(ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk))
        pstab(jl,jk)=max(pstab(jl,jk),gssec)
      endif
  222 continue
  223 continue
!c
!c********************************************************************
!c
!c*     define blocked flow
!c      -------------------
      do 2115 jk=klev,ilevh,-1
      do 2116 jl=kidia,kfdia
      if(jk.ge.kknub(jl).and.jk.le.kknul(jl)) then
        pulow(jl)=pulow(jl)+pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
        pvlow(jl)=pvlow(jl)+pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
      end if
 2116 continue
 2115 continue
      do 2110 jl=kidia,kfdia
      pulow(jl)=pulow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknub(jl)))
      pvlow(jl)=pvlow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknub(jl)))
      znorm(jl)=max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
      pvph(jl,klev+1)=znorm(jl)
 2110 continue
!c
!c*******  setup orography axes and define plane of profiles  *******
!c
      do 2112 jl=kidia,kfdia
      lo=(pulow(jl).lt.gvsec).and.(pulow(jl).ge.-gvsec)
      if(lo) then
        zu=pulow(jl)+2.*gvsec
      else
        zu=pulow(jl)
      endif
      zphi=atan(pvlow(jl)/zu)
      ppsi(jl,klev+1)=ptheta(jl)*rpi/180.-zphi
      zb(jl)=1.-0.18*pgamma(jl)-0.04*pgamma(jl)**2
      zc(jl)=0.48*pgamma(jl)+0.3*pgamma(jl)**2
      pd1(jl)=zb(jl)-(zb(jl)-zc(jl))*(sin(ppsi(jl,klev+1))**2)
      pd2(jl)=(zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))*cos(ppsi(jl,klev+1))
      pdmod(jl)=sqrt(pd1(jl)**2+pd2(jl)**2)
 2112 continue
!c
!c  ************ define flow in plane of lowlevel stress *************
!c
      do 213 jk=1,klev
      do 212 jl=kidia,kfdia
      if(ktest(jl).eq.1)  then
        zvt1       =pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk)
        zvt2       =-pvlow(jl)*pum1(jl,jk)+pulow(jl)*pvm1(jl,jk)
        zvpf(jl,jk)=(zvt1*pd1(jl)+zvt2*pd2(jl))/(znorm(jl)*pdmod(jl))
      endif
      ptau(jl,jk)  =0.0
      pzdep(jl,jk) =0.0
      ppsi(jl,jk)  =0.0
      ll1(jl,jk)   =.false.
  212 continue
  213 continue
      do 215 jk=2,klev
      do 214 jl=kidia,kfdia
      if(ktest(jl).eq.1) then
        zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
        pvph(jl,jk)=((paphm1(jl,jk)-papm1(jl,jk-1))*zvpf(jl,jk)+                     &
       &            (papm1(jl,jk)-paphm1(jl,jk))*zvpf(jl,jk-1))                      &
       &            /zdp(jl,jk)
        if(pvph(jl,jk).lt.gvsec) then
          pvph(jl,jk)=gvsec
          kcrit(jl)=jk
        endif
      endif
  214 continue
  215 continue
!c
!c
!c*         2.2     brunt-vaisala frequency and density at half levels.
!c
  220 continue
!c
      do 2211 jk=ilevh,klev
      do 221 jl=kidia,kfdia
      if(ktest(jl).eq.1) then
      if(jk.ge.(kknub(jl)+1).and.jk.le.kknul(jl)) then
           zst=zcons2/ptm1(jl,jk)*(1.-rcpd*prho(jl,jk)*                              &
       &                   (ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk))
           pstab(jl,klev+1)=pstab(jl,klev+1)+zst*zdp(jl,jk)
           pstab(jl,klev+1)=max(pstab(jl,klev+1),gssec)
           prho(jl,klev+1)=prho(jl,klev+1)+paphm1(jl,jk)*2.*zdp(jl,jk)               &
       &                   *zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
      endif
      endif
  221 continue
 2211 continue
!c
      do 2212 jl=kidia,kfdia
        pstab(jl,klev+1)=pstab(jl,klev+1)/(papm1(jl,kknul(jl))                       &
       &                                          -papm1(jl,kknub(jl)))
        prho(jl,klev+1)=prho(jl,klev+1)/(papm1(jl,kknul(jl))                         &
       &                                          -papm1(jl,kknub(jl)))
        zvar=pstd(jl)
 2212 continue
!c
!c*         2.3     mean flow richardson number.
!c*                 and critical height for froude layer
!c
  230 continue
!c
      do 232 jk=2,klev
      do 231 jl=kidia,kfdia
      if(ktest(jl).eq.1) then
        zdwind=max(abs(zvpf(jl,jk)-zvpf(jl,jk-1)),gvsec)
        pri(jl,jk)=pstab(jl,jk)*(zdp(jl,jk)                                          &
       &          /(rg*prho(jl,jk)*zdwind))**2
        pri(jl,jk)=max(pri(jl,jk),grcrit)
      endif
  231 continue
  232 continue
  
!c
!c
!c*      define top of 'envelope' layer
!c       ----------------------------

      do 233 jl=kidia,kfdia
      pnu (jl)=0.0
      znum(jl)=0.0
 233  continue
      
!      do 234 jk=2,klev-1
!      do 234 jl=kidia,kfdia
! L. Fita. July 2013
      DO jk=2,klev-1
      DO jl=kidia,kfdia

      
      if(ktest(jl).eq.1) then
       
      if (jk.ge.kknub(jl)) then
          
            znum(jl)=pnu(jl)
            zwind=(pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/                     &
       &            max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
            zwind=max(sqrt(zwind**2),gvsec)
            zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
            zstabm=sqrt(max(pstab(jl,jk  ),gssec))
            zstabp=sqrt(max(pstab(jl,jk+1),gssec))
            zrhom=prho(jl,jk  )
            zrhop=prho(jl,jk+1)
            pnu(jl) = pnu(jl) + (zdelp/rg)*                                          &
       &            ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind     
            if((znum(jl).le.gfrcrit).and.(pnu(jl).gt.gfrcrit)                        &
       &                          .and.(kkenvh(jl).eq.klev)) THEN                        
             kkenvh(jl)=jk                                                          
            endif    

      endif
      END IF
      
! 234  continue
      END DO
      END DO
      
!c  calculation of a dynamical mixing height for the breaking
!c  of gravity waves:

              
      do 235 jl=kidia,kfdia
      znup(jl)=0.0
      znum(jl)=0.0
 235  continue

!      do 236 jk=klev-1,2,-1
!      do 236 jl=kidia,kfdia
! L. Fita, July 2013
      DO jk=klev-1,2,-1
      DO jl=kidia,kfdia

      
      if(ktest(jl).eq.1) then

            znum(jl)=znup(jl)
            zwind=(pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/                     &
       &            max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
            zwind=max(sqrt(zwind**2),gvsec)
            zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
            zstabm=sqrt(max(pstab(jl,jk  ),gssec))
            zstabp=sqrt(max(pstab(jl,jk+1),gssec))
            zrhom=prho(jl,jk  )
            zrhop=prho(jl,jk+1)
            znup(jl) = znup(jl) + (zdelp/rg)*                                        &
       &            ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind     
            if((znum(jl).le.rpi/2.).and.(znup(jl).gt.rpi/2.)                         &
       &                          .and.(kkcrith(jl).eq.klev)) THEN                       
              kkcrith(jl)=jk                                                         
            endif
      
      END IF
! 236  continue
      END DO
      END DO
 
      do 237 jl=kidia,kfdia
      kkcrith(jl)=min0(kkcrith(jl),kknu2(jl))
      kkcrith(jl)=max0(kkcrith(jl),ilevh*2)
 237  continue         
!c
!c     directional info for flow blocking ************************* 
!c
      do 251 jk=ilevh,klev    
      do 252 jl=kidia,kfdia
      if(jk.ge.kkenvh(jl)) then
      lo=(pum1(jl,jk).lt.gvsec).and.(pum1(jl,jk).ge.-gvsec)
      if(lo) then
        zu=pum1(jl,jk)+2.*gvsec
      else
        zu=pum1(jl,jk)
      endif
       zphi=atan(pvm1(jl,jk)/zu)
       ppsi(jl,jk)=ptheta(jl)*rpi/180.-zphi
      end if
 252  continue
 251  continue
!c      forms the vertical 'leakiness' **************************

      alpha=3.
      
      do 254  jk=ilevh,klev
      do 253  jl=kidia,kfdia
      if(jk.ge.kkenvh(jl)) then
        zggeenv=amax1(1.,                                                            &
       &          (pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)-1))/2.)      
        zggeom1=amax1(pgeom1(jl,jk),1.)
        zgvar=amax1(pstd(jl)*rg,1.)     
!cmod    pzdep(jl,jk)=sqrt((zggeenv-zggeom1)/(zggeom1+zgvar))      
        pzdep(jl,jk)=(pgeom1(jl,kkenvh(jl)-1)-pgeom1(jl,  jk))/                      &
       &               (pgeom1(jl,kkenvh(jl)-1)-pgeom1(jl,klev))
      end if
 253  continue
 254  continue

 260  continue

      return
      end SUBROUTINE orosetup

      SUBROUTINE gwstress                                                            &
       &         (  nlon  , nlev                                                     &
       &         , ktest, kcrit, kkenvh                                              &
       &         , kknu                                                              &
       &         , prho  , pstab , pvph  , pstd, psig                                &
       &         , pmea , ppic  , ptau                                               &
       &         , pgeom1 , pdmod )
!c
!c**** *gwstress*
!c
!c     purpose.
!c     --------
!c
!c**   interface.
!c     ----------
!c     call *gwstress*  from *gwdrag*
!c
!c        explicit arguments :
!c        --------------------
!c     ==== inputs ===
!c     ==== outputs ===
!c
!c        implicit arguments :   none
!c        --------------------
!c
!c     method.
!c     -------
!c
!c
!c     externals.
!c     ----------
!c
!c
!c     reference.
!c     ----------
!c
!c        see ecmwf research department documentation of the "i.f.s."
!c
!c     author.
!c     -------
!c
!c     modifications.
!c     --------------
!c     f. lott put the new gwd on ifs      22/11/93
!c
!c-----------------------------------------------------------------------
      USE dimphy
      implicit none
!cym#include "dimensions.h"
!cym#include "dimphy.h"
#include "YOMCST.h"
#include "YOEGWD.h"

!c-----------------------------------------------------------------------
!c
!c*       0.1   arguments
!c              ---------
!c
      integer nlon, nlev
      integer kcrit(nlon),                                                           &
       &        ktest(nlon),kkenvh(nlon),kknu(nlon)
!c
      real prho(nlon,nlev+1),pstab(nlon,nlev+1),ptau(nlon,nlev+1),                   &
       &     pvph(nlon,nlev+1),                                                      &
       &     pgeom1(nlon,nlev),pstd(nlon)
!c
      real psig(nlon)
      real pmea(nlon),ppic(nlon)
      real pdmod(nlon)
!c
!c-----------------------------------------------------------------------
!c
!c*       0.2   local arrays
!c              ------------
      integer jl
      real zblock, zvar, zeff
      logical lo
!c
!c-----------------------------------------------------------------------
!c
!c*       0.3   functions
!c              ---------
!c     ------------------------------------------------------------------
!c
!c*         1.    initialization
!c                --------------
!c
 100  continue
!c
!c*         3.1     gravity wave stress.
!c
  300 continue
!c
!c
      do 301 jl=kidia,kfdia
      if(ktest(jl).eq.1) then
      
!c  effective mountain height above the blocked flow
  
         if(kkenvh(jl).eq.klev)then
         zblock=0.0 
         else
         zblock=(pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)+1))/2./rg          
         endif
      
        zvar=ppic(jl)-pmea(jl)
        zeff=amax1(0.,zvar-zblock)

        ptau(jl,klev+1)=prho(jl,klev+1)*gkdrag*psig(jl)*zeff**2                      &
       &    /4./pstd(jl)*pvph(jl,klev+1)*pdmod(jl)*sqrt(pstab(jl,klev+1))

!c  too small value of stress or  low level flow include critical level
!c  or low level flow:  gravity wave stress nul.
                
        lo=(ptau(jl,klev+1).lt.gtsec).or.(kcrit(jl).ge.kknu(jl))                     &
       &      .or.(pvph(jl,klev+1).lt.gvcrit)
!c       if(lo) ptau(jl,klev+1)=0.0
      
      else
      
          ptau(jl,klev+1)=0.0
          
      endif
      
  301 continue
!c
      return
      end SUBROUTINE gwstress

      SUBROUTINE GWPROFIL                                                            &
       &         ( NLON, NLEV                                                        &
       &         , kgwd, kdx , ktest                                                 &
       &         , KKCRITH, KCRIT                                                    &
       &         , PAPHM1, PRHO   , PSTAB  , PVPH , PRI , PTAU                       &
       &         , pdmod   , psig , pvar)

!C**** *GWPROFIL*
!C
!C     PURPOSE.
!C     --------
!C
!C**   INTERFACE.
!C     ----------
!C          FROM *GWDRAG*
!C
!C        EXPLICIT ARGUMENTS :
!C        --------------------
!C     ==== INPUTS ===
!C     ==== OUTPUTS ===
!C
!C        IMPLICIT ARGUMENTS :   NONE
!C        --------------------
!C
!C     METHOD:
!C     -------
!C     THE STRESS PROFILE FOR GRAVITY WAVES IS COMPUTED AS FOLLOWS:
!C     IT IS CONSTANT (NO GWD) AT THE LEVELS BETWEEN THE GROUND
!C     AND THE TOP OF THE BLOCKED LAYER (KKENVH).
!C     IT DECREASES LINEARLY WITH HEIGHTS FROM THE TOP OF THE 
!C     BLOCKED LAYER TO 3*VAROR (kKNU), TO SIMULATES LEE WAVES OR 
!C     NONLINEAR GRAVITY WAVE BREAKING.
!C     ABOVE IT IS CONSTANT, EXCEPT WHEN THE WAVE ENCOUNTERS A CRITICAL
!C     LEVEL (KCRIT) OR WHEN IT BREAKS.
!C     
!C
!C
!C     EXTERNALS.
!C     ----------
!C
!C
!C     REFERENCE.
!C     ----------
!C
!C        SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
!C
!C     AUTHOR.
!C     -------
!C
!C     MODIFICATIONS.
!C     --------------
!C     PASSAGE OF THE NEW GWDRAG TO I.F.S. (F. LOTT, 22/11/93)
!C-----------------------------------------------------------------------
      USE dimphy
      implicit none
!C

!C

!cym#include "dimensions.h"
!cym#include "dimphy.h"
#include "YOMCST.h"
#include "YOEGWD.h"

!C-----------------------------------------------------------------------
!C
!C*       0.1   ARGUMENTS
!C              ---------
!C
      integer nlon,nlev
      INTEGER KKCRITH(NLON),KCRIT(NLON)                                              &
       &       ,kdx(nlon) , ktest(nlon)

!C
      REAL PAPHM1(NLON,NLEV+1), PSTAB(NLON,NLEV+1),                                  &
       &     PRHO  (NLON,NLEV+1), PVPH (NLON,NLEV+1),                                &
       &     PRI   (NLON,NLEV+1), PTAU(NLON,NLEV+1)                                  
      REAL pdmod (NLON) , psig(NLON),                                                &
       &     pvar(NLON)                                                              
!C-----------------------------------------------------------------------
!C
!C*       0.2   LOCAL ARRAYS
!C              ------------
!C
      integer ilevh, ji, kgwd, jl, jk
      real zsqr, zalfa, zriw, zdel, zb, zalpha,zdz2n
      real zdelp, zdelpt 
      REAL ZDZ2 (KLON,KLEV) , ZNORM(KLON) , zoro(KLON)
      REAL ZTAU (KLON,KLEV+1)
!C
!C-----------------------------------------------------------------------
!C
!C*         1.    INITIALIZATION
!C                --------------
!C
!c      print *,' entree gwprofil' 
 100  CONTINUE
!C
!C
!C*    COMPUTATIONAL CONSTANTS.
!C     ------------- ----------
!C
      ilevh=KLEV/3
!C
!c     DO 400 ji=1,kgwd
!c     jl=kdx(ji)
!c  Modif vectorisation 02/04/2004
      DO 400 jl=kidia,kfdia
      if (ktest(jl).eq.1) then
      Zoro(JL)=Psig(JL)*Pdmod(JL)/4./max(pvar(jl),1.0)
      ZTAU(JL,KLEV+1)=PTAU(JL,KLEV+1)
      endif
  400 CONTINUE
  
!C
      DO 430 JK=KLEV,2,-1
!C
!C
!C*         4.1    CONSTANT WAVE STRESS UNTIL TOP OF THE
!C                 BLOCKING LAYER.
  410 CONTINUE
!C
!c     DO 411 ji=1,kgwd
!c     jl=kdx(ji)
!c  Modif vectorisation 02/04/2004
      do 411 jl=kidia,kfdia
      if (ktest(jl).eq.1) then
           IF(JK.GT.KKCRITH(JL)) THEN
           PTAU(JL,JK)=ZTAU(JL,KLEV+1)
!C          ENDIF
!C          IF(JK.EQ.KKCRITH(JL)) THEN
           ELSE                    
           PTAU(JL,JK)=GRAHILO*ZTAU(JL,KLEV+1)
           ENDIF
      endif
 411  CONTINUE             
!C
!C*         4.15   CONSTANT SHEAR STRESS UNTIL THE TOP OF THE
!C                 LOW LEVEL FLOW LAYER.
 415  CONTINUE
!C        
!C
!C*         4.2    WAVE DISPLACEMENT AT NEXT LEVEL.
!C
  420 CONTINUE
!C
!c     DO 421 ji=1,kgwd
!c     jl=kdx(ji)
!c  Modif vectorisation 02/04/2004
      do 421 jl=kidia,kfdia
      if(ktest(jl).eq.1) then
      IF(JK.LT.KKCRITH(JL)) THEN
      ZNORM(JL)=gkdrag*PRHO(JL,JK)*SQRT(PSTAB(JL,JK))*PVPH(JL,JK)                    &
       &                                                    *zoro(jl)
      ZDZ2(JL,JK)=PTAU(JL,JK+1)/max(ZNORM(JL),gssec)
      ENDIF
      endif
  421 CONTINUE
!C
!C*         4.3    WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT
!C*                AND STRESS:  BREAKING EVALUATION AND CRITICAL 
!C                 LEVEL
!C
                          
!c     DO 431 ji=1,kgwd
!c     jl=Kdx(ji)
!c  Modif vectorisation 02/04/2004
      do 431 jl=kidia,kfdia
      if(ktest(jl).eq.1) then

          IF(JK.LT.KKCRITH(JL)) THEN
          IF((PTAU(JL,JK+1).LT.GTSEC).OR.(JK.LE.KCRIT(JL))) THEN
            PTAU(JL,JK)=0.0
          ELSE
               ZSQR=SQRT(PRI(JL,JK))
               ZALFA=SQRT(PSTAB(JL,JK)*ZDZ2(JL,JK))/PVPH(JL,JK)
               ZRIW=PRI(JL,JK)*(1.-ZALFA)/(1+ZALFA*ZSQR)**2
               IF(ZRIW.LT.GRCRIT) THEN
                 ZDEL=4./ZSQR/GRCRIT+1./GRCRIT**2+4./GRCRIT
                 ZB=1./GRCRIT+2./ZSQR
                 ZALPHA=0.5*(-ZB+SQRT(ZDEL))
                 ZDZ2N=(PVPH(JL,JK)*ZALPHA)**2/PSTAB(JL,JK)
                 PTAU(JL,JK)=ZNORM(JL)*ZDZ2N
               ELSE
                 PTAU(JL,JK)=ZNORM(JL)*ZDZ2(JL,JK)
               ENDIF
            PTAU(JL,JK)=MIN(PTAU(JL,JK),PTAU(JL,JK+1))
          ENDIF
          ENDIF
      endif
  431 CONTINUE
  
  430 CONTINUE
  440 CONTINUE
  
!C  REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL

!c     DO 530 ji=1,kgwd
!c     jl=kdx(ji)
!c  Modif vectorisation 02/04/2004
      do 530 jl=kidia,kfdia
      if(ktest(jl).eq.1) then
      ZTAU(JL,KKCRITH(JL))=PTAU(JL,KKCRITH(JL))
      ZTAU(JL,NSTRA)=PTAU(JL,NSTRA)
      endif
 530  CONTINUE      

      DO 531 JK=1,KLEV
      
!c     DO 532 ji=1,kgwd
!c     jl=kdx(ji)
!c  Modif vectorisation 02/04/2004
      do 532 jl=kidia,kfdia
      if(ktest(jl).eq.1) then

                
         IF(JK.GT.KKCRITH(JL))THEN

          ZDELP=PAPHM1(JL,JK)-PAPHM1(JL,KLEV+1    )
          ZDELPT=PAPHM1(JL,KKCRITH(JL))-PAPHM1(JL,KLEV+1    )
          PTAU(JL,JK)=ZTAU(JL,KLEV+1    ) +                                          &
       &                (ZTAU(JL,KKCRITH(JL))-ZTAU(JL,KLEV+1    ) )*                 &
       &                ZDELP/ZDELPT
         ENDIF
            
      endif
 532  CONTINUE    
 
!C  REORGANISATION IN THE STRATOSPHERE

!c     DO 533 ji=1,kgwd
!c     jl=kdx(ji)
!c  Modif vectorisation 02/04/2004
      do 533 jl=kidia,kfdia
      if(ktest(jl).eq.1) then


         IF(JK.LT.NSTRA)THEN

          ZDELP =PAPHM1(JL,NSTRA)
          ZDELPT=PAPHM1(JL,JK)
          PTAU(JL,JK)=ZTAU(JL,NSTRA)*ZDELPT/ZDELP 

        ENDIF

      endif
 533  CONTINUE

!C REORGANISATION IN THE TROPOSPHERE

!c      DO 534 ji=1,kgwd
!c      jl=kdx(ji)
!c  Modif vectorisation 02/04/2004
      do 534 jl=kidia,kfdia
      if(ktest(jl).eq.1) then


         IF(JK.LT.KKCRITH(JL).AND.JK.GT.NSTRA)THEN

           ZDELP=PAPHM1(JL,JK)-PAPHM1(JL,KKCRITH(JL))
           ZDELPT=PAPHM1(JL,NSTRA)-PAPHM1(JL,KKCRITH(JL))
           PTAU(JL,JK)=ZTAU(JL,KKCRITH(JL)) +                                        &
       &                 (ZTAU(JL,NSTRA)-ZTAU(JL,KKCRITH(JL)))*ZDELP                 &
       &                                                     /ZDELPT

       ENDIF
      endif
 534   CONTINUE

 
 531  CONTINUE        


      RETURN
      END SUBROUTINE GWPROFIL

      SUBROUTINE lift_noro (nlon,nlev,dtime,paprs,pplay,                             &
       &                   plat,pmea,pstd, ppic,                                     &
       &                   ktest,                                                    &
       &                   t, u, v,                                                  &
       &                   pulow, pvlow, pustr, pvstr,                               &
       &                   d_t, d_u, d_v)
!c
      USE dimphy
      IMPLICIT none
!c======================================================================
!c Auteur(s): F.Lott (LMD/CNRS) date: 19950201
!c Objet: Frottement de la montagne Interface
!c======================================================================
!c Arguments:
!c dtime---input-R- pas d'integration (s)
!c paprs---input-R-pression pour chaque inter-couche (en Pa)
!c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
!c t-------input-R-temperature (K)
!c u-------input-R-vitesse horizontale (m/s)
!c v-------input-R-vitesse horizontale (m/s)
!c
!c d_t-----output-R-increment de la temperature
!c d_u-----output-R-increment de la vitesse u
!c d_v-----output-R-increment de la vitesse v
!c======================================================================
!cym#include "dimensions.h"
!cym#include "dimphy.h"
#include "YOMCST.h"
!c
!c ARGUMENTS
!c
      INTEGER nlon,nlev
      REAL dtime
      REAL paprs(klon,klev+1)
      REAL pplay(klon,klev)
      REAL plat(nlon),pmea(nlon)
      REAL pstd(nlon)
      REAL ppic(nlon)
      REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon)
      REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev)
      REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev)
!c
      INTEGER i, k, ktest(nlon)
!c
!c Variables locales:
!c
      REAL zgeom(klon,klev)
      REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev)
      REAL pt(klon,klev), pu(klon,klev), pv(klon,klev)
      REAL papmf(klon,klev),papmh(klon,klev+1)
!c
!c initialiser les variables de sortie (pour securite)
!c
      DO i = 1,klon
         pulow(i) = 0.0
         pvlow(i) = 0.0
         pustr(i) = 0.0
         pvstr(i) = 0.0
      ENDDO
      DO k = 1, klev
      DO i = 1, klon
         d_t(i,k) = 0.0
         d_u(i,k) = 0.0
         d_v(i,k) = 0.0
         pdudt(i,k)=0.0
         pdvdt(i,k)=0.0
         pdtdt(i,k)=0.0
      ENDDO
      ENDDO
!c
!c preparer les variables d'entree (attention: l'ordre des niveaux 
!c verticaux augmente du haut vers le bas)
!c
      DO k = 1, klev
      DO i = 1, klon
         pt(i,k) = t(i,klev-k+1) 
         pu(i,k) = u(i,klev-k+1)
         pv(i,k) = v(i,klev-k+1)
         papmf(i,k) = pplay(i,klev-k+1)
      ENDDO
      ENDDO
      DO k = 1, klev+1
      DO i = 1, klon
         papmh(i,k) = paprs(i,klev-k+2)
      ENDDO
      ENDDO
      DO i = 1, klon
         zgeom(i,klev) = RD * pt(i,klev)                                             &
       &                  * LOG(papmh(i,klev+1)/papmf(i,klev))
      ENDDO
      DO k = klev-1, 1, -1
      DO i = 1, klon
         zgeom(i,k) = zgeom(i,k+1) + RD * (pt(i,k)+pt(i,k+1))/2.0                    &
       &               * LOG(papmf(i,k+1)/papmf(i,k))
      ENDDO
      ENDDO
!c
!c appeler la routine principale
!c
      CALL OROLIFT(klon,klev,ktest,                                                  &
       &            dtime,                                                           &
       &            papmh, zgeom,                                                    &
       &            pt, pu, pv,                                                      &
       &            plat,pmea, pstd, ppic,                                           &
       &            pulow,pvlow,                                                     &
       &            pdudt,pdvdt,pdtdt)
!C
      DO k = 1, klev
      DO i = 1, klon
         d_u(i,klev+1-k) = dtime*pdudt(i,k)
         d_v(i,klev+1-k) = dtime*pdvdt(i,k)
         d_t(i,klev+1-k) = dtime*pdtdt(i,k)
         pustr(i)        = pustr(i)                                                  &
!IM BUG .                 +RG*pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))                   &
       &                    +pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))/RG
         pvstr(i)        = pvstr(i)                                                  &
!IM BUG .                 +RG*pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))                   &
       &                    +pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))/RG
      ENDDO
      ENDDO
!c
      RETURN
      END SUBROUTINE lift_noro

      SUBROUTINE OROLIFT( NLON,NLEV                                                  &
       &                 , KTEST                                                     &
       &                 , PTSPHY                                                    &
       &                 , PAPHM1,PGEOM1,PTM1,PUM1,PVM1                              &
       &                 , PLAT                                                      &
       &                 , PMEA, PVAROR, ppic                                        &
!C OUTPUTS                                                                           &
       &                 , PULOW,PVLOW                                               &
       &                 , PVOM,PVOL,PTE )

!C
!C**** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT.
!C
!C     PURPOSE.
!C     --------
!C
!C**   INTERFACE.
!C     ----------
!C          CALLED FROM *lift_noro
!C     ----------
!C
!C     AUTHOR.
!C     -------
!C     F.LOTT  LMD 22/11/95
!C
      USE dimphy
      implicit none
!C
!C
!cym#include "dimensions.h"
!cym#include "dimphy.h"
#include "YOMCST.h"
#include "YOEGWD.h"
!C-----------------------------------------------------------------------
!C
!C*       0.1   ARGUMENTS
!C              ---------
!C
!C
      integer nlon, nlev
      REAL  PTE(NLON,NLEV),                                                          &
       &      PVOL(NLON,NLEV),                                                       &
       &      PVOM(NLON,NLEV),                                                       &
       &      PULOW(NLON),                                                           &
       &      PVLOW(NLON)
      REAL  PUM1(NLON,NLEV),                                                         &
       &      PVM1(NLON,NLEV),                                                       &
       &      PTM1(NLON,NLEV),                                                       &
       &      PLAT(NLON),PMEA(NLON),                                                 &
       &      PVAROR(NLON),                                                          &
       &      ppic(NLON),                                                            &
       &      PGEOM1(NLON,NLEV),                                                     &
       &      PAPHM1(NLON,NLEV+1)
!C
      INTEGER  KTEST(NLON)
      real ptsphy
!C-----------------------------------------------------------------------
!C
!C*       0.2   LOCAL ARRAYS
!C              ------------
      logical lifthigh
!cym      integer klevm1, jl, ilevh, jk
      integer  jl, ilevh, jk
      real zcons1, ztmst, zrtmst,zpi, zhgeo
      real zdelp, zslow, zsqua, zscav, zbet
      INTEGER                                                                        &
       &         IKNUB(klon),                                                        &
       &         IKNUL(klon)
      LOGICAL LL1(KLON,KLEV+1)
!C
      REAL   ZTAU(KLON,KLEV+1),                                                      &
       &       ZTAV(KLON,KLEV+1),                                                    &
       &       ZRHO(KLON,KLEV+1)
      REAL   ZDUDT(KLON),                                                            &
       &       ZDVDT(KLON)
      REAL ZHCRIT(KLON,KLEV)
      CHARACTER (LEN=20) :: modname='orografi'
      CHARACTER (LEN=80) :: abort_message
!C-----------------------------------------------------------------------
!C
!C*         1.1  INITIALIZATIONS
!C               ---------------

      LIFTHIGH=.FALSE.

      IF(NLON.NE.KLON.OR.NLEV.NE.KLEV)THEN
        abort_message = 'pb dimension'
        CALL abort_gcm (modname,abort_message,1)
      ENDIF
      ZCONS1=1./RD
!cym      KLEVM1=KLEV-1
      ZTMST=PTSPHY
      ZRTMST=1./ZTMST
      ZPI=ACOS(-1.)
!C
      DO 1001 JL=kidia,kfdia
      ZRHO(JL,KLEV+1)  =0.0
      PULOW(JL)        =0.0
      PVLOW(JL)        =0.0
      iknub(JL)   =klev
      iknul(JL)   =klev
      ilevh=klev/3
      ll1(jl,klev+1)=.false.
      DO 1000 JK=1,KLEV
      PVOM(JL,JK)=0.0
      PVOL(JL,JK)=0.0
      PTE (JL,JK)=0.0
 1000 CONTINUE
 1001 CONTINUE

!C
!C*         2.1     DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
!C*                 LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
!C*                 THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
!C
!C
!C
      DO 2006 JK=KLEV,1,-1
      DO 2007 JL=kidia,kfdia
      IF(KTEST(JL).EQ.1) THEN
      ZHCRIT(JL,JK)=amax1(Ppic(JL)-pmea(JL),100.)
      ZHGEO=PGEOM1(JL,JK)/RG
      ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
      IF(ll1(JL,JK).neqv.ll1(JL,JK+1)) THEN
        iknub(JL)=JK
      ENDIF
      ENDIF
 2007 CONTINUE
 2006 CONTINUE
!C
      do 2010 jl=kidia,kfdia
      IF(KTEST(JL).EQ.1) THEN
      iknub(jl)=max(iknub(jl),klev/2)
      iknul(jl)=max(iknul(jl),2*klev/3)
      if(iknub(jl).gt.nktopg) iknub(jl)=nktopg
      if(iknub(jl).eq.nktopg) iknul(jl)=klev
      if(iknub(jl).eq.iknul(jl)) iknub(jl)=iknul(jl)-1
      ENDIF
 2010 continue

!C     do 2011 jl=kidia,kfdia
!C     IF(KTEST(JL).EQ.1) THEN
!C       print *,' iknul= ',iknul(jl),'  iknub=',iknub(jl)
!C     ENDIF
!C2011 continue

!C     PRINT *,'  DANS OROLIFT: 2010'

      DO 223 JK=KLEV,2,-1
      DO 222 JL=kidia,kfdia
        ZRHO(JL,JK)=2.*PAPHM1(JL,JK)*ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1))
  222 CONTINUE
  223 CONTINUE
!C     PRINT *,'  DANS OROLIFT: 223'

!C********************************************************************
!C
!C*     DEFINE LOW LEVEL FLOW
!C      -------------------
      DO 2115 JK=klev,1,-1
      DO 2116 JL=kidia,kfdia
      IF(KTEST(JL).EQ.1) THEN
      if(jk.ge.iknub(jl).and.jk.le.iknul(jl)) then
        pulow(JL)=pulow(JL)+PUM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
        pvlow(JL)=pvlow(JL)+PVM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
        zrho(JL,klev+1)=zrho(JL,klev+1)                                              &
       &                 +zrho(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
      end if
      ENDIF
 2116 CONTINUE
 2115 CONTINUE
      DO 2110 JL=kidia,kfdia
      IF(KTEST(JL).EQ.1) THEN
      pulow(JL)=pulow(JL)/(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl)))
      pvlow(JL)=pvlow(JL)/(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl)))
      zrho(JL,klev+1)=zrho(JL,klev+1)                                                &
       &               /(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl)))
      ENDIF
 2110 CONTINUE


200   CONTINUE

!C***********************************************************
!C
!C*         3.      COMPUTE MOUNTAIN LIFT
!C
  300 CONTINUE
!C
      DO 301 JL=kidia,kfdia
      IF(KTEST(JL).EQ.1) THEN
       ZTAU(JL,KLEV+1)= - GKLIFT*ZRHO(JL,KLEV+1)*2.*ROMEGA*                          &
!C    *                 (2*PVAROR(JL)+PMEA(JL))*                                     &
       &                 2*PVAROR(JL)*                                               &
       &                 SIN(ZPI/180.*PLAT(JL))*PVLOW(JL)
       ZTAV(JL,KLEV+1)=   GKLIFT*ZRHO(JL,KLEV+1)*2.*ROMEGA*                          &
!C    *                 (2*PVAROR(JL)+PMEA(JL))*                                     &
       &                 2*PVAROR(JL)*                                               &
       &                 SIN(ZPI/180.*PLAT(JL))*PULOW(JL)
      ELSE
       ZTAU(JL,KLEV+1)=0.0
       ZTAV(JL,KLEV+1)=0.0
      ENDIF
301   CONTINUE

!C
!C*         4.      COMPUTE LIFT PROFILE         
!C*                 --------------------   
!C

  400 CONTINUE

      DO 401 JK=1,KLEV
      DO 401 JL=kidia,kfdia
      IF(KTEST(JL).EQ.1) THEN
      ZTAU(JL,JK)=ZTAU(JL,KLEV+1)*PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1)
      ZTAV(JL,JK)=ZTAV(JL,KLEV+1)*PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1)
      ELSE
      ZTAU(JL,JK)=0.0
      ZTAV(JL,JK)=0.0
      ENDIF
401   CONTINUE
!C
!C
!C*         5.      COMPUTE TENDENCIES.
!C*                 -------------------
      IF(LIFTHIGH)THEN
!C
  500 CONTINUE
!C     PRINT *,'  DANS OROLIFT: 500'
!C
!C  EXPLICIT SOLUTION AT ALL LEVELS
!C
      DO 524 JK=1,klev
      DO 523 JL=KIDIA,KFDIA
      IF(KTEST(JL).EQ.1) THEN
      ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK)
      ZDUDT(JL)=-RG*(ZTAU(JL,JK+1)-ZTAU(JL,JK))/ZDELP
      ZDVDT(JL)=-RG*(ZTAV(JL,JK+1)-ZTAV(JL,JK))/ZDELP
      ENDIF  
  523 CONTINUE
  524 CONTINUE
!C
!C  PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY
!C
      DO 530 JK=1,klev
      DO 530 JL=KIDIA,KFDIA
      IF(KTEST(JL).EQ.1) THEN

        ZSLOW=SQRT(PULOW(JL)**2+PVLOW(JL)**2)
        ZSQUA=AMAX1(SQRT(PUM1(JL,JK)**2+PVM1(JL,JK)**2),GVSEC)
        ZSCAV=-ZDUDT(JL)*PVM1(JL,JK)+ZDVDT(JL)*PUM1(JL,JK)
        IF(ZSQUA.GT.GVSEC)THEN
          PVOM(JL,JK)=-ZSCAV*PVM1(JL,JK)/ZSQUA**2
          PVOL(JL,JK)= ZSCAV*PUM1(JL,JK)/ZSQUA**2
        ELSE
          PVOM(JL,JK)=0.0
          PVOL(JL,JK)=0.0      
        ENDIF  
        ZSQUA=SQRT(PUM1(JL,JK)**2+PUM1(JL,JK)**2)               
        IF(ZSQUA.LT.ZSLOW)THEN
          PVOM(JL,JK)=ZSQUA/ZSLOW*PVOM(JL,JK)
          PVOL(JL,JK)=ZSQUA/ZSLOW*PVOL(JL,JK)
        ENDIF 

      ENDIF  
530   CONTINUE
!C
!C  6.  LOW LEVEL LIFT, SEMI IMPLICIT:
!C  ----------------------------------

      ELSE

        DO 601 JL=KIDIA,KFDIA
        IF(KTEST(JL).EQ.1) THEN
          DO JK=KLEV,IKNUB(JL),-1
          ZBET=GKLIFT*2.*ROMEGA*SIN(ZPI/180.*PLAT(JL))*ztmst*                        &
       &        (PGEOM1(JL,IKNUB(JL)-1)-PGEOM1(JL,  JK))/                            &
       &        (PGEOM1(JL,IKNUB(JL)-1)-PGEOM1(JL,KLEV))
          ZDUDT(JL)=-PUM1(JL,JK)/ztmst/(1+ZBET**2)
          ZDVDT(JL)=-PVM1(JL,JK)/ztmst/(1+ZBET**2)
          PVOM(JL,JK)= ZBET**2*ZDUDT(JL) - ZBET   *ZDVDT(JL)
          PVOL(JL,JK)= ZBET   *ZDUDT(JL) + ZBET**2*ZDVDT(JL)    
          ENDDO
        ENDIF
 601    CONTINUE

      ENDIF

      RETURN
      END SUBROUTINE OROLIFT


      SUBROUTINE SUGWD(NLON,NLEV,paprs,pplay)
      USE dimphy
      USE mod_phys_lmdz_para
      USE mod_grid_phy_lmdz
!c      USE parallel
!C
!C**** *SUGWD* INITIALIZE COMMON YOEGWD CONTROLLING GRAVITY WAVE DRAG
!C
!C     PURPOSE.
!C     --------
!C           INITIALIZE YOEGWD, THE COMMON THAT CONTROLS THE
!C           GRAVITY WAVE DRAG PARAMETRIZATION.
!C
!C**   INTERFACE.
!C     ----------
!C        CALL *SUGWD* FROM *SUPHEC*
!C              -----        ------
!C
!C        EXPLICIT ARGUMENTS :
!C        --------------------
!C        PSIG        : VERTICAL COORDINATE TABLE
!C        NLEV        : NUMBER OF MODEL LEVELS
!C
!C        IMPLICIT ARGUMENTS :
!C        --------------------
!C        COMMON YOEGWD
!C
!C     METHOD.
!C     -------
!C        SEE DOCUMENTATION
!C
!C     EXTERNALS.
!C     ----------
!C        NONE
!C
!C     REFERENCE.
!C     ----------
!C        ECMWF Research Department documentation of the IFS
!C
!C     AUTHOR.
!C     -------
!C        MARTIN MILLER             *ECMWF*
!C
!C     MODIFICATIONS.
!C     --------------
!C        ORIGINAL : 90-01-01
!C     ------------------------------------------------------------------
      implicit none
!C
!C     -----------------------------------------------------------------
#include "YOEGWD.h"
!C      ----------------------------------------------------------------
!C
      integer nlon,nlev, jk
      REAL paprs(nlon,nlev+1)
      REAL pplay(nlon,nlev)
      real zpr,zstra,zsigt,zpm1r
      REAL :: pplay_glo(klon_glo,nlev)
      REAL :: paprs_glo(klon_glo,nlev+1)

!C
!C*       1.    SET THE VALUES OF THE PARAMETERS
!C              --------------------------------
!C
 100  CONTINUE
!C
      PRINT *,' DANS SUGWD NLEV=',NLEV
      GHMAX=10000.
!C
      ZPR=100000.
      ZSTRA=0.1 
      ZSIGT=0.94
!cold  ZPR=80000.
!cold  ZSIGT=0.85
!C
      
      CALL gather(pplay,pplay_glo)
      CALL bcast(pplay_glo)
      CALL gather(paprs,paprs_glo)
      CALL bcast(paprs_glo)
      
            
      DO 110 JK=1,NLEV
      ZPM1R=pplay_glo((klon_glo/2)+1,jk)/paprs_glo((klon_glo/2)+1,1) 
      IF(ZPM1R.GE.ZSIGT)THEN
         nktopg=JK
      ENDIF
      ZPM1R=pplay_glo((klon_glo/2)+1,jk)/paprs_glo((klon_glo/2)+1,1) 
      IF(ZPM1R.GE.ZSTRA)THEN
         NSTRA=JK
      ENDIF
  110 CONTINUE


!c
!c  inversion car dans orodrag on compte les niveaux a l'envers
      nktopg=nlev-nktopg+1
      nstra=nlev-nstra
      print *,' DANS SUGWD nktopg=', nktopg
      print *,' DANS SUGWD nstra=', nstra
!C
      GSIGCR=0.80
!C
      GKDRAG=0.2 
      GRAHILO=1.    
      GRCRIT=0.01
      GFRCRIT=1.0
      GKWAKE=0.50 
!C
      GKLIFT=0.50  
      GVCRIT =0.0
!C
!C
!C      ----------------------------------------------------------------
!C
!C*       2.    SET VALUES OF SECURITY PARAMETERS
!C              ---------------------------------
!C
 200  CONTINUE
!C
      GVSEC=0.10
      GSSEC=1.E-12
!C
      GTSEC=1.E-07
!C
!C      ----------------------------------------------------------------
!C
      RETURN
      END SUBROUTINE SUGWD

END MODULE orografi