Index: /LMDZ6/trunk/libf/phylmd/coef_diff_turb_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/coef_diff_turb_mod.F90	(revision 4229)
+++ /LMDZ6/trunk/libf/phylmd/coef_diff_turb_mod.F90	(revision 4230)
@@ -162,5 +162,5 @@
      
        IF (prt_level > 9) THEN
-          WRITE(lunout,*) 'USTAR = ',yustar
+          WRITE(lunout,*) 'USTAR = ',(yustar(i),i=1,knon)
        ENDIF
          
Index: /LMDZ6/trunk/libf/phylmd/wake.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/wake.F90	(revision 4229)
+++ /LMDZ6/trunk/libf/phylmd/wake.F90	(revision 4230)
@@ -11,6 +11,6 @@
                 dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, &
                 dtke, dqke, omg, dp_deltomg, wkspread, cstar, &
-                d_deltat_gw, &
-                d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2)     ! tendencies
+                d_deltat_gw, &                                                      ! tendencies
+                d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2)             ! tendencies
 
 
@@ -27,5 +27,5 @@
   USE wake_ini_mod , ONLY : wake_ini
   USE wake_ini_mod , ONLY : prt_level,epsim1,RG,RD
-  USE wake_ini_mod , ONLY : stark, wdens_ref, coefgw, alpk, wdens_ref, stark, coefgw, alpk
+  USE wake_ini_mod , ONLY : stark, wdens_ref, coefgw, alpk, pupperbyphs
   USE wake_ini_mod , ONLY : crep_upper, crep_sol, tau_cv, rzero, aa0, flag_wk_check_trgl
   USE wake_ini_mod , ONLY : iflag_wk_act, iflag_wk_check_trgl, iflag_wk_pop_dyn, wdensmin
@@ -119,10 +119,10 @@
   ! 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é
+  ! Tgw    : 1 sur la periode de onde de gravite
+  ! Cgw    : vitesse de propagation de onde de gravite
   ! LL     : distance entre 2 poches
 
   ! -------------------------------------------------------------------------
-  ! Déclaration de variables
+  ! Declaration de variables
   ! -------------------------------------------------------------------------
 
@@ -163,8 +163,9 @@
   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
+  ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields 
+  !                                 computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
+  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_gw
   REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
   REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_awdens2, d_wdens2
@@ -173,5 +174,5 @@
   ! -------------------
 
-  ! Variables à fixer
+  ! Variables a fixer
 
   REAL                                                  :: delta_t_min
@@ -193,9 +194,13 @@
   REAL, DIMENSION(klon)                                 :: f_shear
   REAL, DIMENSION(klon)                                 :: drdt
-  REAL, DIMENSION(klon)                                 :: d_sig_gen, d_sig_death, d_sig_col
+!!  REAL, DIMENSION(klon)                                 :: d_sig_gen, d_sig_death, d_sig_col
   REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
   LOGICAL, DIMENSION (klon)                             :: kill_wake
   REAL                                                  :: drdt_pos
   REAL                                                  :: tau_wk_inv_min
+  ! Some components of the tendencies of state variables  
+  REAL, DIMENSION (klon)                                :: d_sig_gen2, d_sig_death2, d_sig_col2, d_sig_bnd2
+  REAL, DIMENSION (klon)                                :: d_sig_spread2
+  REAL, DIMENSION (klon)                                :: d_dens_gen2, d_dens_death2, d_dens_col2, d_dens_bnd2
 
   ! Variables pour les GW
@@ -219,6 +224,8 @@
   REAL, DIMENSION (klon, klev)                          :: d_deltatw, d_deltaqw
   REAL, DIMENSION (klon, klev)                          :: d_tenv, d_qe
-  REAL, DIMENSION (klon)                                :: d_awdens, d_wdens
-  REAL, DIMENSION (klon)                                :: d_sigmaw, alpha
+  REAL, DIMENSION (klon)                                :: d_awdens, d_wdens, d_sigmaw 
+  REAL, DIMENSION (klon)                                :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd
+  REAL, DIMENSION (klon)                                :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
+  REAL, DIMENSION (klon)                                :: alpha, alpha_tot
   REAL, DIMENSION (klon)                                :: q0_min, q1_min
   LOGICAL, DIMENSION (klon)                             :: wk_adv, ok_qx_qw
@@ -227,6 +234,8 @@
   INTEGER                                               ::isubstep, k, i, igout
 
