Index: LMDZ5/trunk/libf/phylmd/wake.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/wake.F	(revision 1907)
+++ LMDZ5/trunk/libf/phylmd/wake.F90	(revision 1992)
@@ -1,373 +1,1496 @@
-!
+
 ! $Id$
-!
-      Subroutine WAKE (p,ph,pi,dtime,sigd_con
-     :                ,te0,qe0,omgb
-     :                ,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
-      use dimphy
-      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                     wdens   : densite de poches
-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                        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                        pi  : (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_ref: initial number of wakes per unit area (3D) or per
-C                               unit length (2D), at the beginning of each time step
-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 "YOMCST.h"
-      include "cvthermo.h"
-      include "iniprint.h"
-
-c Arguments en entree
-c--------------------
-
-      REAL, dimension(klon,klev) :: p, pi
-      REAL, dimension(klon,klev+1) :: ph, omgb
-      REAL dtime
-      REAL, dimension(klon,klev) :: te0,qe0
-      REAL, dimension(klon,klev) :: dtdwn, dqdwn
-      REAL, dimension(klon,klev) :: wdtPBL,wdqPBL
-      REAL, dimension(klon,klev) :: udtPBL,udqPBL
-      REAL, dimension(klon,klev) :: amdwn, amup
-      REAL, dimension(klon,klev) :: dta, dqa
-      REAL, dimension(klon) :: sigd_con
-
-c Sorties
-c--------
-
-      REAL, dimension(klon,klev) :: deltatw, deltaqw, dth
-      REAL, dimension(klon,klev) :: tu, qu
-      REAL, dimension(klon,klev) :: dtls, dqls
-      REAL, dimension(klon,klev) :: dtKE, dqKE
-      REAL, dimension(klon,klev) :: dtPBL, dqPBL
-      REAL, dimension(klon,klev) :: spread
-      REAL, dimension(klon,klev) :: d_deltatgw
-      REAL, dimension(klon,klev) :: d_deltatw2, d_deltaqw2
-      REAL, dimension(klon,klev+1) :: omgbdth, omg
-      REAL, dimension(klon,klev) :: dp_omgb, dp_deltomg
-      REAL, dimension(klon,klev) :: d_deltat_gw
-      REAL, dimension(klon) :: hw, sigmaw, wape, fip, gfl, Cstar
-      REAL, dimension(klon) :: wdens
-      INTEGER, dimension(klon) :: ktopw
-
-c Variables internes
-c-------------------
-
-c Variables à fixer
-      REAL ALON
-      REAL coefgw
-      REAL :: wdens_ref
-      REAL stark
-      REAL alpk
-      REAL delta_t_min
-      INTEGER nsub
-      REAL dtimesub
-      REAL sigmad, hwmin,wapecut
-      REAL :: sigmaw_max
-      REAL :: dens_rate
-      REAL wdens0
-cIM 080208
-      LOGICAL, dimension(klon) :: gwake
-
-c Variables de sauvegarde
-      REAL, dimension(klon,klev) :: deltatw0
-      REAL, dimension(klon,klev) :: deltaqw0
-      REAL, dimension(klon,klev) :: te, qe
-      REAL, dimension(klon) :: sigmaw0, sigmaw1
-
-c Variables pour les GW
-      REAL, DIMENSION(klon) :: LL
-      REAL, dimension(klon,klev) :: N2
-      REAL, dimension(klon,klev) :: Cgw
-      REAL, dimension(klon,klev) :: Tgw
-
-c 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
-
-c 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/
-
-c 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 Crep_upper, Crep_sol
-
-      REAL, DIMENSION(klon,klev) :: ppi
-
-ccc nrlmd
-      real, dimension(klon) :: death_rate,nat_rate
-      real, dimension(klon,klev) :: entr
-      real, dimension(klon,klev) :: detr
-
-C-------------------------------------------------------------------------
-c         Initialisations
-c-------------------------------------------------------------------------
-
-c      print*, 'wake initialisations'
-
-c   Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
-c-------------------------------------------------------------------------
-
-      DATA wapecut,sigmad, hwmin /5.,.02,10./
-ccc nrlmd
-      DATA sigmaw_max /0.4/
-      DATA dens_rate /0.1/
-ccc
-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-------------------------------------------------------------------------
-
-ccc nrlmd      coefgw=10
-c      coefgw=1
-c      wdens0 = 1.0/(alon**2)
-ccc nrlmd      wdens = 1.0/(alon**2)
-ccc nrlmd      stark = 0.50
-cCRtest
-ccc nrlmd      alpk=0.1
-c      alpk = 1.0
-c      alpk = 0.5
-c      alpk = 0.05
-c
-       stark  = 0.33
-       Alpk   = 0.25
-       wdens_ref  = 8.e-12
-       coefgw = 4.
-      Crep_upper=0.9
-      Crep_sol=1.0
-
-ccc nrlmd Lecture du fichier wake_param.data
-      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
-c
-c   Initialisation de toutes des densites a wdens_ref.
-c   Les densites peuvent evoluer si les poches debordent
-c   (voir au tout debut de la boucle sur les substeps)
-      wdens = wdens_ref
-c
-c      print*,'stark',stark
-c      print*,'alpk',alpk
-c      print*,'wdens',wdens
-c      print*,'coefgw',coefgw
-ccc
-C Minimum value for |T_wake - T_undist|. Used for wake top definition
-c-------------------------------------------------------------------------
-
-      delta_t_min = 0.2
-
-C 1. - Save initial values and initialize tendencies
-C --------------------------------------------------
-
-      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
-      ENDDO
-      ENDDO
-c      sigmaw1=sigmaw
-c      IF (sigd_con.GT.sigmaw1) THEN
-c      print*, 'sigmaw,sigd_con', sigmaw, sigd_con
-c      ENDIF
-      DO i=1, klon
-cc      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)
+
+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
+  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 "dimensions.h"
+  include "YOMCST.h"
+  include "cvthermo.h"
+  include "iniprint.h"
+
+  ! Arguments en entree
+  ! --------------------
+
+  REAL, DIMENSION (klon, klev) :: p, pi
+  REAL, DIMENSION (klon, klev+1) :: ph, omgb
+  REAL dtime
+  REAL, DIMENSION (klon, klev) :: te0, qe0
+  REAL, DIMENSION (klon, klev) :: dtdwn, dqdwn
+  REAL, DIMENSION (klon, klev) :: wdtpbl, wdqpbl
+  REAL, DIMENSION (klon, klev) :: udtpbl, udqpbl
+  REAL, DIMENSION (klon, klev) :: amdwn, amup
+  REAL, DIMENSION (klon, klev) :: dta, dqa
+  REAL, DIMENSION (klon) :: sigd_con
+
+  ! Sorties
+  ! --------
+
+  REAL, DIMENSION (klon, klev) :: deltatw, deltaqw, dth
+  REAL, DIMENSION (klon, klev) :: tu, qu
+  REAL, DIMENSION (klon, klev) :: dtls, dqls
+  REAL, DIMENSION (klon, klev) :: dtke, dqke
+  REAL, DIMENSION (klon, klev) :: dtpbl, dqpbl
+  REAL, DIMENSION (klon, klev) :: spread
+  REAL, DIMENSION (klon, klev) :: d_deltatgw
+  REAL, DIMENSION (klon, klev) :: d_deltatw2, d_deltaqw2
+  REAL, DIMENSION (klon, klev+1) :: omgbdth, omg
+  REAL, DIMENSION (klon, klev) :: dp_omgb, dp_deltomg
+  REAL, DIMENSION (klon, klev) :: d_deltat_gw
+  REAL, DIMENSION (klon) :: hw, sigmaw, wape, fip, gfl, cstar
+  REAL, DIMENSION (klon) :: wdens
+  INTEGER, DIMENSION (klon) :: ktopw
+
+  ! Variables internes
+  ! -------------------
+
+  ! Variables à fixer
+  REAL alon
+  REAL coefgw
+  REAL :: wdens_ref
+  REAL stark
+  REAL alpk
+  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 crep_upper, crep_sol
+
+  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
+
+  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
+  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
+
+  ! 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.
-      wape2(i) = 0.
-      d_sigmaw(i) = 0.
-      ktopw(i) = 0
-      ENDDO
-C
-C
-C 2. - Prognostic part
-C --------------------
-C
-C
-C 2.1 - Undisturbed area and Wake integrals
-C ---------------------------------------------------------
-
-      DO i=1, klon
+      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)
+
+          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)
+
+          ! 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.
-      ktop(i)=0
-      kupper(i) = 0
       sum_thu(i) = 0.
       sum_tu(i) = 0.
