source: LMDZ5/trunk/libf/phylmd/wake.F90 @ 2390

! $Id: wake.F90 2346 2015-08-21 15:13:46Z musat $

SUBROUTINE wake(p, ph, pi, dtime, sigd_con, &
                te0, qe0, omgb, &
                dtdwn, dqdwn, amdwn, amup, dta, dqa, &
                wdtpbl, wdqpbl, udtpbl, udqpbl, &
                deltatw, deltaqw, dth, hw, sigmaw, wape, fip, gfl, &
                dtls, dqls, ktopw, omgbdth, dp_omgb, wdens, tu, qu, &
                dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, spread, cstar, &
                d_deltat_gw, d_deltatw2, d_deltaqw2)


  ! **************************************************************
  ! *
  ! WAKE                                                        *
  ! retour a un Pupper fixe                                *
  ! *
  ! written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
  ! modified by :   ROEHRIG Romain        01/29/2007            *
  ! **************************************************************

  USE dimphy
  use mod_phys_lmdz_para
  USE print_control_mod, ONLY: prt_level
  IMPLICIT NONE
  ! ============================================================================


  ! But : Decrire le comportement des poches froides apparaissant dans les
  ! grands systemes convectifs, et fournir l'energie disponible pour
  ! le declenchement de nouvelles colonnes convectives.

  ! Variables d'etat : deltatw    : ecart de temperature wake-undisturbed
  ! area
  ! deltaqw    : ecart d'humidite wake-undisturbed area
  ! sigmaw     : fraction d'aire occupee par la poche.

  ! Variable de sortie :

  ! wape : WAke Potential Energy
  ! fip  : Front Incident Power (W/m2) - ALP
  ! gfl  : Gust Front Length per unit area (m-1)
  ! dtls : large scale temperature tendency due to wake
  ! dqls : large scale humidity tendency due to wake
  ! hw   : hauteur de la poche
  ! dp_omgb : vertical gradient of large scale omega
  ! wdens   : densite de poches
  ! omgbdth: flux of Delta_Theta transported by LS omega
  ! dtKE   : differential heating (wake - unpertubed)
  ! dqKE   : differential moistening (wake - unpertubed)
  ! omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
  ! dp_deltomg  : vertical gradient of omg (s-1)
  ! spread  : spreading term in dt_wake and dq_wake
  ! deltatw     : updated temperature difference (T_w-T_u).
  ! deltaqw     : updated humidity difference (q_w-q_u).
  ! sigmaw      : updated wake fractional area.
  ! d_deltat_gw : delta T tendency due to GW

  ! Variables d'entree :

  ! aire : aire de la maille
  ! te0  : temperature dans l'environnement  (K)
  ! qe0  : humidite dans l'environnement     (kg/kg)
  ! omgb : vitesse verticale moyenne sur la maille (Pa/s)
  ! dtdwn: source de chaleur due aux descentes (K/s)
  ! dqdwn: source d'humidite due aux descentes (kg/kg/s)
  ! dta  : source de chaleur due courants satures et detrain  (K/s)
  ! dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
  ! amdwn: flux de masse total des descentes, par unite de
  ! surface de la maille (kg/m2/s)
  ! amup : flux de masse total des ascendances, par unite de
  ! surface de la maille (kg/m2/s)
  ! p    : pressions aux milieux des couches (Pa)
  ! ph   : pressions aux interfaces (Pa)
  ! pi  : (p/p_0)**kapa (adim)
  ! dtime: increment temporel (s)

  ! Variables internes :

  ! rhow : masse volumique de la poche froide
  ! rho  : environment density at P levels
  ! rhoh : environment density at Ph levels
  ! te   : environment temperature | may change within
  ! qe   : environment humidity    | sub-time-stepping
  ! the  : environment potential temperature
  ! thu  : potential temperature in undisturbed area
  ! tu   :  temperature  in undisturbed area
  ! qu   : humidity in undisturbed area
  ! dp_omgb: vertical gradient og LS omega
  ! omgbw  : wake average vertical omega
  ! dp_omgbw: vertical gradient of omgbw
  ! omgbdq : flux of Delta_q transported by LS omega
  ! dth  : potential temperature diff. wake-undist.
  ! th1  : first pot. temp. for vertical advection (=thu)
  ! th2  : second pot. temp. for vertical advection (=thw)
  ! q1   : first humidity for vertical advection
  ! q2   : second humidity for vertical advection
  ! d_deltatw   : terme de redistribution pour deltatw
  ! d_deltaqw   : terme de redistribution pour deltaqw
  ! deltatw0   : deltatw initial
  ! deltaqw0   : deltaqw initial
  ! hw0    : hw initial
  ! sigmaw0: sigmaw initial
  ! amflux : horizontal mass flux through wake boundary
  ! wdens_ref: initial number of wakes per unit area (3D) or per
  ! unit length (2D), at the beginning of each time step
  ! Tgw    : 1 sur la période de onde de gravité
  ! Cgw    : vitesse de propagation de onde de gravité
  ! LL     : distance entre 2 poches

  ! -------------------------------------------------------------------------
  ! Déclaration de variables
  ! -------------------------------------------------------------------------

  include "YOMCST.h"
  include "cvthermo.h"

  ! Arguments en entree
  ! --------------------

  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: p, pi
  REAL, DIMENSION (klon, klev+1),   INTENT(IN)          :: ph, omgb
  REAL,                             INTENT(IN)          :: dtime
  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: te0, qe0
  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dtdwn, dqdwn
  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: wdtpbl, wdqpbl, udtpbl, udqpbl ! UNUSED
  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: amdwn, amup
  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dta, dqa
  REAL, DIMENSION (klon),           INTENT(IN)          :: sigd_con

  !
  ! Input/Output
  REAL, DIMENSION (klon, klev),     INTENT(INOUT)       :: deltatw, deltaqw
  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw

  ! Sorties
  ! --------

  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dth
  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: tu, qu
  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtls, dqls
  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke
  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtpbl, dqpbl
  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: spread
  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
  REAL, DIMENSION (klon, klev+1),   INTENT(OUT)         :: omgbdth, omg
  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dp_omgb, dp_deltomg
  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_gw
  REAL, DIMENSION (klon),           INTENT(OUT)         :: hw, wape, fip, gfl, cstar
  REAL, DIMENSION (klon),           INTENT(OUT)         :: wdens
  INTEGER, DIMENSION (klon),        INTENT(OUT)         :: ktopw

  ! Variables internes
  ! -------------------

  ! Variables à fixer
  REAL alon
  LOGICAL, SAVE :: first = .TRUE.
  !$OMP THREADPRIVATE(first)
  REAL, SAVE ::  stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol  
  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
  REAL delta_t_min
  INTEGER nsub
  REAL dtimesub
  REAL sigmad, hwmin, wapecut
  REAL :: sigmaw_max
  REAL :: dens_rate
  REAL wdens0
  ! IM 080208
  LOGICAL, DIMENSION (klon) :: gwake

  ! Variables de sauvegarde
  REAL, DIMENSION (klon, klev) :: deltatw0
  REAL, DIMENSION (klon, klev) :: deltaqw0
  REAL, DIMENSION (klon, klev) :: te, qe
  REAL, DIMENSION (klon) :: sigmaw0, sigmaw1

  ! Variables pour les GW
  REAL, DIMENSION (klon) :: ll
  REAL, DIMENSION (klon, klev) :: n2
  REAL, DIMENSION (klon, klev) :: cgw
  REAL, DIMENSION (klon, klev) :: tgw

  ! Variables liées au calcul de hw
  REAL, DIMENSION (klon) :: ptop_provis, ptop, ptop_new
  REAL, DIMENSION (klon) :: sum_dth
  REAL, DIMENSION (klon) :: dthmin
  REAL, DIMENSION (klon) :: z, dz, hw0
  INTEGER, DIMENSION (klon) :: ktop, kupper

  ! Sub-timestep tendencies and related variables
  REAL d_deltatw(klon, klev), d_deltaqw(klon, klev)
  REAL d_te(klon, klev), d_qe(klon, klev)
  REAL d_sigmaw(klon), alpha(klon)
  REAL q0_min(klon), q1_min(klon)
  LOGICAL wk_adv(klon), ok_qx_qw(klon)
  REAL epsilon
  DATA epsilon/1.E-15/

  ! Autres variables internes
  INTEGER isubstep, k, i

  REAL, DIMENSION (klon) :: sum_thu, sum_tu, sum_qu, sum_thvu
  REAL, DIMENSION (klon) :: sum_dq, sum_rho
  REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn
  REAL, DIMENSION (klon) :: av_thu, av_tu, av_qu, av_thvu
  REAL, DIMENSION (klon) :: av_dth, av_dq, av_rho
  REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn

  REAL, DIMENSION (klon, klev) :: rho, rhow
  REAL, DIMENSION (klon, klev+1) :: rhoh
  REAL, DIMENSION (klon, klev) :: rhow_moyen
  REAL, DIMENSION (klon, klev) :: zh
  REAL, DIMENSION (klon, klev+1) :: zhh
  REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2

  REAL, DIMENSION (klon, klev) :: the, thu

  ! REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw

  REAL, DIMENSION (klon, klev+1) :: omgbw
  REAL, DIMENSION (klon) :: pupper
  REAL, DIMENSION (klon) :: omgtop
  REAL, DIMENSION (klon, klev) :: dp_omgbw
  REAL, DIMENSION (klon) :: ztop, dztop
  REAL, DIMENSION (klon, klev) :: alpha_up

  REAL, DIMENSION (klon) :: rre1, rre2
  REAL :: rrd1, rrd2
  REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2
  REAL, DIMENSION (klon, klev) :: d_th1, d_th2, d_dth
  REAL, DIMENSION (klon, klev) :: d_q1, d_q2, d_dq
  REAL, DIMENSION (klon, klev) :: omgbdq

  REAL, DIMENSION (klon) :: ff, gg
  REAL, DIMENSION (klon) :: wape2, cstar2, heff

  REAL, DIMENSION (klon, klev) :: crep

  REAL, DIMENSION (klon, klev) :: ppi

  ! cc nrlmd
  REAL, DIMENSION (klon) :: death_rate, nat_rate
  REAL, DIMENSION (klon, klev) :: entr
  REAL, DIMENSION (klon, klev) :: detr

  ! -------------------------------------------------------------------------
  ! Initialisations
  ! -------------------------------------------------------------------------

  ! print*, 'wake initialisations'

  ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
  ! -------------------------------------------------------------------------

  DATA wapecut, sigmad, hwmin/5., .02, 10./
  ! cc nrlmd
  DATA sigmaw_max/0.4/
  DATA dens_rate/0.1/
  ! cc
  ! Longueur de maille (en m)
  ! -------------------------------------------------------------------------

  ! ALON = 3.e5
  alon = 1.E6


  ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)

  ! coefgw : Coefficient pour les ondes de gravité
  ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
  ! wdens : Densité de poche froide par maille
  ! -------------------------------------------------------------------------

  ! cc nrlmd      coefgw=10
  ! coefgw=1
  ! wdens0 = 1.0/(alon**2)
  ! cc nrlmd      wdens = 1.0/(alon**2)
  ! cc nrlmd      stark = 0.50
  ! CRtest
  ! cc nrlmd      alpk=0.1
  ! alpk = 1.0
  ! alpk = 0.5
  ! alpk = 0.05

 if (first) then
  stark = 0.33
  alpk = 0.25
  wdens_ref = 8.E-12
  coefgw = 4.
  crep_upper = 0.9
  crep_sol = 1.0

  ! cc nrlmd Lecture du fichier wake_param.data
 !$OMP MASTER
  OPEN (99, FILE='wake_param.data', STATUS='old', FORM='formatted', ERR=9999)
  READ (99, *, END=9998) stark
  READ (99, *, END=9998) alpk
  READ (99, *, END=9998) wdens_ref
  READ (99, *, END=9998) coefgw
