!$gpum horizontal klon
MODULE lmdz_gwd_ogwd
! ============================================================================
! This module contains the routine to compute the orographic gravity wave drag
! (ogwd) and lift effects
!============================================================================

   PRIVATE

   PUBLIC drag_noro_strato, orodrag_strato, orosetup_strato, gwstress_strato, &
        & lift_noro_strato, orolift_strato

CONTAINS
! ============================================================================
   SUBROUTINE drag_noro_strato(klon, klev, dtime, partdrag, 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)

      USE lmdz_gwd_ini, ONLY: rd, rg

      IMPLICIT NONE

      ! ======================================================================
      ! Auteur(s): F.Lott (LMD/CNRS) date: 19950201
      ! Object: Mountain drag interface. Made necessary because:
      ! 1. in the LMD-GCM Layers are from bottom to top,
      ! contrary to most European GCM.
      ! 2. the altitude above ground of each model layers
      ! needs to be known (variable zgeom)
      ! ======================================================================

      ! Declarations
      !====================================================

      ! Input variables
      !-----------------

      INTEGER, INTENT(IN) :: klon, klev ! horizontal and vertical dimensions
      REAL, INTENT(IN) ::  dtime ! timestep [s]
      INTEGER, INTENT(IN) :: partdrag ! controls which part of the drag we consider (total part or GW part)
      REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs ! pressure at bottom layer interface [Pa]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay ! pressure of model layers [Pa]
      REAL, DIMENSION(klon), INTENT(IN) :: pmea ! mean subgrid orography [m]
      REAL, DIMENSION(klon), INTENT(IN) :: pstd ! standard deviation of subgrid orography [m]
      REAL, DIMENSION(klon), INTENT(IN) :: psig ! sso mean slope [-]
      REAL, DIMENSION(klon), INTENT(IN) :: pthe ! sso orientation angle [radian]
      REAL, DIMENSION(klon), INTENT(IN) :: ppic ! elevation of sso peaks [m]
      REAL, DIMENSION(klon), INTENT(IN) :: pval ! elevation of sso valleys [m]
      INTEGER, DIMENSION(klon), INTENT(IN) :: kdx ! locate the physical location of an active point
      INTEGER, DIMENSION(klon), INTENT(IN) :: ktest ! flags to indicate active points
      INTEGER, INTENT(IN) :: kgwd ! total number of grid points where the param is active
      REAL, DIMENSION(klon, klev), INTENT(IN) :: t ! temperature [K]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: u ! zonal wind component [m/s]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: v ! meridional wind component [m/s]

      ! Inout variables
      REAL, DIMENSION(klon), INTENT(INOUT) :: pgam ! sso anisotropy [-]

      ! Output variables
      !------------------

      REAL, DIMENSION(klon), INTENT(OUT) :: pulow ! low level u [m/s]
      REAL, DIMENSION(klon), INTENT(OUT) :: pvlow ! low level v [m/s]
      REAL, DIMENSION(klon), INTENT(OUT) :: pustr ! surface stress due to SSO drag (x component) [Pa]
      REAL, DIMENSION(klon), INTENT(OUT) :: pvstr ! surface stress due to SSO drag (y component) [Pa]

      REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t ! temperature increment [K]
      REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_u  ! u increment [m/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_v  ! v increment [m/s]

      ! Local variables
      !----------------
      ! zgeom-----R: Altitude of layer above ground
      ! pt, pu, pv --R: t u v from top to bottom
      ! pdtdt, pdudt, pdvdt --R: t u v tendencies (from top to bottom)
      ! papmf: pressure at model layer (from top to bottom)
      ! papmh: pressure at model 1/2 layer (from top to bottom)

      INTEGER i, k
      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)
      CHARACTER(LEN=20) :: modname = 'orografi_strato'
      CHARACTER(LEN=80) :: abort_message

      !===========================================================

      ! INITIALIZE OUTPUT VARIABLES

      DO i = 1, klon
         pulow(i) = 0.0
         pvlow(i) = 0.0
         pustr(i) = 0.0
         pvstr(i) = 0.0
      END DO
      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
         END DO
      END DO

      ! PREPARE INPUT VARIABLES FOR ORODRAG (i.e., ORDERED FROM TOP TO BOTTOM)
      ! CALCULATE LAYERS HEIGHT ABOVE GROUND)

      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)
         END DO
      END DO
      DO k = 1, klev + 1
         DO i = 1, klon
            papmh(i, k) = paprs(i, klev - k + 2)
         END DO
      END DO
      DO i = 1, klon
         zgeom(i, klev) = rd*pt(i, klev)*log(papmh(i, klev + 1)/papmf(i, klev))
      END DO
      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))
         END DO
      END DO

      ! CALL SSO DRAG ROUTINES

      CALL orodrag_strato(klon, klev, partdrag, kgwd, kdx, ktest, dtime, papmh, papmf, &
                          zgeom, pt, pu, pv, pmea, pstd, psig, pgam, pthe, ppic, pval, pulow, &
                          pvlow, pdudt, pdvdt, pdtdt)

      ! COMPUTE INCREMENTS AND STRESS FROM TENDENCIES

      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) + pdudt(i, k)*(papmh(i, k + 1) - papmh(i, k))/rg
            pvstr(i) = pvstr(i) + pdvdt(i, k)*(papmh(i, k + 1) - papmh(i, k))/rg
         END DO
      END DO

      RETURN
   END SUBROUTINE drag_noro_strato