+  REAL                                                  :: sigmaw_targ
   REAL                                                  :: wdens_targ
-  REAL                                                  :: sigmaw_targ
+  REAL                                                  :: d_sigmaw_targ
+  REAL                                                  :: d_wdens_targ
 
   REAL, DIMENSION (klon)                                :: sum_thu, sum_tu, sum_qu, sum_thvu
@@ -289,7 +298,7 @@
   ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
 
-  ! coefgw : Coefficient pour les ondes de gravité
+  ! coefgw : Coefficient pour les ondes de gravite
   ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
-  ! wdens : Densité surfacique de poche froide
+  ! wdens : Densite surfacique de poche froide
   ! -------------------------------------------------------------------------
 
@@ -380,33 +389,4 @@
    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
@@ -417,4 +397,43 @@
 !!  ELSE
   ENDIF  ! (iflag_wk_pop_dyn >= 1)
+
+
+  ! 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
+      d_dens_gen2(i)   = 0.
+      d_dens_death2(i) = 0.
+      d_dens_col2(i)   = 0.
+      d_awdens2(i) = 0.
+      wdens_targ = max(wdens(i),wdensmin)
+      d_dens_bnd2(i) = wdens_targ - wdens(i)
+      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)
+    d_sig_gen2(i)   = 0.
+    d_sig_death2(i) = 0.
+    d_sig_col2(i)   = 0.
+    d_sig_spread2(i)= 0.
+    sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
+    d_sig_bnd2(i) = sigmaw_targ - sigmaw(i)
+    d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
+    sigmaw(i) = sigmaw_targ
+!>jyg
+  END DO
 
   wape(:) = 0.
@@ -572,17 +591,7 @@
   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
@@ -657,16 +666,13 @@
   ! -5/ Determination de ktop et kupper
 
+  CALL pkupper (klon, klev, ptop, ph, pupper, 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
+    END DO
+  END DO
+  print*, 'ptop, pupper, ktop, kupper', ptop, pupper, ktop, kupper
+ 
 
 
@@ -808,4 +814,5 @@
 !!      sigmaw(i) = amax1(sigmad, sigd_con(i))
       sigmaw_targ = max(sigmad, sigd_con(i))
+      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
       d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
       sigmaw(i) = sigmaw_targ
@@ -840,4 +847,5 @@
     ok_qx_qw(i) = q0_min(i) >= 0.
     alpha(i) = 1.
+    alpha_tot(i) = 1.
   END DO
 
@@ -855,7 +863,12 @@
   dtimesub = dtime/nsub
 
+
+  
   ! ------------------------------------------------------------
   DO isubstep = 1, nsub
     ! ------------------------------------------------------------
+  CALL pkupper (klon, klev, ptop, ph, pupper, kupper)
+ 
+    print*, 'ptop, pupper, ktop, kupper', ptop, pupper, ktop, kupper
 
     ! wk_adv is the logical flag enabling wake evolution in the time advance
@@ -867,9 +880,10 @@
       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
+    ! negatif de ktop a kupper --------
+    ! cc           On calcule pour cela une densite wdens0 pour laquelle on
     ! aurait un entrainement nul ---
     !jyg<
@@ -878,5 +892,5 @@
     ! 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.
+    ! par les thermiques. La nouvelle valeur de wdens est alors imposee.
 
     DO i = 1, klon
@@ -894,4 +908,5 @@
         IF (wdens(i)<=wdens0*1.1) THEN
           IF (iflag_wk_pop_dyn >= 1) THEN
+             d_dens_bnd2(i) = d_dens_bnd2(i) + wdens0 - wdens(i)
              d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
           ENDIF
@@ -908,4 +923,5 @@
 !!        sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
         sigmaw_targ = min(sigmaw(i), sigmaw_max)
+        d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
         d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
         sigmaw(i) = sigmaw_targ
@@ -914,68 +930,37 @@
     END DO
 
-    IF (iflag_wk_pop_dyn >= 1) THEN
+    IF (iflag_wk_pop_dyn == 1) THEN
+  
+     CALL wake_popdyn_1 (klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
+                  dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
+                  d_awdens, d_wdens, d_sigmaw, &
+                  iflag_wk_act, wk_adv, cin, wape, &
+                  drdt, &
+                  d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
+                  d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
+                  d_wdens_targ, d_sigmaw_targ)
+                      
 !    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
+    
+    ELSEIF (iflag_wk_pop_dyn == 2) THEN
+     CALL wake_popdyn_2 (klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
+                  dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
+                  d_awdens, d_wdens, d_sigmaw, &
+                  iflag_wk_act, wk_adv, cin, wape, &
+                  drdt, &
+                  d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
+                  d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
+                  d_wdens_targ, d_sigmaw_targ)
+     death_rate(:) = 0.
+   
+    ELSEIF (iflag_wk_pop_dyn == 0.) THEN
+    
+    ! cc nrlmd
+
       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
+          ! cc nrlmd          Introduction du taux de mortalite des poches et
           ! test sur sigmaw_max=0.4
           ! cc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
@@ -1008,5 +993,5 @@
     ! 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..?
+    ! IM 060208 au niveau k=1...
     !JYG 161013 Correction : maintenant omg est dimensionne a klev.
     DO k = 1, klev
@@ -1090,5 +1075,5 @@
 
     ! c      DO i=1,klon
-    ! c        print*,'Pente entre 0 et kupper (référence)'
+    ! c        print*,'Pente entre 0 et kupper (reference)'
     ! c     $   	,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
     ! c        print*,'Pente entre ktop et kupper'
@@ -1166,4 +1151,5 @@
         IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
           dth(i, k) = deltatw(i, k)/ppi(i, k)
+  print *, 'VVVVwake k, the(i,k), dth(i,k), sigmaw(i) ', k, the(i,k), dth(i,k), sigmaw(i)
           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
@@ -1214,9 +1200,10 @@
     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)