@@ -381,289 +1504,99 @@
 
       av_thu(i) = 0.
-      av_tu(i) =0.
-      av_qu(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_rho(i) = 0.
+      av_dtdwn(i) = 0.
       av_dqdwn(i) = 0.
-      ENDDO
-c
-c Distance between wakes
-       DO i = 1,klon
-        LL(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
-       ENDDO
-C Potential temperatures and humidity
-c----------------------------------------------------------
-      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 .eq. 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
+    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
-!          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)
-        ENDIF
-!          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)
-       ENDDO
-      ENDDO
-        
-      DO k = 1, klev-1
-      DO i=1, klon 
-        IF(k.eq.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)))
-        ENDIF
-        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)
-      ENDDO
-      ENDDO
-
-      DO i=1, klon
-      N2(i,klev)=0
-      ZH(i,klev)=0
-      Cgw(i,klev)=0
-      Tgw(i,klev)=0
-      ENDDO
-
-c  Calcul de la masse volumique moyenne de la colonne   (bdlmd)
-c-----------------------------------------------------------------
-
-      DO k=1,klev
-       DO i=1, klon
-        epaisseur1(i,k)=0.
-        epaisseur2(i,k)=0.
-       ENDDO
-      ENDDO
-
-      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)
-      ENDDO
-
-      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)
-      ENDDO
-      ENDDO
-
-C
-C Choose an integration bound well above wake top
-c-----------------------------------------------------------------
-c
-C       Pupper = 50000.  ! melting level
-c       Pupper = 60000.
-c       Pupper = 80000.  ! essais pour case_e
-       DO i = 1,klon
-        Pupper(i) = 0.6*ph(i,1)
-        Pupper(i) = max(Pupper(i), 45000.)
-ccc        Pupper(i) = 60000.
-       ENDDO
-
-C
-C    Determine Wake top pressure (Ptop) from buoyancy integral
-C    --------------------------------------------------------
-c
-c-1/ Pressure of the level where dth becomes less than delta_t_min.
-
-      DO i=1,klon
-      ptop_provis(i)=ph(i,1)
-      ENDDO
-      DO k= 2,klev
-      DO i=1,klon
-c
-cIM v3JYG; ptop_provis(i).LT. ph(i,1)
-c
-        IF (dth(i,k) .GT. -delta_t_min .and.
-     $      dth(i,k-1).LT. -delta_t_min .and.
-     $      ptop_provis(i).EQ. 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))
-        ENDIF
-      ENDDO
-      ENDDO
-
-c-2/ dth integral
-
-      DO i=1,klon
-      sum_dth(i) = 0.
-      dthmin(i) = -delta_t_min
-      z(i) = 0.
-      ENDDO
-
-      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) .gt. 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))
-        ENDIF
-      ENDDO
-      ENDDO
-
-c-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))
-      ENDDO
-
-c-4/ now, get Ptop
-
-      DO i=1,klon
-      z(i) = 0.
-      ptop(i) = ph(i,1)
-      ENDDO
-
-      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) .gt. 0) THEN
-         z(i) = z(i)+dz(i)
-         ptop(i) = ph(i,k)-rho(i,k)*rg*dz(i)
-        ENDIF
-      ENDDO
-      ENDDO
-
-
-C-5/ Determination de ktop et kupper
-
-      DO k=klev,1,-1
-      DO i=1,klon
-        IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
-        IF (ph(i,k+1) .lt. pupper(i)) kupper(i)=k
-      ENDDO
-      ENDDO
-
-c      On evite kupper = 1 et kupper = klev
-      DO i=1,klon
-        kupper(i) = max(kupper(i),2)
-        kupper(i) = min(kupper(i),klev-1)
-      ENDDO
-
-
-c-6/ Correct ktop and ptop
-
-      DO i = 1,klon
-        ptop_new(i)=ptop(i)
-      ENDDO
-      DO k= klev,2,-1
-      DO i=1,klon
-        IF (k .LE. ktop(i) .and.
-     $      ptop_new(i) .EQ. ptop(i) .and.
-     $      dth(i,k) .GT. -delta_t_min .and.
-     $      dth(i,k-1).LT. -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))
-        ENDIF
-      ENDDO
-      ENDDO
-
-      DO i=1,klon
-        ptop(i) = ptop_new(i)
-      ENDDO
-
-      DO k=klev,1,-1
-      DO i=1,klon
-        IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
-      ENDDO
-      ENDDO
-c
-c-5/ Set deltatw & deltaqw to 0 above kupper
-c
-      DO k = 1,klev
-      DO i=1,klon
-       IF (k.GE. kupper(i)) THEN
-        deltatw(i,k) = 0.
-        deltaqw(i,k) = 0.
-       ENDIF
-      ENDDO
-      ENDDO
-c
-C
-C Vertical gradient of LS omega
-C
-      DO k = 1,klev
-      DO i=1,klon
-       IF (k.LE. kupper(i)) THEN
-        dp_omgb(i,k) = (omgb(i,k+1) - omgb(i,k))/(ph(i,k+1)-ph(i,k))
-       ENDIF
-      ENDDO
-      ENDDO
-C
-C Integrals (and wake top level number)
-C --------------------------------------
-C
-C Initialize sum_thvu to 1st level virt. pot. temp.
-
-      DO i=1,klon
+          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_thvu(i) = thu(i, 1)*(1.+eps*qu(i,1))*dz(i)
       sum_dth(i) = 0.
