Subroutine WAKE (p,ph,ppi,dtime,sigd_con
     :                ,te0,qe0,omgb,ibas
     :                ,dtdwn,dqdwn,amdwn,amup,dta,dqa
     :                ,wdtPBL,wdqPBL,udtPBL,udqPBL
     o                ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl
     o                ,dtls,dqls
     o                ,ktopw,omgbdth,dp_omgb,wdens
     o                ,tu,qu
     o                ,dtKE,dqKE
     o                ,dtPBL,dqPBL
     o                ,omg,dp_deltomg,spread
     o                ,Cstar,d_deltat_gw
     o                ,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            *
***************************************************************
c
      IMPLICIT none
c============================================================================
C
C
C   But : Decrire le comportement des poches froides apparaissant dans les
C        grands systemes convectifs, et fournir l'energie disponible pour
C        le declenchement de nouvelles colonnes convectives.
C
C   Variables d'etat : deltatw    : ecart de temperature wake-undisturbed area
C                      deltaqw    : ecart d'humidite wake-undisturbed area
C                      sigmaw     : fraction d'aire occupee par la poche.
C
C   Variable de sortie : 
c
c			 wape : WAke Potential Energy
c                        fip  : Front Incident Power (W/m2) - ALP
c                        gfl  : Gust Front Length per unit area (m-1)
C                        dtls : large scale temperature tendency due to wake
C                        dqls : large scale humidity tendency due to wake
C                        hw   : hauteur de la poche
C                     dp_omgb : vertical gradient of large scale omega
C                      omgbdth: flux of Delta_Theta transported by LS omega
C                      dtKE   : differential heating (wake - unpertubed)
C                      dqKE   : differential moistening (wake - unpertubed)
C                      omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
C                 dp_deltomg  : vertical gradient of omg (s-1)
C                     spread  : spreading term in dt_wake and dq_wake
C                 deltatw     : updated temperature difference (T_w-T_u).
C                 deltaqw     : updated humidity difference (q_w-q_u).
C                 sigmaw      : updated wake fractional area.
C                 d_deltat_gw : delta T tendency due to GW
c
C   Variables d'entree : 
c
c		         aire : aire de la maille
c			 te0  : temperature dans l'environnement  (K)
C                        qe0  : humidite dans l'environnement     (kg/kg)
C                        omgb : vitesse verticale moyenne sur la maille (Pa/s)
C                        ibas : cloud base level number
C                        dtdwn: source de chaleur due aux descentes (K/s)
C                        dqdwn: source d'humidite due aux descentes (kg/kg/s)
C			 dta  : source de chaleur due courants satures et detrain  (K/s)
C			 dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
C                        amdwn: flux de masse total des descentes, par unite de
C                                surface de la maille (kg/m2/s)
C                        amup : flux de masse total des ascendances, par unite de
C                                surface de la maille (kg/m2/s)
C                        p    : pressions aux milieux des couches (Pa)
C                        ph   : pressions aux interfaces (Pa)
C                        ppi  : (p/p_0)**kapa (adim)
C                        dtime: increment temporel (s)
c
C   Variables internes :
c
c			 rhow : masse volumique de la poche froide
C                        rho  : environment density at P levels
C                        rhoh : environment density at Ph levels
C                        te   : environment temperature | may change within
C                        qe   : environment humidity    | sub-time-stepping
C                        the  : environment potential temperature
C                        thu  : potential temperature in undisturbed area
C                        tu   :  temperature  in undisturbed area
C                        qu   : humidity in undisturbed area
C                      dp_omgb: vertical gradient og LS omega
C                      omgbw  : wake average vertical omega
C                     dp_omgbw: vertical gradient of omgbw
C                      omgbdq : flux of Delta_q transported by LS omega
C                        dth  : potential temperature diff. wake-undist.
C                        th1  : first pot. temp. for vertical advection (=thu)
C                        th2  : second pot. temp. for vertical advection (=thw)
C                        q1   : first humidity for vertical advection
C                        q2   : second humidity for vertical advection
C                     d_deltatw   : terme de redistribution pour deltatw
C                     d_deltaqw   : terme de redistribution pour deltaqw
C                      deltatw0   : deltatw initial
C                      deltaqw0   : deltaqw initial
C                      hw0    : hw initial
C                      sigmaw0: sigmaw initial
C                      amflux : horizontal mass flux through wake boundary
C                      wdens  : number of wakes per unit area (3D) or per
C                               unit length (2D)
C                      Tgw    : 1 sur la période de onde de gravité
c                      Cgw    : vitesse de propagation de onde de gravité
c                      LL     : distance entre 2 poches

c-------------------------------------------------------------------------
c          Déclaration de variables
c-------------------------------------------------------------------------

#include "dimensions.h"
#include "dimphy.h"
#include "YOMCST.h"
#include "cvthermo.h"

c Arguments en entree
c--------------------

      REAL p(klev),ph(klev+1),ppi(klev)
      REAL dtime
      REAL te0(klev),qe0(klev)
      REAL omgb(klev+1)
      INTEGER ibas
      REAL dtdwn(klev), dqdwn(klev)
      REAL wdtPBL(klev),wdqPBL(klev)
      REAL udtPBL(klev),udqPBL(klev)
      REAL amdwn(klev), amup(klev)
      REAL dta(klev), dqa(klev)
      REAL sigd_con

c Sorties
c--------

      REAL deltatw(klev), deltaqw(klev), dth(klev)
      REAL tu(klev), qu(klev)
      REAL dtls(klev), dqls(klev)
      REAL dtKE(klev), dqKE(klev)
      REAL dtPBL(klev), dqPBL(klev)
      REAL spread(klev)
      REAL d_deltatgw(klev)
      REAL d_deltatw2(klev), d_deltaqw2(klev)
      REAL omgbdth(klev+1), omg(klev+1)
      REAL dp_omgb(klev), dp_deltomg(klev)
      REAL d_deltat_gw(klev)
      REAL hw, sigmaw, wape, fip, gfl, Cstar
      INTEGER ktopw

c Variables internes
c-------------------

c Variables à fixer
      REAL ALON
      REAL coefgw
      REAL wdens0, wdens
      REAL stark
      REAL alpk
      REAL delta_t_min
      REAL Pupper
      INTEGER nsub
      REAL dtimesub
      REAL sigmad, hwmin

c Variables de sauvegarde
      REAL deltatw0(klev)
      REAL deltaqw0(klev)
      REAL te(klev), qe(klev)
      REAL sigmaw0, sigmaw1

c Variables pour les GW
      REAL LL
      REAL N2(klev)
      REAL Cgw(klev)
      REAL Tgw(klev)

c Variables liées au calcul de hw
      REAL ptop_provis, ptop, ptop_new
      REAL sum_dth
      REAL dthmin
      REAL z, dz, hw0
      INTEGER ktop, kupper

c Autres variables internes
      INTEGER isubstep, k

      REAL sum_thu, sum_tu, sum_qu,sum_thvu
      REAL sum_dq, sum_rho
      REAL sum_dtdwn, sum_dqdwn
      REAL av_thu, av_tu, av_qu, av_thvu
      REAL av_dth, av_dq, av_rho
      REAL av_dtdwn, av_dqdwn

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

      REAL pbase
      
      REAL the(klev), thu(klev)

      REAL d_deltatw(klev), d_deltaqw(klev)

      REAL omgbw(klev+1), omgtop
      REAL dp_omgbw(klev)
      REAL ztop, dztop
      REAL alpha_up(klev)
      
      REAL RRe1, RRe2, RRd1, RRd2
      REAL Th1(klev), Th2(klev), q1(klev), q2(klev)
      REAL D_Th1(klev), D_Th2(klev), D_dth(klev)
      REAL D_q1(klev), D_q2(klev), D_dq(klev)
      REAL omgbdq(klev)

      REAL ff, gg
      REAL wape2, Cstar2, heff

      REAL Crep(klev)
      REAL Crep_upper, Crep_sol

C-------------------------------------------------------------------------
c         Initialisations
c-------------------------------------------------------------------------

c      print*, 'wake initialisations'

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

      DATA sigmad, hwmin /.02,10./

C Longueur de maille (en m)
c-------------------------------------------------------------------------

c      ALON = 3.e5
      ALON = 1.e6


C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
c
c      coefgw : Coefficient pour les ondes de gravité
c       stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
c       wdens : Densité de poche froide par maille
c-------------------------------------------------------------------------

      coefgw=10
c      coefgw=1
c      wdens0 = 1.0/(alon**2)   
      wdens = 1.0/(alon**2)       
      stark = 0.50
cCRtest
      alpk=0.1
c      alpk = 1.0 
c      alpk = 0.5
c      alpk = 0.05
      Crep_upper=0.9
      Crep_sol=1.0


C Minimum value for |T_wake - T_undist|. Used for wake top definition
c-------------------------------------------------------------------------

      delta_t_min = 0.2


C Cloud base
c-------------------------------------------------------------------------

      Pbase = P(ibas)


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

      DO k=1,klev
	deltatw0(k) = deltatw(k)
	deltaqw0(k)= deltaqw(k)
	te(k) = te0(k)
	qe(k) = qe0(k)
	dtls(k) = 0.
	dqls(k) = 0.
        d_deltat_gw(k)=0.
        d_deltatw2(k)=0.
        d_deltaqw2(k)=0.
      ENDDO
c      sigmaw1=sigmaw
c      IF (sigd_con.GT.sigmaw1) THEN
c      print*, 'sigmaw,sigd_con', sigmaw, sigd_con
c      ENDIF
      sigmaw = max(sigmaw,sigd_con)
      sigmaw = max(sigmaw,sigmad)
      sigmaw = min(sigmaw,0.99)
      sigmaw0 = sigmaw
c      wdens=wdens0/(10.*sigmaw)
c      IF (sigd_con.GT.sigmaw1) THEN
c      print*, 'sigmaw1,sigd1', sigmaw, sigd_con
c      ENDIF

C 2. - Prognostic part
C =========================================================

c      print *, 'prognostic wake computation'


C 2.1 - Undisturbed area and Wake integrals
C ---------------------------------------------------------

      z = 0.
      ktop=0
      kupper = 0
      sum_thu = 0.
      sum_tu = 0.
      sum_qu = 0.
      sum_thvu = 0.
      sum_dth = 0.
      sum_dq = 0.
      sum_rho = 0.
      sum_dtdwn = 0.
      sum_dqdwn = 0.

      av_thu = 0.
      av_tu =0.
      av_qu =0.
      av_thvu = 0.
      av_dth = 0.
      av_dq = 0.
      av_rho =0.
      av_dtdwn =0.
      av_dqdwn = 0.

C Potential temperatures and humidity
c----------------------------------------------------------

      DO k =1,klev
        rho(k) = p(k)/(rd*te(k))
        IF(k .eq. 1) THEN
          rhoh(k) = ph(k)/(rd*te(k))
          zhh(k)=0
        ELSE
          rhoh(k) = ph(k)*2./(rd*(te(k)+te(k-1)))
          zhh(k)=(ph(k)-ph(k-1))/(-rhoh(k)*RG)+zhh(k-1)
        ENDIF
        the(k) = te(k)/ppi(k)
        thu(k) = (te(k) - deltatw(k)*sigmaw)/ppi(k)
        tu(k) = te(k) - deltatw(k)*sigmaw
        qu(k)  =  qe(k) - deltaqw(k)*sigmaw
        rhow(k) = p(k)/(rd*(te(k)+deltatw(k)))
        dth(k) = deltatw(k)/ppi(k)
        LL = (1-sqrt(sigmaw))/sqrt(wdens)       
      ENDDO
        
      DO k = 1, klev-1
        IF(k.eq.1) THEN
          N2(k)=0
        ELSE
          N2(k)=max(0.,-RG**2/the(k)*rho(k)*(the(k+1)-the(k-1))
     $           /(p(k+1)-p(k-1)))
        ENDIF
        ZH(k)=(zhh(k)+zhh(k+1))/2

        Cgw(k)=sqrt(N2(k))*ZH(k)
        Tgw(k)=coefgw*Cgw(k)/LL
      ENDDO
         
      N2(klev)=0
      ZH(klev)=0
      Cgw(klev)=0
      Tgw(klev)=0

c  Calcul de la masse volumique moyenne de la colonne
c-----------------------------------------------------------------

      DO k=1,klev
        epaisseur1(k)=0.
        epaisseur2(k)=0.
      ENDDO

      epaisseur1(1)= -(Ph(2)-Ph(1))/(rho(1)*rg)+1.
      epaisseur2(1)= -(Ph(2)-Ph(1))/(rho(1)*rg)+1.
      rhow_moyen(1) = rhow(1)

      DO k = 2, klev
        epaisseur1(k)= -(Ph(k+1)-Ph(k))/(rho(k)*rg) +1.
        epaisseur2(k)=epaisseur2(k-1)+epaisseur1(k)
        rhow_moyen(k) = (rhow_moyen(k-1)*epaisseur2(k-1)+
     $                 rhow(k)*epaisseur1(k))/epaisseur2(k)
      ENDDO


C Choose an integration bound well above wake top
c-----------------------------------------------------------------

c       Pupper = 50000.  ! melting level
       Pupper = 60000.
c       Pupper = 70000.


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

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

      Ptop_provis=ph(1)
      DO k= 2,klev
        IF (dth(k) .GT. -delta_t_min .and.
     $      dth(k-1).LT. -delta_t_min) THEN
          Ptop_provis = ((dth(k)+delta_t_min)*p(k-1)
     $          - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
          GO TO 25
        ENDIF
      ENDDO
25    CONTINUE

c-2/ dth integral

      sum_dth = 0.
      dthmin = -delta_t_min
      z = 0.

      DO k = 1,klev
        dz = -(max(ph(k+1),Ptop_provis)-Ph(k))/(rho(k)*rg)
        IF (dz .le. 0) GO TO 40
        z = z+dz
        sum_dth = sum_dth + dth(k)*dz
        dthmin = min(dthmin,dth(k))
      ENDDO
40    CONTINUE

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

      hw0 = 2.*sum_dth/min(dthmin,-0.5)
      hw0 = max(hwmin,hw0)

c-4/ now, get Ptop

      z = 0.
      ptop = ph(1)

      DO k = 1,klev
        dz = min(-(ph(k+1)-Ph(k))/(rho(k)*rg),hw0-z)
        IF (dz .le. 0) GO TO 45
        z = z+dz
        Ptop = Ph(k)-rho(k)*rg*dz
      ENDDO
45    CONTINUE


C-5/ Determination de ktop et kupper

      DO k=klev,1,-1
        IF (ph(k+1) .lt. ptop) ktop=k
        IF (ph(k+1) .lt. pupper) kupper=k
      ENDDO

c-6/ Correct ktop and ptop

      Ptop_new=ptop
      DO k= ktop,2,-1
        IF (dth(k) .GT. -delta_t_min .and.
     $      dth(k-1).LT. -delta_t_min) THEN
          Ptop_new = ((dth(k)+delta_t_min)*p(k-1)
     $          - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
          GO TO 225
        ENDIF
      ENDDO
225   CONTINUE

      ptop = ptop_new

      DO k=klev,1,-1
        IF (ph(k+1) .lt. ptop) ktop=k
      ENDDO

c Set deltatw & deltaqw to 0 above kupper
c-----------------------------------------------------------

      DO k = kupper,klev
        deltatw(k) = 0.
        deltaqw(k) = 0.
      ENDDO


C Vertical gradient of LS omega
C------------------------------------------------------------

      DO k = 1,kupper
        dp_omgb(k) = (omgb(k+1) - omgb(k))/(ph(k+1)-ph(k))
      ENDDO


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

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

      z = 1.
      dz = 1.
      sum_thvu =  thu(1)*(1.+eps*qu(1))*dz
      sum_dth = 0.

      DO k = 1,klev
        dz = -(max(ph(k+1),Ptop)-Ph(k))/(rho(k)*rg)
        IF (dz .LE. 0) GO TO 50
        z = z+dz
        sum_thu = sum_thu + thu(k)*dz
        sum_tu = sum_tu + tu(k)*dz
        sum_qu = sum_qu + qu(k)*dz
        sum_thvu = sum_thvu + thu(k)*(1.+eps*qu(k))*dz
        sum_dth = sum_dth + dth(k)*dz
        sum_dq = sum_dq + deltaqw(k)*dz
        sum_rho = sum_rho + rhow(k)*dz
        sum_dtdwn = sum_dtdwn + dtdwn(k)*dz
        sum_dqdwn = sum_dqdwn + dqdwn(k)*dz
      ENDDO
50    CONTINUE

      hw0 = z

C 2.1 - WAPE and mean forcing computation
C-------------------------------------------------------------

C Means

      av_thu = sum_thu/hw0
      av_tu = sum_tu/hw0
      av_qu = sum_qu/hw0
      av_thvu = sum_thvu/hw0
c      av_thve = sum_thve/hw0
      av_dth = sum_dth/hw0
      av_dq = sum_dq/hw0
      av_rho = sum_rho/hw0
      av_dtdwn = sum_dtdwn/hw0
      av_dqdwn = sum_dqdwn/hw0

      wape = - rg*hw0*(av_dth
     $     + eps*(av_thu*av_dq+av_dth*av_qu+av_dth*av_dq ))/av_thvu

C 2.2 Prognostic variable update
C ------------------------------------------------------------

C Filter out bad wakes

      IF ( wape .LT. 0.) THEN
        print*,'wape<0'
        wape = 0.
        hw = hwmin
        sigmaw = max(sigmad,sigd_con)
        fip = 0.
        DO k = 1,klev
          deltatw(k) = 0.
          deltaqw(k) = 0.
          dth(k) = 0.
        ENDDO
      ELSE
        print*,'wape>0'
        Cstar = stark*sqrt(2.*wape)
      ENDIF

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

c      nsub=36
      nsub=10
      dtimesub=dtime/nsub

c------------------------------------------------------------
      DO isubstep = 1,nsub
c------------------------------------------------------------

c        print*,'---------------','substep=',isubstep,'-------------'

c  Evolution of sigmaw


        gfl = 2.*sqrt(3.14*wdens*sigmaw)            

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

c calcul de la difference de vitesse verticale poche - zone non perturbee

        z= 0.
        dp_deltomg(1:klev)=0.
        omg(1:klev+1)=0.

        omg(1) = 0.
        dp_deltomg(1) = -(gfl*Cstar)/(sigmaw * (1-sigmaw))

        DO k=2,ktop
          dz = -(Ph(k)-Ph(k-1))/(rho(k-1)*rg)
          z = z+dz
          dp_deltomg(k)= dp_deltomg(1)
          omg(k)= dp_deltomg(1)*z
        ENDDO

        dztop=-(Ptop-Ph(ktop))/(rho(ktop)*rg)
        ztop = z+dztop
        omgtop=dp_deltomg(1)*ztop


c Conversion de la vitesse verticale de m/s a Pa/s

        omgtop = -rho(ktop)*rg*omgtop
        dp_deltomg(1) = omgtop/(ptop-ph(1))

        DO k = 1,ktop
          omg(k) = - rho(k)*rg*omg(k)
          dp_deltomg(k) = dp_deltomg(1)
        ENDDO

c   raccordement lineaire de omg de ptop a pupper

      IF (kupper .GT. ktop) THEN
        omg(kupper+1) = - Rg*amdwn(kupper+1)/sigmaw
     $                + Rg*amup(kupper+1)/(1.-sigmaw)
        dp_deltomg(kupper) = (omgtop-omg(kupper+1))/(Ptop-Pupper)
        DO k=ktop+1,kupper
          dp_deltomg(k) = dp_deltomg(kupper)
          omg(k) = omgtop+(ph(k)-Ptop)*dp_deltomg(kupper)
        ENDDO
      ENDIF

c   Compute wake average vertical velocity omgbw

      DO k = 1,klev+1
        omgbw(k) = omgb(k)+(1.-sigmaw)*omg(k)
      ENDDO

c  and its vertical gradient dp_omgbw

      DO k = 1,klev
        dp_omgbw(k) = (omgbw(k+1)-omgbw(k))/(ph(k+1)-ph(k))
      ENDDO


c  Upstream coefficients for omgb velocity
c--    (alpha_up(k) is the coefficient of the value at level k)
c--    (1-alpha_up(k) is the coefficient of the value at level k-1)

      DO k = 1,klev
       alpha_up(k) = 0.
       IF (omgb(k) .GT. 0.) alpha_up(k) = 1.
      ENDDO

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

      RRe1 = 1.-sigmaw
      RRe2 = sigmaw
      RRd1 = -1.
      RRd2 = 1.

c Get [Th1,Th2], dth and [q1,q2]

      DO k = 1,kupper+1
        dth(k) = deltatw(k)/ppi(k)
        Th1(k) = the(k) - sigmaw     *dth(k)   ! undisturbed area
        Th2(k) = the(k) + (1.-sigmaw)*dth(k)   ! wake
        q1(k) = qe(k) - sigmaw     *deltaqw(k) ! undisturbed area
        q2(k) = qe(k) + (1.-sigmaw)*deltaqw(k) ! wake
      ENDDO

      D_Th1(1) = 0.
      D_Th2(1) = 0.
      D_dth(1) = 0.
      D_q1(1) = 0.
      D_q2(1) = 0.
      D_dq(1) = 0.

      DO k = 2,kupper+1 !   loop on interfaces
        D_Th1(k) = Th1(k-1)-Th1(k)
        D_Th2(k) = Th2(k-1)-Th2(k)
        D_dth(k) = dth(k-1)-dth(k)
        D_q1(k) = q1(k-1)-q1(k)
        D_q2(k) = q2(k-1)-q2(k)
        D_dq(k) = deltaqw(k-1)-deltaqw(k)
      ENDDO

      omgbdth(1) = 0.
      omgbdq(1) = 0.

      DO k = 2,kupper+1  !   loop on interfaces
        omgbdth(k) = omgb(k)*(    dth(k-1) -     dth(k))
        omgbdq(k)  = omgb(k)*(deltaqw(k-1) - deltaqw(k))
      ENDDO


c-----------------------------------------------------------------
      DO k=1,kupper-1
c-----------------------------------------------------------------
c
c   Compute redistribution (advective) term
c
         d_deltatw(k) =
     $             dtimesub/(Ph(k)-Ph(k+1))*(
     $       RRd1*omg(k  )*sigmaw     *D_Th1(k)
     $      -RRd2*omg(k+1)*(1.-sigmaw)*D_Th2(k+1)
     $      -(1.-alpha_up(k))*omgbdth(k) - alpha_up(k+1)*omgbdth(k+1)
     $                      )*ppi(k)
c         print*,'d_deltatw=',d_deltatw(k)
c
         d_deltaqw(k) =
     $             dtimesub/(Ph(k)-Ph(k+1))*(
     $       RRd1*omg(k  )*sigmaw     *D_q1(k)
     $      -RRd2*omg(k+1)*(1.-sigmaw)*D_q2(k+1)
     $      -(1.-alpha_up(k))*omgbdq(k) - alpha_up(k+1)*omgbdq(k+1)
     $                      )
c         print*,'d_deltaqw=',d_deltaqw(k)
c
c   and increment large scale tendencies
c
         dtls(k) = dtls(k) +
     $               dtimesub*(
     $        ( RRe1*omg(k  )*sigmaw     *D_Th1(k)
     $         -RRe2*omg(k+1)*(1.-sigmaw)*D_Th2(k+1) )
     $               /(Ph(k)-Ph(k+1))
     $         -sigmaw*(1.-sigmaw)*dth(k)*dp_deltomg(k)
     $                      )*ppi(k)
c         print*,'dtls=',dtls(k)
c
         dqls(k) = dqls(k) +
     $               dtimesub*(
     $        ( RRe1*omg(k  )*sigmaw     *D_q1(k)
     $         -RRe2*omg(k+1)*(1.-sigmaw)*D_q2(k+1) )
     $               /(Ph(k)-Ph(k+1))
     $         -sigmaw*(1.-sigmaw)*deltaqw(k)*dp_deltomg(k)
     $                      )
c         print*,'dqls=',dqls(k)

c-------------------------------------------------------------------
      ENDDO
c------------------------------------------------------------------

C   Increment state variables

      DO k = 1,kupper-1

c Coefficient de répartition

        Crep(k)=Crep_sol*(ph(kupper)-ph(k))/(ph(kupper)-ph(1))
        Crep(k)=Crep(k)+Crep_upper*(ph(1)-ph(k))/(p(1)-ph(kupper))
        

c Reintroduce compensating subsidence term.

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

        dtKE(k)=(dtdwn(k)/sigmaw - dta(k)/(1.-sigmaw))
        dqKE(k)=(dqdwn(k)/sigmaw - dqa(k)/(1.-sigmaw))
c        print*,'dtKE=',dtKE(k)
c        print*,'dqKE=',dqKE(k)
c
        dtPBL(k)=(wdtPBL(k)/sigmaw - udtPBL(k)/(1.-sigmaw))
        dqPBL(k)=(wdqPBL(k)/sigmaw - udqPBL(k)/(1.-sigmaw))
c
        spread(k) = (1.-sigmaw)*dp_deltomg(k)+gfl*Cstar/sigmaw
c        print*,'spread=',spread(k)


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

        d_deltat_gw(k)=d_deltat_gw(k)-Tgw(k)*deltatw(k)* dtimesub
c        print*,'d_delta_gw=',d_deltat_gw(k)
        ff=d_deltatw(k)/dtimesub

c Sans GW
c
c        deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k)) 
c
c GW formule 1
c
c        deltatw(k) = deltatw(k)+dtimesub*
c     $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
c
c GW formule 2

        IF (dtimesub*Tgw(k).lt.1.e-10) THEN
          deltatw(k) = deltatw(k)+dtimesub*
     $          (ff+dtKE(k)+dtPBL(k) 
     $          - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
        ELSE
           deltatw(k) = deltatw(k)+1/Tgw(k)*(1-exp(-dtimesub*Tgw(k)))*
     $          (ff+dtKE(k)+dtPBL(k)
     $          - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
        ENDIF
   
        dth(k) = deltatw(k)/ppi(k)

        gg=d_deltaqw(k)/dtimesub

       deltaqw(k) = deltaqw(k) +
     $         dtimesub*(gg+ dqKE(k)+dqPBL(k) - spread(k)*deltaqw(k))

       d_deltatw2(k)=d_deltatw2(k)+d_deltatw(k)
       d_deltaqw2(k)=d_deltaqw2(k)+d_deltaqw(k)
      ENDDO

C   And update large scale variables

      DO k = 1,kupper
        te(k) = te0(k) + dtls(k)
        qe(k) = qe0(k) + dqls(k)
        the(k) = te(k)/ppi(k)
      ENDDO

c     Determine Ptop from buoyancy integral
c----------------------------------------------------------------------

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

      Ptop_provis=ph(1)

      DO k= 2,klev
        IF (dth(k) .GT. -delta_t_min .and.
     $      dth(k-1).LT. -delta_t_min) THEN
          Ptop_provis = ((dth(k)+delta_t_min)*p(k-1)
     $          - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
        GO TO 65
        ENDIF
      ENDDO
65    CONTINUE

c-2/ dth integral

      sum_dth = 0.
      dthmin = -delta_t_min
      z = 0.

      DO k = 1,klev
        dz = -(max(ph(k+1),Ptop_provis)-Ph(k))/(rho(k)*rg)
        IF (dz .le. 0) GO TO 70
        z = z+dz
        sum_dth = sum_dth + dth(k)*dz
        dthmin = min(dthmin,dth(k))
      ENDDO
70    CONTINUE

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

      hw = 2.*sum_dth/min(dthmin,-0.5)
      hw = max(hwmin,hw)

c-4/ now, get Ptop

      ktop = 0
      z=0.

      DO k = 1,klev
        dz = min(-(ph(k+1)-Ph(k))/(rho(k)*rg),hw-z)
        IF (dz .le. 0) GO TO 75
        z = z+dz
        Ptop = Ph(k)-rho(k)*rg*dz
        ktop = k
      ENDDO
75    CONTINUE

c-5/Correct ktop and ptop

      Ptop_new=ptop

      DO k= ktop,2,-1
        IF (dth(k) .GT. -delta_t_min .and.
     $      dth(k-1).LT. -delta_t_min) THEN
          Ptop_new = ((dth(k)+delta_t_min)*p(k-1)
     $          - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
          GO TO 275
        ENDIF
      ENDDO
275   CONTINUE

      ptop = ptop_new

      DO k=klev,1,-1
        IF (ph(k+1) .LT. ptop) ktop=k
      ENDDO

c-6/ Set deltatw & deltaqw to 0 above kupper

      DO k = kupper,klev
        deltatw(k) = 0.
        deltaqw(k) = 0.
      ENDDO

c------------------------------------------------------------------
      ENDDO      ! end sub-timestep loop
C -----------------------------------------------------------------

c   Get back to tendencies per second

      DO k = 1,kupper-1
        dtls(k) = dtls(k)/dtime
        dqls(k) = dqls(k)/dtime
        d_deltatw2(k)=d_deltatw2(k)/dtime
        d_deltaqw2(k)=d_deltaqw2(k)/dtime
        d_deltat_gw(k) = d_deltat_gw(k)/dtime
      ENDDO

C 2.1 - Undisturbed area and Wake integrals
C ---------------------------------------------------------

      z = 0.
      sum_thu = 0.
      sum_tu = 0.
      sum_qu = 0.
      sum_thvu = 0.
      sum_dth = 0.
      sum_dq = 0.
      sum_rho = 0.
      sum_dtdwn = 0.
      sum_dqdwn = 0.

      av_thu = 0.
      av_tu =0.
      av_qu =0.
      av_thvu = 0.
      av_dth = 0.
      av_dq = 0.
      av_rho =0.
      av_dtdwn =0.
      av_dqdwn = 0.

C Potential temperatures and humidity
c----------------------------------------------------------

      DO k =1,klev
        rho(k) = p(k)/(rd*te(k))
        IF(k .eq. 1) THEN
          rhoh(k) = ph(k)/(rd*te(k))
          zhh(k)=0
        ELSE
          rhoh(k) = ph(k)*2./(rd*(te(k)+te(k-1)))
          zhh(k)=(ph(k)-ph(k-1))/(-rhoh(k)*RG)+zhh(k-1)
        ENDIF
        the(k) = te(k)/ppi(k)
        thu(k) = (te(k) - deltatw(k)*sigmaw)/ppi(k)
        tu(k) = te(k) - deltatw(k)*sigmaw
        qu(k)  =  qe(k) - deltaqw(k)*sigmaw
        rhow(k) = p(k)/(rd*(te(k)+deltatw(k)))
        dth(k) = deltatw(k)/ppi(k)
       
      ENDDO

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

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

      z = 1.
      dz = 1.
      sum_thvu =  thu(1)*(1.+eps*qu(1))*dz
      sum_dth = 0.

      DO k = 1,klev
        dz = -(max(ph(k+1),Ptop)-Ph(k))/(rho(k)*rg)

        IF (dz .LE. 0) GO TO 51
        z = z+dz
        sum_thu = sum_thu + thu(k)*dz
        sum_tu = sum_tu + tu(k)*dz
        sum_qu = sum_qu + qu(k)*dz
        sum_thvu = sum_thvu + thu(k)*(1.+eps*qu(k))*dz
        sum_dth = sum_dth + dth(k)*dz
        sum_dq = sum_dq + deltaqw(k)*dz
        sum_rho = sum_rho + rhow(k)*dz
        sum_dtdwn = sum_dtdwn + dtdwn(k)*dz
        sum_dqdwn = sum_dqdwn + dqdwn(k)*dz
      ENDDO
 51   CONTINUE

      hw0 = z

C 2.1 - WAPE and mean forcing computation
C-------------------------------------------------------------

C Means

      av_thu = sum_thu/hw0
      av_tu = sum_tu/hw0
      av_qu = sum_qu/hw0
      av_thvu = sum_thvu/hw0
      av_dth = sum_dth/hw0
      av_dq = sum_dq/hw0
      av_rho = sum_rho/hw0
      av_dtdwn = sum_dtdwn/hw0
      av_dqdwn = sum_dqdwn/hw0

      wape2 = - rg*hw0*(av_dth
     $     + eps*(av_thu*av_dq+av_dth*av_qu+av_dth*av_dq ))/av_thvu


C 2.2 Prognostic variable update
C ------------------------------------------------------------

C Filter out bad wakes

      IF ( wape2 .LT. 0.) THEN
        print*,'wape2<0'
        wape2 = 0.
        hw = hwmin
        sigmaw = max(sigmad,sigd_con)
        fip = 0.
        DO k = 1,klev
          deltatw(k) = 0.
          deltaqw(k) = 0.
          dth(k) = 0.
        ENDDO
      ELSE
        print*,'wape2>0'
        Cstar2 = stark*sqrt(2.*wape2)

      ENDIF

      ktopw = ktop

      IF (ktopw .gt. 0) then

Cjyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
ccc       heff = 600.
C      Utilisation de la hauteur hw
cc       heff = 0.7*hw
       heff = hw

       FIP = 0.5*rho(ktopw)*Cstar2**3*heff*2*sqrt(sigmaw*wdens*3.14)
       FIP = alpk * FIP
Cjyg2
       ELSE
         FIP = 0.
       ENDIF


C   Limitation de sigmaw
c
C   sécurité : si le wake occuppe plus de 90 % de la surface de la maille,
C              alors il disparait en se mélangeant à la partie undisturbed

      IF ((sigmaw.GT.0.9).or.
     .     ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
c      IF (sigmaw.GT.0.9) THEN
        DO k = 1,klev
          dtls(k) = 0.
          dqls(k) = 0.
          deltatw(k) = 0.
          deltaqw(k) = 0.
        ENDDO
        wape = 0.
        hw = hwmin
        sigmaw = sigmad
        fip = 0.
      ENDIF

      RETURN
      END