! ============================================================================
   SUBROUTINE orodrag_strato(klon, klev, partdrag, kgwd, kdx, ktest, ptsphy, paphm1, &
                             papm1, pgeom1, ptm1, pum1, pvm1, pmea, pstd, psig, pgam, pthe, ppic, pval &
                             , pulow, pvlow, pvom, pvol, pte)

      USE lmdz_gwd_ini, ONLY: nstra, rg, rcpd, gkwake

      IMPLICIT NONE

      ! **** *orodrag* - does the SSO drag  parametrization.

      ! purpose.
      ! --------

      ! this routine computes the physical tendencies of the
      ! prognostic variables u,v  and t due to  vertical transports by
      ! subgridscale orographically excited gravity waves, and to
      ! low level blocked flow drag.

      ! **   interface.
      ! ----------
      ! called from *drag_noro*.

      ! the routine takes its input from the long-term storage:
      ! u,v,t and p at t-1.

      ! method.
      ! -------

      ! externals.
      ! ----------
      INTEGER ismin, ismax
      EXTERNAL ismin, ismax

      ! reference.
      ! ----------

      ! author.
      ! -------
      ! m.miller + b.ritter   e.c.m.w.f.     15/06/86.

      ! f.lott + m. miller    e.c.m.w.f.     22/11/94
      ! -----------------------------------------------------------------------

      ! *       0.1   arguments
      ! ---------

      ! Input variables
      !-----------------

      INTEGER, INTENT(IN) :: klon, klev ! horizontal and vertical dimensions
      REAL, INTENT(IN) ::  ptsphy ! timestep [s]
      INTEGER, INTENT(IN) :: partdrag !control which part of the drag we consider (total part or GW part)
      REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paphm1 ! pressure at bottom layer interface [Pa]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: papm1 ! pressure of model layers [Pa]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: pgeom1 ! geopotential height above ground [m]
      REAL, DIMENSION(klon), INTENT(IN) :: pmea ! mean subgrid orography [m]
      REAL, DIMENSION(klon), INTENT(IN) :: pstd ! standard deviation of subgrid orography [m]
      REAL, DIMENSION(klon), INTENT(IN) :: psig ! sso mean slope [-]
      REAL, DIMENSION(klon), INTENT(IN) :: pthe ! sso orientation angle [radian]
      REAL, DIMENSION(klon), INTENT(IN) :: ppic ! elevation of sso peaks [m]
      REAL, DIMENSION(klon), INTENT(IN) :: pval ! elevation of sso valleys [m]
      INTEGER, INTENT(IN) :: kgwd ! total number of grid points where the param is active
      INTEGER, DIMENSION(klon), INTENT(IN) :: kdx ! locate the physical location of an active point
      INTEGER, DIMENSION(klon), INTENT(IN) :: ktest ! flags to indicate active points
      REAL, DIMENSION(klon, klev), INTENT(IN) :: ptm1 ! temperature [K]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: pum1 ! zonal wind component [m/s]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: pvm1 ! meridional wind component [m/s]

      ! Inout variables
      !-----------------

      REAL, DIMENSION(klon), INTENT(INOUT) :: pgam ! sso anisotropy [-]

      ! Output variables
      !------------------

      REAL, DIMENSION(klon), INTENT(OUT) :: pulow ! low level u [m/s]
      REAL, DIMENSION(klon), INTENT(OUT) :: pvlow ! low level v [m/s]

      REAL, DIMENSION(klon, klev), INTENT(OUT) :: pte ! temperature tendency [m/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT) :: pvol  ! u tendency [m/s/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT) :: pvom  ! v tendency [m/s/s]

      ! -----------------------------------------------------------------------

      ! *       0.2   local arrays
      ! ------------
      INTEGER isect(klon), icrit(klon), ikcrith(klon), ikenvh(klon), iknu(klon), &
         iknu2(klon), ikcrit(klon), ikhlim(klon)

      REAL ztau(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), &
         ztfr(klon), znu(klon), zd1(klon), zd2(klon), zdmod(klon)

      ! local quantities:

      INTEGER jl, jk, ji
      REAL ztmst, zdelp, ztemp, zforc, ztend, rover, facpart
      REAL zb, zc, zconb, zabsv, zzd1, ratio, zbet, zust, zvst, zdis

      ! ------------------------------------------------------------------

      ! *         1.    initialization
      ! --------------

      ! ------------------------------------------------------------------

      ! *         1.1   computational constants
      ! -----------------------

      ! ztmst=twodt
      ! if(nstep.eq.nstart) ztmst=0.5*twodt
      ztmst = ptsphy

      ! ------------------------------------------------------------------

      ! *         1.3   check whether row contains point for printing
      ! ---------------------------------------------

      ! ------------------------------------------------------------------

      ! *         2.     precompute basic state variables.
      ! *                ---------- ----- ----- ----------
      ! *                define low level wind, project winds in plane of
      ! *                low level wind, determine sector in which to take
      ! *                the variance and set indicator for critical levels.

      CALL orosetup_strato(klon, klev, ktest, ikcrit, ikcrith, icrit, isect, &
                           ikhlim, ikenvh, iknu, iknu2, paphm1, papm1, pum1, pvm1, ptm1, pgeom1, &
                           pstd, zrho, zri, zstab, ztau, zvph, zpsi, zzdep, pulow, pvlow, pthe, &
                           pgam, pmea, ppic, pval, znu, zd1, zd2, zdmod)

      ! ***********************************************************

      ! *         3.      compute low level stresses using subcritical and
      ! *                 supercritical forms.computes anisotropy coefficient
      ! *                 as measure of orographic twodimensionality.

      CALL gwstress_strato(klon, klev, ikcrit, isect, ikhlim, ktest, ikcrith, &
                           icrit, ikenvh, iknu, zrho, zstab, zvph, pstd, psig, pmea, ppic, pval, &
                           ztfr, ztau, pgeom1, pgam, zd1, zd2, zdmod, znu)

      ! *         4.      compute stress profile including
      ! trapped waves, wave breaking,
      ! linear decay in stratosphere.

      CALL gwprofil_strato(klon, klev, kgwd, kdx, ktest, ikcrit, ikcrith, icrit, &
                           ikenvh, iknu, iknu2, paphm1, zrho, zstab, ztfr, zvph, zri, ztau &
                           , zdmod, znu, psig, pgam, pstd, ppic, pval)

      ! *         5.      Compute tendencies from waves stress profile.
      ! Compute low level blocked flow drag.
      ! *                 --------------------------------------------

      ! explicit solution at all levels for the gravity wave
      ! implicit solution for the blocked levels

      DO jl = 1, klon
         zvidis(jl) = 0.0
         zdudt(jl) = 0.0
         zdvdt(jl) = 0.0
         zdtdt(jl) = 0.0
      END DO

      DO jk = 1, klev

         ! WAVE STRESS
         ! -------------

         DO ji = 1, klon

            IF (ktest(ji) == 1) THEN

               zdelp = paphm1(ji, jk + 1) - paphm1(ji, jk)
               ztemp = -rg*(ztau(ji, jk + 1) - ztau(ji, jk))/(zvph(ji, klev + 1)*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)

               ! Control Overshoots

               IF (jk >= nstra) THEN
                  rover = 0.10
                  IF (abs(zdudt(ji)) > rover*abs(pum1(ji, jk))/ztmst) zdudt(ji) = rover* &
                                                                         abs(pum1(ji, jk))/ztmst*zdudt(ji)/(abs(zdudt(ji)) + 1.E-10)
                  IF (abs(zdvdt(ji)) > rover*abs(pvm1(ji, jk))/ztmst) zdvdt(ji) = rover* &
                                                                         abs(pvm1(ji, jk))/ztmst*zdvdt(ji)/(abs(zdvdt(ji)) + 1.E-10)
               END IF

               rover = 0.25
               zforc = sqrt(zdudt(ji)**2 + zdvdt(ji)**2)
               ztend = sqrt(pum1(ji, jk)**2 + pvm1(ji, jk)**2)/ztmst

               IF (zforc >= rover*ztend) THEN
                  zdudt(ji) = rover*ztend/zforc*zdudt(ji)
                  zdvdt(ji) = rover*ztend/zforc*zdvdt(ji)
               END IF

               ! BLOCKED FLOW DRAG:
               ! -----------------

               IF (partdrag .GE. 2) THEN
                  facpart = 0.
               ELSE
                  facpart = gkwake
               END IF

               IF (jk > ikenvh(ji)) THEN
                  zb = 1.0 - 0.18*pgam(ji) - 0.04*pgam(ji)**2
                  zc = 0.48*pgam(ji) + 0.3*pgam(ji)**2
                  zconb = 2.*ztmst*facpart*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 + pgam(ji)*sin(zpsi(ji, &
                                                                    jk))**2)/(pgam(ji)*cos(zpsi(ji, jk))**2 + sin(zpsi(ji, jk))**2)
                  zbet = max(0., 2.-1./ratio)*zconb*zzdep(ji, jk)*zzd1*zabsv

                  ! OPPOSED TO THE WIND

                  zdudt(ji) = -pum1(ji, jk)/ztmst
                  zdvdt(ji) = -pvm1(ji, jk)/ztmst

                  ! PERPENDICULAR TO THE SSO MAIN AXIS:

                  ! mod     zdudt(ji)=-(pum1(ji,jk)*cos(pthe(ji)*rpi/180.)
                  ! mod *              +pvm1(ji,jk)*sin(pthe(ji)*rpi/180.))
                  ! mod *              *cos(pthe(ji)*rpi/180.)/ztmst
                  ! mod     zdvdt(ji)=-(pum1(ji,jk)*cos(pthe(ji)*rpi/180.)
                  ! mod *              +pvm1(ji,jk)*sin(pthe(ji)*rpi/180.))
                  ! mod *              *sin(pthe(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

               ! NO TENDENCIES ON TEMPERATURE .....

               ! Instead of, pte(ji,jk)=zdtdt(ji), due to mechanical dissipation

               pte(ji, jk) = 0.0

            END IF

         END DO
      END DO

      RETURN
   END SUBROUTINE orodrag_strato

! ============================================================================
   SUBROUTINE orosetup_strato(klon, klev, ktest, kkcrit, kkcrith, kcrit, ksect, &
                              kkhlim, kkenvh, kknu, kknu2, paphm1, papm1, pum1, pvm1, ptm1, pgeom1, &
                              pstd, prho, pri, pstab, ptau, pvph, ppsi, pzdep, pulow, pvlow, pthe, &
                              pgam, pmea, ppic, pval, pnu, pd1, pd2, pdmod)

      ! **** *gwsetup*

      ! purpose.
      ! --------
      ! SET-UP THE ESSENTIAL PARAMETERS OF THE SSO DRAG SCHEME:
      ! DEPTH OF LOW WBLOCKED LAYER, LOW-LEVEL FLOW, BACKGROUND
      ! STRATIFICATION.....

      ! **   interface.
      ! ----------
      ! from *orodrag*

      ! explicit arguments :
      ! --------------------
      ! ==== inputs ===

      ! klon----input-I-Total number of horizontal points that get into physics
      ! klev----input-I-Number of vertical levels
      ! ktest--input-I: Flags to indicate active points

      ! ptsphy--input-R-Time-step (s)
      ! paphm1--input-R: pressure at model 1/2 layer
      ! papm1---input-R: pressure at model layer
      ! pgeom1--input-R: Altitude of layer above ground
      ! ptm1, pum1, pvm1--R-: t, u and v
      ! pmea----input-R-Mean Orography (m)
      ! pstd----input-R-SSO standard deviation (m)
      ! psig----input-R-SSO slope
      ! pgam----input-R-SSO Anisotropy
      ! pthe----input-R-SSO Angle
      ! ppic----input-R-SSO Peacks elevation (m)
      ! pval----input-R-SSO Valleys elevation (m)

      ! ==== outputs ===
      ! pulow, pvlow -output-R: Low-level wind
      ! kkcrit----I-: Security value for top of low level flow
      ! kcrit-----I-: Critical level
      ! ksect-----I-: Not used
      ! kkhlim----I-: Not used
      ! kkenvh----I-: Top of blocked flow layer
      ! kknu------I-: Layer that sees mountain peacks
      ! kknu2-----I-: Layer that sees mountain peacks above mountain mean
      ! kknub-----I-: Layer that sees mountain mean above valleys
      ! prho------R-: Density at 1/2 layers
      ! pri-------R-: Background Richardson Number, Wind shear measured along GW
      ! stress
      ! pstab-----R-: Brunt-Vaisala freq. at 1/2 layers
      ! pvph------R-: Wind in  plan of GW stress, Half levels.
      ! ppsi------R-: Angle between low level wind and SS0 main axis.
      ! pd1-------R-| Compared the ratio of the stress
      ! pd2-------R-| that is along the wind to that Normal to it.
      ! pdi define the plane of low level stress
      ! compared to the low level wind.
      ! see p. 108 Lott & Miller (1997).
      ! pdmod-----R-: Norme of pdi

      ! === local arrays ===

      ! zvpf------R-: Wind projected in the plan of the low-level stress.

      ! ==== outputs ===

      ! implicit arguments :   none
      ! --------------------

      ! method.
      ! -------

      ! externals.
      ! ----------

      ! reference.
      ! ----------

      ! see ecmwf research department documentation of the "i.f.s."

      ! author.
      ! -------

      ! modifications.
      ! --------------
      ! f.lott  for the new-gwdrag scheme november 1993

      ! -----------------------------------------------------------------------
      USE lmdz_gwd_ini, ONLY: rd, rg, rcpd, rpi, gssec, gtsec, gvsec
      USE lmdz_gwd_ini, ONLY: nktopg, gsigcr, grcrit, gfrcrit

      IMPLICIT NONE

      ! -----------------------------------------------------------------------

      ! *       0.1   arguments
      ! ---------

      INTEGER, INTENT(IN) :: klon, klev
      INTEGER, DIMENSION(klon), INTENT(IN) :: ktest
      REAL, DIMENSION(klon), INTENT(IN) ::  pstd, pmea, ppic, pval, pthe
      REAL, DIMENSION(klon, klev), INTENT(IN) :: papm1, pum1, pvm1, ptm1, pgeom1
      REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paphm1
      REAL, DIMENSION(klon), INTENT(INOUT) :: pgam

      INTEGER, DIMENSION(klon), INTENT(OUT) :: kkcrit, kkcrith, kcrit, ksect, &
                                               kkhlim, kkenvh
      REAL, DIMENSION(klon, klev + 1), INTENT(OUT) :: prho, pri, pstab, ptau, pvph, ppsi
      REAL, DIMENSION(klon, klev), INTENT(OUT) ::  pzdep
      REAL, DIMENSION(klon), INTENT(OUT) :: pulow, pvlow, pnu, pd1, pd2, pdmod

      ! -----------------------------------------------------------------------

      ! *       0.2   local arrays
      ! ------------

      INTEGER ilevh, jl, jk
      REAL zcons1, zcons2, zhgeo, zu, zphi
      REAL zvt1, zvt2, zdwind, zwind, zdelp
      REAL zstabm, zstabp, zrhom, zrhop
      LOGICAL lo
      LOGICAL ll1(klon, klev + 1)
      INTEGER kknu(klon), kknu2(klon), kknub(klon), kknul(klon), kentp(klon), &
         ncount(klon)

      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)

      ! ------------------------------------------------------------------

      ! *         1.    initialization
      ! --------------

      ! ------------------------------------------------------------------

      ! *         1.1   computational constants
      ! -----------------------

      ilevh = klev/3

      zcons1 = 1./rd
      zcons2 = rg**2/rcpd

      ! ------------------------------------------------------------------

      ! *         2.
      ! --------------

      ! ------------------------------------------------------------------

      ! *         2.1     define low level wind, project winds in plane of
      ! *                 low level wind, determine sector in which to take
      ! *                 the variance and set indicator for critical levels.

      DO jl = 1, klon
         kknu(jl) = klev
         kknu2(jl) = klev
         kknub(jl) = klev
         kknul(jl) = klev
         pgam(jl) = max(pgam(jl), gtsec)
         ll1(jl, klev + 1) = .FALSE.
      END DO

      ! Ajouter une initialisation (L. Li, le 23fev99):

      DO jk = klev, ilevh, -1
         DO jl = 1, klon
            ll1(jl, jk) = .FALSE.
         END DO
      END DO

      ! *      define top of low level flow
      ! ----------------------------
      DO jk = klev, ilevh, -1
         DO jl = 1, klon
            IF (ktest(jl) == 1) THEN
               lo = (paphm1(jl, jk)/paphm1(jl, klev + 1)) >= gsigcr
               IF (lo) THEN
                  kkcrit(jl) = jk
               END IF
               zhcrit(jl, jk) = ppic(jl) - pval(jl)
               zhgeo = pgeom1(jl, jk)/rg
               ll1(jl, jk) = (zhgeo > zhcrit(jl, jk))
               IF (ll1(jl, jk) .NEQV. ll1(jl, jk + 1)) THEN
                  kknu(jl) = jk
               END IF
               IF (.NOT. ll1(jl, ilevh)) kknu(jl) = ilevh
            END IF
         END DO
      END DO
      DO jk = klev, ilevh, -1
         DO jl = 1, klon
            IF (ktest(jl) == 1) THEN
               zhcrit(jl, jk) = ppic(jl) - pmea(jl)
               zhgeo = pgeom1(jl, jk)/rg
               ll1(jl, jk) = (zhgeo > zhcrit(jl, jk))
               IF (ll1(jl, jk) .NEQV. ll1(jl, jk + 1)) THEN
                  kknu2(jl) = jk
               END IF
               IF (.NOT. ll1(jl, ilevh)) kknu2(jl) = ilevh
            END IF
         END DO
      END DO
      DO jk = klev, ilevh, -1
         DO jl = 1, klon
            IF (ktest(jl) == 1) THEN
               zhcrit(jl, jk) = amin1(ppic(jl) - pmea(jl), pmea(jl) - pval(jl))
               zhgeo = pgeom1(jl, jk)/rg
               ll1(jl, jk) = (zhgeo > zhcrit(jl, jk))
               IF (ll1(jl, jk) .NEQV. ll1(jl, jk + 1)) THEN
                  kknub(jl) = jk
               END IF
               IF (.NOT. ll1(jl, ilevh)) kknub(jl) = ilevh
            END IF
         END DO
      END DO

      DO jl = 1, klon
         IF (ktest(jl) == 1) THEN
            kknu(jl) = min(kknu(jl), nktopg)
            kknu2(jl) = min(kknu2(jl), nktopg)
            kknub(jl) = min(kknub(jl), nktopg)
            kknul(jl) = klev
         END IF
      END DO

      ! c*     initialize various arrays

      DO jl = 1, klon
         prho(jl, klev + 1) = 0.0
         ! ym correction en attendant mieux
         prho(jl, 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
         pvph(jl, klev + 1) = 0.0
         ! ym correction en attendant mieux
         ! ym      pvph(jl,klev)    =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.
      END DO

      ! *     define flow density and stratification (rho and N2)
      ! at semi layers.
      ! -------------------------------------------------------

      DO jk = klev, 2, -1
         DO jl = 1, klon
            IF (ktest(jl) == 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)
            END IF
         END DO
      END DO

      ! ********************************************************************

      ! *     define Low level flow (between ground and peacks-valleys)
      ! ---------------------------------------------------------
      DO jk = klev, ilevh, -1
         DO jl = 1, klon
            IF (ktest(jl) == 1) THEN
               IF (jk >= kknu2(jl) .AND. jk <= 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) &
                                                        )
                  pstab(jl, klev + 1) = pstab(jl, klev + 1) + pstab(jl, jk)*(paphm1(jl, jk &
                                                                                    + 1) - paphm1(jl, jk))
                  prho(jl, klev + 1) = prho(jl, klev + 1) + prho(jl, jk)*(paphm1(jl, jk + 1) &
                                                                          - paphm1(jl, jk))
               END IF
            END IF
         END DO
      END DO
      DO jl = 1, klon
         IF (ktest(jl) == 1) THEN
            pulow(jl) = pulow(jl)/(paphm1(jl, kknul(jl) + 1) - paphm1(jl, kknu2(jl)))
            pvlow(jl) = pvlow(jl)/(paphm1(jl, kknul(jl) + 1) - paphm1(jl, kknu2(jl)))
            znorm(jl) = max(sqrt(pulow(jl)**2 + pvlow(jl)**2), gvsec)
            pvph(jl, klev + 1) = znorm(jl)
            pstab(jl, klev + 1) = pstab(jl, klev + 1)/(paphm1(jl, kknul(jl) + 1) - paphm1(jl &
                                                                                          , kknu2(jl)))
            prho(jl, klev + 1) = prho(jl, klev + 1)/(paphm1(jl, kknul(jl) + 1) - paphm1(jl, &
                                                                                        kknu2(jl)))
         END IF
      END DO

      ! *******  setup orography orientation relative to the low level
      ! wind and define parameters of the Anisotropic wave stress.

      DO jl = 1, klon
         IF (ktest(jl) == 1) THEN
            lo = (pulow(jl) < gvsec) .AND. (pulow(jl) >= -gvsec)
            IF (lo) THEN
               zu = pulow(jl) + 2.*gvsec
            ELSE
               zu = pulow(jl)
            END IF
            zphi = atan(pvlow(jl)/zu)
            ppsi(jl, klev + 1) = pthe(jl)*rpi/180.-zphi
            zb(jl) = 1.-0.18*pgam(jl) - 0.04*pgam(jl)**2
            zc(jl) = 0.48*pgam(jl) + 0.3*pgam(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)
         END IF
      END DO

      ! ************ projet flow in plane of lowlevel stress *************
      ! ************ Find critical levels...                 *************

      DO jk = 1, klev
         DO jl = 1, klon
            IF (ktest(jl) == 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))
            END IF
            ptau(jl, jk) = 0.0
            pzdep(jl, jk) = 0.0
            ppsi(jl, jk) = 0.0
            ll1(jl, jk) = .FALSE.
         END DO
      END DO
      DO jk = 2, klev
         DO jl = 1, klon
            IF (ktest(jl) == 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) < gvsec) THEN
                  pvph(jl, jk) = gvsec
                  kcrit(jl) = jk
               END IF
            END IF
         END DO
      END DO

      ! *         2.3     mean flow richardson number.

      DO jk = 2, klev
         DO jl = 1, klon
            IF (ktest(jl) == 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)
            END IF
         END DO
      END DO

      ! *      define top of 'envelope' layer
      ! ----------------------------

      DO jl = 1, klon
         pnu(jl) = 0.0
         znum(jl) = 0.0
      END DO

      DO jk = 2, klev - 1
         DO jl = 1, klon

            IF (ktest(jl) == 1) THEN

               IF (jk >= kknu2(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) <= gfrcrit) .AND. (pnu(jl) > gfrcrit) .AND. (kkenvh( &
                                                                             jl) == klev)) kkenvh(jl) = jk

               END IF

            END IF

         END DO
      END DO

      ! calculation of a dynamical mixing height for when the waves
      ! BREAK AT LOW LEVEL: The drag will be repartited over
      ! a depths that depends on waves vertical wavelength,
      ! not just between two adjacent model layers.
      ! of gravity waves:

      DO jl = 1, klon
         znup(jl) = 0.0
         znum(jl) = 0.0
      END DO

      DO jk = klev - 1, 2, -1
         DO jl = 1, klon

            IF (ktest(jl) == 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) <= rpi/4.) .AND. (znup(jl) > rpi/4.) .AND. (kkcrith( &
                                                                         jl) == klev)) kkcrith(jl) = jk

            END IF

         END DO
      END DO

      DO jl = 1, klon
         IF (ktest(jl) == 1) THEN
            kkcrith(jl) = max0(kkcrith(jl), ilevh*2)
            kkcrith(jl) = max0(kkcrith(jl), kknu(jl))
            IF (kcrit(jl) >= kkcrith(jl)) kcrit(jl) = 1
         END IF
      END DO

      ! directional info for flow blocking *************************

      DO jk = 1, klev
         DO jl = 1, klon
            IF (ktest(jl) == 1) THEN
               lo = (pum1(jl, jk) < gvsec) .AND. (pum1(jl, jk) >= -gvsec)
               IF (lo) THEN
                  zu = pum1(jl, jk) + 2.*gvsec
               ELSE
                  zu = pum1(jl, jk)
               END IF
               zphi = atan(pvm1(jl, jk)/zu)
               ppsi(jl, jk) = pthe(jl)*rpi/180.-zphi
            END IF
         END DO
      END DO

      ! forms the vertical 'leakiness' **************************

      DO jk = ilevh, klev
         DO jl = 1, klon
            IF (ktest(jl) == 1) THEN
               pzdep(jl, jk) = 0
               IF (jk >= kkenvh(jl) .AND. kkenvh(jl) /= klev) THEN
                  pzdep(jl, jk) = (pgeom1(jl, kkenvh(jl)) - pgeom1(jl, jk))/ &
                                  (pgeom1(jl, kkenvh(jl)) - pgeom1(jl, klev))
               END IF
            END IF
         END DO
      END DO

      RETURN
   END SUBROUTINE orosetup_strato