-      ENDDO
-
-      DO k = 1,klev
-      DO i=1,klon
+    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) .GT. 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)
-        ENDIF
-      ENDDO
-      ENDDO
-c
-      DO i=1,klon
-        hw0(i) = z(i)
-      ENDDO
-c
-C
-C 2.1 - WAPE and mean forcing computation
-C ---------------------------------------
-C
-C ---------------------------------------
-C
-C Means
-
-      DO i=1,klon
+        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)
-c      av_thve = sum_thve/hw0
       av_dth(i) = sum_dth(i)/hw0(i)
       av_dq(i) = sum_dq(i)/hw0(i)
@@ -672,2225 +1605,1231 @@
       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)
-      ENDDO
-C
-C 2.2 Prognostic variable update
-C ------------------------------
-C
-C Filter out bad wakes
-
-      DO k = 1,klev
-       DO i=1,klon
-        IF ( wape(i) .LT. 0.) THEN
-          deltatw(i,k) = 0.
-          deltaqw(i,k) = 0.
-          dth(i,k) = 0.
-        ENDIF
-       ENDDO
-      ENDDO
-c
-      DO i=1,klon
-      IF ( wape(i) .LT. 0.) THEN
-        wape(i) = 0.
-        Cstar(i) = 0.
+      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))
+        sigmaw(i) = amax1(sigmad, sigd_con(i))
         fip(i) = 0.
         gwake(i) = .FALSE.
       ELSE
-        Cstar(i) = stark*sqrt(2.*wape(i))
+        IF (prt_level>=10) PRINT *, 'wape2>0'
+        cstar2(i) = stark*sqrt(2.*wape2(i))
         gwake(i) = .TRUE.
