Index: LMDZ6/trunk/libf/phylmdiso/calwake.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/calwake.F90	(revision 4581)
+++ LMDZ6/trunk/libf/phylmdiso/calwake.F90	(revision 4594)
@@ -39,4 +39,6 @@
   USE isotopes_mod, ONLY: iso_eau
   USE isotopes_verif_mod
+  USE lmdz_wake, ONLY : wake
+
 #endif
 #endif
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_alim.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_alim.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_alim.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_alim.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_alp.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_alp.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_alp.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_alp.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_closure.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_closure.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_closure.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_closure.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_down.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_down.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_down.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_down.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_dq.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_dq.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_dq.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_dq.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_dry.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_dry.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_dry.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_dry.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_dtke.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_dtke.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_dtke.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_dtke.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_dv2.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_dv2.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_dv2.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_dv2.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_env.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_env.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_env.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_env.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_flux2.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_flux2.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_flux2.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_flux2.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_height.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_height.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_height.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_height.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_ini.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_ini.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_ini.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_ini.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_main.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_main.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_main.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_main.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_old.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_old.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_old.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_old.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_plume.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_plume.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_plume.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_plume.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_plume_6A.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_plume_6A.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_plume_6A.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_plume_6A.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_qsat.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_qsat.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_thermcell_qsat.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_thermcell_qsat.F90
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_wake.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_wake.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_wake.F90	(revision 4594)
@@ -0,0 +1,2629 @@
+MODULE lmdz_wake
+! $Id: wake.F90 3648 2020-03-16 15:36:59Z jghattas $
+CONTAINS
+SUBROUTINE wake(znatsurf, p, ph, pi, dtime, &
+                te0, qe0, omgb, &
+                dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
+                sigd_con, Cin, &
+                deltatw, deltaqw, sigmaw, awdens, wdens, &                  ! state variables
+                dth, hw, wape, fip, gfl, &
+                dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, &
+                dtke, dqke, omg, dp_deltomg, spread, cstar, &
+                d_deltat_gw, &
+                d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2 &     ! tendencies
+
+#ifdef ISO
+                     ,xte0,dxtdwn,dxta,deltaxtw &
+                     ,dxtls,xtu,dxtke,d_deltaxtw2 &
+#endif                  
+                )  
+
+
+  ! **************************************************************
+  ! *
+  ! WAKE                                                        *
+  ! retour a un Pupper fixe                                *
+  ! *
+  ! written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
+  ! modified by :   ROEHRIG Romain        01/29/2007            *
+  ! **************************************************************
+
+  USE ioipsl_getin_p_mod, ONLY : getin_p
+  USE dimphy
+  use mod_phys_lmdz_para
+  USE print_control_mod, ONLY: prt_level
+#ifdef ISO
+  USE infotrac_phy, ONLY : ntraciso=>ntiso
+#ifdef ISOVERIF
+  USE isotopes_verif_mod
+!, ONLY: errmax,errmaxrel
+  USE isotopes_mod, ONLY: iso_eau,iso_hdo,ridicule
+#endif
+#endif
+  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.
+
+  ! State variables : 
+  ! deltatw    : temperature difference between wake and off-wake regions
+  ! deltaqw    : specific humidity difference between wake and off-wake regions
+  ! sigmaw     : fractional area covered by wakes.
+  ! wdens      : number of wakes per unit area
+
+  ! 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   : wake top hight (given by hw*deltatw(1)/2=wape)
+  ! dp_omgb : vertical gradient of large scale omega
+  ! awdens  : densite de poches actives
+  ! 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 d_t_wake and d_q_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)
+  ! wgen : number of wakes generated per unit area and per sec (/m^2/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)
+  ! sigd_con: 
+  ! Cin  : convective inhibition
+  ! 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    : wake top hight (defined as the altitude at which deltatw=0)
+  ! amflux : horizontal mass flux through wake boundary
+  ! wdens_ref: initial number of wakes per unit area (3D) or per
+  ! unit length (2D), at the beginning of each time step
+  ! Tgw    : 1 sur la période de onde de gravité
+  ! Cgw    : vitesse de propagation de onde de gravité
+  ! LL     : distance entre 2 poches
+
+  ! -------------------------------------------------------------------------
+  ! Déclaration de variables
+  ! -------------------------------------------------------------------------
+
+  include "YOMCST.h"
+  include "cvthermo.h"
+
+  ! Arguments en entree
+  ! --------------------
+
+  INTEGER, DIMENSION (klon),        INTENT(IN)          :: znatsurf
+  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: p, pi
+  REAL, DIMENSION (klon, klev+1),   INTENT(IN)          :: ph
+  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: omgb
+  REAL,                             INTENT(IN)          :: dtime
+  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: te0, qe0
+  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dtdwn, dqdwn
+  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: amdwn, amup
+  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dta, dqa
+  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen
+  REAL, DIMENSION (klon),           INTENT(IN)          :: sigd_con
+  REAL, DIMENSION (klon),           INTENT(IN)          :: Cin
+#ifdef ISO
+  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(IN)          :: xte0
+  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(IN)          :: dxtdwn
+  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(IN)          :: dxta
+#endif
+
+  !
+  ! Input/Output
+  ! State variables
+  REAL, DIMENSION (klon, klev),     INTENT(INOUT)       :: deltatw, deltaqw
+  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw
+  REAL, DIMENSION (klon),           INTENT(INOUT)       :: awdens
+  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens
+#ifdef ISO
+  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(INOUT)       :: deltaxtw
+#endif
+
+  ! Sorties
+  ! --------
+
+  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dth
+  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: tu, qu
+  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtls, dqls
+  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke
+  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: spread    !  unused (jyg)
+  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: omgbdth, omg
+  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dp_omgb, dp_deltomg
+  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_gw
+  REAL, DIMENSION (klon),           INTENT(OUT)         :: hw, wape, fip, gfl, cstar
+  INTEGER, DIMENSION (klon),        INTENT(OUT)         :: ktopw
+  ! Tendencies of state variables
+  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
+  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_awdens2, d_wdens2
+#ifdef ISO
+  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(OUT)         :: xtu
+  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(OUT)         :: dxtls
+  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(OUT)         :: dxtke
+  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(OUT)         :: d_deltaxtw2
+#endif
+
+  ! Variables internes
+  ! -------------------
+
+  ! Variables à fixer
+  INTEGER, SAVE                                         :: igout
+  !$OMP THREADPRIVATE(igout)
+  LOGICAL, SAVE                                         :: first = .TRUE.
+  !$OMP THREADPRIVATE(first)
+!jyg<
+!!  REAL, SAVE                                            :: stark, wdens_ref, coefgw, alpk
+  REAL, SAVE, DIMENSION(2)                              :: wdens_ref
+  REAL, SAVE                                            :: stark, coefgw, alpk
+!>jyg
+  REAL, SAVE                                            :: crep_upper, crep_sol  
+  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
+
+  REAL, SAVE                                            :: tau_cv
+  !$OMP THREADPRIVATE(tau_cv)
+  REAL, SAVE                                            :: rzero, aa0 ! minimal wake radius and area
+  !$OMP THREADPRIVATE(rzero, aa0)
+
+  LOGICAL, SAVE                                         :: flag_wk_check_trgl
+  !$OMP THREADPRIVATE(flag_wk_check_trgl)
+  INTEGER, SAVE                                         :: iflag_wk_check_trgl
+  !$OMP THREADPRIVATE(iflag_wk_check_trgl)
+  INTEGER, SAVE                                         :: iflag_wk_pop_dyn
+  !$OMP THREADPRIVATE(iflag_wk_pop_dyn)
+
+  REAL                                                  :: delta_t_min
+  INTEGER                                               :: nsub
+  REAL                                                  :: dtimesub
+  REAL, SAVE                                            :: wdensmin
+  !$OMP THREADPRIVATE(wdensmin)
+  REAL, SAVE                                            :: sigmad, hwmin, wapecut, cstart
+  !$OMP THREADPRIVATE(sigmad, hwmin, wapecut, cstart)
+  REAL, SAVE                                            :: sigmaw_max
+  !$OMP THREADPRIVATE(sigmaw_max)  
+  REAL, SAVE                                            :: dens_rate
+  !$OMP THREADPRIVATE(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)                                :: sigmaw1
+#ifdef ISO
+  REAL, DIMENSION (ntraciso,klon, klev)                          :: deltaxtw0
+  REAL, DIMENSION (ntraciso,klon, klev)                          :: xte
+#endif
+
+  ! Variables liees a la dynamique de population
+  REAL, DIMENSION(klon)                                 :: act
+  REAL, DIMENSION(klon)                                 :: rad_wk, tau_wk_inv
+  REAL, DIMENSION(klon)                                 :: f_shear
+  REAL, DIMENSION(klon)                                 :: drdt
+  REAL, DIMENSION(klon)                                 :: d_sig_gen, d_sig_death, d_sig_col
+  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
+  LOGICAL, DIMENSION (klon)                             :: kill_wake
+  INTEGER, SAVE                                         :: iflag_wk_act
+  !$OMP THREADPRIVATE(iflag_wk_act)
+  REAL                                                  :: drdt_pos
+  REAL                                                  :: tau_wk_inv_min
+
+  ! Variables pour les GW
+  REAL, DIMENSION (klon)                                :: ll
+  REAL, DIMENSION (klon, klev)                          :: n2
+  REAL, DIMENSION (klon, klev)                          :: cgw
+  REAL, DIMENSION (klon, klev)                          :: tgw
+
+  ! Variables liees 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
+
+  ! Variables liees au test de la forme triangulaire du profil de Delta_theta
+  REAL, DIMENSION (klon)                                :: sum_half_dth
+  REAL, DIMENSION (klon)                                :: dz_half
+
+  ! Sub-timestep tendencies and related variables
+  REAL, DIMENSION (klon, klev)                          :: d_deltatw, d_deltaqw
+  REAL, DIMENSION (klon, klev)                          :: d_te, d_qe
+  REAL, DIMENSION (klon)                                :: d_awdens, d_wdens
+  REAL, DIMENSION (klon)                                :: d_sigmaw, alpha
+  REAL, DIMENSION (klon)                                :: q0_min, q1_min
+  LOGICAL, DIMENSION (klon)                             :: wk_adv, ok_qx_qw
+#ifdef ISO
+  REAL, DIMENSION (ntraciso,klon, klev)                          :: d_deltaxtw
+  REAL, DIMENSION (ntraciso,klon, klev)                          :: d_xte
+#endif
+  REAL, SAVE                                            :: epsilon
+  !$OMP THREADPRIVATE(epsilon)
+  DATA epsilon/1.E-15/
+
+  ! Autres variables internes
+  INTEGER                                               ::isubstep, k, i
+
+  REAL                                                  :: wdens_targ
+  REAL                                                  :: sigmaw_targ
+
+  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
+! pas besoin d'isos ici
+
+  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)                          :: 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
+#ifdef ISO
+      REAL, DIMENSION(ntraciso,klon,klev) :: xt1, xt2
+      REAL, DIMENSION(ntraciso,klon,klev) :: D_xt1, D_xt2, D_dxt
+      REAL, DIMENSION(ntraciso,klon,klev) :: omgbdxt
+      integer ixt
+#endif
+
+  REAL, DIMENSION (klon)                                :: ff, gg
+  REAL, DIMENSION (klon)                                :: wape2, cstar2, heff
+                                                        
+  REAL, DIMENSION (klon, klev)                          :: crep
+                                                        
+  REAL, DIMENSION (klon, klev)                          :: ppi
+
+  ! cc nrlmd
+  REAL, DIMENSION (klon)                                :: death_rate
+!!  REAL, DIMENSION (klon)                                :: nat_rate
+  REAL, DIMENSION (klon, klev)                          :: entr
+  REAL, DIMENSION (klon, klev)                          :: detr
+
+  REAL, DIMENSION(klon)                                 :: sigmaw_in             ! pour les prints
+  REAL, DIMENSION(klon)                                 :: awdens_in, wdens_in   ! pour les prints
+
+  ! -------------------------------------------------------------------------
+  ! Initialisations
+  ! -------------------------------------------------------------------------
+
+  ! print*, 'wake initialisations'
+!#ifdef ISOVERIF
+!        write(*,*) 'wake 358: entree'
+!#endif
+
+  ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
+  ! -------------------------------------------------------------------------
+
+!!  DATA wapecut, sigmad, hwmin/5., .02, 10./
+!!  DATA wapecut, sigmad, hwmin/1., .02, 10./
+  DATA sigmad, hwmin/.02, 10./
+!!  DATA wdensmin/1.e-12/
+  DATA wdensmin/1.e-14/
+  ! cc nrlmd
+  DATA sigmaw_max/0.4/
+  DATA dens_rate/0.1/
+  ! cc
+  ! Longueur de maille (en m)
+  ! -------------------------------------------------------------------------
+
+  ! ALON = 3.e5
+  ! alon = 1.E6
+
+  !  Provisionnal; to be suppressed when f_shear is parameterized
+  f_shear(:) = 1.       ! 0. for strong shear, 1. for weak shear
+
+
+  ! 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é surfacique de poche froide
+  ! -------------------------------------------------------------------------
+
+  ! cc nrlmd      coefgw=10
+  ! coefgw=1
+  ! wdens0 = 1.0/(alon**2)
+  ! cc nrlmd      wdens = 1.0/(alon**2)
+  ! cc nrlmd      stark = 0.50
+  ! CRtest
+  ! cc nrlmd      alpk=0.1
+  ! alpk = 1.0
+  ! alpk = 0.5
+  ! alpk = 0.05
+
+ if (first) then
+
+  igout = klon/2+1/klon
+
+  crep_upper = 0.9
+  crep_sol = 1.0
+
+  ! Get wapecut from parameter file
+  wapecut = 1.
+  CALL getin_p('wapecut', wapecut)
+
+  ! cc nrlmd Lecture du fichier wake_param.data
+  stark=0.33
+  CALL getin_p('stark',stark)
+  cstart = stark*sqrt(2.*wapecut)
+
+  alpk=0.25
+  CALL getin_p('alpk',alpk)
+!jyg<
+!!  wdens_ref=8.E-12
+!!  CALL getin_p('wdens_ref',wdens_ref)
+  wdens_ref(1)=8.E-12
+  wdens_ref(2)=8.E-12
+  CALL getin_p('wdens_ref_o',wdens_ref(1))    !wake number per unit area ; ocean
+  CALL getin_p('wdens_ref_l',wdens_ref(2))    !wake number per unit area ; land
+!>jyg
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!  Population dynamics parameters    !!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!------------------------------------------------------------------------ 
+
+  iflag_wk_pop_dyn = 0
+  CALL getin_p('iflag_wk_pop_dyn',iflag_wk_pop_dyn) ! switch between wdens prescribed 
+                                                    ! and wdens prognostic
+  iflag_wk_act = 0
+  CALL getin_p('iflag_wk_act',iflag_wk_act) ! 0: act(:)=0.
+                                            ! 1: act(:)=1.
+                                            ! 2: act(:)=f(Wape)
+
+  rzero = 5000.
+  CALL getin_p('rzero_wk', rzero)
+  aa0 = 3.14*rzero*rzero
+!
+  tau_cv = 4000.
+  CALL getin_p('tau_cv', tau_cv)
+
+!------------------------------------------------------------------------ 
+
+  coefgw=4.
+  CALL getin_p('coefgw',coefgw)
+
+  WRITE(*,*) 'stark=', stark
+  WRITE(*,*) 'alpk=', alpk
+!jyg<
+!!  WRITE(*,*) 'wdens_ref=', wdens_ref
+  WRITE(*,*) 'wdens_ref_o=', wdens_ref(1)
+  WRITE(*,*) 'wdens_ref_l=', wdens_ref(2)
+!>jyg
+  WRITE(*,*) 'iflag_wk_pop_dyn=',iflag_wk_pop_dyn
+  WRITE(*,*) 'iflag_wk_act',iflag_wk_act
+  WRITE(*,*) 'coefgw=', coefgw 
+
+  flag_wk_check_trgl=.false.
+  CALL getin_p('flag_wk_check_trgl ', flag_wk_check_trgl)
+  WRITE(*,*) 'flag_wk_check_trgl=', flag_wk_check_trgl
+  WRITE(*,*) 'flag_wk_check_trgl OBSOLETE. Utilisr iflag_wk_check_trgl plutot'
+  iflag_wk_check_trgl=0 ; IF (flag_wk_check_trgl) iflag_wk_check_trgl=1
+  CALL getin_p('iflag_wk_check_trgl ', iflag_wk_check_trgl)
+  WRITE(*,*) 'iflag_wk_check_trgl=', iflag_wk_check_trgl
+
+  first=.false.
+ endif
+
+ IF (iflag_wk_pop_dyn == 0) THEN
+  ! 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)
+  !jyg<
+  !!  wdens(:) = wdens_ref
+   DO i = 1,klon
+     wdens(i) = wdens_ref(znatsurf(i)+1)
+   ENDDO 
+  !>jyg
+ ENDIF  ! (iflag_wk_pop_dyn == 0)
+
+  ! 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, initialize tendencies, initialize output fields
+  ! ------------------------------------------------------------------------
+
+!jyg<
+!!  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
+      ppi(:,:) = pi(:,:)
+      deltatw0(:,:) = deltatw(:,:)
+      deltaqw0(:,:) = deltaqw(:,:)
+      te(:,:) = te0(:,:)
+      qe(:,:) = qe0(:,:)
+      dtls(:,:) = 0.
+      dqls(:,:) = 0.
+      d_deltat_gw(:,:) = 0.
+      d_te(:,:) = 0.
+      d_qe(:,:) = 0.
+      d_deltatw(:,:) = 0.
+      d_deltaqw(:,:) = 0.
+      d_deltatw2(:,:) = 0.
+      d_deltaqw2(:,:) = 0. 
+#ifdef ISO
+      deltaxtw0(:,:,:)= deltaxtw(:,:,:)
+      xte(:,:,:) = xte0(:,:,:)
+      dxtls(:,:,:) = 0.
+      d_xte(:,:,:) = 0.
+      d_deltaxtw(:,:,:) = 0. 
+      d_deltaxtw2(:,:,:)=0. 
+      xt1(:,:,:) = 0. 
+      xt2(:,:,:)=0.   
+      ! init non indispensable mais facilite verif iso
+      q1(:,:)=0.0
+      q2(:,:)=0.0
+#endif
+
+      IF (iflag_wk_act == 0) THEN
+        act(:) = 0.
+      ELSEIF (iflag_wk_act == 1) THEN
+        act(:) = 1.
+      ENDIF
+
+!!  DO i = 1, klon
+!!   sigmaw_in(i) = sigmaw(i)
+!!  END DO
+   sigmaw_in(:) = sigmaw(:)
+!>jyg
+
+  ! sigmaw1=sigmaw
+  ! IF (sigd_con.GT.sigmaw1) THEN
+  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
+  ! ENDIF
+  IF (iflag_wk_pop_dyn >=1) THEN
+    DO i = 1, klon
+      wdens_targ = max(wdens(i),wdensmin)
+      d_wdens2(i) = wdens_targ - wdens(i)
+      wdens(i) = wdens_targ
+    END DO
+  ELSE
+    DO i = 1, klon
+      d_awdens2(i) = 0.
+      d_wdens2(i) = 0.
+    END DO
+  ENDIF  ! (iflag_wk_pop_dyn >=1)
+!
+  DO i = 1, klon
+    ! c      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
+!jyg<
+!!    sigmaw(i) = amax1(sigmaw(i), sigmad)
+!!    sigmaw(i) = amin1(sigmaw(i), 0.99)
+    sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
+    d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
+    sigmaw(i) = sigmaw_targ
+!>jyg
+  END DO
+
+!
+  IF (iflag_wk_pop_dyn >= 1) THEN
+    awdens_in(:) = awdens(:)
+    wdens_in(:) = wdens(:)
+!!    wdens(:) = wdens(:) + wgen(:)*dtime
+!!    d_wdens2(:) = wgen(:)*dtime
+!!  ELSE
+  ENDIF  ! (iflag_wk_pop_dyn >= 1)
+
+  wape(:) = 0.
+  wape2(:) = 0.
+  d_sigmaw(:) = 0.
+  ktopw(:) = 0
+!
+!<jyg
+dth(:,:) = 0.
+tu(:,:) = 0.
+qu(:,:) = 0.
+dtke(:,:) = 0.
+dqke(:,:) = 0.
+spread(:,:) = 0.
+omgbdth(:,:) = 0.
+omg(:,:) = 0.
+dp_omgb(:,:) = 0.
+dp_deltomg(:,:) = 0.
+hw(:) = 0.
+wape(:) = 0.
+fip(:) = 0.
+gfl(:) = 0.
+cstar(:) = 0.
+ktopw(:) = 0
+!
+!  Vertical advection local variables
+omgbw(:,:) = 0.
+omgtop(:) = 0
+dp_omgbw(:,:) = 0.
+omgbdq(:,:) = 0.
+#ifdef ISO
+xtu(:,:,:) = 0.
+dxtke(:,:,:) = 0.
+omgbdxt(:,:,:) = 0.
+#endif
+!>jyg
+!
+  IF (prt_level>=10) THEN
+    PRINT *, 'wake-1, sigmaw(igout) ', sigmaw(igout)
+    PRINT *, 'wake-1, deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev)
+    PRINT *, 'wake-1, deltaqw(igout,k) ', (k,deltaqw(igout,k), k=1,klev)
+    PRINT *, 'wake-1, dowwdraughts, amdwn(igout,k) ', (k,amdwn(igout,k), k=1,klev)
+    PRINT *, 'wake-1, dowwdraughts, dtdwn(igout,k) ', (k,dtdwn(igout,k), k=1,klev)
+    PRINT *, 'wake-1, dowwdraughts, dqdwn(igout,k) ', (k,dqdwn(igout,k), k=1,klev)
+    PRINT *, 'wake-1, updraughts, amup(igout,k) ', (k,amup(igout,k), k=1,klev)
+    PRINT *, 'wake-1, updraughts, dta(igout,k) ', (k,dta(igout,k), k=1,klev)
+    PRINT *, 'wake-1, updraughts, dqa(igout,k) ', (k,dqa(igout,k), k=1,klev)
+  ENDIF
+
+  ! 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.
+    ! pas besoin d'isos ici
+  END DO
+
+
+#ifdef ISOVERIF
+      ! TODO      
+#endif 
+
+  ! 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)
+#ifdef ISO
+        do ixt=1,ntraciso
+          xtu(ixt,i,k)  =  xte(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)
+        enddo !do ixt=1,ntraciso
+#ifdef ISOVERIF
+        if (iso_eau.gt.0) then
+              call iso_verif_egalite(deltaqw(i,k),deltaxtw(iso_eau,i,k),'wake 723a')
+              call iso_verif_egalite(qu(i,k),xtu(iso_eau,i,k),'wake 723b')
+        endif
+        if (iso_HDO.gt.0) then
+               if (iso_verif_aberrant_enc_choix_nostop(xtu(iso_hdo,i,k),qu(i,k),ridicule,deltalim, &
+                        'wake 723c xtu').eq.1) then
+                     stop
+               endif
+         endif  
+#endif
+#endif
+    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
+
+  IF (prt_level>=10) THEN
+    PRINT *, 'wake-2, ptop_provis(igout), ptop(igout) ', ptop_provis(igout), ptop(igout)
+  ENDIF
+
+
+  ! -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
+
+  IF (prt_level>=10) THEN
+    PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
+  ENDIF
+
+  ! -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.
+        d_deltatw2(i,k) = -deltatw0(i,k)
+#ifdef ISO
+        do ixt=1,ntraciso
+        deltaxtw(ixt,i, k) = 0.
+        d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
+        enddo
+#endif
+      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.+epsim1*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.+epsim1*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)+ &
+        epsim1*(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.
+        d_deltatw2(i,k) = -deltatw0(i,k)
+        d_deltaqw2(i,k) = -deltaqw0(i,k)
+#ifdef ISO
+        do ixt=1,ntraciso
+        deltaxtw(ixt,i, k) = 0.
+        d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
+        enddo
+#endif
+      END IF
+    END DO
+  END DO
+
+  DO i = 1, klon
+    IF (wape(i)<0.) THEN
+      wape(i) = 0.
+      cstar(i) = 0.
+      hw(i) = hwmin
+!jyg<
+!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
+      sigmaw_targ = max(sigmad, sigd_con(i))
+      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
+      sigmaw(i) = sigmaw_targ
+!>jyg
+      fip(i) = 0.
+      gwake(i) = .FALSE.
+    ELSE
+      hw(i) = hw0(i)
+      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
+
+  IF (prt_level>=10) THEN
+    PRINT *, 'wake-4, sigmaw(igout), cstar(igout), wape(igout), ktop(igout) ', &
+                      sigmaw(igout), cstar(igout), wape(igout), ktop(igout)
+  ENDIF
+
+
+  ! 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
+    IF (prt_level>=10) THEN
+      PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', &
+                          isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout)
+    ENDIF
+
+    ! 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 ---
+    !jyg<
+    ! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou 
+    ! les poches sont insuffisantes pour accueillir tout le flux de masse 
+    ! des descentes unsaturees. Nous faisons alors l'hypothese que la 
+    ! convection profonde cree directement de nouvelles poches, sans passer 
+    ! par les thermiques. La nouvelle valeur de wdens est alors imposée.
+
+    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 (prt_level >= 10) THEN
+             print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)', &
+                     omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)
+        ENDIF
+        IF (wdens(i)<=wdens0*1.1) THEN
+          IF (iflag_wk_pop_dyn >= 1) THEN
+             d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
+          ENDIF
+          wdens(i) = wdens0
+        END IF
+      END IF
+    END DO
+
+    DO i = 1, klon
+      IF (wk_adv(i)) THEN
+        gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
+        rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
+!jyg<
+!!        sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
+        sigmaw_targ = min(sigmaw(i), sigmaw_max)
+        d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
+        sigmaw(i) = sigmaw_targ
+!>jyg
+      END IF
+    END DO
+
+    IF (iflag_wk_pop_dyn >= 1) THEN
+!    The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.
+!    Here, it has to be set to zero.
+      death_rate(:) = 0.
+
+      IF (iflag_wk_act ==2) THEN
+      DO i = 1, klon
+        IF (wk_adv(i)) THEN
+          wape1_act(i) = abs(cin(i))
+          wape2_act(i) = 2.*wape1_act(i) + 1.
+          act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
+        ENDIF  ! (wk_adv(i))
+      ENDDO
+      ENDIF  ! (iflag_wk_act ==2)
+
+
+      DO i = 1, klon
+        IF (wk_adv(i)) THEN
+!!          tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.)
+          tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
+          tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
+          drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
+                    (1 + 2*f_shear(i)*(2.*sigmaw(i)-aa0*wdens(i)) - 2.*sigmaw(i))
+!!                    (1 - 2*sigmaw(i)*(1.-f_shear(i)))
+          drdt_pos=max(drdt(i),0.)
+
+!!          d_wdens(i) = ( wgen(i)*(1.+2.*(sigmaw(i)-sigmad)) &
+!!                     - wdens(i)*tau_wk_inv_min &
+!!                     - 2.*gfl(i)*wdens(i)*Cstar(i) )*dtimesub
+          d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub
+          d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min -  &
+                         2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub
+          d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
+
+!!          d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) &
+!!                      + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) &
+!!                      - sigmaw(i)*tau_wk_inv_min )*dtimesub
+          d_sig_gen(i) = wgen(i)*aa0
+          d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min
+!!          d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos
+          d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos
+          d_sigmaw(i) = ( d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + gfl(i)*cstar(i) )*dtimesub
+          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
+        ENDIF
+      ENDDO
+
+      IF (prt_level >= 10) THEN
+        print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
+                       cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
+        print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
+                       wdens(1), awdens(1), act(1), d_awdens(1)
+        print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', &
+                       wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1)
+        print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
+                       d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
+      ENDIF
+    
+    ELSE  !  (iflag_wk_pop_dyn >= 1)
+
+    ! cc nrlmd
+
+      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
+
+    ENDIF   !  (iflag_wk_pop_dyn >= 1)
+
+
+    ! 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..?
+    !JYG 161013 Correction : maintenant omg est dimensionne a klev.
+    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
+      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
+
+    IF (prt_level>=10) THEN
+      PRINT *, 'wake-4.2, omg(igout,k) ', (k,omg(igout,k), k=1,klev)
+      PRINT *, 'wake-4.2, omgtop(igout), ptop(igout), ktop(igout) ', &
+                          omgtop(igout), ptop(igout), ktop(igout)
+    ENDIF
+
+    ! -----------------
+    ! 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
+!!    print *,'omg(igout,k) ', (k,omg(igout,k),k=1,klev)
+    ! 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
+      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-1
+      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
+    DO i = 1, klon
+      IF (wk_adv(i)) THEN
+          dp_omgbw(i, klev) = 0.
+      END IF
+    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
+#ifdef ISO
+          do ixt=1,ntraciso
+            xt1(ixt,i,k) = xte(ixt,i,k) -sigmaw(i)     *deltaxtw(ixt,i,k) ! undisturbed area
+            xt2(ixt,i,k) = xte(ixt,i,k) +(1.-sigmaw(i))*deltaxtw(ixt,i,k) ! wake
+          enddo
+#endif
+        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.
+#ifdef ISO
+        do ixt=1,ntraciso
+          d_xt1(ixt,i,1) = 0.
+          d_xt2(ixt,i,1) = 0.
+          d_dxt(ixt,i,1) = 0.
+        enddo !do ixt=1,ntraciso
+#endif
+      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)
+#ifdef ISO
+        do ixt=1,ntraciso
+          d_xt1(ixt,i,k) = xt1(ixt,i,k-1)-xt1(ixt,i,k)
+          d_xt2(ixt,i,k) = xt2(ixt,i,k-1)-xt2(ixt,i,k)
+          d_dxt(ixt,i,k) = deltaxtw(ixt,i,k-1)-deltaxtw(ixt,i,k)
+        enddo !do ixt=1,ntraciso
+#endif
+        END IF
+      END DO
+    END DO
+
+    DO i = 1, klon
+      IF (wk_adv(i)) THEN
+        omgbdth(i, 1) = 0.
+        omgbdq(i, 1) = 0.
+#ifdef ISO
+        do ixt=1,ntraciso
+         omgbdxt(ixt,i,1) = 0.
+        enddo !do ixt=1,ntraciso
+#endif
+      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))
+#ifdef ISO
+        do ixt=1,ntraciso
+          omgbdxt(ixt,i,k)  = omgb(i,k)*(deltaxtw(ixt,i,k-1) - deltaxtw(ixt,i,k))
+        enddo !do ixt=1,ntraciso
+#ifdef ISOVERIF
+      if (iso_eau.gt.0) then
+        call iso_verif_egalite(deltaqw(i,k-1),deltaxtw(iso_eau,i,k-1),'wake 1460a')
+        call iso_verif_egalite(deltaqw(i,k),deltaxtw(iso_eau,i,k),'wake 1460b')
+        call iso_verif_egalite(omgbdq(i,k),omgbdxt(iso_eau,i,k),'wake 1460c')
+      endif
+#endif
+#endif
+        END IF
+      END DO
+    END DO
+
+    IF (prt_level>=10) THEN
+      PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,klev)
+      PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,klev)
+      PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,klev)
+      PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,klev)
+    ENDIF
+
+    ! -----------------------------------------------------------------
+    DO k = 1, klev-1
+      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=', k, 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=', k, d_deltaqw(i,k)
+#ifdef ISO
+        do ixt=1,ntraciso
+          d_deltaxtw(ixt,i,k) = dtimesub/(Ph(i,k)-Ph(i,k+1))* &
+            (rrd1*omg(i,k  )*sigmaw(i)     *d_xt1(ixt,i,k)- &
+            rrd2*omg(i,k+1)*(1.-sigmaw(i))*d_xt2(ixt,i,k+1)- &
+            (1.-alpha_up(i,k))*omgbdxt(ixt,i,k)- &
+            alpha_up(i,k+1)*omgbdxt(ixt,i,k+1))
+        enddo !do ixt=1,ntraciso
+#ifdef ISOVERIF
+      if (iso_eau.gt.0) then
+        call iso_verif_egalite(d_q1(i,k),d_xt1(iso_eau,i,k),'wake 1502a')
+        call iso_verif_egalite(d_q2(i,k),d_xt2(iso_eau,i,k),'wake 1502b')
+        call iso_verif_egalite(omgbdq(i,k),omgbdxt(iso_eau,i,k),'wake 1502c')
+        call iso_verif_egalite(omgbdq(i,k+1),omgbdxt(iso_eau,i,k+1),'wake 1502d')
+        call iso_verif_egalite(d_deltaqw(i,k),d_deltaxtw(iso_eau,i,k),'wake 1502e')
+      endif
+#endif
+#endif
+
+          ! 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)) &
+                                 -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/ &
+                                 (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)- &
+                                  rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/ &
+                                 (ph(i,k)-ph(i,k+1)) &
+                                 -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i,k+1))/ &
+                                 (ph(i,k)-ph(i,k+1)) )
+#ifdef ISO
+        do ixt=1,ntraciso
+         d_xte(ixt,i,k) =  dtimesub*( &
+             ( rre1(i)*omg(i,k  )*sigmaw(i)     *d_xt1(ixt,i,k) &
+              -rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_xt2(ixt,i,k+1) ) &
+                    /(Ph(i,k)-Ph(i,k+1)) &
+              -sigmaw(i)*(1.-sigmaw(i))*deltaxtw(ixt,i,k) &
+              *(omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1)) &
+              )
+        enddo !do ixt=1,ntraciso
+#endif 
+        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)))
+
+#ifdef ISO
+        do ixt=1,ntraciso
+         d_xte(ixt,i,k) =  dtimesub*( &
+            ( rre1(i)*omg(i,k  )*sigmaw(i)  *d_xt1(ixt,i,k)     &   
+             /(Ph(i,k)-Ph(i,k+1))) &
+                            )
+        enddo !do ixt=1,ntraciso
+#endif
+        END IF
+        ! cc
+      END DO
+    END DO
+    ! ------------------------------------------------------------------
+
+    IF (prt_level>=10) THEN
+      PRINT *, 'wake-4.3, d_deltatw(igout,k) ', (k,d_deltatw(igout,k), k=1,klev)
+      PRINT *, 'wake-4.3, d_deltaqw(igout,k) ', (k,d_deltaqw(igout,k), k=1,klev)
+    ENDIF
+
+    ! Increment state variables
+!jyg<
+    IF (iflag_wk_pop_dyn >= 1) THEN
+      DO k = 1, klev
+        DO i = 1, klon
+          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
+            detr(i,k) = - d_sig_death(i) - d_sig_col(i)      
+            entr(i,k) = d_sig_gen(i)
+          ENDIF
+        ENDDO
+      ENDDO
+      ELSE  ! (iflag_wk_pop_dyn >= 1)
+      DO k = 1, klev
+        DO i = 1, klon
+          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
+            detr(i, k) = 0.
+   
+            entr(i, k) = 0.
+          ENDIF
+        ENDDO
+      ENDDO
+    ENDIF  ! (iflag_wk_pop_dyn >= 1)
+
+    
+
+    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)))
+#ifdef ISO
+        do ixt=1,ntraciso
+          dxtke(ixt,i,k)=(dxtdwn(ixt,i,k)/sigmaw(i) - dxta(ixt,i,k) &
+               /(1.-sigmaw(i)))
+        enddo !do ixt=1,ntraciso
+#ifdef ISOVERIF
+      if (iso_eau.gt.0) then
+        call iso_verif_egalite(dqke(i,k),dxtKE(iso_eau,i,k),'wake 1621a')
+        call iso_verif_egalite(dqdwn(i,k),dxtdwn(iso_eau,i,k),'wake 1621b')
+        call iso_verif_egalite(dqa(i,k),dxta(iso_eau,i,k),'wake 1621c')
+        call iso_verif_egalite(d_deltaqw(i,k),d_deltaxtw(iso_eau,i,k),'wake 1621d')
+      endif
+#endif
+#endif
+          ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
+
+!
+
+          ! cc nrlmd          Prise en compte du taux de mortalité
+          ! cc               Définitions de entr, detr
+!jyg<
+!!            detr(i, k) = 0.
+!!   
+!!            entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
+!!              sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)
+!!
+            entr(i, k) = entr(i,k) + gfl(i)*cstar(i) + &
+                         sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)   
+!>jyg
+            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) - & 
+               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) - &
+                entr(i,k)*deltatw(i,k)/sigmaw(i) - &
+                (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) - &
+                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) - & 
+            entr(i,k)*deltaqw(i,k)/sigmaw(i) - &
+            (death_rate(i)*sigmaw(i)+detr(i,k))*deltaqw(i,k)/(1.-sigmaw(i)))
+#ifdef ISO
+        do ixt=1,ntraciso
+          gg(i)=d_deltaxtw(ixt,i,k)/dtimesub
+          d_deltaxtw(ixt,i,k) = dtimesub*(gg(i) + dxtKE(ixt,i,k) - &
+             entr(i,k)*deltaxtw(ixt,i,k)/sigmaw(i) - &
+             (death_rate(i)*sigmaw(i)+detr(i,k))*deltaxtw(ixt,i,k)/(1.-sigmaw(i)))
+        enddo !do ixt=1,ntraciso
+#ifdef ISOVERIF
+      if (iso_eau.gt.0) then
+        call iso_verif_egalite(dqke(i,k),dxtKE(iso_eau,i,k),'wake 1692a')
+        call iso_verif_egalite(deltaqw(i,k),deltaxtw(iso_eau,i,k),'wake 1692b')
+        call iso_verif_egalite(d_deltaqw(i,k),d_deltaxtw(iso_eau,i,k),'wake 1692c')
+      endif
+#endif
+#endif
+          ! 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
+
+#ifdef ISO
+#ifdef ISOVERIF
+      if (iso_eau.gt.0) then
+             call iso_verif_egalite_vect2D(d_deltaxtw,d_deltaqw, &
+                 'wake 1359',ntraciso,klon,klev)
+      endif      
+#endif         
+#endif  
+
+    ! 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)
+#ifdef ISO
+        do ixt=1,ntraciso
+          d_xte(ixt,i,k)=alpha(i)*d_xte(ixt,i,k)
+          d_deltaxtw(ixt,i,k)=alpha(i)*d_deltaxtw(ixt,i,k)
+        enddo !do ixt=1,ntraciso
+#endif
+        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)
+#ifdef ISO
+        do ixt=1,ntraciso
+          dxtls(ixt,i,k)=dxtls(ixt,i,k)+d_xte(ixt,i,k)
+        enddo !do ixt=1,ntraciso
+#endif
+          ! 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)
+#ifdef ISO
+        do ixt=1,ntraciso
+          d_deltaxtw2(ixt,i,k)=d_deltaxtw2(ixt,i,k)+d_deltaxtw(ixt,i,k)
+        enddo !do ixt=1,ntraciso
+#endif
+          ! cc
+        END IF
+      END DO
+    END DO
+
+
+#ifdef ISO
+#ifdef ISOVERIF
+      if (iso_eau.gt.0) then
+        DO k= 1,klev 
+          DO i = 1,klon 
+              call iso_verif_egalite_choix(dxtls(iso_eau,i,k), &
+                dqls(i,k),'wake 1379',errmax,errmaxrel)
+          enddo ! DO i = 1,klon     
+        enddo ! DO k= 1,klev 
+      endif !if (iso_eau.gt.0) then
+#endif
+#endif
+
+
+    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)
+#ifdef ISO
+        do ixt=1,ntraciso
+          xte(ixt,i,k) = xte0(ixt,i,k) + dxtls(ixt,i,k)
+          deltaxtw(ixt,i,k) = deltaxtw(ixt,i,k)+d_deltaxtw(ixt,i,k)
+        enddo !do ixt=1,ntraciso
+#endif
+        END IF
+      END DO
+    END DO
+!
+    DO i = 1, klon
+      IF (wk_adv(i)) THEN
+        sigmaw(i) = sigmaw(i) + d_sigmaw(i)
+        d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
+      END IF
+    END DO
+
+#ifdef ISO
+#ifdef ISOVERIF
+      write(*,*) 'wake 1859'
+      if (iso_eau.gt.0) then
+        DO k= 1,klev 
+          DO i = 1,klon 
+              call iso_verif_egalite_choix(xte(iso_eau,i,k), &
+                qe(i,k),'wake 1379',errmax,errmaxrel)
+          enddo ! DO i = 1,klon     
+        enddo ! DO k= 1,klev 
+      endif !if (iso_eau.gt.0) then     
+      if (iso_hdo.gt.0) then
+      call iso_verif_aberrant_enc_vect2D( &
+            xte,qe, &
+            'wake 1456, xte apres modifs',ntraciso,klon,klev)
+!      call iso_verif_aberrant_enc_vect2D_ns(
+!     :       deltaxtw,deltaqw,
+!     :       'wake 1518, deltaqw apres modifs',ntraciso,klon,klev)
+      endif
+#endif
+#endif 
+
+!jyg<
+    IF (iflag_wk_pop_dyn >= 1) THEN
+      DO i = 1, klon
+        IF (wk_adv(i)) THEN
+          awdens(i) = awdens(i) + d_awdens(i)
+          wdens(i) = wdens(i) + d_wdens(i)
+          d_awdens2(i) = d_awdens2(i) + d_awdens(i)
+          d_wdens2(i) = d_wdens2(i) + d_wdens(i)
+        END IF
+      END DO
+      DO i = 1, klon
+        IF (wk_adv(i)) THEN
+          wdens_targ = max(wdens(i),wdensmin)
+          d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
+          wdens(i) = wdens_targ
+!
+          wdens_targ = min( max(awdens(i),0.), wdens(i) )
+          d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i)
+          awdens(i) = wdens_targ
+        END IF
+      END DO
+      DO i = 1, klon
+        IF (wk_adv(i)) THEN
+          sigmaw_targ = max(sigmaw(i),sigmad)
+          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
+          sigmaw(i) = sigmaw_targ
+        END IF
+      END DO
+    ENDIF  ! (iflag_wk_pop_dyn >= 1)
+!>jyg
+
+
+    ! 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.
+          d_deltatw2(i,k) = -deltatw0(i,k)
+          d_deltaqw2(i,k) = -deltaqw0(i,k)
+#ifdef ISO
+        do ixt=1,ntraciso
+          deltaxtw(ixt,i,k) = 0.
+          d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
+        enddo !do ixt=1,ntraciso
+#endif
+        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.+epsim1*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.+epsim1*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)+epsim1*(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.
+            d_deltatw2(i,k) = -deltatw0(i,k)
+            d_deltaqw2(i,k) = -deltaqw0(i,k)
+#ifdef ISO
+        do ixt=1,ntraciso
+          deltaxtw(ixt,i,k) = 0.
+          d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
+        enddo !do ixt=1,ntraciso
+#endif
+          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
+!jyg<
+!!          sigmaw(i) = max(sigmad, sigd_con(i))
+          sigmaw_targ = max(sigmad, sigd_con(i))
+          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
+          sigmaw(i) = sigmaw_targ
+!>jyg
+          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
+
+  IF (prt_level>=10) THEN
+    PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', &
+                      sigmaw(igout), cstar(igout), wape(igout), ptop(igout)
+  ENDIF
+
+
+  ! ----------------------------------------------------------
+  ! Determine wake final state; recompute wape, cstar, ktop;
+  ! filter out bad wakes.
+  ! ----------------------------------------------------------
+
+  ! 2.1 - Undisturbed area and Wake integrals
+  ! ---------------------------------------------------------
+
+  DO i = 1, klon
+    ! cc nrlmd       if (wk_adv(i)) then !!! nrlmd
+    IF (ok_qx_qw(i)) THEN
+      ! cc
+      z(i) = 0.
+      sum_thu(i) = 0.
+      sum_tu(i) = 0.
+      sum_qu(i) = 0.
+      sum_thvu(i) = 0.
+      sum_dth(i) = 0.
+      sum_half_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.
+
+      dthmin(i) = -delta_t_min
+    END IF
+  END DO
+  ! Potential temperatures and humidity
+  ! ----------------------------------------------------------
+
+  DO k = 1, klev
+    DO i = 1, klon
+      ! cc nrlmd       IF ( wk_adv(i)) THEN
+      IF (ok_qx_qw(i)) THEN
+        ! cc
+        rho(i, k) = p(i, k)/(rd*te(i,k))
+        IF (k==1) THEN
+          rhoh(i, k) = ph(i, k)/(rd*te(i,k))
+          zhh(i, k) = 0
+        ELSE
+          rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
+          zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
+        END IF
+        the(i, k) = te(i, k)/ppi(i, k)
+        thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
+        tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
+        qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
+        rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
+        dth(i, k) = deltatw(i, k)/ppi(i, k)
+#ifdef ISO
+        do ixt=1,ntraciso
+          xtu(ixt,i,k)  =  xte(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)
+        enddo !do ixt=1,ntraciso
+#endif
+      END IF
+    END DO
+  END DO
+
+#ifdef ISO
+#ifdef ISOVERIF
+      if (iso_hdo.gt.0) then
+      call iso_verif_aberrant_enc_vect2D( &
+              xtu,qu, &
+              'wake 1834, apres modifs',ntraciso,klon,klev)
+      endif   
+#endif
+#endif   
+
+  ! 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.
+      dz_half(i) = 1.
+      sum_thvu(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)
+      sum_dth(i) = 0.
+    END IF
+  END DO
+
+  DO k = 1, klev
+    DO i = 1, klon
+      ! cc nrlmd       IF ( wk_adv(i)) THEN
+      IF (ok_qx_qw(i)) THEN
+        ! cc
+        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
+        dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-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.+epsim1*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)
+!
+          dthmin(i) = min(dthmin(i), dth(i,k))
+        END IF
+        IF (dz_half(i)>0) THEN
+          sum_half_dth(i) = sum_half_dth(i) + dth(i, k)*dz_half(i)
+        END IF
+      END IF
+    END DO
+  END DO
+
+  DO i = 1, klon
+    ! cc nrlmd       IF ( wk_adv(i)) THEN
+    IF (ok_qx_qw(i)) THEN
+      ! cc
+      hw0(i) = z(i)
+    END IF
+  END DO
+
+  ! - WAPE and mean forcing computation
+  ! -------------------------------------------------------------
+
+  ! Means
+
+  DO i = 1, klon
+    ! cc nrlmd       IF ( wk_adv(i)) THEN
+    IF (ok_qx_qw(i)) THEN
+      ! cc
+      av_thu(i) = sum_thu(i)/hw0(i)
+      av_tu(i) = sum_tu(i)/hw0(i)
+      av_qu(i) = sum_qu(i)/hw0(i)
+      av_thvu(i) = sum_thvu(i)/hw0(i)
+      av_dth(i) = sum_dth(i)/hw0(i)
+      av_dq(i) = sum_dq(i)/hw0(i)
+      av_rho(i) = sum_rho(i)/hw0(i)
+      av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
+      av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
+
+      wape2(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(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
+
+  IF (iflag_wk_check_trgl>=1) THEN
+    ! Check triangular shape of dth profile
+    DO i = 1, klon
+      IF (ok_qx_qw(i)) THEN
+        !! print *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i)
+        !! print *,'wake, 2.*sum_dth(i)/(hw0(i)*dthmin(i)) ', &
+        !!                2.*sum_dth(i)/(hw0(i)*dthmin(i))
+        !! print *,'wake, sum_half_dth(i), sum_dth(i) ', &
+        !!                sum_half_dth(i), sum_dth(i)
+        IF ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min) ) THEN
+          wape2(i) = -1.
+          !! print *,'wake, rej 1'
+        ELSE IF (iflag_wk_check_trgl==1.AND.abs(2.*sum_dth(i)/(hw0(i)*dthmin(i)) - 1.) > 0.5) THEN
+          wape2(i) = -1.
+          !! print *,'wake, rej 2'
+        ELSE IF (abs(sum_half_dth(i)) < 0.5*abs(sum_dth(i)) ) THEN
+          wape2(i) = -1.
+          !! print *,'wake, rej 3'
+        END IF
+      END IF
+    END DO
+  END IF
+
+
+  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.
+        d_deltatw2(i,k) = -deltatw0(i,k)
+        d_deltaqw2(i,k) = -deltaqw0(i,k)
+#ifdef ISO
+        do ixt=1,ntraciso
+          deltaxtw(ixt,i,k) = 0.
+          d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
+        enddo !do ixt=1,ntraciso
+#endif
+      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
+!jyg<
+!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
+      sigmaw_targ = max(sigmad, sigd_con(i))
+      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
+      sigmaw(i) = sigmaw_targ
+!>jyg
+        fip(i) = 0.
+        gwake(i) = .FALSE.
+      ELSE
+        IF (prt_level>=10) PRINT *, 'wape2>0'
+        cstar2(i) = stark*sqrt(2.*wape2(i))
+        gwake(i) = .TRUE.
+      END IF
+    END IF
+  END DO
+
+  DO i = 1, klon
+    ! cc nrlmd       IF ( wk_adv(i)) THEN
+    IF (ok_qx_qw(i)) THEN
+      ! cc
+      ktopw(i) = ktop(i)
+    END IF
+  END DO
+
+  DO i = 1, klon
+    ! cc nrlmd       IF ( wk_adv(i)) THEN
+    IF (ok_qx_qw(i)) THEN
+      ! cc
+      IF (ktopw(i)>0 .AND. gwake(i)) THEN
+
+        ! jyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
+        ! cc       heff = 600.
+        ! Utilisation de la hauteur hw
+        ! c       heff = 0.7*hw
+        heff(i) = hw(i)
+
+        fip(i) = 0.5*rho(i, ktopw(i))*cstar2(i)**3*heff(i)*2* &
+          sqrt(sigmaw(i)*wdens(i)*3.14)
+        fip(i) = alpk*fip(i)
+        ! jyg2
+      ELSE
+        fip(i) = 0.
+      END IF
+    END IF
+  END DO
+
+  ! Limitation de sigmaw
+
+  ! cc nrlmd
+  ! DO i=1,klon
+  ! IF (OK_qx_qw(i)) THEN
+  ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
+  ! ENDIF
+  ! ENDDO
+  ! cc
+
+  !jyg<
+  IF (iflag_wk_pop_dyn >= 1) THEN
+    DO i = 1, klon
+      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
+          .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
+    ENDDO
+  ELSE  ! (iflag_wk_pop_dyn >= 1)
+    DO i = 1, klon
+      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
+          .NOT. ok_qx_qw(i)
+    ENDDO
+  ENDIF  ! (iflag_wk_pop_dyn >= 1)
+  !>jyg
+
+  DO k = 1, klev
+    DO i = 1, klon
+!!jyg      IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
+!!jyg          .NOT. ok_qx_qw(i)) THEN
+      IF (kill_wake(i)) THEN
+        ! cc
+        dtls(i, k) = 0.
+        dqls(i, k) = 0.
+        deltatw(i, k) = 0.
+        deltaqw(i, k) = 0.
+        d_deltatw2(i,k) = -deltatw0(i,k)
+        d_deltaqw2(i,k) = -deltaqw0(i,k)
+#ifdef ISO
+        do ixt=1,ntraciso
+          dxtls(ixt,i,k) = 0.
+          deltaxtw(ixt,i,k) = 0.
+          d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
+        enddo !do ixt=1,ntraciso
+#endif
+      END IF  ! (kill_wake(i))
+    END DO
+  END DO
+
+  DO i = 1, klon
+!!jyg    IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
+!!jyg        .NOT. ok_qx_qw(i)) THEN
+      IF (kill_wake(i)) THEN
+      ktopw(i) = 0
+      wape(i) = 0.
+      cstar(i) = 0.
+!!jyg   Outside subroutine "Wake" hw, wdens and sigmaw are zero when there are no wakes
+!!      hw(i) = hwmin                       !jyg
+!!      sigmaw(i) = sigmad                  !jyg
+      hw(i) = 0.                            !jyg
+      fip(i) = 0.
+!!      sigmaw(i) = 0.                        !jyg
+      sigmaw_targ = 0.
+      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
+      sigmaw(i) = sigmaw_targ
+      IF (iflag_wk_pop_dyn >= 1) THEN
+!!        awdens(i) = 0.
+!!        wdens(i) = 0.
+        wdens_targ = 0.
+        d_wdens2(i) = wdens_targ - wdens(i)
+        wdens(i) = wdens_targ
+        wdens_targ = 0.
+        d_awdens2(i) = wdens_targ - awdens(i)
+        awdens(i) = wdens_targ
+      ENDIF  ! (iflag_wk_pop_dyn >= 1)
+    ELSE  ! (kill_wake(i))
+      wape(i) = wape2(i)
+      cstar(i) = cstar2(i)
+    END IF  ! (kill_wake(i))
+    ! c        print*,'wape wape2 ktopw OK_qx_qw =',
+    ! c     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
+  END DO
+
+  IF (prt_level>=10) THEN
+    PRINT *, 'wake-6, wape wape2 ktopw OK_qx_qw =', &
+                      wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout)
+  ENDIF
+
+
+  ! -----------------------------------------------------------------
+  ! 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
+!jyg<
+!!      IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
+      IF (ok_qx_qw(i)) THEN
+!>jyg
+        ! 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)
+#ifdef ISO
+        do ixt=1,ntraciso
+          dxtls(ixt,i, k) = dxtls(ixt,i, k)/dtime
+          d_deltaxtw2(ixt,i, k) = d_deltaxtw2(ixt,i, k)/dtime
+        enddo
+#endif
+      END IF
+    END DO
+  END DO
+!jyg<
+  DO i = 1, klon
+    d_sigmaw2(i) = d_sigmaw2(i)/dtime
+    d_awdens2(i) = d_awdens2(i)/dtime
+    d_wdens2(i) = d_wdens2(i)/dtime
+  ENDDO
+!>jyg
+
+
+
+  RETURN
+END SUBROUTINE wake
+
+SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, &
+    d_deltaqw, sigmaw, d_sigmaw, alpha)
+  ! ------------------------------------------------------
+  ! D\'etermination du coefficient alpha tel que les tendances
+  ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
+  ! a une humidite positive dans la zone (x) et dans la zone (w).
+  ! ------------------------------------------------------
+  IMPLICIT NONE
+
+  ! Input
+  REAL qe(nlon, nl), d_qe(nlon, nl)
+  REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
+  REAL sigmaw(nlon), d_sigmaw(nlon)
+  LOGICAL wk_adv(nlon)
+  INTEGER nl, nlon
+  ! Output
+  REAL alpha(nlon)
+  ! Internal variables
+  REAL zeta(nlon, nl)
+  REAL alpha1(nlon)
+  REAL x, a, b, c, discrim
+  REAL epsilon
+  ! DATA epsilon/1.e-15/
+  INTEGER i,k
+
+  DO k = 1, nl
+    DO i = 1, nlon
+      IF (wk_adv(i)) THEN
+        IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
+          zeta(i, k) = 0.
+        ELSE
+          zeta(i, k) = 1.
+        END IF
+      END IF
+    END DO
+    DO i = 1, nlon
+      IF (wk_adv(i)) THEN
+        x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
+          (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i) * &
+          (deltaqw(i,k)+d_deltaqw(i,k))
+        a = -d_sigmaw(i)*d_deltaqw(i, k)
+        b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
+          deltaqw(i, k)*d_sigmaw(i)
+        c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon
+        discrim = b*b - 4.*a*c
+        ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
+        IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap
+          alpha1(i) = 1.
+        ELSE
+          IF (x>=0.) THEN
+            alpha1(i) = 1.
+          ELSE
+            IF (a>0.) THEN
+              alpha1(i) = 0.9*min( (2.*c)/(-b+sqrt(discrim)),  &
+                                   (-b+sqrt(discrim))/(2.*a) )
+            ELSE IF (a==0.) THEN
+              alpha1(i) = 0.9*(-c/b)
+            ELSE
+              ! print*,'a,b,c discrim',a,b,c discrim
+              alpha1(i) = 0.9*max( (2.*c)/(-b+sqrt(discrim)),  &
+                                   (-b+sqrt(discrim))/(2.*a))
+            END IF
+          END IF
+        END IF
+        alpha(i) = min(alpha(i), alpha1(i))
+      END IF
+    END DO
+  END DO
+
+  RETURN
+END SUBROUTINE wake_vec_modulation
+
+END MODULE lmdz_wake
Index: LMDZ6/trunk/libf/phylmdiso/lmdz_wake_ini.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/lmdz_wake_ini.F90	(revision 4594)
+++ LMDZ6/trunk/libf/phylmdiso/lmdz_wake_ini.F90	(revision 4594)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_wake_ini.F90
Index: LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4581)
+++ LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4594)
@@ -80,8 +80,9 @@
 #endif
     USE lscp_mod, ONLY : lscp