9998 CONTINUE
  CLOSE (99)
9999 CONTINUE
 !$OMP END MASTER
  CALL bcast(stark)
  CALL bcast(alpk)
  CALL bcast(wdens_ref)
  CALL bcast(coefgw)

  first=.false.
 endif

  ! Initialisation de toutes des densites a wdens_ref.
  ! Les densites peuvent evoluer si les poches debordent
  ! (voir au tout debut de la boucle sur les substeps)
  wdens = wdens_ref

  ! print*,'stark',stark
  ! print*,'alpk',alpk
  ! print*,'wdens',wdens
  ! print*,'coefgw',coefgw
  ! cc
  ! Minimum value for |T_wake - T_undist|. Used for wake top definition
  ! -------------------------------------------------------------------------

  delta_t_min = 0.2

  ! 1. - Save initial values and initialize tendencies
  ! --------------------------------------------------

  DO k = 1, klev
    DO i = 1, klon
      ppi(i, k) = pi(i, k)
      deltatw0(i, k) = deltatw(i, k)
      deltaqw0(i, k) = deltaqw(i, k)
      te(i, k) = te0(i, k)
      qe(i, k) = qe0(i, k)
      dtls(i, k) = 0.
      dqls(i, k) = 0.
      d_deltat_gw(i, k) = 0.
      d_te(i, k) = 0.
      d_qe(i, k) = 0.
      d_deltatw(i, k) = 0.
      d_deltaqw(i, k) = 0.
      ! IM 060508 beg
      d_deltatw2(i, k) = 0.
      d_deltaqw2(i, k) = 0.
      ! IM 060508 end
    END DO
  END DO
  ! sigmaw1=sigmaw
  ! IF (sigd_con.GT.sigmaw1) THEN
  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
  ! ENDIF
  DO i = 1, klon
    ! c      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
    sigmaw(i) = amax1(sigmaw(i), sigmad)
    sigmaw(i) = amin1(sigmaw(i), 0.99)
    sigmaw0(i) = sigmaw(i)
    wape(i) = 0.
    wape2(i) = 0.
    d_sigmaw(i) = 0.
    ktopw(i) = 0
  END DO


  ! 2. - Prognostic part
  ! --------------------


  ! 2.1 - Undisturbed area and Wake integrals
  ! ---------------------------------------------------------

  DO i = 1, klon
    z(i) = 0.
    ktop(i) = 0
    kupper(i) = 0
    sum_thu(i) = 0.
    sum_tu(i) = 0.
    sum_qu(i) = 0.
    sum_thvu(i) = 0.
    sum_dth(i) = 0.
    sum_dq(i) = 0.
    sum_rho(i) = 0.
    sum_dtdwn(i) = 0.
    sum_dqdwn(i) = 0.

    av_thu(i) = 0.
    av_tu(i) = 0.
    av_qu(i) = 0.
    av_thvu(i) = 0.
    av_dth(i) = 0.
    av_dq(i) = 0.
    av_rho(i) = 0.
    av_dtdwn(i) = 0.
    av_dqdwn(i) = 0.
  END DO

  ! Distance between wakes
  DO i = 1, klon
    ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
  END DO
  ! Potential temperatures and humidity
  ! ----------------------------------------------------------
  DO k = 1, klev
    DO i = 1, klon
      ! write(*,*)'wake 1',i,k,rd,te(i,k)
      rho(i, k) = p(i, k)/(rd*te(i,k))
      ! write(*,*)'wake 2',rho(i,k)
      IF (k==1) THEN
        ! write(*,*)'wake 3',i,k,rd,te(i,k)
        rhoh(i, k) = ph(i, k)/(rd*te(i,k))
        ! write(*,*)'wake 4',i,k,rd,te(i,k)
        zhh(i, k) = 0
      ELSE
        ! write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1))
        rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
        ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
        zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
      END IF
      ! write(*,*)'wake 7',ppi(i,k)
      the(i, k) = te(i, k)/ppi(i, k)
      thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
      tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
      qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
      ! write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k)))
      rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
      dth(i, k) = deltatw(i, k)/ppi(i, k)
    END DO
  END DO

  DO k = 1, klev - 1
    DO i = 1, klon
      IF (k==1) THEN
        n2(i, k) = 0
      ELSE
        n2(i, k) = amax1(0., -rg**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i, &
          k-1))/(p(i,k+1)-p(i,k-1)))
      END IF
      zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2

      cgw(i, k) = sqrt(n2(i,k))*zh(i, k)
      tgw(i, k) = coefgw*cgw(i, k)/ll(i)
    END DO
  END DO

  DO i = 1, klon
    n2(i, klev) = 0
    zh(i, klev) = 0
    cgw(i, klev) = 0
    tgw(i, klev) = 0
  END DO

  ! Calcul de la masse volumique moyenne de la colonne   (bdlmd)
  ! -----------------------------------------------------------------

  DO k = 1, klev
    DO i = 1, klon
      epaisseur1(i, k) = 0.
      epaisseur2(i, k) = 0.
    END DO
  END DO

  DO i = 1, klon
    epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
    epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
    rhow_moyen(i, 1) = rhow(i, 1)
  END DO

  DO k = 2, klev
    DO i = 1, klon
      epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1.
      epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k)
      rhow_moyen(i, k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+rhow(i,k)* &
        epaisseur1(i,k))/epaisseur2(i, k)
    END DO
  END DO


  ! Choose an integration bound well above wake top
  ! -----------------------------------------------------------------

  ! Pupper = 50000.  ! melting level
  ! Pupper = 60000.
  ! Pupper = 80000.  ! essais pour case_e
  DO i = 1, klon
    pupper(i) = 0.6*ph(i, 1)
    pupper(i) = max(pupper(i), 45000.)
    ! cc        Pupper(i) = 60000.
  END DO


  ! Determine Wake top pressure (Ptop) from buoyancy integral
  ! --------------------------------------------------------

  ! -1/ Pressure of the level where dth becomes less than delta_t_min.

  DO i = 1, klon
    ptop_provis(i) = ph(i, 1)
  END DO
  DO k = 2, klev
    DO i = 1, klon

      ! IM v3JYG; ptop_provis(i).LT. ph(i,1)

      IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. &
          ptop_provis(i)==ph(i,1)) THEN
        ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
          k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
      END IF
    END DO
  END DO

  ! -2/ dth integral

  DO i = 1, klon
    sum_dth(i) = 0.
    dthmin(i) = -delta_t_min
    z(i) = 0.
  END DO

  DO k = 1, klev
    DO i = 1, klon
      dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
      IF (dz(i)>0) THEN
        z(i) = z(i) + dz(i)
        sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
        dthmin(i) = amin1(dthmin(i), dth(i,k))
      END IF
    END DO
  END DO

  ! -3/ height of triangle with area= sum_dth and base = dthmin

  DO i = 1, klon
    hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
    hw0(i) = amax1(hwmin, hw0(i))
  END DO

  ! -4/ now, get Ptop

  DO i = 1, klon
    z(i) = 0.
    ptop(i) = ph(i, 1)
  END DO

  DO k = 1, klev
    DO i = 1, klon
      dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i))
      IF (dz(i)>0) THEN
        z(i) = z(i) + dz(i)
        ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i)
      END IF
    END DO
  END DO


  ! -5/ Determination de ktop et kupper

  DO k = klev, 1, -1
    DO i = 1, klon
      IF (ph(i,k+1)<ptop(i)) ktop(i) = k
      IF (ph(i,k+1)<pupper(i)) kupper(i) = k
    END DO
  END DO

  ! On evite kupper = 1 et kupper = klev
  DO i = 1, klon
    kupper(i) = max(kupper(i), 2)
    kupper(i) = min(kupper(i), klev-1)
  END DO


  ! -6/ Correct ktop and ptop

  DO i = 1, klon
    ptop_new(i) = ptop(i)
  END DO
  DO k = klev, 2, -1
    DO i = 1, klon
      IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
          dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
        ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
          k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
      END IF
    END DO
  END DO

  DO i = 1, klon
    ptop(i) = ptop_new(i)
  END DO

  DO k = klev, 1, -1
    DO i = 1, klon
      IF (ph(i,k+1)<ptop(i)) ktop(i) = k
    END DO
  END DO

  ! -5/ Set deltatw & deltaqw to 0 above kupper

  DO k = 1, klev
    DO i = 1, klon
      IF (k>=kupper(i)) THEN
        deltatw(i, k) = 0.
        deltaqw(i, k) = 0.
      END IF
    END DO
  END DO


  ! Vertical gradient of LS omega

  DO k = 1, klev
    DO i = 1, klon
      IF (k<=kupper(i)) THEN
        dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k))
      END IF
    END DO
  END DO

  ! Integrals (and wake top level number)
  ! --------------------------------------

  ! Initialize sum_thvu to 1st level virt. pot. temp.

  DO i = 1, klon
    z(i) = 1.
    dz(i) = 1.
    sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i)
    sum_dth(i) = 0.
  END DO

  DO k = 1, klev
    DO i = 1, klon
      dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
      IF (dz(i)>0) THEN
        z(i) = z(i) + dz(i)
        sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
        sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
        sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
        sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i)
        sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
        sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
        sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
        sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
        sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
      END IF
    END DO
  END DO

  DO i = 1, klon
    hw0(i) = z(i)
  END DO


  ! 2.1 - WAPE and mean forcing computation
  ! ---------------------------------------

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

  ! Means

  DO i = 1, klon
    av_thu(i) = sum_thu(i)/hw0(i)
    av_tu(i) = sum_tu(i)/hw0(i)
    av_qu(i) = sum_qu(i)/hw0(i)
    av_thvu(i) = sum_thvu(i)/hw0(i)
    ! av_thve = sum_thve/hw0
    av_dth(i) = sum_dth(i)/hw0(i)
    av_dq(i) = sum_dq(i)/hw0(i)
    av_rho(i) = sum_rho(i)/hw0(i)
    av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
    av_dqdwn(i) = sum_dqdwn(i)/hw0(i)

    wape(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i &
      )+av_dth(i)*av_dq(i)))/av_thvu(i)
  END DO

  ! 2.2 Prognostic variable update
  ! ------------------------------

  ! Filter out bad wakes

  DO k = 1, klev
    DO i = 1, klon
      IF (wape(i)<0.) THEN
        deltatw(i, k) = 0.
        deltaqw(i, k) = 0.
        dth(i, k) = 0.
      END IF
    END DO
  END DO

  DO i = 1, klon
    IF (wape(i)<0.) THEN
      wape(i) = 0.
      cstar(i) = 0.
      hw(i) = hwmin
      sigmaw(i) = amax1(sigmad, sigd_con(i))
      fip(i) = 0.
      gwake(i) = .FALSE.
    ELSE
      cstar(i) = stark*sqrt(2.*wape(i))
      gwake(i) = .TRUE.
    END IF
  END DO


  ! Check qx and qw positivity
  ! --------------------------
  DO i = 1, klon
    q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)), (qe(i, &
      1)+(1.-sigmaw(i))*deltaqw(i,1)))
  END DO
  DO k = 2, klev
    DO i = 1, klon
      q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), (qe(i, &
        k)+(1.-sigmaw(i))*deltaqw(i,k)))
      IF (q1_min(i)<=q0_min(i)) THEN
        q0_min(i) = q1_min(i)
      END IF
    END DO
  END DO

  DO i = 1, klon
    ok_qx_qw(i) = q0_min(i) >= 0.
    alpha(i) = 1.
  END DO

  ! C -----------------------------------------------------------------
  ! Sub-time-stepping
  ! -----------------

  nsub = 10
  dtimesub = dtime/nsub

  ! ------------------------------------------------------------
  DO isubstep = 1, nsub
    ! ------------------------------------------------------------

    ! wk_adv is the logical flag enabling wake evolution in the time advance
    ! loop
    DO i = 1, klon
      wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
    END DO

    ! cc nrlmd   Ajout d'un recalcul de wdens dans le cas d'un entrainement
    ! négatif de ktop à kupper --------
    ! cc           On calcule pour cela une densité wdens0 pour laquelle on
    ! aurait un entrainement nul ---
    DO i = 1, klon
      ! c       print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
      ! c     $           isubstep,wk_adv(i),cstar(i),wape(i)
      IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
        omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
          rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
        wdens0 = (sigmaw(i)/(4.*3.14))*((1.-sigmaw(i))*omg(i,kupper(i)+1)/(( &
          ph(i,1)-pupper(i))*cstar(i)))**(2)
        IF (wdens(i)<=wdens0*1.1) THEN
          wdens(i) = wdens0
        END IF
        ! c        print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
        ! c     $     ,ph(i,1)-pupper(i)',
        ! c     $             omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
        ! c     $     ,ph(i,1)-pupper(i)
      END IF
    END DO

    ! cc nrlmd

    DO i = 1, klon
      IF (wk_adv(i)) THEN
        gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
        sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
      END IF
    END DO
    DO i = 1, klon
      IF (wk_adv(i)) THEN
        ! cc nrlmd          Introduction du taux de mortalité des poches et
        ! test sur sigmaw_max=0.4
        ! cc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
        IF (sigmaw(i)>=sigmaw_max) THEN
          death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
        ELSE
          death_rate(i) = 0.
        END IF
        d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* &
          dtimesub
        ! $              - nat_rate(i)*sigmaw(i)*dtimesub
        ! c        print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
        ! c     $  death_rate(i),ktop(i),kupper(i)',
        ! c     $                d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
        ! c     $  death_rate(i),ktop(i),kupper(i)

        ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
        ! sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
        ! wdens = wdens0/(10.*sigmaw)
        ! sigmaw =max(sigmaw,sigd_con)
        ! sigmaw =max(sigmaw,sigmad)
      END IF
    END DO


    ! calcul de la difference de vitesse verticale poche - zone non perturbee
    ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
    ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on
    ! definit
    ! IM 060208 au niveau k=1..?
    DO k = 1, klev
      DO i = 1, klon
        IF (wk_adv(i)) THEN !!! nrlmd
          dp_deltomg(i, k) = 0.
        END IF
      END DO
    END DO
    DO k = 1, klev + 1
      DO i = 1, klon
        IF (wk_adv(i)) THEN !!! nrlmd
          omg(i, k) = 0.
        END IF
      END DO
    END DO

    DO i = 1, klon
      IF (wk_adv(i)) THEN
        z(i) = 0.
        omg(i, 1) = 0.
        dp_deltomg(i, 1) = -(gfl(i)*cstar(i))/(sigmaw(i)*(1-sigmaw(i)))
      END IF
    END DO

    DO k = 2, klev
      DO i = 1, klon
        IF (wk_adv(i) .AND. k<=ktop(i)) THEN
          dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
          z(i) = z(i) + dz(i)
          dp_deltomg(i, k) = dp_deltomg(i, 1)
          omg(i, k) = dp_deltomg(i, 1)*z(i)
        END IF
      END DO
    END DO

    DO i = 1, klon
      IF (wk_adv(i)) THEN
        dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
        ztop(i) = z(i) + dztop(i)
        omgtop(i) = dp_deltomg(i, 1)*ztop(i)
      END IF
    END DO

    ! -----------------
    ! From m/s to Pa/s
    ! -----------------

    DO i = 1, klon
      IF (wk_adv(i)) THEN
        omgtop(i) = -rho(i, ktop(i))*rg*omgtop(i)
        dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
      END IF
    END DO

    DO k = 1, klev
      DO i = 1, klon
        IF (wk_adv(i) .AND. k<=ktop(i)) THEN
          omg(i, k) = -rho(i, k)*rg*omg(i, k)
          dp_deltomg(i, k) = dp_deltomg(i, 1)
        END IF
      END DO
    END DO

    ! raccordement lineaire de omg de ptop a pupper

    DO i = 1, klon
      IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
        omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
          rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
        dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
          (ptop(i)-pupper(i))
      END IF
    END DO

    ! c      DO i=1,klon
    ! c        print*,'Pente entre 0 et kupper (référence)'
    ! c     $           ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
    ! c        print*,'Pente entre ktop et kupper'
    ! c     $   ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
    ! c      ENDDO
    ! c
    DO k = 1, klev
      DO i = 1, klon
        IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN
          dp_deltomg(i, k) = dp_deltomg(i, kupper(i))
          omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i))
        END IF
      END DO
    END DO
    ! cc nrlmd
    ! c      DO i=1,klon
    ! c      print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
    ! c      END DO
    ! cc


    ! --    Compute wake average vertical velocity omgbw


    DO k = 1, klev + 1
      DO i = 1, klon
        IF (wk_adv(i)) THEN
          omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k)
        END IF
      END DO
    END DO
    ! --    and its vertical gradient dp_omgbw

    DO k = 1, klev
      DO i = 1, klon
        IF (wk_adv(i)) THEN
          dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
        END IF
      END DO
    END DO

    ! --    Upstream coefficients for omgb velocity
    ! --    (alpha_up(k) is the coefficient of the value at level k)
    ! --    (1-alpha_up(k) is the coefficient of the value at level k-1)
    DO k = 1, klev
      DO i = 1, klon
        IF (wk_adv(i)) THEN
          alpha_up(i, k) = 0.
          IF (omgb(i,k)>0.) alpha_up(i, k) = 1.
        END IF
      END DO
    END DO

    ! Matrix expressing [The,deltatw] from  [Th1,Th2]

    DO i = 1, klon
      IF (wk_adv(i)) THEN
        rre1(i) = 1. - sigmaw(i)
        rre2(i) = sigmaw(i)
      END IF
    END DO
    rrd1 = -1.
    rrd2 = 1.

    ! --    Get [Th1,Th2], dth and [q1,q2]

    DO k = 1, klev
      DO i = 1, klon
        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
          dth(i, k) = deltatw(i, k)/ppi(i, k)
          th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
          th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
          q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
          q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
        END IF
      END DO
    END DO

    DO i = 1, klon
      IF (wk_adv(i)) THEN !!! nrlmd
        d_th1(i, 1) = 0.
        d_th2(i, 1) = 0.
        d_dth(i, 1) = 0.
        d_q1(i, 1) = 0.
        d_q2(i, 1) = 0.
        d_dq(i, 1) = 0.
      END IF
    END DO

    DO k = 2, klev
      DO i = 1, klon
        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
          d_th1(i, k) = th1(i, k-1) - th1(i, k)
          d_th2(i, k) = th2(i, k-1) - th2(i, k)
          d_dth(i, k) = dth(i, k-1) - dth(i, k)
          d_q1(i, k) = q1(i, k-1) - q1(i, k)
          d_q2(i, k) = q2(i, k-1) - q2(i, k)
          d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k)
        END IF
      END DO
    END DO

    DO i = 1, klon
      IF (wk_adv(i)) THEN
        omgbdth(i, 1) = 0.
        omgbdq(i, 1) = 0.
      END IF
    END DO

    DO k = 2, klev
      DO i = 1, klon
        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN !   loop on interfaces
          omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k))
          omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k))
        END IF
      END DO
    END DO

    ! -----------------------------------------------------------------
    DO k = 1, klev
      DO i = 1, klon
        IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
          ! -----------------------------------------------------------------

          ! Compute redistribution (advective) term

          d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
            (rrd1*omg(i,k)*sigmaw(i)*d_th1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( &
            i))*d_th2(i,k+1)-(1.-alpha_up(i,k))*omgbdth(i,k)-alpha_up(i,k+1)* &
            omgbdth(i,k+1))*ppi(i, k)
          ! print*,'d_deltatw=',d_deltatw(i,k)

          d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
            (rrd1*omg(i,k)*sigmaw(i)*d_q1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( &
            i))*d_q2(i,k+1)-(1.-alpha_up(i,k))*omgbdq(i,k)-alpha_up(i,k+1)* &
            omgbdq(i,k+1))
          ! print*,'d_deltaqw=',d_deltaqw(i,k)

          ! and increment large scale tendencies




          ! C
          ! -----------------------------------------------------------------
          d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, &
            k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/(ph(i,k)-ph(i, &
            k+1)) &                ! cc nrlmd     $
                                   ! -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k)
            -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/(ph(i, &
            k)-ph(i,k+1)) &        ! cc
            )*ppi(i, k)

          d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, &
            k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/(ph(i,k)-ph(i, &
            k+1)) &                ! cc nrlmd     $
                                   ! -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k)
            -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i, &
            k+1))/(ph(i,k)-ph(i,k+1)) & ! cc
            )
          ! cc nrlmd
        ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
          d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, &
            k)/(ph(i,k)-ph(i,k+1))))*ppi(i, k)

          d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, &
            k)/(ph(i,k)-ph(i,k+1))))

        END IF
        ! cc
      END DO
    END DO
    ! ------------------------------------------------------------------

    ! Increment state variables

    DO k = 1, klev
      DO i = 1, klon
        ! cc nrlmd       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
          ! cc



          ! Coefficient de répartition

          crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
            (ph(i,kupper(i))-ph(i,1))
          crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-ph(i &
            ,kupper(i)))


          ! Reintroduce compensating subsidence term.

          ! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
          ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
          ! .                   /(1-sigmaw)
          ! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
          ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
          ! .                   /(1-sigmaw)

          ! dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
          ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
          ! .                   /(1-sigmaw)
          ! dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
          ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
          ! .                   /(1-sigmaw)

          dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))
          dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))
          ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)