-      ENDIF
-      ENDDO
-
-c
-c Check qx and qw positivity
-c --------------------------
-      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))  )
-      ENDDO
-      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).le.q0_min(i)) THEN
-          q0_min(i)=q1_min(i)
-        ENDIF
-      ENDDO
-      ENDDO
-c
-      DO i = 1,klon
-       OK_qx_qw(i) = q0_min(i) .GE. 0.
-       alpha(i) = 1.
-      ENDDO
-c
-CC -----------------------------------------------------------------
-C    Sub-time-stepping
-C    -----------------
-C
-      nsub=10
-      dtimesub=dtime/nsub
-c
-c------------------------------------------------------------
-      DO isubstep = 1,nsub
-c------------------------------------------------------------
-c
-c 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) .GE. 1.
-      ENDDO
-c
-ccc nrlmd   Ajout d'un recalcul de wdens dans le cas d'un entrainement négatif de ktop à kupper --------
-ccc           On calcule pour cela une densité wdens0 pour laquelle on aurait un entrainement nul ---
-      DO i=1,klon
-cc       print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
-cc     $           isubstep,wk_adv(i),cstar(i),wape(i)
-        IF (wk_adv(i) .AND. cstar(i).GT.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) .LE. wdens0*1.1 ) THEN
-            wdens(i) = wdens0
-         ENDIF
-cc	   print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
-cc     $     ,ph(i,1)-pupper(i)',
-cc     $             omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
-cc     $     ,ph(i,1)-pupper(i)
-        ENDIF
-      ENDDO
-
-ccc 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)
-       ENDIF
-      ENDDO
-      DO i=1,klon
-        IF (wk_adv(i)) THEN
-ccc nrlmd          Introduction du taux de mortalité des poches et test sur sigmaw_max=0.4
-ccc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
-           IF (sigmaw(i).ge.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
-c     $              - nat_rate(i)*sigmaw(i)*dtimesub
-cc        print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
-cc     $  death_rate(i),ktop(i),kupper(i)',
-cc     $	         d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
-cc     $  death_rate(i),ktop(i),kupper(i)
-
-c        sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
-c        sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
-c        wdens = wdens0/(10.*sigmaw)
-c        sigmaw =max(sigmaw,sigd_con)
-c        sigmaw =max(sigmaw,sigmad)
-        ENDIF
-      ENDDO
-C
-C
-c calcul de la difference de vitesse verticale poche - zone non perturbee
-cIM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
-cIM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
-cIM 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
-      ENDDO
-      ENDDO
-      DO k= 1,klev+1
-      DO i = 1,klon
-      if (wk_adv(i)) THEN !!! nrlmd
-        omg(i,k)=0.
-      end if
-      ENDDO
-      ENDDO
-c
-      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)))
-        ENDIF
-      ENDDO
-c
-      DO k= 2,klev
-      DO i = 1,klon
-       IF( wk_adv(i) .AND. k .LE. 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)
-       ENDIF
-      ENDDO
-      ENDDO
-c
-      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)
-        ENDIF
-      ENDDO
-c
-c        -----------------
-c        From m/s to Pa/s
-c        -----------------
-c
-       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))
-        ENDIF
-       ENDDO
-c
-      DO k= 1,klev
-      DO i = 1,klon
-       IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN
-          omg(i,k) = - rho(i,k)*rg*omg(i,k)
-          dp_deltomg(i,k) = dp_deltomg(i,1)
-       ENDIF
-      ENDDO
-      ENDDO
-c
-c   raccordement lineaire de omg de ptop a pupper
-
-      DO i=1,klon
-      IF ( wk_adv(i) .AND. kupper(i) .GT. 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))
-      ENDIF
-      ENDDO
-c
-cc      DO i=1,klon
-cc        print*,'Pente entre 0 et kupper (référence)'
-cc     $   	,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
-cc        print*,'Pente entre ktop et kupper'
-cc     $  	,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
-cc      ENDDO
-cc
-      DO k= 1,klev
-      DO i = 1,klon
-       IF( wk_adv(i) .AND. k .GT. ktop(i) .AND. k .LE. 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))
-       ENDIF
-      ENDDO
-      ENDDO
-ccc nrlmd
-cc      DO i=1,klon
-cc      print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
-cc      END DO
-ccc
-c
-c
-c--    Compute wake average vertical velocity omgbw
-c
-c
-      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)
-        ENDIF
-      ENDDO
-      ENDDO
-c--    and its vertical gradient dp_omgbw
-c
-      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))
-        ENDIF
-      ENDDO
-      ENDDO
-C
-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
-      DO i=1,klon
-        IF ( wk_adv(i)) THEN
-         alpha_up(i,k) = 0.
-         IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1.
-        ENDIF
-      ENDDO
-      ENDDO
-
-c  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)
-        ENDIF
-      ENDDO
-      RRd1 = -1.
-      RRd2 = 1.
-c
-c--    Get [Th1,Th2], dth and [q1,q2]
-c
-      DO k= 1,klev
-      DO i = 1,klon
-       IF( wk_adv(i) .AND. k .LE. 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
-       ENDIF
-      ENDDO
-      ENDDO
-
-      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
-      ENDDO
-
-      DO k= 2,klev
-      DO i = 1,klon
-       IF( wk_adv(i) .AND. k .LE. 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)
-       ENDIF
-      ENDDO
-      ENDDO
-
-      DO i=1,klon
-        IF( wk_adv(i)) THEN
-         omgbdth(i,1) = 0.
-         omgbdq(i,1) = 0.
-        ENDIF
-      ENDDO
-
-      DO k= 2,klev
-      DO i = 1,klon
-       IF( wk_adv(i) .AND. k .LE. 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))
-       ENDIF
-      ENDDO
-      ENDDO
-c
-c-----------------------------------------------------------------
-      DO k= 1,klev
-      DO i = 1,klon
-       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
-c-----------------------------------------------------------------
-c
-c   Compute redistribution (advective) term
-c
-         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)
-c         print*,'d_deltatw=',d_deltatw(i,k)
-c
-         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))
-c         print*,'d_deltaqw=',d_deltaqw(i,k)
-c
-c   and increment large scale tendencies
-c
-
-c
-C
-CC -----------------------------------------------------------------
-         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))
-ccc 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))
-ccc
-     $                      )*ppi(i,k)
-c
-         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))
-ccc 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))
-ccc
-     $                      )
-ccc nrlmd
-       ELSE IF(wk_adv(i) .AND. k .EQ. 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)))
-     $                       )
-
-       ENDIF
-ccc
-      ENDDO
-      ENDDO
-c------------------------------------------------------------------
-C
-C   Increment state variables
-
-      DO k= 1,klev
-      DO i = 1,klon
-ccc nrlmd       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
-        IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
-ccc
-
-
-c
-c 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)))
-        
-
-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(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)))
-c        print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
-c
-        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)))
-c        print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k)
-c
-ccc nrlmd          Prise en compte du taux de mortalité
-ccc               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)
-ccc        spread(i,k) = (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
-ccc     $  sigmaw(i)
-
-
-c 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
-
-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(i,k).lt.1.e-10) THEN
-          d_deltatw(i,k) = dtimesub*
-     $       (ff(i)+dtKE(i,k)+dtPBL(i,k)
-ccc     $       -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))
-ccc
-     $       -Tgw(i,k)*deltatw(i,k))
+      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.
+      hw(i) = hwmin
+      sigmaw(i) = sigmad
+      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).
+  ! ------------------------------------------------------
+
+
+  ! 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/
+
+  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
-           d_deltatw(i,k) = 1/Tgw(i,k)*(1-exp(-dtimesub*
-     $       Tgw(i,k)))*
-     $       (ff(i)+dtKE(i,k)+dtPBL(i,k)
-ccc     $       -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))
-ccc
-     $       -Tgw(i,k)*deltatw(i,k))
-        ENDIF
-
-        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)
-ccc     $     -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)))
-ccc
-
-ccc nrlmd
-ccc       d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
-ccc       d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
-ccc
-       ENDIF
-      ENDDO
-      ENDDO
-
-C
-C   Scale tendencies so that water vapour remains positive in w and x.
-C
-      call wake_vec_modulation(klon,klev,wk_adv,epsilon,qe,d_qe,deltaqw,
-     $                d_deltaqw,sigmaw,d_sigmaw,alpha)
-c
-ccc nrlmd
-cc      print*,'alpha'
-cc      do i=1,klon
-cc         print*,alpha(i)
-cc      end do
-ccc
-      DO k = 1,klev
-      DO i = 1,klon
-       IF( wk_adv(i) .AND. k .LE. 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)
-       ENDIF
-      ENDDO
-      ENDDO
-      DO i = 1,klon
-       IF( wk_adv(i)) THEN
-        d_sigmaw(i)=alpha(i)*d_sigmaw(i)
-       ENDIF
-      ENDDO
-
-C   Update large scale variables and wake variables
-cIM 060208 manque DO i + remplace DO k=1,kupper(i)
-cIM 060208     DO k = 1,kupper(i)
-      DO k= 1,klev
-      DO i = 1,klon
-       IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
-        dtls(i,k)=dtls(i,k)+d_te(i,k)
-        dqls(i,k)=dqls(i,k)+d_qe(i,k)
-ccc 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)
-ccc
-       ENDIF
-      ENDDO
-      ENDDO
-      DO k= 1,klev
-      DO i = 1,klon
-       IF( wk_adv(i) .AND. k .LE. 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)
-cc      print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k)
-cc     $        ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
-       ENDIF
-      ENDDO
-      ENDDO
-      DO i = 1,klon
-       IF( wk_adv(i)) THEN
-        sigmaw(i) = sigmaw(i)+d_sigmaw(i)
-       ENDIF
-      ENDDO
-c
-C
-c     Determine Ptop from buoyancy integral
-c     ---------------------------------------
-c
-c-     1/ Pressure of the level where dth changes sign.
-c
-      DO i=1,klon
-       IF ( wk_adv(i)) THEN
-        Ptop_provis(i)=ph(i,1)
-       ENDIF
-      ENDDO
-c
-      DO k= 2,klev
-      DO i=1,klon
-        IF ( wk_adv(i) .AND.
-     $       Ptop_provis(i) .EQ. ph(i,1) .AND.
-     $      dth(i,k) .GT. -delta_t_min .and.
-     $      dth(i,k-1).LT. -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))
-        ENDIF
-      ENDDO
-      ENDDO
-c
-c-     2/ dth integral
-c
-      DO i=1,klon
-       if (wk_adv(i)) then !!! nrlmd
-       sum_dth(i) = 0.
-       dthmin(i) = -delta_t_min
-       z(i) = 0.
-       end if
-      ENDDO
-
-      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) .gt. 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))
-        ENDIF
-       ENDIF
-      ENDDO
-      ENDDO
-c
-c-     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))
-       ENDIF
-      ENDDO
-c
-c-     4/ now, get Ptop
-c
-      DO i=1,klon
-       if (wk_adv(i)) then !!! nrlmd
-       ktop(i) = 0
-       z(i)=0.
-       end if
-      ENDDO
-c
-      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) .gt. 0) THEN
-         z(i) = z(i)+dz(i)
-         Ptop(i) = Ph(i,k)-rho(i,k)*rg*dz(i)
-         ktop(i) = k
-        ENDIF
-       ENDIF
-      ENDDO
-      ENDDO
-c
-c      4.5/Correct ktop and ptop
-c
-      DO i=1,klon
-       IF ( wk_adv(i)) THEN
-        Ptop_new(i)=ptop(i)
-       ENDIF
-      ENDDO
-c
-      DO k= klev,2,-1
-      DO i=1,klon
-cIM v3JYG; IF (k .GE. ktop(i)
-       IF ( wk_adv(i) .AND.
-     $      k .LE. ktop(i) .AND.
-     $      ptop_new(i) .EQ. ptop(i) .AND.
-     $      dth(i,k) .GT. -delta_t_min .and.
-     $      dth(i,k-1).LT. -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))
-        ENDIF
-      ENDDO
-      ENDDO
-c
-c
-      DO i=1,klon
-       IF ( wk_adv(i)) THEN
-        ptop(i) = ptop_new(i)
-       ENDIF
-      ENDDO
-
-      DO k=klev,1,-1
-      DO i=1,klon
-      if (wk_adv(i)) then !!! nrlmd
-        IF (ph(i,k+1) .LT. ptop(i)) ktop(i)=k
-      end if
-      ENDDO
-      ENDDO
-c
-c      5/ Set deltatw & deltaqw to 0 above kupper
-c
-      DO k = 1,klev
-      DO i=1,klon
-        IF ( wk_adv(i) .AND. k .GE. kupper(i)) THEN
-         deltatw(i,k) = 0.
-         deltaqw(i,k) = 0.
-        ENDIF
-      ENDDO
-      ENDDO
-c
-C
-c-------------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
-      ENDDO
-C
-C Integrals (and wake top level number)
-C --------------------------------------
-C
-C 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
-      ENDDO
-
-      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) .GT. 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)
-        ENDIF
-       end if
-      ENDDO
-      ENDDO
-c
-      DO i=1,klon
-       if (wk_adv(i)) then !!! nrlmd
-        hw0(i) = z(i)
-       end if
-      ENDDO
-c
-C
-C - WAPE and mean forcing computation
-C ---------------------------------------
-C
-C ---------------------------------------
-C
-C 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)
-c
-       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
-      ENDDO
-C
-C Filter out bad wakes
-
-      DO k = 1,klev
-       DO i=1,klon
-        if (wk_adv(i)) then !!! nrlmd
-        IF ( wape(i) .LT. 0.) THEN
-          deltatw(i,k) = 0.
-          deltaqw(i,k) = 0.
-          dth(i,k) = 0.
-        ENDIF
-	end if
-       ENDDO
-      ENDDO
-c
-      DO i=1,klon
-      if (wk_adv(i)) then !!! nrlmd
-      IF ( wape(i) .LT. 0.) THEN
-        wape(i) = 0.
-        Cstar(i) = 0.
-        hw(i) = hwmin
-        sigmaw(i) = max(sigmad,sigd_con(i))
-        fip(i) = 0.
-        gwake(i) = .FALSE.
+          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
+
+SUBROUTINE wake_scal(p, ph, ppi, 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
+  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
+  ! 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)
+  ! ppi  : (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  : number of wakes per unit area (3D) or per
+  ! unit length (2D)
+  ! 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 "dimensions.h"
+  ! ccc      include "dimphy.h"
+  include "YOMCST.h"
+  include "cvthermo.h"
+  include "iniprint.h"
+
+  ! Arguments en entree
+  ! --------------------
+
+  REAL p(klev), ph(klev+1), ppi(klev)
+  REAL dtime
+  REAL te0(klev), qe0(klev)
+  REAL omgb(klev+1)
+  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
+
+  ! Sorties
+  ! --------
+
+  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
+
+  ! Variables internes
+  ! -------------------
+
+  ! 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
+
+  ! Variables de sauvegarde
+  REAL deltatw0(klev)
+  REAL deltaqw0(klev)
+  REAL te(klev), qe(klev)
+  REAL sigmaw0, sigmaw1
+
+  ! Variables pour les GW
+  REAL ll
+  REAL n2(klev)
+  REAL cgw(klev)
+  REAL tgw(klev)
+
+  ! 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
+
+  ! 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 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
+
+  ! -------------------------------------------------------------------------
+  ! Initialisations
+  ! -------------------------------------------------------------------------
+
+  ! print*, 'wake initialisations'
+
+  ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
+  ! -------------------------------------------------------------------------
+
+  DATA sigmad, hwmin/.02, 10./
+
+  ! 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
+  ! -------------------------------------------------------------------------
+
+  coefgw = 10
+  ! coefgw=1
+  ! wdens0 = 1.0/(alon**2)
+  wdens = 1.0/(alon**2)
+  stark = 0.50
+  ! CRtest
+  alpk = 0.1
+  ! alpk = 1.0
+  ! alpk = 0.5
+  ! alpk = 0.05
+  crep_upper = 0.9
+  crep_sol = 1.0
+
+
+  ! 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
+    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.
+  END DO
+  ! sigmaw1=sigmaw
+  ! IF (sigd_con.GT.sigmaw1) THEN
+  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
+  ! ENDIF
+  sigmaw = max(sigmaw, sigd_con)
+  sigmaw = max(sigmaw, sigmad)
+  sigmaw = min(sigmaw, 0.99)
+  sigmaw0 = sigmaw
+  ! wdens=wdens0/(10.*sigmaw)
+  ! IF (sigd_con.GT.sigmaw1) THEN
+  ! print*, 'sigmaw1,sigd1', sigmaw, sigd_con
+  ! ENDIF
+
+  ! 2. - Prognostic part
+  ! =========================================================
+
+  ! print *, 'prognostic wake computation'
+
+
+  ! 2.1 - Undisturbed area and Wake integrals
+  ! ---------------------------------------------------------
+
+  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.
+
+  ! Potential temperatures and humidity
+  ! ----------------------------------------------------------
+
+  DO k = 1, klev
+    rho(k) = p(k)/(rd*te(k))
+    IF (k==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)
+    END IF
+    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)
+  END DO
+
+  DO k = 1, klev - 1
+    IF (k==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)))
+    END IF
+    zh(k) = (zhh(k)+zhh(k+1))/2
+
+    cgw(k) = sqrt(n2(k))*zh(k)
+    tgw(k) = coefgw*cgw(k)/ll
+  END DO
+
+  n2(klev) = 0
+  zh(klev) = 0
+  cgw(klev) = 0
+  tgw(klev) = 0
+
+  ! Calcul de la masse volumique moyenne de la colonne
+  ! -----------------------------------------------------------------
+
+  DO k = 1, klev
+    epaisseur1(k) = 0.
+    epaisseur2(k) = 0.
+  END DO
+
+  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)
+  END DO
+
+
+  ! Choose an integration bound well above wake top
+  ! -----------------------------------------------------------------
+
+  ! Pupper = 50000.  ! melting level
+  pupper = 60000.
+  ! Pupper = 70000.
+
+
+  ! Determine Wake top pressure (Ptop) from buoyancy integral
+  ! -----------------------------------------------------------------
+
+  ! -1/ Pressure of the level where dth becomes less than delta_t_min.
+
+  ptop_provis = ph(1)
+  DO k = 2, klev
+    IF (dth(k)>-delta_t_min .AND. dth(k-1)<-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
+    END IF
+  END DO
+25 CONTINUE
+
+  ! -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<=0) GO TO 40
+    z = z + dz
+    sum_dth = sum_dth + dth(k)*dz
+    dthmin = min(dthmin, dth(k))
+  END DO
+40 CONTINUE
+
+  ! -3/ height of triangle with area= sum_dth and base = dthmin
+
+  hw0 = 2.*sum_dth/min(dthmin, -0.5)
+  hw0 = max(hwmin, hw0)
+
+  ! -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<=0) GO TO 45
+    z = z + dz
+    ptop = ph(k) - rho(k)*rg*dz
+  END DO
+45 CONTINUE
+
+
+  ! -5/ Determination de ktop et kupper
+
+  DO k = klev, 1, -1
+    IF (ph(k+1)<ptop) ktop = k
+    IF (ph(k+1)<pupper) kupper = k
+  END DO
+
+  ! -6/ Correct ktop and ptop
+
+  ptop_new = ptop
+  DO k = ktop, 2, -1
+    IF (dth(k)>-delta_t_min .AND. dth(k-1)<-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
+    END IF
+  END DO
+225 CONTINUE
+
+  ptop = ptop_new
+
+  DO k = klev, 1, -1
+    IF (ph(k+1)<ptop) ktop = k
+  END DO
+
+  ! Set deltatw & deltaqw to 0 above kupper
+  ! -----------------------------------------------------------
+
+  DO k = kupper, klev
+    deltatw(k) = 0.
+    deltaqw(k) = 0.
+  END DO
+
+
+  ! Vertical gradient of LS omega
+  ! ------------------------------------------------------------
+
+  DO k = 1, kupper
+    dp_omgb(k) = (omgb(k+1)-omgb(k))/(ph(k+1)-ph(k))
+  END DO
+
+
+  ! Integrals (and wake top level number)
+  ! -----------------------------------------------------------
+
+  ! 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<=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
+  END DO
+50 CONTINUE
+
+  hw0 = z
+
+  ! 2.1 - WAPE and mean forcing computation
+  ! -------------------------------------------------------------
+
+  ! Means
+
+  av_thu = sum_thu/hw0
+  av_tu = sum_tu/hw0
+  av_qu = sum_qu/hw0
+  av_thvu = sum_thvu/hw0
+  ! 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
+
+  ! 2.2 Prognostic variable update
+  ! ------------------------------------------------------------
+
+  ! Filter out bad wakes
+
+  IF (wape<0.) THEN
+    IF (prt_level>=10) 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.
+    END DO
+  ELSE
+    IF (prt_level>=10) PRINT *, 'wape>0'
+    cstar = stark*sqrt(2.*wape)
+  END IF
+
+  ! ------------------------------------------------------------------
+  ! Sub-time-stepping
+  ! ------------------------------------------------------------------
+
+  ! nsub=36
+  nsub = 10
+  dtimesub = dtime/nsub
+
+  ! ------------------------------------------------------------
+  DO isubstep = 1, nsub
+    ! ------------------------------------------------------------
+
+    ! print*,'---------------','substep=',isubstep,'-------------'
+
+    ! Evolution of sigmaw
+
+
+    gfl = 2.*sqrt(3.14*wdens*sigmaw)
+
+    sigmaw = sigmaw + gfl*cstar*dtimesub
+    sigmaw = min(sigmaw, 0.99) !!!!!!!!
+    ! wdens = wdens0/(10.*sigmaw)
+    ! sigmaw =max(sigmaw,sigd_con)
+    ! sigmaw =max(sigmaw,sigmad)
+
+    ! 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
+    END DO
+
+    dztop = -(ptop-ph(ktop))/(rho(ktop)*rg)
+    ztop = z + dztop
+    omgtop = dp_deltomg(1)*ztop
+
+
+    ! 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)
+    END DO
+
+    ! raccordement lineaire de omg de ptop a pupper
+
+    IF (kupper>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)
+      END DO
+    END IF
+
+    ! Compute wake average vertical velocity omgbw
+
+    DO k = 1, klev + 1
+      omgbw(k) = omgb(k) + (1.-sigmaw)*omg(k)
+    END DO
+
+    ! and its vertical gradient dp_omgbw
+
+    DO k = 1, klev
+      dp_omgbw(k) = (omgbw(k+1)-omgbw(k))/(ph(k+1)-ph(k))
+    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
+      alpha_up(k) = 0.
+      IF (omgb(k)>0.) alpha_up(k) = 1.
+    END DO
+
+    ! Matrix expressing [The,deltatw] from  [Th1,Th2]
+
+    rre1 = 1. - sigmaw
+    rre2 = sigmaw
+    rrd1 = -1.
+    rrd2 = 1.
+
+    ! 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
+    END DO
+
+    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)
+    END DO
+
+    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))
+    END DO
+
+
+    ! -----------------------------------------------------------------
+    DO k = 1, kupper - 1
+      ! -----------------------------------------------------------------
+
+      ! Compute redistribution (advective) term
+
+      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)
+      ! print*,'d_deltatw=',d_deltatw(k)
+
+      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))
+      ! print*,'d_deltaqw=',d_deltaqw(k)
+
+      ! and increment large scale tendencies
+
+      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)
+      ! print*,'dtls=',dtls(k)
+
+      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))
+      ! print*,'dqls=',dqls(k)
+
+      ! -------------------------------------------------------------------
+    END DO
+    ! ------------------------------------------------------------------
+
+    ! Increment state variables
+
+    DO k = 1, kupper - 1
+
+      ! 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))
+
+
+      ! 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(k) = (dtdwn(k)/sigmaw-dta(k)/(1.-sigmaw))
+      dqke(k) = (dqdwn(k)/sigmaw-dqa(k)/(1.-sigmaw))
+      ! print*,'dtKE=',dtKE(k)
+      ! print*,'dqKE=',dqKE(k)
+
+      dtpbl(k) = (wdtpbl(k)/sigmaw-udtpbl(k)/(1.-sigmaw))
+      dqpbl(k) = (wdqpbl(k)/sigmaw-udqpbl(k)/(1.-sigmaw))
+
+      spread(k) = (1.-sigmaw)*dp_deltomg(k) + gfl*cstar/sigmaw
+      ! print*,'spread=',spread(k)
+
+
+      ! 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
+      ! print*,'d_delta_gw=',d_deltat_gw(k)
+      ff = d_deltatw(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(k)<1.E-10) THEN
+        deltatw(k) = deltatw(k) + dtimesub*(ff+dtke(k)+dtpbl(k)-spread(k)* &
+          deltatw(k)-tgw(k)*deltatw(k))
       ELSE