-    USE wake_ini_mod, ONLY : wake_ini
+    USE lmdz_wake_ini, ONLY : wake_ini
     USE yamada_ini_mod, ONLY : yamada_ini
     USE atke_turbulence_ini_mod, ONLY : atke_ini
-    USE thermcell_ini_mod, ONLY : thermcell_ini
+    USE lmdz_thermcell_ini, ONLY : thermcell_ini
+    USE lmdz_thermcell_dtke, ONLY : thermcell_dtke
     USE blowing_snow_ini_mod, ONLY : blowing_snow_ini , qbst_bs 
     USE lscp_ini_mod, ONLY : lscp_ini
Index: LMDZ6/trunk/libf/phylmdiso/thermcell.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_alim.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_alim.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_alim.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_alp.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_alp.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_alp.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_closure.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_closure.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_closure.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_down.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_down.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_down.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_dq.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_dq.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_dq.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_dry.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_dry.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_dry.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_dtke.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_dtke.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_dtke.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_dv2.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_dv2.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_dv2.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_env.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_env.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_env.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_flux2.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_flux2.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_flux2.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_height.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_height.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_height.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_ini_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_ini_mod.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_ini_mod.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_main.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_main.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_main.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_old.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_old.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_old.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_plume.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_plume.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_plume.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_plume_6A.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_plume_6A.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_plume_6A.F90
Index: LMDZ6/trunk/libf/phylmdiso/thermcell_qsat.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/thermcell_qsat.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/thermcell_qsat.F90
Index: LMDZ6/trunk/libf/phylmdiso/wake.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/wake.F90	(revision 4581)
+++ 	(revision )
@@ -1,2631 +1,0 @@
-
-! $Id: wake.F90 3648 2020-03-16 15:36:59Z jghattas $
-
-SUBROUTINE wake(znatsurf, p, ph, pi, dtime, &
-                te0, qe0, omgb, &
-                dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
-                sigd_con, Cin, &
-                deltatw, deltaqw, sigmaw, awdens, wdens, &                  ! state variables
-                dth, hw, wape, fip, gfl, &
-                dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, &
-                dtke, dqke, omg, dp_deltomg, spread, cstar, &
-                d_deltat_gw, &
-                d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2 &     ! tendencies
-
-#ifdef ISO
-                     ,xte0,dxtdwn,dxta,deltaxtw &
-                     ,dxtls,xtu,dxtke,d_deltaxtw2 &
-#endif                  
-                )  
-
-
-  ! **************************************************************
-  ! *
-  ! WAKE                                                        *
-  ! retour a un Pupper fixe                                *
-  ! *
-  ! written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
-  ! modified by :   ROEHRIG Romain        01/29/2007            *
-  ! **************************************************************
-
-  USE ioipsl_getin_p_mod, ONLY : getin_p
-  USE dimphy
-  use mod_phys_lmdz_para
-  USE print_control_mod, ONLY: prt_level
-#ifdef ISO
-  USE infotrac_phy, ONLY : ntraciso=>ntiso
-#ifdef ISOVERIF
-  USE isotopes_verif_mod
-!, ONLY: errmax,errmaxrel
-  USE isotopes_mod, ONLY: iso_eau,iso_hdo,ridicule
-#endif
-#endif
-  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.
-
-  ! State variables : 
-  ! deltatw    : temperature difference between wake and off-wake regions
-  ! deltaqw    : specific humidity difference between wake and off-wake regions
-  ! sigmaw     : fractional area covered by wakes.
-  ! wdens      : number of wakes per unit area
-
-  ! 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   : wake top hight (given by hw*deltatw(1)/2=wape)
-  ! dp_omgb : vertical gradient of large scale omega
-  ! awdens  : densite de poches actives
-  ! 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 d_t_wake and d_q_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)
-  ! wgen : number of wakes generated per unit area and per sec (/m^2/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)
-  ! sigd_con: 
-  ! Cin  : convective inhibition
-  ! 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    : wake top hight (defined as the altitude at which deltatw=0)
-  ! amflux : horizontal mass flux through wake boundary
-  ! wdens_ref: initial number of wakes per unit area (3D) or per
-  ! unit length (2D), at the beginning of each time step
-  ! Tgw    : 1 sur la période de onde de gravité
-  ! Cgw    : vitesse de propagation de onde de gravité
-  ! LL     : distance entre 2 poches
-
-  ! -------------------------------------------------------------------------
-  ! Déclaration de variables
-  ! -------------------------------------------------------------------------
-
-  include "YOMCST.h"
-  include "cvthermo.h"
-
-  ! Arguments en entree
-  ! --------------------
-
-  INTEGER, DIMENSION (klon),        INTENT(IN)          :: znatsurf
-  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: p, pi
-  REAL, DIMENSION (klon, klev+1),   INTENT(IN)          :: ph
-  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: omgb
-  REAL,                             INTENT(IN)          :: dtime
-  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: te0, qe0
-  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dtdwn, dqdwn
-  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: amdwn, amup
-  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dta, dqa
-  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen
-  REAL, DIMENSION (klon),           INTENT(IN)          :: sigd_con
-  REAL, DIMENSION (klon),           INTENT(IN)          :: Cin
-#ifdef ISO
-  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(IN)          :: xte0
-  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(IN)          :: dxtdwn
-  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(IN)          :: dxta
-#endif
-
-  !
-  ! Input/Output
-  ! State variables
-  REAL, DIMENSION (klon, klev),     INTENT(INOUT)       :: deltatw, deltaqw
-  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw
-  REAL, DIMENSION (klon),           INTENT(INOUT)       :: awdens
-  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens
-#ifdef ISO
-  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(INOUT)       :: deltaxtw
-#endif
-
-  ! Sorties
-  ! --------
-
-  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dth
-  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: tu, qu
-  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtls, dqls
-  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke
-  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: spread    !  unused (jyg)
-  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: omgbdth, omg
-  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dp_omgb, dp_deltomg
-  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_gw
-  REAL, DIMENSION (klon),           INTENT(OUT)         :: hw, wape, fip, gfl, cstar
-  INTEGER, DIMENSION (klon),        INTENT(OUT)         :: ktopw
-  ! Tendencies of state variables
-  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
-  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_awdens2, d_wdens2
-#ifdef ISO
-  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(OUT)         :: xtu
-  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(OUT)         :: dxtls
-  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(OUT)         :: dxtke
-  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(OUT)         :: d_deltaxtw2
-#endif
-
-  ! Variables internes
-  ! -------------------
-
-  ! Variables à fixer
-  INTEGER, SAVE                                         :: igout
-  !$OMP THREADPRIVATE(igout)
-  LOGICAL, SAVE                                         :: first = .TRUE.
-  !$OMP THREADPRIVATE(first)
-!jyg<
-!!  REAL, SAVE                                            :: stark, wdens_ref, coefgw, alpk
-  REAL, SAVE, DIMENSION(2)                              :: wdens_ref
-  REAL, SAVE                                            :: stark, coefgw, alpk
-!>jyg
-  REAL, SAVE                                            :: crep_upper, crep_sol  
-  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
-
-  REAL, SAVE                                            :: tau_cv
-  !$OMP THREADPRIVATE(tau_cv)
-  REAL, SAVE                                            :: rzero, aa0 ! minimal wake radius and area
-  !$OMP THREADPRIVATE(rzero, aa0)
-
-  LOGICAL, SAVE                                         :: flag_wk_check_trgl
-  !$OMP THREADPRIVATE(flag_wk_check_trgl)
-  INTEGER, SAVE                                         :: iflag_wk_check_trgl
-  !$OMP THREADPRIVATE(iflag_wk_check_trgl)
-  INTEGER, SAVE                                         :: iflag_wk_pop_dyn
-  !$OMP THREADPRIVATE(iflag_wk_pop_dyn)
-
-  REAL                                                  :: delta_t_min
-  INTEGER                                               :: nsub
-  REAL                                                  :: dtimesub
-  REAL, SAVE                                            :: wdensmin
-  !$OMP THREADPRIVATE(wdensmin)
-  REAL, SAVE                                            :: sigmad, hwmin, wapecut, cstart
-  !$OMP THREADPRIVATE(sigmad, hwmin, wapecut, cstart)
-  REAL, SAVE                                            :: sigmaw_max
-  !$OMP THREADPRIVATE(sigmaw_max)  
-  REAL, SAVE                                            :: dens_rate
-  !$OMP THREADPRIVATE(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)                                :: sigmaw1
-#ifdef ISO
-  REAL, DIMENSION (ntraciso,klon, klev)                          :: deltaxtw0
-  REAL, DIMENSION (ntraciso,klon, klev)                          :: xte
-#endif
-
-  ! Variables liees a la dynamique de population
-  REAL, DIMENSION(klon)                                 :: act
-  REAL, DIMENSION(klon)                                 :: rad_wk, tau_wk_inv
-  REAL, DIMENSION(klon)                                 :: f_shear
-  REAL, DIMENSION(klon)                                 :: drdt
-  REAL, DIMENSION(klon)                                 :: d_sig_gen, d_sig_death, d_sig_col
-  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
-  LOGICAL, DIMENSION (klon)                             :: kill_wake
-  INTEGER, SAVE                                         :: iflag_wk_act
-  !$OMP THREADPRIVATE(iflag_wk_act)
-  REAL                                                  :: drdt_pos
-  REAL                                                  :: tau_wk_inv_min
-
-  ! Variables pour les GW
-  REAL, DIMENSION (klon)                                :: ll
-  REAL, DIMENSION (klon, klev)                          :: n2
-  REAL, DIMENSION (klon, klev)                          :: cgw
-  REAL, DIMENSION (klon, klev)                          :: tgw
-
-  ! Variables liees 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
-
-  ! Variables liees au test de la forme triangulaire du profil de Delta_theta
-  REAL, DIMENSION (klon)                                :: sum_half_dth
-  REAL, DIMENSION (klon)                                :: dz_half
-
-  ! Sub-timestep tendencies and related variables
-  REAL, DIMENSION (klon, klev)                          :: d_deltatw, d_deltaqw
-  REAL, DIMENSION (klon, klev)                          :: d_te, d_qe
-  REAL, DIMENSION (klon)                                :: d_awdens, d_wdens
-  REAL, DIMENSION (klon)                                :: d_sigmaw, alpha
-  REAL, DIMENSION (klon)                                :: q0_min, q1_min
-  LOGICAL, DIMENSION (klon)                             :: wk_adv, ok_qx_qw
-#ifdef ISO
-  REAL, DIMENSION (ntraciso,klon, klev)                          :: d_deltaxtw
-  REAL, DIMENSION (ntraciso,klon, klev)                          :: d_xte
-#endif
-  REAL, SAVE                                            :: epsilon
-  !$OMP THREADPRIVATE(epsilon)
-  DATA epsilon/1.E-15/
-
-  ! Autres variables internes
-  INTEGER                                               ::isubstep, k, i
-
-  REAL                                                  :: wdens_targ
-  REAL                                                  :: sigmaw_targ
-
-  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
-! pas besoin d'isos ici
-
-  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)                          :: 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
-#ifdef ISO
-      REAL, DIMENSION(ntraciso,klon,klev) :: xt1, xt2
-      REAL, DIMENSION(ntraciso,klon,klev) :: D_xt1, D_xt2, D_dxt
-      REAL, DIMENSION(ntraciso,klon,klev) :: omgbdxt
-      integer ixt
-#endif
-
-  REAL, DIMENSION (klon)                                :: ff, gg
-  REAL, DIMENSION (klon)                                :: wape2, cstar2, heff
-                                                        
-  REAL, DIMENSION (klon, klev)                          :: crep
-                                                        
-  REAL, DIMENSION (klon, klev)                          :: ppi
-
-  ! cc nrlmd
-  REAL, DIMENSION (klon)                                :: death_rate
-!!  REAL, DIMENSION (klon)                                :: nat_rate
-  REAL, DIMENSION (klon, klev)                          :: entr
-  REAL, DIMENSION (klon, klev)                          :: detr
-
-  REAL, DIMENSION(klon)                                 :: sigmaw_in             ! pour les prints
-  REAL, DIMENSION(klon)                                 :: awdens_in, wdens_in   ! pour les prints
-
-  ! -------------------------------------------------------------------------
-  ! Initialisations
-  ! -------------------------------------------------------------------------
-
-  ! print*, 'wake initialisations'
-!#ifdef ISOVERIF
-!        write(*,*) 'wake 358: entree'
-!#endif
-
-  ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
-  ! -------------------------------------------------------------------------
-
-!!  DATA wapecut, sigmad, hwmin/5., .02, 10./
-!!  DATA wapecut, sigmad, hwmin/1., .02, 10./
-  DATA sigmad, hwmin/.02, 10./
-!!  DATA wdensmin/1.e-12/
-  DATA wdensmin/1.e-14/
-  ! cc nrlmd
-  DATA sigmaw_max/0.4/
-  DATA dens_rate/0.1/
-  ! cc
-  ! Longueur de maille (en m)
-  ! -------------------------------------------------------------------------
-
-  ! ALON = 3.e5
-  ! alon = 1.E6
-
-  !  Provisionnal; to be suppressed when f_shear is parameterized
-  f_shear(:) = 1.       ! 0. for strong shear, 1. for weak shear
-
-
-  ! 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é surfacique de poche froide
-  ! -------------------------------------------------------------------------
-
-  ! cc nrlmd      coefgw=10
-  ! coefgw=1
-  ! wdens0 = 1.0/(alon**2)
-  ! cc nrlmd      wdens = 1.0/(alon**2)
-  ! cc nrlmd      stark = 0.50
-  ! CRtest
-  ! cc nrlmd      alpk=0.1
-  ! alpk = 1.0
-  ! alpk = 0.5
-  ! alpk = 0.05
-
- if (first) then
-
-  igout = klon/2+1/klon
-
-  crep_upper = 0.9
-  crep_sol = 1.0
-
-  ! Get wapecut from parameter file
-  wapecut = 1.
-  CALL getin_p('wapecut', wapecut)
-
-  ! cc nrlmd Lecture du fichier wake_param.data
-  stark=0.33
-  CALL getin_p('stark',stark)
-  cstart = stark*sqrt(2.*wapecut)
-
-  alpk=0.25
-  CALL getin_p('alpk',alpk)
-!jyg<
-!!  wdens_ref=8.E-12
-!!  CALL getin_p('wdens_ref',wdens_ref)
-  wdens_ref(1)=8.E-12
-  wdens_ref(2)=8.E-12
-  CALL getin_p('wdens_ref_o',wdens_ref(1))    !wake number per unit area ; ocean
-  CALL getin_p('wdens_ref_l',wdens_ref(2))    !wake number per unit area ; land
-!>jyg
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!  Population dynamics parameters    !!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!------------------------------------------------------------------------ 
-
-  iflag_wk_pop_dyn = 0
-  CALL getin_p('iflag_wk_pop_dyn',iflag_wk_pop_dyn) ! switch between wdens prescribed 
-                                                    ! and wdens prognostic
-  iflag_wk_act = 0
-  CALL getin_p('iflag_wk_act',iflag_wk_act) ! 0: act(:)=0.
-                                            ! 1: act(:)=1.
-                                            ! 2: act(:)=f(Wape)
-
-  rzero = 5000.
-  CALL getin_p('rzero_wk', rzero)
-  aa0 = 3.14*rzero*rzero
-!
-  tau_cv = 4000.
-  CALL getin_p('tau_cv', tau_cv)
-
-!------------------------------------------------------------------------ 
-
-  coefgw=4.
-  CALL getin_p('coefgw',coefgw)
-
-  WRITE(*,*) 'stark=', stark
-  WRITE(*,*) 'alpk=', alpk
-!jyg<
-!!  WRITE(*,*) 'wdens_ref=', wdens_ref
-  WRITE(*,*) 'wdens_ref_o=', wdens_ref(1)
-  WRITE(*,*) 'wdens_ref_l=', wdens_ref(2)
-!>jyg
-  WRITE(*,*) 'iflag_wk_pop_dyn=',iflag_wk_pop_dyn
-  WRITE(*,*) 'iflag_wk_act',iflag_wk_act
-  WRITE(*,*) 'coefgw=', coefgw 
-
-  flag_wk_check_trgl=.false.
-  CALL getin_p('flag_wk_check_trgl ', flag_wk_check_trgl)
-  WRITE(*,*) 'flag_wk_check_trgl=', flag_wk_check_trgl
-  WRITE(*,*) 'flag_wk_check_trgl OBSOLETE. Utilisr iflag_wk_check_trgl plutot'
-  iflag_wk_check_trgl=0 ; IF (flag_wk_check_trgl) iflag_wk_check_trgl=1
-  CALL getin_p('iflag_wk_check_trgl ', iflag_wk_check_trgl)
-  WRITE(*,*) 'iflag_wk_check_trgl=', iflag_wk_check_trgl
-
-  first=.false.
- endif
-
- IF (iflag_wk_pop_dyn == 0) THEN
-  ! 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)
-  !jyg<
-  !!  wdens(:) = wdens_ref
-   DO i = 1,klon
-     wdens(i) = wdens_ref(znatsurf(i)+1)
-   ENDDO 
-  !>jyg
- ENDIF  ! (iflag_wk_pop_dyn == 0)
-
-  ! 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, initialize tendencies, initialize output fields
-  ! ------------------------------------------------------------------------
-
-!jyg<
-!!  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
-      ppi(:,:) = pi(:,:)
-      deltatw0(:,:) = deltatw(:,:)
-      deltaqw0(:,:) = deltaqw(:,:)
-      te(:,:) = te0(:,:)
-      qe(:,:) = qe0(:,:)
-      dtls(:,:) = 0.
-      dqls(:,:) = 0.
-      d_deltat_gw(:,:) = 0.
-      d_te(:,:) = 0.
-      d_qe(:,:) = 0.
-      d_deltatw(:,:) = 0.
-      d_deltaqw(:,:) = 0.
-      d_deltatw2(:,:) = 0.
-      d_deltaqw2(:,:) = 0. 
-#ifdef ISO
-      deltaxtw0(:,:,:)= deltaxtw(:,:,:)
-      xte(:,:,:) = xte0(:,:,:)
-      dxtls(:,:,:) = 0.
-      d_xte(:,:,:) = 0.
-      d_deltaxtw(:,:,:) = 0. 
-      d_deltaxtw2(:,:,:)=0. 
-      xt1(:,:,:) = 0. 
-      xt2(:,:,:)=0.   
-      ! init non indispensable mais facilite verif iso
-      q1(:,:)=0.0
-      q2(:,:)=0.0
-#endif
-
-      IF (iflag_wk_act == 0) THEN
-        act(:) = 0.
-      ELSEIF (iflag_wk_act == 1) THEN
-        act(:) = 1.
-      ENDIF
-
-!!  DO i = 1, klon
-!!   sigmaw_in(i) = sigmaw(i)
-!!  END DO
-   sigmaw_in(:) = sigmaw(:)
-!>jyg
-
-  ! sigmaw1=sigmaw
-  ! IF (sigd_con.GT.sigmaw1) THEN
-  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
-  ! ENDIF
-  IF (iflag_wk_pop_dyn >=1) THEN
-    DO i = 1, klon
-      wdens_targ = max(wdens(i),wdensmin)
-      d_wdens2(i) = wdens_targ - wdens(i)
-      wdens(i) = wdens_targ
-    END DO
-  ELSE
-    DO i = 1, klon
-      d_awdens2(i) = 0.
-      d_wdens2(i) = 0.
-    END DO
-  ENDIF  ! (iflag_wk_pop_dyn >=1)
-!
-  DO i = 1, klon
-    ! c      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
-!jyg<
-!!    sigmaw(i) = amax1(sigmaw(i), sigmad)
-!!    sigmaw(i) = amin1(sigmaw(i), 0.99)
-    sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
-    d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
-    sigmaw(i) = sigmaw_targ
-!>jyg
-  END DO
-
-!
-  IF (iflag_wk_pop_dyn >= 1) THEN
-    awdens_in(:) = awdens(:)
-    wdens_in(:) = wdens(:)
-!!    wdens(:) = wdens(:) + wgen(:)*dtime
-!!    d_wdens2(:) = wgen(:)*dtime
-!!  ELSE
-  ENDIF  ! (iflag_wk_pop_dyn >= 1)
-
-  wape(:) = 0.
-  wape2(:) = 0.
-  d_sigmaw(:) = 0.
-  ktopw(:) = 0
-!
-!<jyg
-dth(:,:) = 0.
-tu(:,:) = 0.
-qu(:,:) = 0.
-dtke(:,:) = 0.
-dqke(:,:) = 0.
-spread(:,:) = 0.
-omgbdth(:,:) = 0.
-omg(:,:) = 0.
-dp_omgb(:,:) = 0.
-dp_deltomg(:,:) = 0.
-hw(:) = 0.
-wape(:) = 0.
-fip(:) = 0.
-gfl(:) = 0.
-cstar(:) = 0.
-ktopw(:) = 0
-!
-!  Vertical advection local variables
-omgbw(:,:) = 0.
-omgtop(:) = 0
-dp_omgbw(:,:) = 0.
-omgbdq(:,:) = 0.
-#ifdef ISO
-xtu(:,:,:) = 0.
-dxtke(:,:,:) = 0.
-omgbdxt(:,:,:) = 0.
-#endif
-!>jyg
-!
-  IF (prt_level>=10) THEN
-    PRINT *, 'wake-1, sigmaw(igout) ', sigmaw(igout)
-    PRINT *, 'wake-1, deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev)
-    PRINT *, 'wake-1, deltaqw(igout,k) ', (k,deltaqw(igout,k), k=1,klev)
-    PRINT *, 'wake-1, dowwdraughts, amdwn(igout,k) ', (k,amdwn(igout,k), k=1,klev)
-    PRINT *, 'wake-1, dowwdraughts, dtdwn(igout,k) ', (k,dtdwn(igout,k), k=1,klev)
-    PRINT *, 'wake-1, dowwdraughts, dqdwn(igout,k) ', (k,dqdwn(igout,k), k=1,klev)
-    PRINT *, 'wake-1, updraughts, amup(igout,k) ', (k,amup(igout,k), k=1,klev)
-    PRINT *, 'wake-1, updraughts, dta(igout,k) ', (k,dta(igout,k), k=1,klev)
-    PRINT *, 'wake-1, updraughts, dqa(igout,k) ', (k,dqa(igout,k), k=1,klev)
-  ENDIF
-
-  ! 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.
-    ! pas besoin d'isos ici
-  END DO
-
-
-#ifdef ISOVERIF
-      ! TODO      
-#endif 
-
-  ! 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)
-#ifdef ISO
-        do ixt=1,ntraciso
-          xtu(ixt,i,k)  =  xte(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)
-        enddo !do ixt=1,ntraciso
-#ifdef ISOVERIF
-        if (iso_eau.gt.0) then
-              call iso_verif_egalite(deltaqw(i,k),deltaxtw(iso_eau,i,k),'wake 723a')
-              call iso_verif_egalite(qu(i,k),xtu(iso_eau,i,k),'wake 723b')
-        endif
-        if (iso_HDO.gt.0) then
-               if (iso_verif_aberrant_enc_choix_nostop(xtu(iso_hdo,i,k),qu(i,k),ridicule,deltalim, &
-                        'wake 723c xtu').eq.1) then
-                     stop
-               endif
-         endif  
-#endif
-#endif
-    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
-
-  IF (prt_level>=10) THEN
-    PRINT *, 'wake-2, ptop_provis(igout), ptop(igout) ', ptop_provis(igout), ptop(igout)
-  ENDIF
-
-
-  ! -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
-
-  IF (prt_level>=10) THEN
-    PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
-  ENDIF
-
-  ! -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.
-        d_deltatw2(i,k) = -deltatw0(i,k)
-#ifdef ISO
-        do ixt=1,ntraciso
-        deltaxtw(ixt,i, k) = 0.
-        d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
-        enddo
-#endif
-      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.+epsim1*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.+epsim1*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)+ &
-        epsim1*(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.
-        d_deltatw2(i,k) = -deltatw0(i,k)
-        d_deltaqw2(i,k) = -deltaqw0(i,k)
-#ifdef ISO
-        do ixt=1,ntraciso
-        deltaxtw(ixt,i, k) = 0.
-        d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
-        enddo
-#endif
-      END IF
-    END DO
-  END DO
-
-  DO i = 1, klon
-    IF (wape(i)<0.) THEN
-      wape(i) = 0.
-      cstar(i) = 0.
-      hw(i) = hwmin
-!jyg<
-!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
-      sigmaw_targ = max(sigmad, sigd_con(i))
-      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
-      sigmaw(i) = sigmaw_targ
-!>jyg
-      fip(i) = 0.
-      gwake(i) = .FALSE.
-    ELSE
-      hw(i) = hw0(i)
-      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
-
-  IF (prt_level>=10) THEN
-    PRINT *, 'wake-4, sigmaw(igout), cstar(igout), wape(igout), ktop(igout) ', &
-                      sigmaw(igout), cstar(igout), wape(igout), ktop(igout)
-  ENDIF
-
-
-  ! 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
-    IF (prt_level>=10) THEN
-      PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', &
-                          isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout)
-    ENDIF
-
-    ! 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 ---
-    !jyg<
-    ! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou 
-    ! les poches sont insuffisantes pour accueillir tout le flux de masse 
-    ! des descentes unsaturees. Nous faisons alors l'hypothese que la 
-    ! convection profonde cree directement de nouvelles poches, sans passer 
-    ! par les thermiques. La nouvelle valeur de wdens est alors imposée.
-
-    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 (prt_level >= 10) THEN
-             print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)', &
-                     omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)
-        ENDIF
-        IF (wdens(i)<=wdens0*1.1) THEN
-          IF (iflag_wk_pop_dyn >= 1) THEN
-             d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
-          ENDIF
-          wdens(i) = wdens0
-        END IF
-      END IF
-    END DO
-
-    DO i = 1, klon
-      IF (wk_adv(i)) THEN
-        gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
-        rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
-!jyg<
-!!        sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
-        sigmaw_targ = min(sigmaw(i), sigmaw_max)
-        d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
-        sigmaw(i) = sigmaw_targ
-!>jyg
-      END IF
-    END DO
-
-    IF (iflag_wk_pop_dyn >= 1) THEN
-!    The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.
-!    Here, it has to be set to zero.
-      death_rate(:) = 0.
-
-      IF (iflag_wk_act ==2) THEN
-      DO i = 1, klon
-        IF (wk_adv(i)) THEN
-          wape1_act(i) = abs(cin(i))
-          wape2_act(i) = 2.*wape1_act(i) + 1.
-          act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
-        ENDIF  ! (wk_adv(i))
-      ENDDO
-      ENDIF  ! (iflag_wk_act ==2)
-
-
-      DO i = 1, klon
-        IF (wk_adv(i)) THEN
-!!          tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.)
-          tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
-          tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
-          drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
-                    (1 + 2*f_shear(i)*(2.*sigmaw(i)-aa0*wdens(i)) - 2.*sigmaw(i))
-!!                    (1 - 2*sigmaw(i)*(1.-f_shear(i)))
-          drdt_pos=max(drdt(i),0.)
-
-!!          d_wdens(i) = ( wgen(i)*(1.+2.*(sigmaw(i)-sigmad)) &
-!!                     - wdens(i)*tau_wk_inv_min &
-!!                     - 2.*gfl(i)*wdens(i)*Cstar(i) )*dtimesub
-          d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub
-          d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min -  &
-                         2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub
-          d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
-
-!!          d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) &
-!!                      + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) &
-!!                      - sigmaw(i)*tau_wk_inv_min )*dtimesub
-          d_sig_gen(i) = wgen(i)*aa0
-          d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min
-!!          d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos
-          d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos
-          d_sigmaw(i) = ( d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + gfl(i)*cstar(i) )*dtimesub
-          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
-        ENDIF
-      ENDDO
-
-      IF (prt_level >= 10) THEN
-        print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
-                       cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
-        print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
-                       wdens(1), awdens(1), act(1), d_awdens(1)
-        print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', &
-                       wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1)
-        print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
-                       d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
-      ENDIF
-    
-    ELSE  !  (iflag_wk_pop_dyn >= 1)
-
-    ! cc nrlmd
-
-      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
-
-    ENDIF   !  (iflag_wk_pop_dyn >= 1)
-
-
-    ! 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..?
-    !JYG 161013 Correction : maintenant omg est dimensionne a klev.
-    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
-      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
-
-    IF (prt_level>=10) THEN
-      PRINT *, 'wake-4.2, omg(igout,k) ', (k,omg(igout,k), k=1,klev)
-      PRINT *, 'wake-4.2, omgtop(igout), ptop(igout), ktop(igout) ', &
-                          omgtop(igout), ptop(igout), ktop(igout)
-    ENDIF
-
-    ! -----------------
-    ! 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
-!!    print *,'omg(igout,k) ', (k,omg(igout,k),k=1,klev)
-    ! 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
-      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-1
-      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
-    DO i = 1, klon
-      IF (wk_adv(i)) THEN
-          dp_omgbw(i, klev) = 0.
-      END IF
-    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
-#ifdef ISO
-          do ixt=1,ntraciso
-            xt1(ixt,i,k) = xte(ixt,i,k) -sigmaw(i)     *deltaxtw(ixt,i,k) ! undisturbed area
-            xt2(ixt,i,k) = xte(ixt,i,k) +(1.-sigmaw(i))*deltaxtw(ixt,i,k) ! wake
-          enddo
-#endif
-        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.
-#ifdef ISO
-        do ixt=1,ntraciso
-          d_xt1(ixt,i,1) = 0.
-          d_xt2(ixt,i,1) = 0.
-          d_dxt(ixt,i,1) = 0.
-        enddo !do ixt=1,ntraciso
-#endif
-      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)
-#ifdef ISO
-        do ixt=1,ntraciso
-          d_xt1(ixt,i,k) = xt1(ixt,i,k-1)-xt1(ixt,i,k)
-          d_xt2(ixt,i,k) = xt2(ixt,i,k-1)-xt2(ixt,i,k)
-          d_dxt(ixt,i,k) = deltaxtw(ixt,i,k-1)-deltaxtw(ixt,i,k)
-        enddo !do ixt=1,ntraciso
-#endif
-        END IF
-      END DO
-    END DO
-
-    DO i = 1, klon
-      IF (wk_adv(i)) THEN
-        omgbdth(i, 1) = 0.
-        omgbdq(i, 1) = 0.
-#ifdef ISO
-        do ixt=1,ntraciso
-         omgbdxt(ixt,i,1) = 0.
-        enddo !do ixt=1,ntraciso
-#endif
-      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))
-#ifdef ISO
-        do ixt=1,ntraciso
-          omgbdxt(ixt,i,k)  = omgb(i,k)*(deltaxtw(ixt,i,k-1) - deltaxtw(ixt,i,k))
-        enddo !do ixt=1,ntraciso
-#ifdef ISOVERIF
-      if (iso_eau.gt.0) then
-        call iso_verif_egalite(deltaqw(i,k-1),deltaxtw(iso_eau,i,k-1),'wake 1460a')
-        call iso_verif_egalite(deltaqw(i,k),deltaxtw(iso_eau,i,k),'wake 1460b')
-        call iso_verif_egalite(omgbdq(i,k),omgbdxt(iso_eau,i,k),'wake 1460c')
-      endif
-#endif
-#endif
-        END IF
-      END DO
-    END DO
-
-    IF (prt_level>=10) THEN
-      PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,klev)
-      PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,klev)
-      PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,klev)
-      PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,klev)
-    ENDIF
-
-    ! -----------------------------------------------------------------
-    DO k = 1, klev-1
-      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=', k, 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=', k, d_deltaqw(i,k)
-#ifdef ISO
-        do ixt=1,ntraciso
-          d_deltaxtw(ixt,i,k) = dtimesub/(Ph(i,k)-Ph(i,k+1))* &
-            (rrd1*omg(i,k  )*sigmaw(i)     *d_xt1(ixt,i,k)- &
-            rrd2*omg(i,k+1)*(1.-sigmaw(i))*d_xt2(ixt,i,k+1)- &
-            (1.-alpha_up(i,k))*omgbdxt(ixt,i,k)- &
-            alpha_up(i,k+1)*omgbdxt(ixt,i,k+1))
-        enddo !do ixt=1,ntraciso
-#ifdef ISOVERIF
-      if (iso_eau.gt.0) then
-        call iso_verif_egalite(d_q1(i,k),d_xt1(iso_eau,i,k),'wake 1502a')
-        call iso_verif_egalite(d_q2(i,k),d_xt2(iso_eau,i,k),'wake 1502b')
-        call iso_verif_egalite(omgbdq(i,k),omgbdxt(iso_eau,i,k),'wake 1502c')
-        call iso_verif_egalite(omgbdq(i,k+1),omgbdxt(iso_eau,i,k+1),'wake 1502d')
-        call iso_verif_egalite(d_deltaqw(i,k),d_deltaxtw(iso_eau,i,k),'wake 1502e')
-      endif
-#endif
-#endif
-
-          ! 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)) &
-                                 -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/ &
-                                 (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)- &
-                                  rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/ &
-                                 (ph(i,k)-ph(i,k+1)) &
-                                 -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i,k+1))/ &
-                                 (ph(i,k)-ph(i,k+1)) )
-#ifdef ISO
-        do ixt=1,ntraciso
-         d_xte(ixt,i,k) =  dtimesub*( &
-             ( rre1(i)*omg(i,k  )*sigmaw(i)     *d_xt1(ixt,i,k) &
-              -rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_xt2(ixt,i,k+1) ) &
-                    /(Ph(i,k)-Ph(i,k+1)) &
-              -sigmaw(i)*(1.-sigmaw(i))*deltaxtw(ixt,i,k) &
-              *(omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1)) &
-              )
-        enddo !do ixt=1,ntraciso
-#endif 
-        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)))
-
-#ifdef ISO
-        do ixt=1,ntraciso
-         d_xte(ixt,i,k) =  dtimesub*( &
-            ( rre1(i)*omg(i,k  )*sigmaw(i)  *d_xt1(ixt,i,k)     &   
-             /(Ph(i,k)-Ph(i,k+1))) &
-                            )
-        enddo !do ixt=1,ntraciso
-#endif
-        END IF
-        ! cc
-      END DO
-    END DO
-    ! ------------------------------------------------------------------
-
-    IF (prt_level>=10) THEN
-      PRINT *, 'wake-4.3, d_deltatw(igout,k) ', (k,d_deltatw(igout,k), k=1,klev)
-      PRINT *, 'wake-4.3, d_deltaqw(igout,k) ', (k,d_deltaqw(igout,k), k=1,klev)
-    ENDIF
-
-    ! Increment state variables
-!jyg<
-    IF (iflag_wk_pop_dyn >= 1) THEN
-      DO k = 1, klev
-        DO i = 1, klon
-          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
-            detr(i,k) = - d_sig_death(i) - d_sig_col(i)      
-            entr(i,k) = d_sig_gen(i)
-          ENDIF
-        ENDDO
-      ENDDO
-      ELSE  ! (iflag_wk_pop_dyn >= 1)
-      DO k = 1, klev
-        DO i = 1, klon
-          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
-            detr(i, k) = 0.
-   
-            entr(i, k) = 0.
-          ENDIF
-        ENDDO
-      ENDDO
-    ENDIF  ! (iflag_wk_pop_dyn >= 1)
-
-    
-
-    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)))
-#ifdef ISO
-        do ixt=1,ntraciso
-          dxtke(ixt,i,k)=(dxtdwn(ixt,i,k)/sigmaw(i) - dxta(ixt,i,k) &
-               /(1.-sigmaw(i)))
-        enddo !do ixt=1,ntraciso
-#ifdef ISOVERIF
-      if (iso_eau.gt.0) then
-        call iso_verif_egalite(dqke(i,k),dxtKE(iso_eau,i,k),'wake 1621a')
-        call iso_verif_egalite(dqdwn(i,k),dxtdwn(iso_eau,i,k),'wake 1621b')
-        call iso_verif_egalite(dqa(i,k),dxta(iso_eau,i,k),'wake 1621c')
-        call iso_verif_egalite(d_deltaqw(i,k),d_deltaxtw(iso_eau,i,k),'wake 1621d')
-      endif
-#endif
-#endif
-          ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
-
-!
-
-          ! cc nrlmd          Prise en compte du taux de mortalité
-          ! cc               Définitions de entr, detr
-!jyg<
-!!            detr(i, k) = 0.
-!!   
-!!            entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
-!!              sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)
-!!
-            entr(i, k) = entr(i,k) + gfl(i)*cstar(i) + &
-                         sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)   
-!>jyg
-            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) - & 
-               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) - &
-                entr(i,k)*deltatw(i,k)/sigmaw(i) - &
-                (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) - &
-                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) - & 
-            entr(i,k)*deltaqw(i,k)/sigmaw(i) - &
-            (death_rate(i)*sigmaw(i)+detr(i,k))*deltaqw(i,k)/(1.-sigmaw(i)))
-#ifdef ISO
-        do ixt=1,ntraciso
-          gg(i)=d_deltaxtw(ixt,i,k)/dtimesub
-          d_deltaxtw(ixt,i,k) = dtimesub*(gg(i) + dxtKE(ixt,i,k) - &
-             entr(i,k)*deltaxtw(ixt,i,k)/sigmaw(i) - &
-             (death_rate(i)*sigmaw(i)+detr(i,k))*deltaxtw(ixt,i,k)/(1.-sigmaw(i)))
-        enddo !do ixt=1,ntraciso
-#ifdef ISOVERIF
-      if (iso_eau.gt.0) then
-        call iso_verif_egalite(dqke(i,k),dxtKE(iso_eau,i,k),'wake 1692a')
-        call iso_verif_egalite(deltaqw(i,k),deltaxtw(iso_eau,i,k),'wake 1692b')
-        call iso_verif_egalite(d_deltaqw(i,k),d_deltaxtw(iso_eau,i,k),'wake 1692c')
-      endif
-#endif
-#endif
-          ! 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
-
-#ifdef ISO
-#ifdef ISOVERIF
-      if (iso_eau.gt.0) then
-             call iso_verif_egalite_vect2D(d_deltaxtw,d_deltaqw, &
-                 'wake 1359',ntraciso,klon,klev)
-      endif      
-#endif         
-#endif  
-
-    ! 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)
-#ifdef ISO
-        do ixt=1,ntraciso
-          d_xte(ixt,i,k)=alpha(i)*d_xte(ixt,i,k)
-          d_deltaxtw(ixt,i,k)=alpha(i)*d_deltaxtw(ixt,i,k)
-        enddo !do ixt=1,ntraciso
-#endif
-        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)
-#ifdef ISO
-        do ixt=1,ntraciso
-          dxtls(ixt,i,k)=dxtls(ixt,i,k)+d_xte(ixt,i,k)
-        enddo !do ixt=1,ntraciso
-#endif
-          ! 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)
-#ifdef ISO
-        do ixt=1,ntraciso
-          d_deltaxtw2(ixt,i,k)=d_deltaxtw2(ixt,i,k)+d_deltaxtw(ixt,i,k)
-        enddo !do ixt=1,ntraciso
-#endif
-          ! cc
-        END IF
-      END DO
-    END DO
-
-
-#ifdef ISO
-#ifdef ISOVERIF
-      if (iso_eau.gt.0) then
-        DO k= 1,klev 
-          DO i = 1,klon 
-              call iso_verif_egalite_choix(dxtls(iso_eau,i,k), &
-                dqls(i,k),'wake 1379',errmax,errmaxrel)
-          enddo ! DO i = 1,klon     
-        enddo ! DO k= 1,klev 
-      endif !if (iso_eau.gt.0) then
-#endif
-#endif
-
-
-    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)
-#ifdef ISO
-        do ixt=1,ntraciso
-          xte(ixt,i,k) = xte0(ixt,i,k) + dxtls(ixt,i,k)
-          deltaxtw(ixt,i,k) = deltaxtw(ixt,i,k)+d_deltaxtw(ixt,i,k)
-        enddo !do ixt=1,ntraciso
-#endif
-        END IF
-      END DO
-    END DO
-!
-    DO i = 1, klon
-      IF (wk_adv(i)) THEN
-        sigmaw(i) = sigmaw(i) + d_sigmaw(i)
-        d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
-      END IF
-    END DO
-
-#ifdef ISO
-#ifdef ISOVERIF
-      write(*,*) 'wake 1859'
-      if (iso_eau.gt.0) then
-        DO k= 1,klev 
-          DO i = 1,klon 
-              call iso_verif_egalite_choix(xte(iso_eau,i,k), &
-                qe(i,k),'wake 1379',errmax,errmaxrel)
-          enddo ! DO i = 1,klon     
-        enddo ! DO k= 1,klev 
-      endif !if (iso_eau.gt.0) then     
-      if (iso_hdo.gt.0) then
-      call iso_verif_aberrant_enc_vect2D( &
-            xte,qe, &
-            'wake 1456, xte apres modifs',ntraciso,klon,klev)
-!      call iso_verif_aberrant_enc_vect2D_ns(
-!     :       deltaxtw,deltaqw,
-!     :       'wake 1518, deltaqw apres modifs',ntraciso,klon,klev)
-      endif
-#endif
-#endif 
-
-!jyg<
-    IF (iflag_wk_pop_dyn >= 1) THEN
-      DO i = 1, klon
-        IF (wk_adv(i)) THEN
-          awdens(i) = awdens(i) + d_awdens(i)
-          wdens(i) = wdens(i) + d_wdens(i)
-          d_awdens2(i) = d_awdens2(i) + d_awdens(i)
-          d_wdens2(i) = d_wdens2(i) + d_wdens(i)
-        END IF
-      END DO
-      DO i = 1, klon
-        IF (wk_adv(i)) THEN
-          wdens_targ = max(wdens(i),wdensmin)
-          d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
-          wdens(i) = wdens_targ
-!
-          wdens_targ = min( max(awdens(i),0.), wdens(i) )
-          d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i)
-          awdens(i) = wdens_targ
-        END IF
-      END DO
-      DO i = 1, klon
-        IF (wk_adv(i)) THEN
-          sigmaw_targ = max(sigmaw(i),sigmad)
-          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
-          sigmaw(i) = sigmaw_targ
-        END IF
-      END DO
-    ENDIF  ! (iflag_wk_pop_dyn >= 1)
-!>jyg
-
-
-    ! 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.
-          d_deltatw2(i,k) = -deltatw0(i,k)
-          d_deltaqw2(i,k) = -deltaqw0(i,k)
-#ifdef ISO
-        do ixt=1,ntraciso
-          deltaxtw(ixt,i,k) = 0.
-          d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
-        enddo !do ixt=1,ntraciso
-#endif
-        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.+epsim1*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.+epsim1*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)+epsim1*(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.
-            d_deltatw2(i,k) = -deltatw0(i,k)
-            d_deltaqw2(i,k) = -deltaqw0(i,k)
-#ifdef ISO
-        do ixt=1,ntraciso
-          deltaxtw(ixt,i,k) = 0.
-          d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
-        enddo !do ixt=1,ntraciso
-#endif
-          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
-!jyg<
-!!          sigmaw(i) = max(sigmad, sigd_con(i))
-          sigmaw_targ = max(sigmad, sigd_con(i))
-          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
-          sigmaw(i) = sigmaw_targ
-!>jyg
-          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
-
-  IF (prt_level>=10) THEN
-    PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', &
-                      sigmaw(igout), cstar(igout), wape(igout), ptop(igout)
-  ENDIF
-
-
-  ! ----------------------------------------------------------
-  ! Determine wake final state; recompute wape, cstar, ktop;
-  ! filter out bad wakes.
-  ! ----------------------------------------------------------
-
-  ! 2.1 - Undisturbed area and Wake integrals
-  ! ---------------------------------------------------------
-
-  DO i = 1, klon
-    ! cc nrlmd       if (wk_adv(i)) then !!! nrlmd
-    IF (ok_qx_qw(i)) THEN
-      ! cc
-      z(i) = 0.
-      sum_thu(i) = 0.
-      sum_tu(i) = 0.
-      sum_qu(i) = 0.
-      sum_thvu(i) = 0.
-      sum_dth(i) = 0.
-      sum_half_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.
-
-      dthmin(i) = -delta_t_min
-    END IF
-  END DO
-  ! Potential temperatures and humidity
-  ! ----------------------------------------------------------
-
-  DO k = 1, klev
-    DO i = 1, klon
-      ! cc nrlmd       IF ( wk_adv(i)) THEN
-      IF (ok_qx_qw(i)) THEN
-        ! cc
-        rho(i, k) = p(i, k)/(rd*te(i,k))
-        IF (k==1) THEN
-          rhoh(i, k) = ph(i, k)/(rd*te(i,k))
-          zhh(i, k) = 0
-        ELSE
-          rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
-          zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
-        END IF
-        the(i, k) = te(i, k)/ppi(i, k)
-        thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
-        tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
-        qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
-        rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
-        dth(i, k) = deltatw(i, k)/ppi(i, k)
-#ifdef ISO
-        do ixt=1,ntraciso
-          xtu(ixt,i,k)  =  xte(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)
-        enddo !do ixt=1,ntraciso
-#endif
-      END IF
-    END DO
-  END DO
-
-#ifdef ISO
-#ifdef ISOVERIF
-      if (iso_hdo.gt.0) then
-      call iso_verif_aberrant_enc_vect2D( &
-              xtu,qu, &
-              'wake 1834, apres modifs',ntraciso,klon,klev)
-      endif   
-#endif
-#endif   
-
-  ! 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.
-      dz_half(i) = 1.
-      sum_thvu(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)
-      sum_dth(i) = 0.
-    END IF
-  END DO
-
-  DO k = 1, klev
-    DO i = 1, klon
-      ! cc nrlmd       IF ( wk_adv(i)) THEN
-      IF (ok_qx_qw(i)) THEN
-        ! cc
-        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
-        dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-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.+epsim1*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)
-!
-          dthmin(i) = min(dthmin(i), dth(i,k))
-        END IF
-        IF (dz_half(i)>0) THEN
-          sum_half_dth(i) = sum_half_dth(i) + dth(i, k)*dz_half(i)
-        END IF
-      END IF
-    END DO
-  END DO
-
-  DO i = 1, klon
-    ! cc nrlmd       IF ( wk_adv(i)) THEN
-    IF (ok_qx_qw(i)) THEN
-      ! cc
-      hw0(i) = z(i)
-    END IF
-  END DO
-
-  ! - WAPE and mean forcing computation
-  ! -------------------------------------------------------------
-
-  ! Means
-
-  DO i = 1, klon
-    ! cc nrlmd       IF ( wk_adv(i)) THEN
-    IF (ok_qx_qw(i)) THEN
-      ! cc
-      av_thu(i) = sum_thu(i)/hw0(i)
-      av_tu(i) = sum_tu(i)/hw0(i)
-      av_qu(i) = sum_qu(i)/hw0(i)
-      av_thvu(i) = sum_thvu(i)/hw0(i)
-      av_dth(i) = sum_dth(i)/hw0(i)
-      av_dq(i) = sum_dq(i)/hw0(i)
-      av_rho(i) = sum_rho(i)/hw0(i)
-      av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
-      av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
-
-      wape2(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(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
-
-  IF (iflag_wk_check_trgl>=1) THEN
-    ! Check triangular shape of dth profile
-    DO i = 1, klon
-      IF (ok_qx_qw(i)) THEN
-        !! print *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i)
-        !! print *,'wake, 2.*sum_dth(i)/(hw0(i)*dthmin(i)) ', &
-        !!                2.*sum_dth(i)/(hw0(i)*dthmin(i))
-        !! print *,'wake, sum_half_dth(i), sum_dth(i) ', &
-        !!                sum_half_dth(i), sum_dth(i)
-        IF ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min) ) THEN
-          wape2(i) = -1.
-          !! print *,'wake, rej 1'
-        ELSE IF (iflag_wk_check_trgl==1.AND.abs(2.*sum_dth(i)/(hw0(i)*dthmin(i)) - 1.) > 0.5) THEN
-          wape2(i) = -1.
-          !! print *,'wake, rej 2'
-        ELSE IF (abs(sum_half_dth(i)) < 0.5*abs(sum_dth(i)) ) THEN
-          wape2(i) = -1.
-          !! print *,'wake, rej 3'
-        END IF
-      END IF
-    END DO
-  END IF
-
-
-  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.
-        d_deltatw2(i,k) = -deltatw0(i,k)
-        d_deltaqw2(i,k) = -deltaqw0(i,k)
-#ifdef ISO
-        do ixt=1,ntraciso
-          deltaxtw(ixt,i,k) = 0.
-          d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
-        enddo !do ixt=1,ntraciso
-#endif
-      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
-!jyg<
-!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
-      sigmaw_targ = max(sigmad, sigd_con(i))
-      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
-      sigmaw(i) = sigmaw_targ
-!>jyg
-        fip(i) = 0.
-        gwake(i) = .FALSE.
-      ELSE
-        IF (prt_level>=10) PRINT *, 'wape2>0'
-        cstar2(i) = stark*sqrt(2.*wape2(i))
-        gwake(i) = .TRUE.
-      END IF
-    END IF
-  END DO
-
-  DO i = 1, klon
-    ! cc nrlmd       IF ( wk_adv(i)) THEN
-    IF (ok_qx_qw(i)) THEN
-      ! cc
-      ktopw(i) = ktop(i)
-    END IF
-  END DO
-
-  DO i = 1, klon
-    ! cc nrlmd       IF ( wk_adv(i)) THEN
-    IF (ok_qx_qw(i)) THEN
-      ! cc
-      IF (ktopw(i)>0 .AND. gwake(i)) THEN
-
-        ! jyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
-        ! cc       heff = 600.
-        ! Utilisation de la hauteur hw
-        ! c       heff = 0.7*hw
-        heff(i) = hw(i)
-
-        fip(i) = 0.5*rho(i, ktopw(i))*cstar2(i)**3*heff(i)*2* &
-          sqrt(sigmaw(i)*wdens(i)*3.14)
-        fip(i) = alpk*fip(i)
-        ! jyg2
-      ELSE
-        fip(i) = 0.
-      END IF
-    END IF
-  END DO
-
-  ! Limitation de sigmaw
-
-  ! cc nrlmd
-  ! DO i=1,klon
-  ! IF (OK_qx_qw(i)) THEN
-  ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
-  ! ENDIF
-  ! ENDDO
-  ! cc
-
-  !jyg<
-  IF (iflag_wk_pop_dyn >= 1) THEN
-    DO i = 1, klon
-      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
-          .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
-    ENDDO
-  ELSE  ! (iflag_wk_pop_dyn >= 1)
-    DO i = 1, klon
-      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
-          .NOT. ok_qx_qw(i)
-    ENDDO
-  ENDIF  ! (iflag_wk_pop_dyn >= 1)
-  !>jyg
-
-  DO k = 1, klev
-    DO i = 1, klon
-!!jyg      IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
-!!jyg          .NOT. ok_qx_qw(i)) THEN
-      IF (kill_wake(i)) THEN
-        ! cc
-        dtls(i, k) = 0.
-        dqls(i, k) = 0.
-        deltatw(i, k) = 0.
-        deltaqw(i, k) = 0.
-        d_deltatw2(i,k) = -deltatw0(i,k)
-        d_deltaqw2(i,k) = -deltaqw0(i,k)
-#ifdef ISO
-        do ixt=1,ntraciso
-          dxtls(ixt,i,k) = 0.
-          deltaxtw(ixt,i,k) = 0.
-          d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
-        enddo !do ixt=1,ntraciso
-#endif
-      END IF  ! (kill_wake(i))
-    END DO
-  END DO
-
-  DO i = 1, klon
-!!jyg    IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
-!!jyg        .NOT. ok_qx_qw(i)) THEN
-      IF (kill_wake(i)) THEN
-      ktopw(i) = 0
-      wape(i) = 0.
-      cstar(i) = 0.
-!!jyg   Outside subroutine "Wake" hw, wdens and sigmaw are zero when there are no wakes
-!!      hw(i) = hwmin                       !jyg
-!!      sigmaw(i) = sigmad                  !jyg
-      hw(i) = 0.                            !jyg
-      fip(i) = 0.
-!!      sigmaw(i) = 0.                        !jyg
-      sigmaw_targ = 0.
-      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
-      sigmaw(i) = sigmaw_targ
-      IF (iflag_wk_pop_dyn >= 1) THEN
-!!        awdens(i) = 0.
-!!        wdens(i) = 0.
-        wdens_targ = 0.
-        d_wdens2(i) = wdens_targ - wdens(i)
-        wdens(i) = wdens_targ
-        wdens_targ = 0.
-        d_awdens2(i) = wdens_targ - awdens(i)
-        awdens(i) = wdens_targ
-      ENDIF  ! (iflag_wk_pop_dyn >= 1)
-    ELSE  ! (kill_wake(i))
-      wape(i) = wape2(i)
-      cstar(i) = cstar2(i)
-    END IF  ! (kill_wake(i))
-    ! c        print*,'wape wape2 ktopw OK_qx_qw =',
-    ! c     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
-  END DO
-
-  IF (prt_level>=10) THEN
-    PRINT *, 'wake-6, wape wape2 ktopw OK_qx_qw =', &
-                      wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout)
-  ENDIF
-
-
-  ! -----------------------------------------------------------------
-  ! 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
-!jyg<
-!!      IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
-      IF (ok_qx_qw(i)) THEN
-!>jyg
-        ! 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)
-#ifdef ISO
-        do ixt=1,ntraciso
-          dxtls(ixt,i, k) = dxtls(ixt,i, k)/dtime
-          d_deltaxtw2(ixt,i, k) = d_deltaxtw2(ixt,i, k)/dtime
-        enddo
-#endif
-      END IF
-    END DO
-  END DO
-!jyg<
-  DO i = 1, klon
-    d_sigmaw2(i) = d_sigmaw2(i)/dtime
-    d_awdens2(i) = d_awdens2(i)/dtime
-    d_wdens2(i) = d_wdens2(i)/dtime
-  ENDDO
-!>jyg
-
-
-
-  RETURN
-END SUBROUTINE wake
-
-SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, &
-    d_deltaqw, sigmaw, d_sigmaw, alpha)
-  ! ------------------------------------------------------
-  ! D\'etermination du coefficient alpha tel que les tendances
-  ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
-  ! a une humidite positive dans la zone (x) et dans la zone (w).
-  ! ------------------------------------------------------
-  IMPLICIT NONE
-
-  ! Input
-  REAL qe(nlon, nl), d_qe(nlon, nl)
-  REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
-  REAL sigmaw(nlon), d_sigmaw(nlon)
-  LOGICAL wk_adv(nlon)
-  INTEGER nl, nlon
-  ! Output
-  REAL alpha(nlon)
-  ! Internal variables
-  REAL zeta(nlon, nl)
-  REAL alpha1(nlon)
-  REAL x, a, b, c, discrim
-  REAL epsilon
-  ! DATA epsilon/1.e-15/
-  INTEGER i,k
-
-  DO k = 1, nl
-    DO i = 1, nlon
-      IF (wk_adv(i)) THEN
-        IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
-          zeta(i, k) = 0.
-        ELSE
-          zeta(i, k) = 1.
-        END IF
-      END IF
-    END DO
-    DO i = 1, nlon
-      IF (wk_adv(i)) THEN
-        x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
-          (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i) * &
-          (deltaqw(i,k)+d_deltaqw(i,k))
-        a = -d_sigmaw(i)*d_deltaqw(i, k)
-        b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
-          deltaqw(i, k)*d_sigmaw(i)
-        c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon
-        discrim = b*b - 4.*a*c
-        ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
-        IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap
-          alpha1(i) = 1.
-        ELSE
-          IF (x>=0.) THEN
-            alpha1(i) = 1.
-          ELSE
-            IF (a>0.) THEN
-              alpha1(i) = 0.9*min( (2.*c)/(-b+sqrt(discrim)),  &
-                                   (-b+sqrt(discrim))/(2.*a) )
-            ELSE IF (a==0.) THEN
-              alpha1(i) = 0.9*(-c/b)
-            ELSE
-              ! print*,'a,b,c discrim',a,b,c discrim
-              alpha1(i) = 0.9*max( (2.*c)/(-b+sqrt(discrim)),  &
-                                   (-b+sqrt(discrim))/(2.*a))
-            END IF
-          END IF
-        END IF
-        alpha(i) = min(alpha(i), alpha1(i))
-      END IF
-    END DO
-  END DO
-
-  RETURN
-END SUBROUTINE wake_vec_modulation
-
-
-
-
Index: LMDZ6/trunk/libf/phylmdiso/wake_ini_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmdiso/wake_ini_mod.F90	(revision 4581)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../phylmd/wake_ini_mod.F90