! ============================================================================
   SUBROUTINE gwstress_strato(klon, klev, kkcrit, ksect, kkhlim, ktest, kkcrith, &
                              kcrit, kkenvh, kknu, prho, pstab, pvph, pstd, psig, pmea, ppic, pval, &
                              ptfr, ptau, pgeom1, pgamma, pd1, pd2, pdmod, pnu)

      ! **** *gwstress*

      ! purpose.
      ! --------
      ! Compute the surface stress due to Gravity Waves, according
      ! to the Phillips (1979) theory of 3-D flow above
      ! anisotropic elliptic ridges.

      ! The stress is reduced two account for cut-off flow over
      ! hill.  The flow only see that part of the ridge located
      ! above the blocked layer (see zeff).

      ! **   interface.
      ! ----------
      ! call *gwstress*  from *gwdrag*

      ! explicit arguments :
      ! --------------------
      ! ==== inputs ===
      ! ==== outputs ===

      ! implicit arguments :   none
      ! --------------------

      ! method.
      ! -------

      ! externals.
      ! ----------

      ! reference.
      ! ----------

      ! LOTT and MILLER (1997)  &  LOTT (1999)

      ! author.
      ! -------

      ! modifications.
      ! --------------
      ! f. lott put the new gwd on ifs      22/11/93

      ! -----------------------------------------------------------------------
      USE lmdz_gwd_ini, ONLY: gkdrag, gfrcrit

      IMPLICIT NONE

      ! -----------------------------------------------------------------------

      ! *       0.1   arguments
      ! ---------

      INTEGER, INTENT(IN) :: klon, klev
      INTEGER, DIMENSION(klon) :: kkcrit, kkcrith, kcrit, ksect, &
                                  kkhlim, ktest, kkenvh, kknu

      REAL prho(klon, klev + 1), pstab(klon, klev + 1), ptau(klon, klev + 1), &
         pvph(klon, klev + 1), ptfr(klon), pgeom1(klon, klev), pstd(klon)

      REAL pd1(klon), pd2(klon), pnu(klon), psig(klon), pgamma(klon)
      REAL pmea(klon), ppic(klon), pval(klon)
      REAL pdmod(klon)

      ! -----------------------------------------------------------------------

      ! *       0.2   local arrays
      ! ------------
      ! zeff--real: effective height seen by the flow when there is blocking

      INTEGER jl
      REAL zeff

      ! -----------------------------------------------------------------------

      ! *       0.3   functions
      ! ---------
      ! ------------------------------------------------------------------

      ! *         1.    initialization
      ! --------------

      ! *         3.1     gravity wave stress.

      DO jl = 1, klon
         IF (ktest(jl) == 1) THEN

            ! effective mountain height above the blocked flow

            zeff = ppic(jl) - pval(jl)
            IF (kkenvh(jl) < klev) THEN
               zeff = amin1(gfrcrit*pvph(jl, klev + 1)/sqrt(pstab(jl, klev + 1)), zeff)
            END IF

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

            ! too small value of stress or  low level flow include critical level
            ! 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)
            ! if(lo) ptau(jl,klev+1)=0.0

         ELSE

            ptau(jl, klev + 1) = 0.0

         END IF

      END DO

      RETURN
   END SUBROUTINE gwstress_strato