-        Cstar(i) = stark*sqrt(2.*wape(i))
-        gwake(i) = .TRUE.
-      ENDIF
-      end if
-      ENDDO
-
-       ENDDO      ! end sub-timestep loop
-C
-C -----------------------------------------------------------------
-c   Get back to tendencies per second
-c
-      DO k = 1,klev
-      DO i=1,klon
-
-ccc nrlmd        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
-        IF ( OK_qx_qw(i) .AND. k .LE. kupper(i)) THEN
-ccc
-         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
-cc      print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
-cc     $         ,death_rate(i)*sigmaw(i)
-        ENDIF
-      ENDDO
-      ENDDO
-
-c
-c----------------------------------------------------------
-c   Determine wake final state; recompute wape, cstar, ktop;
-c   filter out bad wakes.
-c----------------------------------------------------------
-c
-C 2.1 - Undisturbed area and Wake integrals
-C ---------------------------------------------------------
-
-      DO i=1,klon
-ccc nrlmd       if (wk_adv(i)) then !!! nrlmd
-      if (OK_qx_qw(i)) then
-ccc
-        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	
-      ENDDO
-C Potential temperatures and humidity
-c----------------------------------------------------------
-
-      DO k =1,klev
-      DO i=1,klon
-ccc nrlmd       IF ( wk_adv(i)) THEN
-       if (OK_qx_qw(i)) then
-ccc
-        rho(i,k) = p(i,k)/(rd*te(i,k))
-        IF(k .eq. 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)
-        ENDIF
-        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)
-       ENDIF
-      ENDDO
-      ENDDO
-
-C Integrals (and wake top level number)
-C -----------------------------------------------------------
-
-C Initialize sum_thvu to 1st level virt. pot. temp.
-
-      DO i=1,klon
-ccc nrlmd       IF ( wk_adv(i)) THEN
-      if (OK_qx_qw(i)) then
-ccc
-        z(i) = 1.
-        dz(i) = 1.
-        sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
-        sum_dth(i) = 0.
-      ENDIF
-      ENDDO
-
-      DO k = 1,klev
-      DO i=1,klon
-ccc nrlmd       IF ( wk_adv(i)) THEN
-       if (OK_qx_qw(i)) then
-ccc
-        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
-        IF (dz(i) .GT. 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)
-        ENDIF
-       ENDIF
-      ENDDO
-      ENDDO
-c
-      DO i=1,klon
-ccc nrlmd       IF ( wk_adv(i)) THEN
-       if (OK_qx_qw(i)) then
-ccc
-        hw0(i) = z(i)
-       ENDIF
-      ENDDO
-c
-C - WAPE and mean forcing computation
-C-------------------------------------------------------------
-
-C Means
-
-      DO i=1, klon
-ccc nrlmd       IF ( wk_adv(i)) THEN
-      if (OK_qx_qw(i)) then
-ccc
-        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)
-       ENDIF
-      ENDDO
-
-C Prognostic variable update
-C ------------------------------------------------------------
-
-C Filter out bad wakes
-c
-      DO k = 1,klev
-      DO i=1,klon
-ccc nrlmd        IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
-      if (OK_qx_qw(i) .AND. wape2(i) .LT. 0.) then
-ccc
-          deltatw(i,k) = 0.
-          deltaqw(i,k) = 0.
-          dth(i,k) = 0.
-        ENDIF
-      ENDDO
-      ENDDO
-c
-
-      DO i=1, klon
-ccc nrlmd       IF ( wk_adv(i)) THEN
-      if (OK_qx_qw(i)) then
-ccc
-       IF ( wape2(i) .LT. 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.ge.10) print*,'wape2>0'
-        Cstar2(i) = stark*sqrt(2.*wape2(i))
-        gwake(i) = .TRUE.
-      ENDIF
-      ENDIF
-      ENDDO
-c
-      DO i=1, klon
-ccc nrlmd       IF ( wk_adv(i)) THEN
-       if (OK_qx_qw(i)) then
-ccc
-        ktopw(i) = ktop(i)
-       ENDIF
-      ENDDO
-c
-      DO i=1, klon
-ccc nrlmd       IF ( wk_adv(i)) THEN
-       if (OK_qx_qw(i)) then
-ccc
-       IF (ktopw(i) .gt. 0 .and. gwake(i)) 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(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)
-Cjyg2
-       ELSE
-         FIP(i) = 0.
-       ENDIF
-       ENDIF
-      ENDDO
-c
-C   Limitation de sigmaw
-
-ccc nrlmd
-c       DO i=1,klon
-c         IF (OK_qx_qw(i)) THEN
-c	   IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
-c	 ENDIF
-c       ENDDO
-ccc
-      DO k = 1,klev
-       DO i=1, klon
-
-ccc nrlmd      On maintient désormais constant sigmaw en régime permanent
-ccc      IF ((sigmaw(i).GT.sigmaw_max).or.
-        IF     ( ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.
-     $         (ktopw(i).le.2) .OR.
-     $         .not. OK_qx_qw(i) ) THEN
-ccc
-          dtls(i,k) = 0.
-          dqls(i,k) = 0.
-          deltatw(i,k) = 0.
-          deltaqw(i,k) = 0.
-        ENDIF
-       ENDDO
-      ENDDO
-c
-ccc nrlmd      On maintient désormais constant sigmaw en régime permanent
-      DO i=1, klon
-        IF  ( ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.
-     $        (ktopw(i).le.2) .OR.
-     $        .not. OK_qx_qw(i)   ) THEN
-         wape(i) = 0.
-         cstar(i)=0.
-         hw(i) = hwmin
-         sigmaw(i) = sigmad
-         fip(i) = 0.
-        ELSE
-         wape(i) = wape2(i)
-         cstar(i)=cstar2(i)
-        ENDIF
-cc        print*,'wape wape2 ktopw OK_qx_qw =',
-cc     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
-      ENDDO
-c
-c
-      RETURN
-      END
-
-      SUBROUTINE wake_vec_modulation(nlon,nl,wk_adv,epsilon,qe,d_qe,
-     $           deltaqw,d_deltaqw,sigmaw,d_sigmaw,alpha)
-c------------------------------------------------------
-cDtermination du coefficient alpha tel que les tendances
-c corriges alpha*d_G, pour toutes les grandeurs G, correspondent
-c a une humidite positive dans la zone (x) et dans la zone (w).
-c------------------------------------------------------
-c
- 
-c  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
-c  Output
-      REAL alpha(nlon)
-c  Internal variables
-      REAL zeta(nlon,nl)
-      REAL alpha1(nlon)
-      REAL x,a,b,c,discrim
-      REAL epsilon
-!      DATA epsilon/1.e-15/
-c
-      DO k=1,nl
-      DO i = 1,nlon
-       IF (wk_adv(i)) THEN
-        IF ((deltaqw(i,k)+d_deltaqw(i,k)).ge.0.) then
-         zeta(i,k)=0.
-        ELSE
-         zeta(i,k)=1.
-        END IF
-       ENDIF
-      ENDDO
-      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
-c      print*, 'x, a, b, c, discrim', x, a, b, c, discrim
-        IF (a+b .GE. 0.) THEN !! Condition suffisante pour la positivité de ovap
-         alpha1(i)=1.
-        ELSE
-         IF (x .GE. 0.) THEN
-            alpha1(i)=1.
-         ELSE
-              IF (a .GT. 0.) THEN
-                 alpha1(i)=0.9*min(   (2.*c)/(-b+sqrt(discrim)),
-     $                        (-b+sqrt(discrim))/(2.*a)   )
-              ELSE IF (a .eq. 0.) then
-                 alpha1(i)=0.9*(-c/b)
-              ELSE
-c         print*,'a,b,c discrim',a,b,c discrim
-                 alpha1(i)=0.9*max(   (2.*c)/(-b+sqrt(discrim)),
-     $                        (-b+sqrt(discrim))/(2.*a)   )
-              ENDIF
-         ENDIF
-        ENDIF
-       alpha(i) = min(alpha(i),alpha1(i))
-       ENDIF
-      ENDDO
-      ENDDO
-!
-      return
-      end
-
-      Subroutine WAKE_scal (p,ph,ppi,dtime,sigd_con
-     :                ,te0,qe0,omgb
-     :                ,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
-      USE dimphy
-      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                        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"
-cccc      include "dimphy.h"
-      include "YOMCST.h"
-      include "cvthermo.h"
-      include "iniprint.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)
-      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 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 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
-        if(prt_level.ge.10) 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
-        if(prt_level.ge.10) 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))
+        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))
+      END IF
+
+      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)
+    END DO
+
+    ! 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)
+    END DO
+
+    ! Determine Ptop from buoyancy integral
+    ! ----------------------------------------------------------------------
+
+    ! -1/ Pressure of the level where dth changes sign.
+
+    ptop_provis = ph(1)
+
+    DO k = 2, klev
+      IF (dth(k)>-delta_t_min .AND. dth(k-1)<-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
-        if(prt_level.ge.10) 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
-        if(prt_level.ge.10) 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
-
-! correction NICOLAS     .     ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
-      IF ((sigmaw.GT.0.9).or.
-     .     ((wape.ge.wape2).and.(wape2.le.1.0)).or.(ktopw.le.2)) THEN
-cIM cf NR/JYG 251108    .     ((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
-
-
-
+      END IF
+    END DO
+65  CONTINUE
+
+    ! -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<=0) GO TO 70
+      z = z + dz
+      sum_dth = sum_dth + dth(k)*dz
+      dthmin = min(dthmin, dth(k))
+    END DO
+70  CONTINUE
+
+    ! -3/ height of triangle with area= sum_dth and base = dthmin
+
+    hw = 2.*sum_dth/min(dthmin, -0.5)
+    hw = max(hwmin, hw)
+
+    ! -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<=0) GO TO 75
+      z = z + dz
+      ptop = ph(k) - rho(k)*rg*dz
+      ktop = k
+    END DO
+75  CONTINUE
+
+    ! -5/Correct ktop and ptop
+
+    ptop_new = ptop
+
+    DO k = ktop, 2, -1
+      IF (dth(k)>-delta_t_min .AND. dth(k-1)<-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
+      END IF
+    END DO
+275 CONTINUE
+
+    ptop = ptop_new
+
+    DO k = klev, 1, -1
+      IF (ph(k+1)<ptop) ktop = k
+    END DO
+
+    ! -6/ Set deltatw & deltaqw to 0 above kupper
+
+    DO k = kupper, klev
+      deltatw(k) = 0.
+      deltaqw(k) = 0.
+    END DO
+
+    ! ------------------------------------------------------------------
+  END DO ! end sub-timestep loop
+  ! -----------------------------------------------------------------
+
+  ! 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
+  END DO
+
+  ! 2.1 - Undisturbed area and Wake integrals
+  ! ---------------------------------------------------------
+
+  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.
+
+  ! Potential temperatures and humidity
+  ! ----------------------------------------------------------
+
+  DO k = 1, klev
+    rho(k) = p(k)/(rd*te(k))
+    IF (k==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)
+    END IF
+    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)
+
+  END DO
+
+  ! Integrals (and wake top level number)
+  ! -----------------------------------------------------------
+
+  ! 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<=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
+  END DO
+51 CONTINUE
+
+  hw0 = z
+
+  ! 2.1 - WAPE and mean forcing computation
+  ! -------------------------------------------------------------
+
+  ! 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
+
+
+  ! 2.2 Prognostic variable update
+  ! ------------------------------------------------------------
+
+  ! Filter out bad wakes
+
+  IF (wape2<0.) THEN
+    IF (prt_level>=10) 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.
+    END DO
+  ELSE
+    IF (prt_level>=10) PRINT *, 'wape2>0'
+    cstar2 = stark*sqrt(2.*wape2)
+
+  END IF
+
+  ktopw = ktop
+
+  IF (ktopw>0) THEN
+
+    ! jyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
+    ! cc       heff = 600.
+    ! Utilisation de la hauteur hw
+    ! c       heff = 0.7*hw
+    heff = hw
+
+    fip = 0.5*rho(ktopw)*cstar2**3*heff*2*sqrt(sigmaw*wdens*3.14)
+    fip = alpk*fip
+    ! jyg2
+  ELSE
+    fip = 0.
+  END IF
+
+
+  ! Limitation de sigmaw
+
+  ! sécurité : si le wake occuppe plus de 90 % de la surface de la maille,
+  ! alors il disparait en se mélangeant à la partie undisturbed
+
+  ! correction NICOLAS     .     ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
+  IF ((sigmaw>0.9) .OR. ((wape>=wape2) .AND. (wape2<= &
+      1.0)) .OR. (ktopw<=2)) THEN
+    ! IM cf NR/JYG 251108    .     ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
+    ! IF (sigmaw.GT.0.9) THEN
+    DO k = 1, klev
+      dtls(k) = 0.
+      dqls(k) = 0.
+      deltatw(k) = 0.
+      deltaqw(k) = 0.
+    END DO
+    wape = 0.
+    hw = hwmin
+    sigmaw = sigmad
+    fip = 0.
+  END IF
+
+  RETURN
+END SUBROUTINE wake_scal
+
+
+