+!!    IF (prt_level>=10) THEN
+    IF (prt_level>=10 .and. wk_adv(igout)) THEN
+      PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,kupper(igout))
+      PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,kupper(igout))
+      PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,kupper(igout))
+      PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,kupper(igout))
     ENDIF
 
@@ -1310,5 +1297,5 @@
 
 
-          ! Coefficient de répartition
+          ! Coefficient de repartition
 
           crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
@@ -1340,6 +1327,6 @@
 !
 
-          ! cc nrlmd          Prise en compte du taux de mortalité
-          ! cc               Définitions de entr, detr
+          ! cc nrlmd          Prise en compte du taux de mortalite
+          ! cc               Definitions de entr, detr
 !jyg<
 !!            detr(i, k) = 0.
@@ -1358,5 +1345,5 @@
 
 
-          ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU
+          ! ajout d'un effet onde de gravite -Tgw(k)*deltatw(k) 03/02/06 YU
           ! Jingmei
 
@@ -1414,4 +1401,11 @@
     CALL wake_vec_modulation(klon, klev, wk_adv, epsilon_loc, qe, d_qe, deltaqw, &
       d_deltaqw, sigmaw, d_sigmaw, alpha)
+    !
+    ! Alpha_tot = Product of all the alpha's
+    DO i = 1, klon
+      IF (wk_adv(i)) THEN
+        alpha_tot(i) = alpha_tot(i)*alpha(i)    
+      END IF
+    END DO
 
     ! cc nrlmd
@@ -1478,8 +1472,21 @@
       DO i = 1, klon
         IF (wk_adv(i)) THEN
+          d_sig_gen2(i)   = d_sig_gen2(i)   + d_sig_gen(i)
+          d_sig_death2(i) = d_sig_death2(i) + d_sig_death(i)
+          d_sig_col2(i)   = d_sig_col2(i)   + d_sig_col(i)
+          d_sig_spread2(i)= d_sig_spread2(i)+ d_sig_spread(i)
+          d_sig_bnd2(i)   = d_sig_bnd2(i)   + d_sig_bnd(i)
+        END IF
+      END DO
+      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)
+          d_dens_gen2(i)   = d_dens_gen2(i)   + d_dens_gen(i)
+          d_dens_death2(i) = d_dens_death2(i) + d_dens_death(i)
+          d_dens_col2(i)   = d_dens_col2(i)   + d_dens_col(i)
+          d_dens_bnd2(i)   = d_dens_bnd2(i)   + d_dens_bnd(i)
         END IF
       END DO
@@ -1487,4 +1494,5 @@
         IF (wk_adv(i)) THEN
           wdens_targ = max(wdens(i),wdensmin)
+          d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
           d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
           wdens(i) = wdens_targ
@@ -1498,4 +1506,5 @@
         IF (wk_adv(i)) THEN
           sigmaw_targ = max(sigmaw(i),sigmad)