! ============================================================================
   SUBROUTINE gwprofil_strato(klon, klev, kgwd, kdx, ktest, kkcrit, kkcrith, &
                              kcrit, kkenvh, kknu, kknu2, paphm1, prho, pstab, ptfr, pvph, pri, ptau, &
                              pdmod, pnu, psig, pgamma, pstd, ppic, pval)

      ! **** *gwprofil*

      ! purpose.
      ! --------

      ! **   interface.
      ! ----------
      ! from *gwdrag*

      ! explicit arguments :
      ! --------------------
      ! ==== inputs ===

      ! ==== outputs ===

      ! implicit arguments :   none
      ! --------------------

      ! method:
      ! -------
      ! the stress profile for gravity waves is computed as follows:
      ! it decreases linearly with heights from the ground
      ! to the low-level indicated by kkcrith,
      ! to simulates lee waves or
      ! low-level gravity wave breaking.
      ! above it is constant, except when the waves encounter a critical
      ! level (kcrit) or when they break.
      ! The stress is also uniformly distributed above the level
      ! nstra.

      USE lmdz_gwd_ini, ONLY: nstra, grahilo, gssec, gvsec, gtsec, grcrit

      IMPLICIT NONE

      ! -----------------------------------------------------------------------

      ! *       0.1   ARGUMENTS
      ! ---------

      INTEGER klon, klev, kgwd
      INTEGER kkcrit(klon), kkcrith(klon), kcrit(klon), kdx(klon), ktest(klon), &
         kkenvh(klon), kknu(klon), kknu2(klon)

      REAL paphm1(klon, klev + 1), pstab(klon, klev + 1), prho(klon, klev + 1), &
         pvph(klon, klev + 1), pri(klon, klev + 1), ptfr(klon), ptau(klon, klev + 1)

      REAL pdmod(klon), pnu(klon), psig(klon), pgamma(klon), pstd(klon), &
         ppic(klon), pval(klon)

      ! -----------------------------------------------------------------------

      ! *       0.2   local arrays
      ! ------------

      INTEGER jl, jk
      REAL zsqr, zalfa, zriw, zdel, zb, zalpha, zdz2n, zdelp, zdelpt

      REAL zdz2(klon, klev), znorm(klon), zoro(klon)
      REAL ztau(klon, klev + 1)

      ! -----------------------------------------------------------------------

      ! *         1.    INITIALIZATION
      ! --------------

      ! *    COMPUTATIONAL CONSTANTS.
      ! ------------- ----------

      DO jl = 1, klon
         IF (ktest(jl) == 1) THEN
            zoro(jl) = psig(jl)*pdmod(jl)/4./pstd(jl)
            ztau(jl, klev + 1) = ptau(jl, klev + 1)
            ztau(jl, kkcrith(jl)) = grahilo*ptau(jl, klev + 1)
         END IF
      END DO

      DO jk = klev + 1, 1, -1
         ! *         4.1    constant shear stress until top of the
         ! low-level breaking/trapped layer

         DO jl = 1, klon
            IF (ktest(jl) == 1) THEN
               IF (jk > 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) + zdelp/zdelpt*(ztau(jl, kkcrith(jl) &
                                                                         ) - ztau(jl, klev + 1))
               ELSE
                  ptau(jl, jk) = ztau(jl, kkcrith(jl))
               END IF
            END IF
         END DO

         ! *         4.15   constant shear stress until the top of the
         ! low level flow layer.

         ! *         4.2    wave displacement at next level.

      END DO

      ! *         4.4    wave richardson number, new wave displacement
      ! *                and stress:  breaking evaluation and critical
      ! level

      DO jk = klev, 1, -1

         DO jl = 1, klon
            IF (ktest(jl) == 1) THEN
               znorm(jl) = prho(jl, jk)*sqrt(pstab(jl, jk))*pvph(jl, jk)
               zdz2(jl, jk) = ptau(jl, jk)/amax1(znorm(jl), gssec)/zoro(jl)
            END IF
         END DO

         DO jl = 1, klon
            IF (ktest(jl) == 1) THEN
               IF (jk < kkcrith(jl)) THEN
                  IF ((ptau(jl, jk + 1) < gtsec) .OR. (jk <= 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 < 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*zoro(jl)
                     END IF

                     ptau(jl, jk) = amin1(ptau(jl, jk), ptau(jl, jk + 1))

                  END IF
               END IF
            END IF
         END DO
      END DO

      ! REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL

      DO jl = 1, klon
         IF (ktest(jl) == 1) THEN
            ztau(jl, kkcrith(jl) - 1) = ptau(jl, kkcrith(jl) - 1)
            ztau(jl, nstra) = ptau(jl, nstra)
         END IF
      END DO

      DO jk = 1, klev

         DO jl = 1, klon
            IF (ktest(jl) == 1) THEN

               IF (jk > kkcrith(jl) - 1) THEN

                  zdelp = paphm1(jl, jk) - paphm1(jl, klev + 1)
                  zdelpt = paphm1(jl, kkcrith(jl) - 1) - paphm1(jl, klev + 1)
                  ptau(jl, jk) = ztau(jl, klev + 1) + (ztau(jl, kkcrith(jl) - 1) - ztau(jl, &
                                                                                        klev + 1))*zdelp/zdelpt

               END IF
            END IF

         END DO

         ! REORGANISATION AT THE MODEL TOP....

         DO jl = 1, klon
            IF (ktest(jl) == 1) THEN

               IF (jk < nstra) THEN

                  zdelp = paphm1(jl, nstra)
                  zdelpt = paphm1(jl, jk)
                  ptau(jl, jk) = ztau(jl, nstra)*zdelpt/zdelp

               END IF

            END IF

         END DO

      END DO

123   FORMAT(I4, 1X, 20(F6.3, 1X))

      RETURN
   END SUBROUTINE gwprofil_strato

!==========================================================================
   SUBROUTINE lift_noro_strato(klon, klev, dtime, paprs, pplay, plat, pmea, &
                               pstd, psig, pgam, pthe, ppic, pval, kgwd, kdx, ktest, t, u, v, pulow, &
                               pvlow, pustr, pvstr, d_t, d_u, d_v)

      USE lmdz_gwd_ini, ONLY: rd, rg

      IMPLICIT NONE
      ! ======================================================================
      ! Auteur(s): F.Lott (LMD/CNRS) date: 19950201
      ! Object: Mountain lift interface (enhanced vortex stretching).
      ! Made necessary because:
      ! 1. in the LMD-GCM Layers are from bottom to top,
      ! contrary to most European GCM.
      ! 2. the altitude above ground of each model layers
      ! needs to be known (variable zgeom)
      ! ======================================================================

      ! Declarations
      !====================================================

      ! Input variables
      !-----------------

      INTEGER, INTENT(IN) :: klon, klev ! horizontal and vertical dimensions
      REAL, INTENT(IN) ::  dtime ! timestep [s]
      REAL, DIMENSION(klon), INTENT(IN) :: plat !latitude [degrees]
      REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs ! pressure at bottom layer interface [Pa]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay ! pressure of model layers [Pa]
      REAL, DIMENSION(klon), INTENT(IN) :: pmea ! mean subgrid orography [m]
      REAL, DIMENSION(klon), INTENT(IN) :: pstd ! standard deviation of subgrid orography [m]
      REAL, DIMENSION(klon), INTENT(IN) :: psig ! sso mean slope [-]
      REAL, DIMENSION(klon), INTENT(IN) :: pgam ! sso anisotropy [-]
      REAL, DIMENSION(klon), INTENT(IN) :: pthe ! sso orientation angle [radian]
      REAL, DIMENSION(klon), INTENT(IN) :: ppic ! elevation of sso peaks [m]
      REAL, DIMENSION(klon), INTENT(IN) :: pval ! elevation of sso valleys [m]
      INTEGER, DIMENSION(klon), INTENT(IN) :: kdx ! locate the physical location of an active point
      INTEGER, DIMENSION(klon), INTENT(IN) :: ktest ! flags to indicate active points
      INTEGER, INTENT(IN) :: kgwd ! total number of grid points where the param is active
      REAL, DIMENSION(klon, klev), INTENT(IN) :: t ! temperature [K]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: u ! zonal wind component [m/s]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: v ! meridional wind component [m/s]

      ! Output variables
      !------------------

      REAL, DIMENSION(klon), INTENT(OUT) :: pulow ! low level u [m/s]
      REAL, DIMENSION(klon), INTENT(OUT) :: pvlow ! low level v [m/s]
      REAL, DIMENSION(klon), INTENT(OUT) :: pustr ! surface stress due to SSO drag (x component) [Pa]
      REAL, DIMENSION(klon), INTENT(OUT) :: pvstr ! surface stress due to SSO drag (y component) [Pa]

      REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t ! temperature increment [K]
      REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_u  ! u increment [m/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_v  ! v increment [m/s]

      ! Local variables
      !-----------------
      ! zgeom-----R: Altitude of layer above ground
      ! pt, pu, pv --R: t u v from top to bottom
      ! pdtdt, pdudt, pdvdt --R: t u v tendencies (from top to bottom)
      ! papmf: pressure at model layer (from top to bottom)
      ! papmh: pressure at model 1/2 layer (from top to bottom)
      INTEGER i, k
      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)

      ! initialise output variables (for security)

      DO i = 1, klon
         pulow(i) = 0.0
         pvlow(i) = 0.0
         pustr(i) = 0.0
         pvstr(i) = 0.0
      END DO
      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
         END DO
      END DO

      ! prepare intput variables (beware that vertical levels' order increase from top top bottom)

      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)
         END DO
      END DO
      DO k = 1, klev + 1
         DO i = 1, klon
            papmh(i, k) = paprs(i, klev - k + 2)
         END DO
      END DO
      DO i = 1, klon
         zgeom(i, klev) = rd*pt(i, klev)*log(papmh(i, klev + 1)/papmf(i, klev))
      END DO
      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))
         END DO
      END DO

      ! appeler la routine principale

      CALL orolift_strato(klon, klev, kgwd, kdx, ktest, dtime, papmh, papmf, &
                          zgeom, pt, pu, pv, plat, pmea, pstd, psig, pgam, pthe, ppic, pval, pulow, &
                          pvlow, pdudt, pdvdt, pdtdt)

      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) + pdudt(i, k)*(papmh(i, k + 1) - papmh(i, k))/rg
            pvstr(i) = pvstr(i) + pdvdt(i, k)*(papmh(i, k + 1) - papmh(i, k))/rg
         END DO
      END DO

      RETURN
   END SUBROUTINE lift_noro_strato

