MODULE lmdz_wake

! $Id: lmdz_wake.F90 4908 2024-04-15 17:30:55Z jyg $

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
  USE yomcst_mod_h
  USE cvthermo_mod_h
  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

  ! 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