+          d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
           d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
           sigmaw(i) = sigmaw_targ
@@ -1744,4 +1753,5 @@
 !!          sigmaw(i) = max(sigmad, sigd_con(i))
           sigmaw_targ = max(sigmad, sigd_con(i))
+          d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
           d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
           sigmaw(i) = sigmaw_targ
@@ -1961,4 +1971,5 @@
 !!      sigmaw(i) = amax1(sigmad, sigd_con(i))
       sigmaw_targ = max(sigmad, sigd_con(i))
+      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
       d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
       sigmaw(i) = sigmaw_targ
@@ -2058,5 +2069,7 @@
 !!      sigmaw(i) = 0.                        !jyg
       sigmaw_targ = 0.
-      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
+      d_sig_bnd2(i) = d_sig_bnd2(i) + sigmaw_targ - sigmaw(i)
+!!      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
+      d_sigmaw2(i) = sigmaw_targ - sigmaw_in(i)      ! _in = correction jyg 20220124
       sigmaw(i) = sigmaw_targ
       IF (iflag_wk_pop_dyn >= 1) THEN
@@ -2064,4 +2077,5 @@
 !!        wdens(i) = 0.
         wdens_targ = 0.
+        d_dens_bnd2(i) = d_dens_bnd2(i) + wdens_targ - wdens(i)
         d_wdens2(i) = wdens_targ - wdens(i)
         wdens(i) = wdens_targ
@@ -2106,10 +2120,25 @@
     END DO
   END DO
-
-  DO i = 1, klon
+!jyg<
+  IF (iflag_wk_pop_dyn >= 1) THEN
+  DO i = 1, klon
+      IF (ok_qx_qw(i)) THEN
+    d_sig_gen2(i) = d_sig_gen2(i)/dtime
+    d_sig_death2(i) = d_sig_death2(i)/dtime
+    d_sig_col2(i) = d_sig_col2(i)/dtime
+    d_sig_spread2(i) = d_sig_spread2(i)/dtime
+    d_sig_bnd2(i) = d_sig_bnd2(i)/dtime
     d_sigmaw2(i) = d_sigmaw2(i)/dtime
+!
+    d_dens_gen2(i) = d_dens_gen2(i)/dtime
+    d_dens_death2(i) = d_dens_death2(i)/dtime
+    d_dens_col2(i) = d_dens_col2(i)/dtime
+    d_dens_bnd2(i) = d_dens_bnd2(i)/dtime
     d_awdens2(i) = d_awdens2(i)/dtime
     d_wdens2(i) = d_wdens2(i)/dtime
+      ENDIF
   ENDDO
+  ENDIF
+!>jyg
 
  RETURN
@@ -2161,5 +2190,5 @@
         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
+        IF (a+b>=0.) THEN !! Condition suffisante pour la positivite de ovap
           alpha1(i) = 1.
         ELSE
@@ -2186,2 +2215,380 @@
   RETURN
 END SUBROUTINE wake_vec_modulation