! ========================================================================
   SUBROUTINE orolift_strato(klon, klev, kgwd, kdx, ktest, ptsphy, paphm1, &
                             papm1, pgeom1, ptm1, pum1, pvm1, plat, pmea, pstd, psig, pgam, pthe, &
                             ppic, pval, pulow, pvlow, pvom, pvol, pte)

      ! **** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT.

      ! PURPOSE.
      ! --------
      ! this routine computes the physical tendencies of the
      ! prognostic variables u,v  when enhanced vortex stretching
      ! is needed.

      ! **   INTERFACE.
      ! ----------
      ! CALLED FROM *lift_noro

      ! Implicit Arguments:
      ! ===================

      ! klon-common-I: Number of points seen by the physics
      ! klev-common-I: Number of vertical layers

      ! ----------

      ! AUTHOR.
      ! -------
      ! F.LOTT  LMD 22/11/95

      USE lmdz_gwd_ini, ONLY: rpi, rd, rg, romega, gklift, gvsec, nktopg

      IMPLICIT NONE

      ! -----------------------------------------------------------------------

      ! *       0.1   ARGUMENTS
      ! ---------

      ! Input variables
      !-----------------

      INTEGER, INTENT(IN) :: klon, klev ! horizontal and vertical dimensions
      REAL, INTENT(IN) ::  ptsphy ! timestep [s]
      REAL, DIMENSION(klon), INTENT(IN) :: plat !latitude [degrees]
      INTEGER, DIMENSION(klon), INTENT(IN) :: kdx ! locate the physical location of an active point
      INTEGER, DIMENSION(klon), INTENT(IN) :: ktest ! flags to indicate active points
      REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paphm1 ! pressure at bottom layer interface [Pa]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: papm1 ! pressure of model layers [Pa]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: pgeom1 ! geopotential height above ground [m]
      REAL, DIMENSION(klon), INTENT(IN) :: pmea ! mean subgrid orography [m]
      REAL, DIMENSION(klon), INTENT(IN) :: pstd ! standard deviation of subgrid orography [m]
      REAL, DIMENSION(klon), INTENT(IN) :: psig ! sso mean slope [-]
      REAL, DIMENSION(klon), INTENT(IN) :: pgam ! sso anisotropy [-]
      REAL, DIMENSION(klon), INTENT(IN) :: pthe ! sso orientation angle [radian]
      REAL, DIMENSION(klon), INTENT(IN) :: ppic ! elevation of sso peaks [m]
      REAL, DIMENSION(klon), INTENT(IN) :: pval ! elevation of sso valleys [m]
      INTEGER, INTENT(IN) :: kgwd ! total number of grid points where the param is active
      REAL, DIMENSION(klon, klev), INTENT(IN) :: ptm1 ! temperature [K]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: pum1 ! zonal wind component [m/s]
      REAL, DIMENSION(klon, klev), INTENT(IN) :: pvm1 ! meridional wind component [m/s]

      ! Output variables
      !------------------

      REAL, DIMENSION(klon), INTENT(OUT) :: pulow ! low level u [m/s]
      REAL, DIMENSION(klon), INTENT(OUT) :: pvlow ! low level v [m/s]

      REAL, DIMENSION(klon, klev), INTENT(OUT) :: pte ! temperature tendency [K/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT) :: pvom  ! u tendency [m/s/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT) :: pvol  ! v tendency [m/s/s]

      ! -----------------------------------------------------------------------

      ! *       0.2   local arrays

      INTEGER jl, ilevh, jk
      REAL zhgeo, zdelp, zslow, zsqua, zscav, zbet
      ! ------------
      INTEGER iknub(klon), iknul(klon)
      LOGICAL ll1(klon, klev + 1)

      REAL ztau(klon, klev + 1), ztav(klon, klev + 1), zrho(klon, klev + 1)
      REAL zdudt(klon), zdvdt(klon)
      REAL zhcrit(klon, klev)

      LOGICAL lifthigh
      REAL zcons1, ztmst
      CHARACTER(LEN=20) :: modname = 'orolift_strato'
      CHARACTER(LEN=80) :: abort_message

      ! -----------------------------------------------------------------------

      ! *         1.1  initialisations
      ! ---------------

      lifthigh = .FALSE.

      IF (klon /= klon .OR. klev /= klev) THEN
         abort_message = 'pb dimension'
         CALL abort_physic(modname, abort_message, 1)
      END IF
      zcons1 = 1./rd
      ztmst = ptsphy

      DO jl = 1, klon
         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 jk = 1, klev
            pvom(jl, jk) = 0.0
            pvol(jl, jk) = 0.0
            pte(jl, jk) = 0.0
         END DO
      END DO

      ! *         2.1     DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
      ! *                 LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
      ! *                 THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.

      DO jk = klev, 1, -1
         DO jl = 1, klon
            IF (ktest(jl) == 1) THEN
               zhcrit(jl, jk) = amax1(ppic(jl) - pval(jl), 100.)
               zhgeo = pgeom1(jl, jk)/rg
               ll1(jl, jk) = (zhgeo > zhcrit(jl, jk))
               IF (ll1(jl, jk) .NEQV. ll1(jl, jk + 1)) THEN
                  iknub(jl) = jk
               END IF
            END IF
         END DO
      END DO

      DO jl = 1, klon
         IF (ktest(jl) == 1) THEN
            iknub(jl) = max(iknub(jl), klev/2)
            iknul(jl) = max(iknul(jl), 2*klev/3)
            IF (iknub(jl) > nktopg) iknub(jl) = nktopg
            IF (iknub(jl) == nktopg) iknul(jl) = klev
            IF (iknub(jl) == iknul(jl)) iknub(jl) = iknul(jl) - 1
         END IF
      END DO

      DO jk = klev, 2, -1
         DO jl = 1, klon
            zrho(jl, jk) = 2.*paphm1(jl, jk)*zcons1/(ptm1(jl, jk) + ptm1(jl, jk - 1))
         END DO
      END DO

      ! ********************************************************************

      ! *     define low level flow
      ! -------------------
      DO jk = klev, 1, -1
         DO jl = 1, klon
            IF (ktest(jl) == 1) THEN
               IF (jk >= iknub(jl) .AND. jk <= 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
            END IF
         END DO
      END DO
      DO jl = 1, klon
         IF (ktest(jl) == 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)))
         END IF
      END DO

      ! ***********************************************************

      ! *         3.      COMPUTE MOUNTAIN LIFT

      DO jl = 1, klon
         IF (ktest(jl) == 1) THEN
            ztau(jl, klev + 1) = -gklift*zrho(jl, klev + 1)*2.*romega* & ! *
                                 ! (2*pstd(jl)+pmea(jl))*
                                 2*pstd(jl)*sin(rpi/180.*plat(jl))*pvlow(jl)
            ztav(jl, klev + 1) = gklift*zrho(jl, klev + 1)*2.*romega* & ! *
                                 ! (2*pstd(jl)+pmea(jl))*
                                 2*pstd(jl)*sin(rpi/180.*plat(jl))*pulow(jl)
         ELSE
            ztau(jl, klev + 1) = 0.0
            ztav(jl, klev + 1) = 0.0
         END IF
      END DO

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

      DO jk = 1, klev
         DO jl = 1, klon
            IF (ktest(jl) == 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
            END IF
         END DO
      END DO

      ! *         5.      COMPUTE TENDENCIES.
      ! *                 -------------------
      IF (lifthigh) THEN
         ! EXPLICIT SOLUTION AT ALL LEVELS

         DO jk = 1, klev
            DO jl = 1, klon
               IF (ktest(jl) == 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
               END IF
            END DO
         END DO

         ! PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY

         DO jk = 1, klev
            DO jl = 1, klon
               IF (ktest(jl) == 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 > 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
                  END IF
                  zsqua = sqrt(pum1(jl, jk)**2 + pum1(jl, jk)**2)
                  IF (zsqua < zslow) THEN
                     pvom(jl, jk) = zsqua/zslow*pvom(jl, jk)
                     pvol(jl, jk) = zsqua/zslow*pvol(jl, jk)
                  END IF

               END IF
            END DO
         END DO

         ! 6.  LOW LEVEL LIFT, SEMI IMPLICIT:
         ! ----------------------------------

      ELSE

         DO jl = 1, klon
            IF (ktest(jl) == 1) THEN
               DO jk = klev, iknub(jl), -1
                  zbet = gklift*2.*romega*sin(rpi/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)
               END DO
            END IF
         END DO

      END IF

      RETURN
   END SUBROUTINE orolift_strato

END MODULE lmdz_gwd_ogwd