!jyg<
!!
!!---------------------------------------------------------------
!! The change of delta_T due to PBL (vertical diffusion plus thermal plumes)
!! is accounted for by the PBL and the Thermals schemes. It is now set to zero 
!! within the Wake scheme.
!!---------------------------------------------------------------
          dtPBL(i,k) = 0.
          dqPBL(i,k) = 0.
! 
!!          dtPBL(i,k)=wdtPBL(i,k) - udtPBL(i,k)
!!          dqPBL(i,k)=wdqPBL(i,k) - udqPBL(i,k)
! 
!!          dtpbl(i, k) = (wdtpbl(i,k)/sigmaw(i)-udtpbl(i,k)/(1.-sigmaw(i)))
!!          dqpbl(i, k) = (wdqpbl(i,k)/sigmaw(i)-udqpbl(i,k)/(1.-sigmaw(i)))
          ! print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k)
!>jyg
!

          ! cc nrlmd          Prise en compte du taux de mortalité
          ! cc               Définitions de entr, detr
          detr(i, k) = 0.

          entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
            sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)

          spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
          ! cc        spread(i,k) =
          ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
          ! cc     $  sigmaw(i)


          ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU
          ! Jingmei

          ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
          ! &  Tgw(i,k),deltatw(i,k)
          d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)*deltatw(i, k)* &
            dtimesub
          ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
          ff(i) = d_deltatw(i, k)/dtimesub

          ! Sans GW

          ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))

          ! GW formule 1

          ! deltatw(k) = deltatw(k)+dtimesub*
          ! $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))

          ! GW formule 2

          IF (dtimesub*tgw(i,k)<1.E-10) THEN
            d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc
                                                                     ! $
                                                                     ! -spread(i,k)*deltatw(i,k)
              -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( &
              i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
              -tgw(i,k)*deltatw(i,k))
          ELSE
            d_deltatw(i, k) = 1/tgw(i, k)*(1-exp(-dtimesub*tgw(i, &
              k)))*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc     $
                                                 ! -spread(i,k)*deltatw(i,k)
              -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( &
              i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
              -tgw(i,k)*deltatw(i,k))
          END IF

          dth(i, k) = deltatw(i, k)/ppi(i, k)

          gg(i) = d_deltaqw(i, k)/dtimesub

          d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k)+dqpbl(i,k) & ! cc     $
                                                                   ! -spread(i,k)*deltaqw(i,k))
            -entr(i,k)*deltaqw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw(i)+detr( &
            i,k))*deltaqw(i,k)/(1.-sigmaw(i)))
          ! cc

          ! cc nrlmd
          ! cc       d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
          ! cc       d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
          ! cc
        END IF
      END DO
    END DO


    ! Scale tendencies so that water vapour remains positive in w and x.

    CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, &
      d_deltaqw, sigmaw, d_sigmaw, alpha)

    ! cc nrlmd
    ! c      print*,'alpha'
    ! c      do i=1,klon
    ! c         print*,alpha(i)
    ! c      end do
    ! cc
    DO k = 1, klev
      DO i = 1, klon
        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
          d_te(i, k) = alpha(i)*d_te(i, k)
          d_qe(i, k) = alpha(i)*d_qe(i, k)
          d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
          d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
          d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
        END IF
      END DO
    END DO
    DO i = 1, klon
      IF (wk_adv(i)) THEN
        d_sigmaw(i) = alpha(i)*d_sigmaw(i)
      END IF
    END DO

    ! Update large scale variables and wake variables
    ! IM 060208 manque DO i + remplace DO k=1,kupper(i)
    ! IM 060208     DO k = 1,kupper(i)
    DO k = 1, klev
      DO i = 1, klon
        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
          dtls(i, k) = dtls(i, k) + d_te(i, k)
          dqls(i, k) = dqls(i, k) + d_qe(i, k)
          ! cc nrlmd
          d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
          d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
          ! cc
        END IF
      END DO
    END DO
    DO k = 1, klev
      DO i = 1, klon
        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
          te(i, k) = te0(i, k) + dtls(i, k)
          qe(i, k) = qe0(i, k) + dqls(i, k)
          the(i, k) = te(i, k)/ppi(i, k)
          deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
          deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
          dth(i, k) = deltatw(i, k)/ppi(i, k)
          ! c      print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k)
          ! c     $        ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
        END IF
      END DO
    END DO
    DO i = 1, klon
      IF (wk_adv(i)) THEN
        sigmaw(i) = sigmaw(i) + d_sigmaw(i)
      END IF
    END DO


    ! Determine Ptop from buoyancy integral
    ! ---------------------------------------

    ! -     1/ Pressure of the level where dth changes sign.

    DO i = 1, klon
      IF (wk_adv(i)) THEN
        ptop_provis(i) = ph(i, 1)
      END IF
    END DO

    DO k = 2, klev
      DO i = 1, klon
        IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
          ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
            k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
        END IF
      END DO
    END DO

    ! -     2/ dth integral

    DO i = 1, klon
      IF (wk_adv(i)) THEN !!! nrlmd
        sum_dth(i) = 0.
        dthmin(i) = -delta_t_min
        z(i) = 0.
      END IF
    END DO

    DO k = 1, klev
      DO i = 1, klon
        IF (wk_adv(i)) THEN
          dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
          IF (dz(i)>0) THEN
            z(i) = z(i) + dz(i)
            sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
            dthmin(i) = amin1(dthmin(i), dth(i,k))
          END IF
        END IF
      END DO
    END DO

    ! -     3/ height of triangle with area= sum_dth and base = dthmin

    DO i = 1, klon
      IF (wk_adv(i)) THEN
        hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
        hw(i) = amax1(hwmin, hw(i))
      END IF
    END DO

    ! -     4/ now, get Ptop

    DO i = 1, klon
      IF (wk_adv(i)) THEN !!! nrlmd
        ktop(i) = 0
        z(i) = 0.
      END IF
    END DO

    DO k = 1, klev
      DO i = 1, klon
        IF (wk_adv(i)) THEN
          dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i))
          IF (dz(i)>0) THEN
            z(i) = z(i) + dz(i)
            ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i)
            ktop(i) = k
          END IF
        END IF
      END DO
    END DO

    ! 4.5/Correct ktop and ptop

    DO i = 1, klon
      IF (wk_adv(i)) THEN
        ptop_new(i) = ptop(i)
      END IF
    END DO

    DO k = klev, 2, -1
      DO i = 1, klon
        ! IM v3JYG; IF (k .GE. ktop(i)
        IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
          ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
            k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
        END IF
      END DO
    END DO


    DO i = 1, klon
      IF (wk_adv(i)) THEN
        ptop(i) = ptop_new(i)
      END IF
    END DO

    DO k = klev, 1, -1
      DO i = 1, klon
        IF (wk_adv(i)) THEN !!! nrlmd
          IF (ph(i,k+1)<ptop(i)) ktop(i) = k
        END IF
      END DO
    END DO

    ! 5/ Set deltatw & deltaqw to 0 above kupper

    DO k = 1, klev
      DO i = 1, klon
        IF (wk_adv(i) .AND. k>=kupper(i)) THEN
          deltatw(i, k) = 0.
          deltaqw(i, k) = 0.
        END IF
      END DO
    END DO


    ! -------------Cstar computation---------------------------------
    DO i = 1, klon
      IF (wk_adv(i)) THEN !!! nrlmd
        sum_thu(i) = 0.
        sum_tu(i) = 0.
        sum_qu(i) = 0.
        sum_thvu(i) = 0.
        sum_dth(i) = 0.
        sum_dq(i) = 0.
        sum_rho(i) = 0.
        sum_dtdwn(i) = 0.
        sum_dqdwn(i) = 0.

        av_thu(i) = 0.
        av_tu(i) = 0.
        av_qu(i) = 0.
        av_thvu(i) = 0.
        av_dth(i) = 0.
        av_dq(i) = 0.
        av_rho(i) = 0.
        av_dtdwn(i) = 0.
        av_dqdwn(i) = 0.
      END IF
    END DO

    ! Integrals (and wake top level number)
    ! --------------------------------------

    ! Initialize sum_thvu to 1st level virt. pot. temp.

    DO i = 1, klon
      IF (wk_adv(i)) THEN !!! nrlmd
        z(i) = 1.
        dz(i) = 1.
        sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i)
        sum_dth(i) = 0.
      END IF
    END DO

    DO k = 1, klev
      DO i = 1, klon
        IF (wk_adv(i)) THEN !!! nrlmd
          dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
          IF (dz(i)>0) THEN
            z(i) = z(i) + dz(i)
            sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
            sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
            sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
            sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i)
            sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
            sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
            sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
            sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
            sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
          END IF
        END IF
      END DO
    END DO

    DO i = 1, klon
      IF (wk_adv(i)) THEN !!! nrlmd
        hw0(i) = z(i)
      END IF
    END DO


    ! - WAPE and mean forcing computation
    ! ---------------------------------------

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

    ! Means

    DO i = 1, klon
      IF (wk_adv(i)) THEN !!! nrlmd
        av_thu(i) = sum_thu(i)/hw0(i)
        av_tu(i) = sum_tu(i)/hw0(i)
        av_qu(i) = sum_qu(i)/hw0(i)
        av_thvu(i) = sum_thvu(i)/hw0(i)
        av_dth(i) = sum_dth(i)/hw0(i)
        av_dq(i) = sum_dq(i)/hw0(i)
        av_rho(i) = sum_rho(i)/hw0(i)
        av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
        av_dqdwn(i) = sum_dqdwn(i)/hw0(i)

        wape(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* &
          av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
      END IF
    END DO

    ! Filter out bad wakes

    DO k = 1, klev
      DO i = 1, klon
        IF (wk_adv(i)) THEN !!! nrlmd
          IF (wape(i)<0.) THEN
            deltatw(i, k) = 0.
            deltaqw(i, k) = 0.
            dth(i, k) = 0.
          END IF
        END IF
      END DO
    END DO

    DO i = 1, klon
      IF (wk_adv(i)) THEN !!! nrlmd
        IF (wape(i)<0.) THEN
          wape(i) = 0.
          cstar(i) = 0.
          hw(i) = hwmin
          sigmaw(i) = max(sigmad, sigd_con(i))
          fip(i) = 0.
          gwake(i) = .FALSE.
        ELSE
          cstar(i) = stark*sqrt(2.*wape(i))
          gwake(i) = .TRUE.
        END IF
      END IF
    END DO

  END DO ! end sub-timestep loop

  ! -----------------------------------------------------------------
  ! Get back to tendencies per second

  DO k = 1, klev
    DO i = 1, klon

      ! cc nrlmd        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
      IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
        ! cc
        dtls(i, k) = dtls(i, k)/dtime
        dqls(i, k) = dqls(i, k)/dtime
        d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
        d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
        d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
        ! c      print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
        ! c     $         ,death_rate(i)*sigmaw(i)
      END IF
    END DO
  END DO


  ! ----------------------------------------------------------
  ! Determine wake final state; recompute wape, cstar, ktop;
  ! filter out bad wakes.
  ! ----------------------------------------------------------

  ! 2.1 - Undisturbed area and Wake integrals
  ! ---------------------------------------------------------

  DO i = 1, klon
    ! cc nrlmd       if (wk_adv(i)) then !!! nrlmd
    IF (ok_qx_qw(i)) THEN
      ! cc
      z(i) = 0.
      sum_thu(i) = 0.
      sum_tu(i) = 0.
      sum_qu(i) = 0.
      sum_thvu(i) = 0.
      sum_dth(i) = 0.
      sum_dq(i) = 0.
      sum_rho(i) = 0.
      sum_dtdwn(i) = 0.
      sum_dqdwn(i) = 0.

      av_thu(i) = 0.
      av_tu(i) = 0.
      av_qu(i) = 0.
      av_thvu(i) = 0.
      av_dth(i) = 0.
      av_dq(i) = 0.
      av_rho(i) = 0.
      av_dtdwn(i) = 0.
      av_dqdwn(i) = 0.
    END IF
  END DO
  ! Potential temperatures and humidity
  ! ----------------------------------------------------------

  DO k = 1, klev
    DO i = 1, klon
      ! cc nrlmd       IF ( wk_adv(i)) THEN
      IF (ok_qx_qw(i)) THEN
        ! cc
        rho(i, k) = p(i, k)/(rd*te(i,k))
        IF (k==1) THEN
          rhoh(i, k) = ph(i, k)/(rd*te(i,k))
          zhh(i, k) = 0
        ELSE
          rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
          zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
        END IF
        the(i, k) = te(i, k)/ppi(i, k)
        thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
        tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
        qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
        rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
        dth(i, k) = deltatw(i, k)/ppi(i, k)
      END IF
    END DO
  END DO

  ! Integrals (and wake top level number)
  ! -----------------------------------------------------------

  ! Initialize sum_thvu to 1st level virt. pot. temp.

  DO i = 1, klon
    ! cc nrlmd       IF ( wk_adv(i)) THEN
    IF (ok_qx_qw(i)) THEN
      ! cc
      z(i) = 1.
      dz(i) = 1.
      sum_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i)
      sum_dth(i) = 0.
    END IF
  END DO

  DO k = 1, klev
    DO i = 1, klon
      ! cc nrlmd       IF ( wk_adv(i)) THEN
      IF (ok_qx_qw(i)) THEN
        ! cc
        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
        IF (dz(i)>0) THEN
          z(i) = z(i) + dz(i)
          sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
          sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
          sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
          sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+eps*qu(i,k))*dz(i)
          sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
          sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
          sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
          sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
          sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
        END IF
      END IF
    END DO
  END DO

  DO i = 1, klon
    ! cc nrlmd       IF ( wk_adv(i)) THEN
    IF (ok_qx_qw(i)) THEN
      ! cc
      hw0(i) = z(i)
    END IF
  END DO

  ! - WAPE and mean forcing computation
  ! -------------------------------------------------------------

  ! Means

  DO i = 1, klon
    ! cc nrlmd       IF ( wk_adv(i)) THEN
    IF (ok_qx_qw(i)) THEN
      ! cc
      av_thu(i) = sum_thu(i)/hw0(i)
      av_tu(i) = sum_tu(i)/hw0(i)
      av_qu(i) = sum_qu(i)/hw0(i)
      av_thvu(i) = sum_thvu(i)/hw0(i)
      av_dth(i) = sum_dth(i)/hw0(i)
      av_dq(i) = sum_dq(i)/hw0(i)
      av_rho(i) = sum_rho(i)/hw0(i)
      av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
      av_dqdwn(i) = sum_dqdwn(i)/hw0(i)

      wape2(i) = -rg*hw0(i)*(av_dth(i)+eps*(av_thu(i)*av_dq(i)+av_dth(i)* &
        av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
    END IF
  END DO

  ! Prognostic variable update
  ! ------------------------------------------------------------

  ! Filter out bad wakes

  DO k = 1, klev
    DO i = 1, klon
      ! cc nrlmd        IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
      IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
        ! cc
        deltatw(i, k) = 0.
        deltaqw(i, k) = 0.
        dth(i, k) = 0.
      END IF
    END DO
  END DO


  DO i = 1, klon
    ! cc nrlmd       IF ( wk_adv(i)) THEN
    IF (ok_qx_qw(i)) THEN
      ! cc
      IF (wape2(i)<0.) THEN
        wape2(i) = 0.
        cstar2(i) = 0.
        hw(i) = hwmin
        sigmaw(i) = amax1(sigmad, sigd_con(i))
        fip(i) = 0.
        gwake(i) = .FALSE.
      ELSE
        IF (prt_level>=10) PRINT *, 'wape2>0'
        cstar2(i) = stark*sqrt(2.*wape2(i))
        gwake(i) = .TRUE.
      END IF
    END IF
  END DO

  DO i = 1, klon
    ! cc nrlmd       IF ( wk_adv(i)) THEN
    IF (ok_qx_qw(i)) THEN
      ! cc
      ktopw(i) = ktop(i)
    END IF
  END DO

  DO i = 1, klon
    ! cc nrlmd       IF ( wk_adv(i)) THEN
    IF (ok_qx_qw(i)) THEN
      ! cc
      IF (ktopw(i)>0 .AND. gwake(i)) THEN

        ! jyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
        ! cc       heff = 600.
        ! Utilisation de la hauteur hw
        ! c       heff = 0.7*hw
        heff(i) = hw(i)

        fip(i) = 0.5*rho(i, ktopw(i))*cstar2(i)**3*heff(i)*2* &
          sqrt(sigmaw(i)*wdens(i)*3.14)
        fip(i) = alpk*fip(i)
        ! jyg2
      ELSE
        fip(i) = 0.
      END IF
    END IF
  END DO

  ! Limitation de sigmaw

  ! cc nrlmd
  ! DO i=1,klon
  ! IF (OK_qx_qw(i)) THEN
  ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
  ! ENDIF
  ! ENDDO
  ! cc
  DO k = 1, klev
    DO i = 1, klon

      ! cc nrlmd      On maintient désormais constant sigmaw en régime
      ! permanent
      ! cc      IF ((sigmaw(i).GT.sigmaw_max).or.
      IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
          .NOT. ok_qx_qw(i)) THEN
        ! cc
        dtls(i, k) = 0.
        dqls(i, k) = 0.
        deltatw(i, k) = 0.
        deltaqw(i, k) = 0.
      END IF
    END DO
  END DO

  ! cc nrlmd      On maintient désormais constant sigmaw en régime permanent
  DO i = 1, klon
    IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
        .NOT. ok_qx_qw(i)) THEN
      wape(i) = 0.
      cstar(i) = 0.
!!jyg   Outside subroutine "Wake" hw and sigmaw are zero when there are no wakes
!!      hw(i) = hwmin                       !jyg
!!      sigmaw(i) = sigmad                  !jyg
      hw(i) = 0.                            !jyg
      sigmaw(i) = 0.                        !jyg
      fip(i) = 0.
    ELSE
      wape(i) = wape2(i)
      cstar(i) = cstar2(i)
    END IF
    ! c        print*,'wape wape2 ktopw OK_qx_qw =',
    ! c     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
  END DO


  RETURN
END SUBROUTINE wake

SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, &
    d_deltaqw, sigmaw, d_sigmaw, alpha)
  ! ------------------------------------------------------
  ! Dtermination du coefficient alpha tel que les tendances
  ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
  ! a une humidite positive dans la zone (x) et dans la zone (w).
  ! ------------------------------------------------------
  IMPLICIT NONE

  ! Input
  REAL qe(nlon, nl), d_qe(nlon, nl)
  REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
  REAL sigmaw(nlon), d_sigmaw(nlon)
  LOGICAL wk_adv(nlon)
  INTEGER nl, nlon
  ! Output
  REAL alpha(nlon)
  ! Internal variables
  REAL zeta(nlon, nl)
  REAL alpha1(nlon)
  REAL x, a, b, c, discrim
  REAL epsilon
  ! DATA epsilon/1.e-15/
  INTEGER i,k

  DO k = 1, nl
    DO i = 1, nlon
      IF (wk_adv(i)) THEN
        IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
          zeta(i, k) = 0.
        ELSE
          zeta(i, k) = 1.
        END IF
      END IF
    END DO
    DO i = 1, nlon
      IF (wk_adv(i)) THEN
        x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
          (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i)*(deltaqw(i,k)+ &
          d_deltaqw(i,k))
        a = -d_sigmaw(i)*d_deltaqw(i, k)
        b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
          deltaqw(i, k)*d_sigmaw(i)
        c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon
        discrim = b*b - 4.*a*c
        ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
        IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap
          alpha1(i) = 1.
        ELSE
          IF (x>=0.) THEN
            alpha1(i) = 1.
          ELSE
            IF (a>0.) THEN
              alpha1(i) = 0.9*min((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
                ))/(2.*a))
            ELSE IF (a==0.) THEN
              alpha1(i) = 0.9*(-c/b)
            ELSE
              ! print*,'a,b,c discrim',a,b,c discrim
              alpha1(i) = 0.9*max((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
                ))/(2.*a))
            END IF
          END IF
        END IF
        alpha(i) = min(alpha(i), alpha1(i))
      END IF
    END DO
  END DO

  RETURN
END SUBROUTINE wake_vec_modulation