+
+
+
+SUBROUTINE pkupper (klon, klev, ptop, ph, pupper, kupper)
+
+USE wake_ini_mod , ONLY : pupperbyphs
+IMPLICIT NONE
+
+INTEGER,  INTENT(IN)                              :: klon,klev
+REAL,     INTENT(IN),   DIMENSION (klon,klev+1)   :: ph
+REAL,     INTENT(IN),   DIMENSION (klon)          :: ptop
+REAL,     INTENT(OUT),  DIMENSION (klon)          :: pupper
+INTEGER,  INTENT(OUT),  DIMENSION (klon)          :: kupper
+INTEGER                                           :: i,k
+
+
+ kupper = 0
+ 
+IF (pupperbyphs<1.) THEN
+ ! 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) = pupperbyphs*ph(i, 1)
+    pupper(i) = max(pupper(i), 45000.)
+    ! cc        Pupper(i) = 60000.
+  END DO
+
+ELSE
+  
+  DO i=1, klon
+     ! pupper(i) = pupperbyphs*ptop(i)+(1.-pupperbyphs)*ph(i, 1) 
+     pupper(i) = min( pupperbyphs*ptop(i)+(1.-pupperbyphs)*ph(i, 1) , ptop(i)-5000.) 
+  END DO
+END IF
+ 
+  ! -5/ Determination de kupper
+
+  DO k = klev, 1, -1
+    DO i = 1, klon
+      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 
+    RETURN
+END SUBROUTINE pkupper
+
+
+SUBROUTINE wake_popdyn_1(klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
+                  dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
+                  d_awdens, d_wdens, d_sigmaw, &
+                  iflag_wk_act, wk_adv, cin, wape, &
+                  drdt, &
+                  d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
+                  d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
+                  d_wdens_targ, d_sigmaw_targ)
+                
+
+  USE wake_ini_mod , ONLY : wake_ini
+  USE wake_ini_mod , ONLY : prt_level,RG
+  USE wake_ini_mod , ONLY : stark, wdens_ref
+  USE wake_ini_mod , ONLY : tau_cv, rzero, aa0
+  USE wake_ini_mod , ONLY : iflag_wk_pop_dyn, wdensmin
+  USE wake_ini_mod , ONLY : sigmad, cstart, sigmaw_max
+  
+IMPLICIT NONE
+
+  INTEGER, INTENT(IN)                                   :: klon,klev
+  LOGICAL, DIMENSION (klon),        INTENT(IN)          :: wk_adv
+  REAL,                             INTENT(IN)          :: dtime
+  REAL,                             INTENT(IN)          :: dtimesub
+  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen
+  REAL, DIMENSION (klon),           INTENT(IN)          :: wdens
+  REAL, DIMENSION (klon),           INTENT(IN)          :: awdens
+  REAL, DIMENSION (klon),           INTENT(IN)          :: sigmaw
+  REAL, DIMENSION (klon),           INTENT(IN)          :: gfl, cstar
+  REAL, DIMENSION (klon),           INTENT(IN)          :: cin, wape
+  REAL, DIMENSION (klon),           INTENT(IN)          :: rad_wk
+  REAL, DIMENSION (klon),           INTENT(IN)          :: f_shear
+  INTEGER,                          INTENT(IN)          :: iflag_wk_act
+
+  
+  !
+ 
+  ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields 
+  !                                 computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
+  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw, d_awdens, d_wdens
+  REAL, DIMENSION (klon),           INTENT(OUT)         :: drdt
+  ! Some components of the tendencies of state variables  
+  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sig_gen, d_sig_death, d_sig_col, d_sig_bnd
+  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sig_spread
+  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
+  REAL,                             INTENT(OUT)         :: d_wdens_targ, d_sigmaw_targ
+ 
+  
+  REAL                                                  :: delta_t_min
+  INTEGER                                               :: nsub
+  INTEGER                                               :: i, k
+  REAL                                                  :: wdens0
+  ! IM 080208
+  LOGICAL, DIMENSION (klon)                             :: gwake
+  
+   ! Variables liees a la dynamique de population
+  REAL, DIMENSION(klon)                                 :: act
+  REAL, DIMENSION(klon)                                 :: tau_wk_inv
+  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
+  LOGICAL, DIMENSION (klon)                             :: kill_wake
+  REAL                                                  :: drdt_pos
+  REAL                                                  :: tau_wk_inv_min
+  
+     
+
+      IF (iflag_wk_act == 0) THEN
+        act(:) = 0.
+      ELSEIF (iflag_wk_act == 1) THEN
+        act(:) = 1.
+      ELSEIF (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
+       ! print*, 'XXX wk_adv(i)', wk_adv(i)
+        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
+!jyg+mlt<
+          d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub
+          d_dens_gen(i) = wgen(i)
+          d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min 
+          d_dens_col(i) =  -2.*wdens(i)*gfl(i)*drdt_pos 
+          d_dens_gen(i) =  d_dens_gen(i)*dtimesub
+          d_dens_death(i) = d_dens_death(i)*dtimesub
+          d_dens_col(i) =  d_dens_col(i)*dtimesub 
+
+          d_wdens(i) = d_dens_gen(i)+d_dens_death(i)+d_dens_col(i)
+!!          d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min -  &
+!!                         2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub
+!>jyg+mlt
+!
+!jyg<
+          d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
+!!          d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
+          d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
+          d_wdens(i) = d_wdens_targ
+!!          d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
+!>jyg
+
+!jyg+mlt<
+!!          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
+          print*, 'XXX sigmaw(i), awdens(i), wdens(i), tau_wk_inv_min', &
+                   sigmaw(i), awdens(i), wdens(i), tau_wk_inv_min
+          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_sig_spread(i) = gfl(i)*cstar(i)
+          d_sig_gen(i) =  d_sig_gen(i)*dtimesub
+          d_sig_death(i) = d_sig_death(i)*dtimesub
+          d_sig_col(i) =  d_sig_col(i)*dtimesub 
+          d_sig_spread(i) =  d_sig_spread(i)*dtimesub 
+          d_sigmaw(i) =  d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
+!>jyg+mlt
+!
+!jyg<
+          d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
+!!          d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
+!!          d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
+          d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
+          d_sigmaw(i) = d_sigmaw_targ
+!!          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
+!>jyg
+        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
+    
+    RETURN 
+    END SUBROUTINE wake_popdyn_1
+    
+    SUBROUTINE wake_popdyn_2(klon, klev, dtime, cstar, tau_wk_inv, wgen, wdens, awdens, sigmaw, &
+                  dtimesub, gfl, rad_wk, f_shear, drdt_pos, &
+                  d_awdens, d_wdens, d_sigmaw, &
+                  iflag_wk_act, wk_adv, cin, wape, &
+                  drdt, &
+                  d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd, &
+                  d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd, &
+                  d_wdens_targ, d_sigmaw_targ)
+                
+
+  USE wake_ini_mod , ONLY : wake_ini
+  USE wake_ini_mod , ONLY : prt_level,RG
+  USE wake_ini_mod , ONLY : stark, wdens_ref
+  USE wake_ini_mod , ONLY : tau_cv, rzero, aa0
+  USE wake_ini_mod , ONLY : iflag_wk_pop_dyn, wdensmin
+  USE wake_ini_mod , ONLY : sigmad, cstart, sigmaw_max
+  
+IMPLICIT NONE
+
+  INTEGER, INTENT(IN)                                   :: klon,klev
+  LOGICAL, DIMENSION (klon),        INTENT(IN)          :: wk_adv
+  REAL,                             INTENT(IN)          :: dtime
+  REAL,                             INTENT(IN)          :: dtimesub
+  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen
+  REAL, DIMENSION (klon),           INTENT(IN)          :: wdens
+  REAL, DIMENSION (klon),           INTENT(IN)          :: awdens
+  REAL, DIMENSION (klon),           INTENT(IN)          :: sigmaw
+  REAL, DIMENSION (klon),           INTENT(IN)          :: gfl, cstar
+  REAL, DIMENSION (klon),           INTENT(IN)          :: cin, wape
+  REAL, DIMENSION (klon),           INTENT(IN)          :: rad_wk
+  REAL, DIMENSION (klon),           INTENT(IN)          :: f_shear
+  INTEGER,                          INTENT(IN)          :: iflag_wk_act
+
+  
+  !
+ 
+  ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields 
+  !                                 computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)
+  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw, d_awdens, d_wdens
+  REAL, DIMENSION (klon),           INTENT(OUT)         :: drdt
+  ! Some components of the tendencies of state variables  
+  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sig_gen, d_sig_death, d_sig_col, d_sig_bnd
+  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sig_spread
+  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_dens_gen, d_dens_death, d_dens_col, d_dens_bnd
+  REAL,                             INTENT(OUT)         :: d_wdens_targ, d_sigmaw_targ
+ 
+  
+  REAL                                                  :: delta_t_min
+  INTEGER                                               :: nsub
+  INTEGER                                               :: i, k
+  REAL                                                  :: wdens0
+  ! IM 080208
+  LOGICAL, DIMENSION (klon)                             :: gwake
+  
+   ! Variables liees a la dynamique de population
+  REAL, DIMENSION(klon)                                 :: act
+  REAL, DIMENSION(klon)                                 :: tau_wk_inv
+  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
+  LOGICAL, DIMENSION (klon)                             :: kill_wake
+  REAL                                                  :: drdt_pos
+  REAL                                                  :: tau_wk_inv_min
+  
+     
+
+      IF (iflag_wk_act == 0) THEN
+        act(:) = 0.
+      ELSEIF (iflag_wk_act == 1) THEN
+        act(:) = 1.
+      ELSEIF (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
+       ! print*, 'XXX wk_adv(i)', wk_adv(i)
+        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
+!jyg+mlt<
+          d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub
+          d_dens_gen(i) = wgen(i)
+          d_dens_death(i) = - (wdens(i)-awdens(i))*tau_wk_inv_min 
+          d_dens_col(i) =  -2.*wdens(i)*gfl(i)*drdt_pos 
+          d_dens_gen(i) =  d_dens_gen(i)*dtimesub
+          d_dens_death(i) = d_dens_death(i)*dtimesub
+          d_dens_col(i) =  d_dens_col(i)*dtimesub 
+
+          d_wdens(i) = d_dens_gen(i)+d_dens_death(i)+d_dens_col(i)
+!!          d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min -  &
+!!                         2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub
+!>jyg+mlt
+!
+!jyg<
+          d_wdens_targ = max(d_wdens(i), wdensmin-wdens(i))
+!!          d_dens_bnd(i) = d_dens_bnd(i) + d_wdens_targ - d_wdens(i)
+          d_dens_bnd(i) = d_wdens_targ - d_wdens(i)
+          d_wdens(i) = d_wdens_targ
+!!          d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
+!>jyg
+
+!jyg+mlt<
+!!          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
+          print*, 'XXX sigmaw(i), awdens(i), wdens(i), tau_wk_inv_min', &
+                   sigmaw(i), awdens(i), wdens(i), tau_wk_inv_min
+          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_sig_spread(i) = gfl(i)*cstar(i)
+          d_sig_gen(i) =  d_sig_gen(i)*dtimesub
+          d_sig_death(i) = d_sig_death(i)*dtimesub
+          d_sig_col(i) =  d_sig_col(i)*dtimesub 
+          d_sig_spread(i) =  d_sig_spread(i)*dtimesub 
+          d_sigmaw(i) =  d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + d_sig_spread(i)
+!>jyg+mlt
+!
+!jyg<
+          d_sigmaw_targ = max(d_sigmaw(i), sigmad-sigmaw(i))
+!!          d_sig_bnd(i) = d_sig_bnd(i) + d_sigmaw_targ - d_sigmaw(i)
+!!          d_sig_bnd_provis(i) = d_sigmaw_targ - d_sigmaw(i)
+          d_sig_bnd(i) = d_sigmaw_targ - d_sigmaw(i)
+          d_sigmaw(i) = d_sigmaw_targ
+!!          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
+!>jyg
+        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
+    
+    RETURN 
+    END SUBROUTINE wake_popdyn_2  
+  
+  
Index: /LMDZ6/trunk/libf/phylmd/wake_ini_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/wake_ini_mod.F90	(revision 4229)
+++ /LMDZ6/trunk/libf/phylmd/wake_ini_mod.F90	(revision 4230)
@@ -17,16 +17,16 @@
 
   ! -------------------------------------------------------------------------
-  ! Déclaration de variables
-  ! -------------------------------------------------------------------------
-
-  ! Variables à fixer
+  ! Declaration de variables
+  ! -------------------------------------------------------------------------
+
+  ! Variables a fixer
 !jyg<
 !!  REAL, SAVE                                            :: stark, wdens_ref, coefgw, alpk
   INTEGER, SAVE, PROTECTED                                         :: prt_level
   REAL, SAVE, PROTECTED, DIMENSION(2)                              :: wdens_ref
-  REAL, SAVE, PROTECTED                                            :: stark, coefgw, alpk
+  REAL, SAVE, PROTECTED                                            :: stark, coefgw, alpk, pupperbyphs
 !>jyg
   REAL, SAVE, PROTECTED                                            :: crep_upper, crep_sol  
-  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
+  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, pupperbyphs, crep_upper, crep_sol)
 
   REAL, SAVE, PROTECTED                                            :: tau_cv
@@ -111,7 +111,7 @@
   ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
 
-  ! coefgw : Coefficient pour les ondes de gravité
+  ! coefgw : Coefficient pour les ondes de gravite
   ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
-  ! wdens : Densité surfacique de poche froide
+  ! wdens : Densite surfacique de poche froide
   ! -------------------------------------------------------------------------
 
@@ -149,4 +149,9 @@
   alpk=0.25
   CALL getin_p('alpk',alpk)
+  
+  pupperbyphs=0.6
+  CALL getin_p('pupperbyphs',pupperbyphs)
+
+
 !jyg<
 !!  wdens_ref=8.E-12
@@ -184,4 +189,5 @@
   WRITE(*,*) 'stark=', stark
   WRITE(*,*) 'alpk=', alpk
+  WRITE(*,*) 'pupperbyphs=', pupperbyphs
 !jyg<
 !!  WRITE(*,*) 'wdens_ref=', wdens_ref
