Index: /LMDZ5/trunk/libf/phylmd/YOETHF.h
===================================================================
--- /LMDZ5/trunk/libf/phylmd/YOETHF.h	(revision 2634)
+++ /LMDZ5/trunk/libf/phylmd/YOETHF.h	(revision 2635)
@@ -24,5 +24,5 @@
      &               RALFDCP,RTWAT,RTBER,RTBERCU,                       &
      &               RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,&
-     &	  	     RKOOP2
+     &               RKOOP2
 
 !$OMP THREADPRIVATE(/YOETHF/)
Index: /LMDZ5/trunk/libf/phylmd/add_wake_tend.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/add_wake_tend.F90	(revision 2635)
+++ /LMDZ5/trunk/libf/phylmd/add_wake_tend.F90	(revision 2635)
@@ -0,0 +1,64 @@
+SUBROUTINE add_wake_tend(zddeltat, zddeltaq, zds, zddens, zoccur, text, abortphy)
+!===================================================================
+! Ajoute les tendances liées aux diverses parametrisations physiques aux
+! variables d'etat des poches froides.
+!===================================================================
+!======================================================================
+! Declarations
+!======================================================================
+
+USE dimphy, ONLY: klon, klev
+USE phys_state_var_mod, ONLY: wake_deltat, wake_deltaq, wake_s, wake_dens
+
+USE print_control_mod, ONLY: prt_level
+IMPLICIT none
+
+! Arguments :
+!------------
+  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: zddeltat, zddeltaq
+  REAL, DIMENSION(klon),         INTENT (IN)         :: zds, zddens
+  INTEGER, DIMENSION(klon),      INTENT (IN)         :: zoccur
+  CHARACTER*(*),                 INTENT (IN)         :: text
+  INTEGER,                       INTENT (IN)         :: abortphy
+
+! Local :
+!--------
+
+INTEGER                                              :: i, l
+
+
+
+     IF (prt_level >= 5) then
+        write (*,*) "In add_wake_tend, after ",text
+        call flush
+     end if
+
+     IF (abortphy==1) RETURN ! on n ajoute pas les tendance si le modele
+                              ! a deja plante.
+
+!======================================================================
+!    Add tendencies to wake state variables
+!======================================================================
+         DO l = 1, klev
+           DO i = 1, klon
+             IF (zoccur(i) .GE. 1) THEN
+               wake_deltat(i, l) = wake_deltat(i, l) + zddeltat(i,l)
+               wake_deltaq(i, l) = wake_deltaq(i, l) + zddeltaq(i,l)
+             ELSE
+               wake_deltat(i, l) = 0.
+               wake_deltaq(i, l) = 0.
+             ENDIF   ! (zoccur(i) .GE. 1)
+           END DO
+         END DO
+         DO i = 1, klon
+           IF (zoccur(i) .GE. 1) THEN
+             wake_s(i)    = wake_s(i)    + zds(i)
+             wake_dens(i) = wake_dens(i) + zddens(i)
+           ELSE
+             wake_s(i)    = 0.
+             wake_dens(i) = 0.
+           ENDIF   ! (zoccur(i) .GE. 1)
+         END DO
+
+RETURN
+END
Index: /LMDZ5/trunk/libf/phylmd/calwake.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/calwake.F90	(revision 2634)
+++ /LMDZ5/trunk/libf/phylmd/calwake.F90	(revision 2635)
@@ -2,11 +2,17 @@
 ! $Id$
 
-SUBROUTINE calwake(paprs, pplay, dtime, t, q, omgb, dt_dwn, dq_dwn, m_dwn, &
-    m_up, dt_a, dq_a, sigd, wdt_pbl, wdq_pbl, udt_pbl, udq_pbl, wake_deltat, &
-    wake_deltaq, wake_dth, wake_h, wake_s, wake_dens, wake_pe, wake_fip, &
-    wake_gfl, dt_wake, dq_wake, wake_k, undi_t, undi_q, wake_omgbdth, &
-    wake_dp_omgb, wake_dtke, wake_dqke, wake_dtpbl, wake_dqpbl, wake_omg, &
-    wake_dp_deltomg, wake_spread, wake_cstar, wake_d_deltat_gw, wake_ddeltat, &
-    wake_ddeltaq)
+SUBROUTINE calwake(iflag_wake_tend, paprs, pplay, dtime, &
+    t, q, omgb, &
+    dt_dwn, dq_dwn, m_dwn, m_up, dt_a, dq_a, &
+    sigd, &
+    wake_deltat, wake_deltaq, wake_s, wake_dens, &
+    wake_dth, wake_h, &
+    wake_pe, wake_fip, wake_gfl, &
+    dt_wake, dq_wake, wake_k, t_x, q_x, wake_omgbdth, &
+    wake_dp_omgb, &
+    wake_dtke, wake_dqke, &
+    wake_omg, wake_dp_deltomg, &
+    wake_spread, wake_cstar, wake_d_deltat_gw, &
+    wake_ddeltat, wake_ddeltaq, wake_ds, wake_ddens)
   ! **************************************************************
   ! *
@@ -26,57 +32,67 @@
   ! Arguments
   ! ----------
-
-  INTEGER i, l, ktopw(klon)
-  REAL dtime
-
-  REAL paprs(klon, klev+1), pplay(klon, klev)
-  REAL t(klon, klev), q(klon, klev), omgb(klon, klev)
-  REAL dt_dwn(klon, klev), dq_dwn(klon, klev), m_dwn(klon, klev)
-  REAL m_up(klon, klev)
-  REAL dt_a(klon, klev), dq_a(klon, klev)
-  REAL wdt_pbl(klon, klev), wdq_pbl(klon, klev)
-  REAL udt_pbl(klon, klev), udq_pbl(klon, klev)
-  REAL wake_deltat(klon, klev), wake_deltaq(klon, klev)
-  REAL dt_wake(klon, klev), dq_wake(klon, klev)
-  REAL wake_d_deltat_gw(klon, klev)
-  REAL wake_h(klon), wake_s(klon)
-  REAL wake_dth(klon, klev)
-  REAL wake_pe(klon), wake_fip(klon), wake_gfl(klon)
-  REAL undi_t(klon, klev), undi_q(klon, klev)
-  REAL wake_omgbdth(klon, klev), wake_dp_omgb(klon, klev)
-  REAL wake_dtke(klon, klev), wake_dqke(klon, klev)
-  REAL wake_dtpbl(klon, klev), wake_dqpbl(klon, klev)
-  REAL wake_omg(klon, klev), wake_dp_deltomg(klon, klev)
-  REAL wake_spread(klon, klev), wake_cstar(klon)
-  REAL wake_ddeltat(klon, klev), wake_ddeltaq(klon, klev)
-  REAL d_deltatw(klon, klev), d_deltaqw(klon, klev)
-  INTEGER wake_k(klon)
-  REAL sigd(klon)
-  REAL wake_dens(klon)
+  ! Input
+  ! ----
+  INTEGER,                       INTENT (IN)         :: iflag_wake_tend
+  REAL,                          INTENT (IN)         :: dtime
+  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: pplay
+  REAL, DIMENSION(klon, klev+1), INTENT (IN)         :: paprs
+  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: t, q, omgb
+  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: dt_dwn, dq_dwn
+  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: m_up, m_dwn
+  REAL, DIMENSION(klon, klev),   INTENT (IN)         :: dt_a, dq_a
+  REAL, DIMENSION(klon),         INTENT (IN)         :: sigd
+  ! Input/Output
+  ! ------------
+  REAL, DIMENSION(klon, klev),   INTENT (INOUT)      :: wake_deltat, wake_deltaq
+  REAL, DIMENSION(klon),         INTENT (INOUT)      :: wake_s
+  REAL, DIMENSION(klon),         INTENT (INOUT)      :: wake_dens
+  ! Output
+  ! ------
+  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: dt_wake, dq_wake
+  INTEGER, DIMENSION(klon),      INTENT (OUT)        :: wake_k
+  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_d_deltat_gw
+  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_h
+  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_dth
+  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_pe, wake_fip, wake_gfl
+  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: t_x, q_x
+  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_omgbdth, wake_dp_omgb
+  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_dtke, wake_dqke
+  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_omg, wake_dp_deltomg
+  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_spread
+  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_cstar
+  REAL, DIMENSION(klon, klev),   INTENT (OUT)        :: wake_ddeltat, wake_ddeltaq
+  REAL, DIMENSION(klon),         INTENT (OUT)        :: wake_ds, wake_ddens
+
 
   ! Variable internes
   ! -----------------
-
-  REAL aire
-  REAL p(klon, klev), ph(klon, klev+1), pi(klon, klev)
-  REAL te(klon, klev), qe(klon, klev), omgbe(klon, klev+1)
-  REAL dtdwn(klon, klev), dqdwn(klon, klev)
-  REAL dta(klon, klev), dqa(klon, klev)
-  REAL wdtpbl(klon, klev), wdqpbl(klon, klev)
-  REAL udtpbl(klon, klev), udqpbl(klon, klev)
-  REAL amdwn(klon, klev), amup(klon, klev)
-  REAL dtw(klon, klev), dqw(klon, klev), dth(klon, klev)
-  REAL d_deltat_gw(klon, klev)
-  REAL dtls(klon, klev), dqls(klon, klev)
-  REAL tu(klon, klev), qu(klon, klev)
-  REAL hw(klon), sigmaw(klon), wape(klon), fip(klon), gfl(klon)
-  REAL omgbdth(klon, klev+1), dp_omgb(klon, klev)
-  REAL dtke(klon, klev), dqke(klon, klev)
-  REAL dtpbl(klon, klev), dqpbl(klon, klev)
-  REAL omg(klon, klev+1), dp_deltomg(klon, klev), spread(klon, klev)
-  REAL cstar(klon)
-  REAL sigd0(klon), wdens(klon)
-
-  REAL rdcp
+  INTEGER                                            :: i, l
+  REAL                                               :: aire
+  REAL, DIMENSION(klon, klev)                        :: p,  pi
+  REAL, DIMENSION(klon, klev+1)                      ::  ph, omgbe
+  REAL, DIMENSION(klon, klev)                        :: te, qe
+  REAL, DIMENSION(klon, klev)                        :: dtdwn, dqdwn
+  REAL, DIMENSION(klon, klev)                        :: dta, dqa
+  REAL, DIMENSION(klon, klev)                        :: amdwn, amup
+  REAL, DIMENSION(klon, klev)                        :: dtw, dqw, dth
+  REAL, DIMENSION(klon, klev)                        :: dtls, dqls
+  REAL, DIMENSION(klon, klev)                        :: tx, qx
+  REAL, DIMENSION(klon)                              :: hw, wape, fip, gfl
+  REAL, DIMENSION(klon)                              :: sigmaw, wdens
+  REAL, DIMENSION(klon, klev+1)                      :: omgbdth
+  REAL, DIMENSION(klon, klev)                        :: dp_omgb
+  REAL, DIMENSION(klon, klev)                        :: dtke, dqke
+  REAL, DIMENSION(klon, klev+1)                      :: omg
+  REAL, DIMENSION(klon, klev)                        :: dp_deltomg, spread
+  REAL, DIMENSION(klon)                              :: cstar
+  REAL, DIMENSION(klon)                              :: sigd0
+  INTEGER, DIMENSION(klon)                           :: ktopw
+  REAL, DIMENSION(klon, klev)                        :: d_deltat_gw
+  REAL, DIMENSION(klon, klev)                        :: d_deltatw, d_deltaqw
+  REAL, DIMENSION(klon)                              :: d_sigmaw, d_wdens
+
+  REAL                                               :: rdcp
+
 
   ! print *, '-> calwake, wake_s ', wake_s(1)
@@ -104,8 +120,4 @@
       dta(i, l) = dt_a(i, l)
       dqa(i, l) = dq_a(i, l)
-      wdtpbl(i, l) = wdt_pbl(i, l)
-      wdqpbl(i, l) = wdq_pbl(i, l)
-      udtpbl(i, l) = udt_pbl(i, l)
-      udqpbl(i, l) = udq_pbl(i, l)
     END DO
   END DO
@@ -125,4 +137,9 @@
   END DO
 
+  DO i = 1, klon
+    hw(i) = wake_h(i)
+  END DO
+!
+!    Make a copy of state variables
   DO l = 1, klev
     DO i = 1, klon
@@ -132,14 +149,10 @@
   END DO
 
-  DO l = 1, klev
-    DO i = 1, klon
-      dtls(i, l) = dt_wake(i, l)
-      dqls(i, l) = dq_wake(i, l)
-    END DO
-  END DO
-
-  DO i = 1, klon
-    hw(i) = wake_h(i)
+  DO i = 1, klon
     sigmaw(i) = wake_s(i)
+  END DO
+
+  DO i = 1, klon
+    wdens(i) = wake_dens(i)
   END DO
 
@@ -166,15 +179,19 @@
   END DO
 
-  CALL wake(p, ph, pi, dtime, sigd0, te, qe, omgbe, dtdwn, dqdwn, amdwn, &
-    amup, dta, dqa, wdtpbl, wdqpbl, udtpbl, udqpbl, dtw, dqw, dth, hw, &
-    sigmaw, wape, fip, gfl, dtls, dqls, ktopw, omgbdth, dp_omgb, wdens, tu, &
-    qu, dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, spread, cstar, &
-    d_deltat_gw, d_deltatw, d_deltaqw)
-
+  CALL wake(p, ph, pi, dtime, &
+    te, qe, omgbe, &
+    dtdwn, dqdwn, amdwn, amup, dta, dqa, &
+    sigd0, &
+    dtw, dqw, sigmaw, wdens, &                                   ! state variables
+    dth, hw, wape, fip, gfl, &
+    dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, &
+    dtke, dqke, omg, dp_deltomg, spread, cstar, &
+    d_deltat_gw, &
+    d_deltatw, d_deltaqw, d_sigmaw, d_wdens)                     ! tendencies
+
+!
   DO l = 1, klev
     DO i = 1, klon
       IF (ktopw(i)>0) THEN
-        wake_deltat(i, l) = dtw(i, l)
-        wake_deltaq(i, l) = dqw(i, l)
         wake_d_deltat_gw(i, l) = d_deltat_gw(i, l)
         wake_omgbdth(i, l) = omgbdth(i, l)
@@ -182,19 +199,13 @@
         wake_dtke(i, l) = dtke(i, l)
         wake_dqke(i, l) = dqke(i, l)
-        wake_dtpbl(i, l) = dtpbl(i, l)
-        wake_dqpbl(i, l) = dqpbl(i, l)
         wake_omg(i, l) = omg(i, l)
         wake_dp_deltomg(i, l) = dp_deltomg(i, l)
         wake_spread(i, l) = spread(i, l)
         wake_dth(i, l) = dth(i, l)
-        dt_wake(i, l) = dtls(i, l)
-        dq_wake(i, l) = dqls(i, l)
-        undi_t(i, l) = tu(i, l)
-        undi_q(i, l) = qu(i, l)
-        wake_ddeltat(i, l) = d_deltatw(i, l)
-        wake_ddeltaq(i, l) = d_deltaqw(i, l)
+        dt_wake(i, l) = dtls(i, l)*dtime         ! derivative -> tendency
+        dq_wake(i, l) = dqls(i, l)*dtime         ! derivative -> tendency
+        t_x(i, l) = tx(i, l)
+        q_x(i, l) = qx(i, l)
       ELSE
-        wake_deltat(i, l) = 0.
-        wake_deltaq(i, l) = 0.
         wake_d_deltat_gw(i, l) = 0.
         wake_omgbdth(i, l) = 0.
@@ -202,6 +213,4 @@
         wake_dtke(i, l) = 0.
         wake_dqke(i, l) = 0.
-        wake_dtpbl(i, l) = 0.
-        wake_dqpbl(i, l) = 0.
         wake_omg(i, l) = 0.
         wake_dp_deltomg(i, l) = 0.
@@ -210,8 +219,6 @@
         dt_wake(i, l) = 0.
         dq_wake(i, l) = 0.
-        undi_t(i, l) = te(i, l)
-        undi_q(i, l) = qe(i, l)
-        wake_ddeltat(i, l) = 0.
-        wake_ddeltaq(i, l) = 0.
+        t_x(i, l) = te(i, l)
+        q_x(i, l) = qe(i, l)
       END IF
     END DO
@@ -220,5 +227,4 @@
   DO i = 1, klon
     wake_h(i) = hw(i)
-    wake_s(i) = sigmaw(i)
     wake_pe(i) = wape(i)
     wake_fip(i) = fip(i)
@@ -226,7 +232,50 @@
     wake_k(i) = ktopw(i)
     wake_cstar(i) = cstar(i)
-    wake_dens(i) = wdens(i)
-  END DO
+  END DO
+
+!  Tendencies of state variables
+  DO l = 1, klev
+    DO i = 1, klon
+      IF (ktopw(i)>0) THEN
+        wake_ddeltat(i, l) = d_deltatw(i, l)*dtime
+        wake_ddeltaq(i, l) = d_deltaqw(i, l)*dtime
+      ELSE
+        wake_ddeltat(i, l) = -wake_deltat(i, l)
+        wake_ddeltaq(i, l) = -wake_deltaq(i, l)
+      END IF
+    END DO
+  END DO
+  DO i = 1, klon
+    IF (ktopw(i)>0) THEN
+      wake_ds(i) = d_sigmaw(i)*dtime
+      wake_ddens(i) = d_wdens(i)*dtime
+    ELSE
+      wake_ds(i)   = -wake_s(i)
+      wake_ddens(i)= -wake_dens(i)
+    END IF
+  END DO
+
+!jyg<  
+  IF (iflag_wake_tend .EQ. 0) THEN
+!  Update State variables
+    DO l = 1, klev
+      DO i = 1, klon
+        IF (ktopw(i)>0) THEN
+          wake_deltat(i, l) = dtw(i, l)
+          wake_deltaq(i, l) = dqw(i, l)
+        ELSE
+          wake_deltat(i, l) = 0.
+          wake_deltaq(i, l) = 0.
+        END IF
+      END DO
+    END DO
+    DO i = 1, klon
+      wake_s(i) = sigmaw(i)
+      wake_dens(i) = wdens(i)
+    END DO
+  ENDIF
+!>jyg
 
   RETURN
 END SUBROUTINE calwake
+
Index: /LMDZ5/trunk/libf/phylmd/dyn1d/lmdz1d.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/dyn1d/lmdz1d.F90	(revision 2634)
+++ /LMDZ5/trunk/libf/phylmd/dyn1d/lmdz1d.F90	(revision 2635)
@@ -16,17 +16,17 @@
        ftsol, pbl_tke, pctsrf, radsol, rain_fall, snow_fall, ratqs, &
        rnebcon, rugoro, sig1, w01, solaire_etat0, sollw, sollwdown, &
-       solsw, t_ancien, q_ancien, u_ancien, v_ancien, wake_cstar, wake_deltaq, &
-       wake_deltat, wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
-       wake_s, zgam, &
-       zmax0, zmea, zpic, zsig, &
+       solsw, t_ancien, q_ancien, u_ancien, v_ancien, wake_cstar, &
+       wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
+       wake_deltaq, wake_deltat, wake_s, wake_dens, &
+       zgam, zmax0, zmea, zpic, zsig, &
        zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl
-   use dimphy
-   use surface_data, only : type_ocean,ok_veget
-   use pbl_surface_mod, only : ftsoil, pbl_surface_init,                     &
-     &                            pbl_surface_final
-      use fonte_neige_mod, only : fonte_neige_init, fonte_neige_final 
-
-   use infotrac ! new
-   use control_mod
+   USE dimphy
+   USE surface_data, only : type_ocean,ok_veget
+   USE pbl_surface_mod, only : ftsoil, pbl_surface_init, &
+                                 pbl_surface_final
+   USE fonte_neige_mod, only : fonte_neige_init, fonte_neige_final 
+
+   USE infotrac ! new
+   USE control_mod
    USE indice_sol_mod
    USE phyaqua_mod
@@ -787,4 +787,5 @@
         wake_pe = 0.
         wake_s = 0.
+        wake_dens = 0.
         ale_bl = 0.
         ale_bl_trig = 0.
@@ -810,5 +811,6 @@
 ! t_ancien,q_ancien,,frugs(:,is_oce),clwcon(:,1),rnebcon(:,1),ratqs(:,1)
 ! run_off_lic_0,pbl_tke(:,1:klev,nsrf), zmax0,f0,sig1,w01
-! wake_deltat,wake_deltaq,wake_s,wake_cstar,wake_fip,wake_delta_pbl_tke(:,1:klev,nsrf)
+! wake_deltat,wake_deltaq,wake_s,wake_dens,wake_cstar,
+! wake_fip,wake_delta_pbl_tke(:,1:klev,nsrf)
 !
 ! NB2: The content of the startphy.nc file depends on some flags defined in
Index: /LMDZ5/trunk/libf/phylmd/phyetat0.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phyetat0.F90	(revision 2634)
+++ /LMDZ5/trunk/libf/phylmd/phyetat0.F90	(revision 2635)
@@ -17,5 +17,5 @@
        solsw, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
        wake_deltat, wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
-       wake_s, zgam, zmax0, zmea, zpic, zsig, &
+       wake_s, wake_dens, zgam, zmax0, zmea, zpic, zsig, &
        zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, u10m, v10m
   USE geometry_mod, ONLY : longitude_deg, latitude_deg
@@ -377,5 +377,6 @@
   found=phyetat0_get(klev,wake_deltat,"WAKE_DELTAT","Delta T wake/env",0.)
   found=phyetat0_get(klev,wake_deltaq,"WAKE_DELTAQ","Delta hum. wake/env",0.)
-  found=phyetat0_get(1,wake_s,"WAKE_S","WAKE_S",0.)
+  found=phyetat0_get(1,wake_s,"WAKE_S","Wake frac. area",0.)
+  found=phyetat0_get(1,wake_dens,"WAKE_DENS","Wake num. /unit area",0.)
   found=phyetat0_get(1,wake_cstar,"WAKE_CSTAR","WAKE_CSTAR",0.)
   found=phyetat0_get(1,wake_pe,"WAKE_PE","WAKE_PE",0.)
Index: /LMDZ5/trunk/libf/phylmd/phyredem.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phyredem.F90	(revision 2634)
+++ /LMDZ5/trunk/libf/phylmd/phyredem.F90	(revision 2635)
@@ -18,5 +18,6 @@
                                 v_ancien, clwcon, rnebcon, ratqs, pbl_tke,   &
                                 wake_delta_pbl_tke, zmax0, f0, sig1, w01,    &
-                                wake_deltat, wake_deltaq, wake_s, wake_cstar,&
+                                wake_deltat, wake_deltaq, wake_s, wake_dens, &
+                                wake_cstar,                                  &
                                 wake_pe, wake_fip, fm_therm, entr_therm,     &
                                 detr_therm, Ale_bl, Ale_bl_trig, Alp_bl,     &
@@ -258,5 +259,7 @@
   CALL put_field("WAKE_DELTAQ", "WAKE_DELTAQ", wake_deltaq)
 
-  CALL put_field("WAKE_S", "WAKE_S", wake_s)
+  CALL put_field("WAKE_S", "Wake frac. area", wake_s)
+
+  CALL put_field("WAKE_DENS", "Wake num. /unit area", wake_dens)
 
   CALL put_field("WAKE_CSTAR", "WAKE_CSTAR", wake_cstar)
Index: /LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 2634)
+++ /LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 2635)
@@ -253,6 +253,32 @@
       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: sens, flwp, fiwp 
 !$OMP THREADPRIVATE(sens, flwp, fiwp)
-      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: ale_wake, alp_wake, bils
-!$OMP THREADPRIVATE(ale_wake, alp_wake, bils)
+!!
+!!         Wake variables
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:)            :: ale_wake, alp_wake
+!$OMP THREADPRIVATE(ale_wake, alp_wake)           
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:)            :: wake_h
+!$OMP THREADPRIVATE(wake_h)                       
+    REAL,ALLOCATABLE,SAVE,DIMENSION(:,:)            :: wake_omg
+!$OMP THREADPRIVATE(wake_omg)                     
+    REAL, SAVE, ALLOCATABLE,DIMENSION(:,:)          :: d_deltat_wk, d_deltaq_wk
+!$OMP THREADPRIVATE(d_deltat_wk, d_deltaq_wk)
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:)            :: d_s_wk, d_dens_wk
+!$OMP THREADPRIVATE(d_s_wk, d_dens_wk)
+    REAL, SAVE, ALLOCATABLE,DIMENSION(:,:)          :: d_deltat_wk_gw, d_deltaq_wk_gw
+!$OMP THREADPRIVATE(d_deltat_wk_gw, d_deltaq_wk_gw)
+    REAL, SAVE, ALLOCATABLE,DIMENSION(:,:)          :: d_deltat_vdf, d_deltaq_vdf
+!$OMP THREADPRIVATE(d_deltat_vdf, d_deltaq_vdf)
+!!!      REAL,ALLOCATABLE,SAVE,DIMENSION(:)          :: d_s_vdf, d_dens_vdf
+!!!OMP THREADPRIVATE(d_s_vdf, d_dens_vdf)
+    REAL, SAVE, ALLOCATABLE,DIMENSION(:,:)          :: d_deltat_the, d_deltaq_the
+!$OMP THREADPRIVATE(d_deltat_the, d_deltaq_the)
+!!!      REAL,ALLOCATABLE,SAVE,DIMENSION(:)          :: d_s_the, d_dens_the
+!!!OMP THREADPRIVATE(d_s_the, d_dens_the)
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:)           :: d_deltat_ajs_cv, d_deltaq_ajs_cv
+!$OMP THREADPRIVATE(d_deltat_ajs_cv, d_deltaq_ajs_cv)                       
+!!         End of Wake variables
+!!
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: bils
+!$OMP THREADPRIVATE(bils)
       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: cdragm, cdragh
 !$OMP THREADPRIVATE(cdragm, cdragh)
@@ -312,6 +338,4 @@
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: dqvdf_x, dqvdf_w
 !$OMP THREADPRIVATE(dqvdf_x, dqvdf_w)
-      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: undi_tke, wake_tke
-!$OMP THREADPRIVATE(undi_tke, wake_tke)
 ! Variables supplémentaires dans physiq.F relative au splitting de la surface
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: pbl_tke_input
@@ -330,6 +354,6 @@
 !>jyg+nrlmd
   !
-      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: wake_h, wbeff, zmax_th, zq2m, zt2m
-!$OMP THREADPRIVATE(wake_h, wbeff, zmax_th, zq2m, zt2m)
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: wbeff, zmax_th, zq2m, zt2m
+!$OMP THREADPRIVATE(wbeff, zmax_th, zq2m, zt2m)
       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zt2m_min_mon, zt2m_max_mon
 !$OMP THREADPRIVATE(zt2m_min_mon, zt2m_max_mon)
@@ -365,6 +389,6 @@
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: ref_liq_pi, ref_ice_pi
 !$OMP THREADPRIVATE(ref_liq_pi, ref_ice_pi)
-      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: wake_omg, zx_rh
-!$OMP THREADPRIVATE(wake_omg, zx_rh)
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: zx_rh
+!$OMP THREADPRIVATE(zx_rh)
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: pmflxr, pmflxs, prfl, psfl, fraca
 !$OMP THREADPRIVATE(pmflxr, pmflxs, prfl, psfl, fraca)
@@ -530,5 +554,20 @@
       ALLOCATE(tal1(klon), pal1(klon), pab1(klon), pab2(klon))
       ALLOCATE(ptstar(klon),pt0(klon),slp(klon))
-      ALLOCATE(ale_wake(klon), alp_wake(klon), bils(klon))
+!!
+!!          Wake variables
+      ALLOCATE(ale_wake(klon), alp_wake(klon))
+      ALLOCATE(wake_h(klon))
+      ALLOCATE(wake_omg(klon, klev))
+      ALLOCATE(d_deltat_wk(klon, klev), d_deltaq_wk(klon, klev))
+      ALLOCATE(d_s_wk(klon), d_dens_wk(klon))
+      ALLOCATE(d_deltat_wk_gw(klon, klev), d_deltaq_wk_gw(klon, klev))
+      ALLOCATE(d_deltat_vdf(klon, klev), d_deltaq_vdf(klon, klev))
+!!      ALLOCATE( d_s_vdf(klon), d_dens_vdf(klon))
+      ALLOCATE(d_deltat_the(klon, klev), d_deltaq_the(klon, klev))
+!!      ALLOCATE( d_s_the(klon), d_dens_the(klon))
+      ALLOCATE(d_deltat_ajs_cv(klon, klev), d_deltaq_ajs_cv(klon, klev))
+!!         End of wake variables
+!!
+      ALLOCATE(bils(klon))
       ALLOCATE(cdragm(klon), cdragh(klon), cldh(klon), cldl(klon))
       ALLOCATE(cldm(klon), cldq(klon), cldt(klon), qsat2m(klon))
@@ -561,5 +600,4 @@
       ALLOCATE(dtvdf_x(klon,klev), dtvdf_w(klon,klev))
       ALLOCATE(dqvdf_x(klon,klev), dqvdf_w(klon,klev))
-      ALLOCATE(undi_tke(klon,klev), wake_tke(klon,klev))
       ALLOCATE(pbl_tke_input(klon,klev+1,nbsrf))
       ALLOCATE(t_therm(klon,klev), q_therm(klon,klev),u_therm(klon,klev), v_therm(klon,klev))
@@ -568,5 +606,5 @@
       ALLOCATE(kh(klon), kh_x(klon), kh_w(klon))
 !
-      ALLOCATE(wake_h(klon), wbeff(klon), zmax_th(klon))
+      ALLOCATE(wbeff(klon), zmax_th(klon))
       ALLOCATE(zq2m(klon), zt2m(klon), weak_inversion(klon))
       ALLOCATE(zt2m_min_mon(klon), zt2m_max_mon(klon))
@@ -589,5 +627,5 @@
       ALLOCATE(ref_liq(klon, klev), ref_ice(klon, klev), theta(klon, klev))
       ALLOCATE(ref_liq_pi(klon, klev), ref_ice_pi(klon, klev))
-      ALLOCATE(zphi(klon, klev), wake_omg(klon, klev), zx_rh(klon, klev))
+      ALLOCATE(zphi(klon, klev), zx_rh(klon, klev))
       ALLOCATE(pmfd(klon, klev), pmfu(klon, klev))
 
@@ -741,5 +779,18 @@
       DEALLOCATE(tal1, pal1, pab1, pab2)
       DEALLOCATE(ptstar, pt0, slp)
-      DEALLOCATE(ale_wake, alp_wake, bils)
+!
+      DEALLOCATE(ale_wake, alp_wake)
+      DEALLOCATE(wake_h)
+      DEALLOCATE(wake_omg)
+      DEALLOCATE(d_deltat_wk, d_deltaq_wk)
+      DEALLOCATE(d_s_wk, d_dens_wk)
+      DEALLOCATE(d_deltat_wk_gw, d_deltaq_wk_gw)
+      DEALLOCATE(d_deltat_vdf, d_deltaq_vdf)
+!!      DEALLOCATE( d_s_vdf, d_dens_vdf)
+      DEALLOCATE(d_deltat_the, d_deltaq_the)
+!!      DEALLOCATE( d_s_the, d_dens_the)
+      DEALLOCATE(d_deltat_ajs_cv, d_deltaq_ajs_cv)
+!
+      DEALLOCATE(bils)
       DEALLOCATE(cdragm, cdragh, cldh, cldl)
       DEALLOCATE(cldm, cldq, cldt, qsat2m)
@@ -770,5 +821,4 @@
       DEALLOCATE(dtvdf_x, dtvdf_w)
       DEALLOCATE(dqvdf_x, dqvdf_w)
-      DEALLOCATE(undi_tke, wake_tke)
       DEALLOCATE(pbl_tke_input)
       DEALLOCATE(t_therm, q_therm, u_therm, v_therm)
@@ -777,5 +827,5 @@
       DEALLOCATE(kh, kh_x, kh_w)
 !
-      DEALLOCATE(wake_h, wbeff, zmax_th)
+      DEALLOCATE(wbeff, zmax_th)
       DEALLOCATE(zq2m, zt2m, weak_inversion)
       DEALLOCATE(zt2m_min_mon, zt2m_max_mon)
@@ -798,5 +848,5 @@
       DEALLOCATE(ref_liq, ref_ice, theta)
       DEALLOCATE(ref_liq_pi, ref_ice_pi)
-      DEALLOCATE(zphi, wake_omg, zx_rh)
+      DEALLOCATE(zphi, zx_rh)
       DEALLOCATE(pmfd, pmfu)
 
Index: /LMDZ5/trunk/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 2634)
+++ /LMDZ5/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 2635)
@@ -203,8 +203,8 @@
       REAL,ALLOCATABLE,SAVE :: cin(:)
 !$OMP THREADPRIVATE(cin)
-! ftd : differential heating between wake and environment
+! ftd : convective heating due to unsaturated downdraughts
       REAL,ALLOCATABLE,SAVE :: ftd(:,:)
 !$OMP THREADPRIVATE(ftd)
-! fqd : differential moistening between wake and environment
+! fqd : convective moistening due to unsaturated downdraughts
       REAL,ALLOCATABLE,SAVE :: fqd(:,:)     
 !$OMP THREADPRIVATE(fqd)
@@ -232,25 +232,24 @@
 ! wake_deltat : ecart de temperature avec la zone non perturbee
 ! wake_deltaq : ecart d'humidite avec la zone non perturbee
+! wake_s      : fraction surfacique occupee par la poche froide
+! wake_dens   : number of wakes per unit area
+! wake_occ    : occurence of wakes (= 1 if wakes occur, =0 otherwise)
 ! wake_Cstar  : vitesse d'etalement de la poche
-! wake_s      : fraction surfacique occupee par la poche froide
 ! wake_pe     : wake potential energy - WAPE
 ! wake_fip    : Gust Front Impinging power - ALP
-! dt_wake, dq_wake: LS tendencies due to wake
       REAL,ALLOCATABLE,SAVE :: wake_deltat(:,:)
 !$OMP THREADPRIVATE(wake_deltat)
       REAL,ALLOCATABLE,SAVE :: wake_deltaq(:,:)
 !$OMP THREADPRIVATE(wake_deltaq)
+      REAL,ALLOCATABLE,SAVE :: wake_s(:)
+!$OMP THREADPRIVATE(wake_s)
+      REAL,ALLOCATABLE,SAVE :: wake_dens(:)
+!$OMP THREADPRIVATE(wake_dens)
       REAL,ALLOCATABLE,SAVE :: wake_Cstar(:)
 !$OMP THREADPRIVATE(wake_Cstar)
-      REAL,ALLOCATABLE,SAVE :: wake_s(:)
-!$OMP THREADPRIVATE(wake_s)
       REAL,ALLOCATABLE,SAVE :: wake_pe(:)
 !$OMP THREADPRIVATE(wake_pe)
       REAL,ALLOCATABLE,SAVE :: wake_fip(:)
 !$OMP THREADPRIVATE(wake_fip)
-      REAL,ALLOCATABLE,SAVE :: dt_wake(:,:)
-!$OMP THREADPRIVATE(dt_wake)
-      REAL,ALLOCATABLE,SAVE :: dq_wake(:,:)
-!$OMP THREADPRIVATE(dq_wake)
 !
 !jyg<
@@ -524,7 +523,7 @@
       ALLOCATE(wght_th(klon,klev))
       ALLOCATE(wake_deltat(klon,klev), wake_deltaq(klon,klev))
-      ALLOCATE(wake_Cstar(klon), wake_s(klon))
+      ALLOCATE(wake_s(klon), wake_dens(klon))
+      ALLOCATE(wake_Cstar(klon))
       ALLOCATE(wake_pe(klon), wake_fip(klon))
-      ALLOCATE(dt_wake(klon,klev), dq_wake(klon,klev))
 !jyg<
       ALLOCATE(wake_delta_pbl_TKE(klon,klev+1,nbsrf+1))
@@ -664,6 +663,6 @@
       deallocate(lalim_conv, wght_th)
       deallocate(wake_deltat, wake_deltaq)
-      deallocate(wake_Cstar, wake_s, wake_pe, wake_fip)
-      deallocate(dt_wake, dq_wake)
+      deallocate(wake_s, wake_dens)
+      deallocate(wake_Cstar, wake_pe, wake_fip)
 !jyg<
       deallocate(wake_delta_pbl_TKE)
Index: /LMDZ5/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/physiq_mod.F90	(revision 2634)
+++ /LMDZ5/trunk/libf/phylmd/physiq_mod.F90	(revision 2635)
@@ -133,5 +133,4 @@
        dtvdf_x, dtvdf_w, &
        dqvdf_x, dqvdf_w, &
-       undi_tke, wake_tke, &
        pbl_tke_input, &
        t_therm, q_therm, u_therm, v_therm, &
@@ -141,5 +140,17 @@
        !
        ale_wake, alp_wake, & 
-       wake_h, wbeff, zmax_th, &
+       wake_h, wake_omg, &
+                       ! tendencies of delta T and delta q:
+       d_deltat_wk, d_deltaq_wk, &         ! due to wakes
+       d_deltat_wk_gw, d_deltaq_wk_gw, &   ! due to wake induced gravity waves
+       d_deltat_vdf, d_deltaq_vdf, &       ! due to vertical diffusion
+       d_deltat_the, d_deltaq_the, &       ! due to thermals
+       d_deltat_ajs_cv, d_deltaq_ajs_cv, & ! due to dry adjustment of (w) before convection
+                       ! tendencies of wake fractional area and wake number per unit area:
+       d_s_wk,  d_dens_wk, &             ! due to wakes
+!!!       d_s_vdf, d_dens_vdf, &            ! due to vertical diffusion
+!!!       d_s_the, d_dens_the, &            ! due to thermals
+       !                                  
+       wbeff, zmax_th, &
        sens, flwp, fiwp,  &
        ale_bl_stat,alp_bl_conv,alp_bl_det,  &
@@ -157,5 +168,5 @@
        ref_liq, ref_ice, theta,  &
        ref_liq_pi, ref_ice_pi,  &
-       zphi, wake_omg, zx_rh,  &
+       zphi, zx_rh,  &
        pmfd, pmfu,  &
        !
@@ -529,13 +540,13 @@
     !RC
     ! Variables li\'ees \`a la poche froide (jyg et rr)
-    ! Version diagnostique pour l'instant : pas de r\'etroaction sur
-    ! la convection
-
-    REAL t_wake(klon,klev),q_wake(klon,klev) ! wake pour la convection
+
+    INTEGER,  SAVE               :: iflag_wake_tend  ! wake: if =0, then wake state variables are
+                                                     ! updated within calwake
+    !$OMP THREADPRIVATE(iflag_wake_tend)
+    REAL t_w(klon,klev),q_w(klon,klev) ! temperature and moisture profiles in the wake region
+    REAL t_x(klon,klev),q_x(klon,klev) ! temperature and moisture profiles in the off-wake region
 
     REAL wake_dth(klon,klev)        ! wake : temp pot difference
 
-    REAL wake_d_deltat_gw(klon,klev)! wake : delta T tendency due to
-    ! Gravity Wave (/s)
     REAL wake_omgbdth(klon,klev)    ! Wake : flux of Delta_Theta
     ! transported by LS omega
@@ -546,9 +557,4 @@
     REAL wake_dqKE(klon,klev)       ! Wake : differential moistening
     ! (wake - unpertubed) CONV
-    REAL wake_dtPBL(klon,klev)      ! Wake : differential heating
-    ! (wake - unpertubed) PBL
-    REAL wake_dqPBL(klon,klev)      ! Wake : differential moistening
-    ! (wake - unpertubed) PBL
-    REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
     REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
     REAL wake_spread(klon,klev)     ! spreading term in wake_delt
@@ -557,9 +563,4 @@
     !
     INTEGER wake_k(klon)            ! Wake sommet
-    !
-    REAL t_undi(klon,klev)          ! temperature moyenne dans la zone
-    ! non perturbee
-    REAL q_undi(klon,klev)          ! humidite moyenne dans la zone
-    ! non perturbee
     !
     !jyg<
@@ -568,13 +569,9 @@
 
     REAL wake_gfl(klon)             ! Gust Front Length
-    REAL wake_dens(klon)
+!!!    REAL wake_dens(klon)         ! moved to phys_state_var_mod
     !
     !
     REAL dt_dwn(klon,klev)
     REAL dq_dwn(klon,klev)
-    REAL wdt_PBL(klon,klev)
-    REAL udt_PBL(klon,klev)
-    REAL wdq_PBL(klon,klev)
-    REAL udq_PBL(klon,klev)
     REAL M_dwn(klon,klev)
     REAL M_up(klon,klev)
@@ -589,12 +586,4 @@
     REAL, SAVE :: alp_offset
     !$OMP THREADPRIVATE(alp_offset)
-
-    ! !!
-    !=================================================================
-    !         PROVISOIRE : DECOUPLAGE PBL/WAKE
-    !         --------------------------------
-    REAL wake_deltat_sav(klon,klev)
-    REAL wake_deltaq_sav(klon,klev)
-    !=================================================================
 
     !
@@ -815,4 +804,6 @@
     ! tendance nulles
     REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0
+    REAL, dimension(klon)     :: dsig0, ddens0
+    INTEGER, dimension(klon)  :: wkoccur1
     !
     ! Flag pour pouvoir ne pas ajouter les tendances. 
@@ -1209,4 +1200,6 @@
        CALL getin_p('ratqsp0',ratqsp0)
        CALL getin_p('ratqsdp',ratqsdp)
+       iflag_wake_tend = 0
+       CALL getin_p('iflag_wake_tend',iflag_wake_tend)
     ENDIF
 
@@ -1698,4 +1691,7 @@
     dql0(:,:)=0.
     dqi0(:,:)=0.
+    dsig0(:) = 0.
+    ddens0(:) = 0.
+    wkoccur1(:)=1
     !
     ! Mettre a zero des variables de sortie (pour securite)
@@ -2135,13 +2131,4 @@
        ENDIF
        ! !!
-       !=================================================================
-       !         PROVISOIRE : DECOUPLAGE PBL/WAKE
-       !         --------------------------------
-       !
-       !!      wake_deltat_sav(:,:)=wake_deltat(:,:)
-       !!      wake_deltaq_sav(:,:)=wake_deltaq(:,:)
-       !!      wake_deltat(:,:)=0.
-       !!      wake_deltaq(:,:)=0.
-       !=================================================================
        !>jyg+nrlmd
        !
@@ -2218,16 +2205,15 @@
             )
        !
-       !=================================================================
-       !         PROVISOIRE : DECOUPLAGE PBL/WAKE
-       !         --------------------------------
-       !
-       !!      wake_deltat(:,:)=wake_deltat_sav(:,:)
-       !!      wake_deltaq(:,:)=wake_deltaq_sav(:,:)
-       !=================================================================
-       !
        !  Add turbulent diffusion tendency to the wake difference variables
        IF (mod(iflag_pbl_split,2) .NE. 0) THEN
-          wake_deltat(:,:) = wake_deltat(:,:) + (d_t_vdf_w(:,:)-d_t_vdf_x(:,:))
-          wake_deltaq(:,:) = wake_deltaq(:,:) + (d_q_vdf_w(:,:)-d_q_vdf_x(:,:))
+!jyg<
+          d_deltat_vdf(:,:) = d_t_vdf_w(:,:)-d_t_vdf_x(:,:)
+          d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
+          CALL add_wake_tend &
+             (d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, wkoccur1, 'vdf', abortphy) 
+       ELSE
+          d_deltat_vdf(:,:) = 0.
+          d_deltaq_vdf(:,:) = 0.
+!>jyg
        ENDIF
 
@@ -2404,24 +2390,24 @@
        !=======================================================================
        !ajout pour la parametrisation des poches froides: calcul de
-       !t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
-       do k=1,klev
-          do i=1,klon
-             if (iflag_wake>=1) then
-                t_wake(i,k) = t_seri(i,k) &
+       !t_w et t_x: si pas de poches froides, t_w=t_x=t_seri
+       if (iflag_wake>=1) then
+         do k=1,klev
+            do i=1,klon
+                t_w(i,k) = t_seri(i,k) &
                      +(1-wake_s(i))*wake_deltat(i,k)
-                q_wake(i,k) = q_seri(i,k) &
+                q_w(i,k) = q_seri(i,k) &
                      +(1-wake_s(i))*wake_deltaq(i,k)
-                t_undi(i,k) = t_seri(i,k) &
+                t_x(i,k) = t_seri(i,k) &
                      -wake_s(i)*wake_deltat(i,k)
-                q_undi(i,k) = q_seri(i,k) &
+                q_x(i,k) = q_seri(i,k) &
                      -wake_s(i)*wake_deltaq(i,k)
-             else
-                t_wake(i,k) = t_seri(i,k)
-                q_wake(i,k) = q_seri(i,k)
-                t_undi(i,k) = t_seri(i,k)
-                q_undi(i,k) = q_seri(i,k)
-             endif
-          enddo
-       enddo
+            enddo
+         enddo
+       else
+                t_w(:,:) = t_seri(:,:)
+                q_w(:,:) = q_seri(:,:)
+                t_x(:,:) = t_seri(:,:)
+                q_x(:,:) = q_seri(:,:)
+       endif
        !
        !jyg<
@@ -2432,5 +2418,5 @@
           IF (ok_adjwk) THEN
              limbas(:) = 1
-             CALL ajsec(paprs, pplay, t_wake, q_wake, limbas, &
+             CALL ajsec(paprs, pplay, t_w, q_w, limbas, &
                   d_t_adjwk, d_q_adjwk)
           ENDIF
@@ -2439,11 +2425,16 @@
              DO i=1,klon
                 IF (wake_s(i) .GT. 1.e-3) THEN
-                   t_wake(i,k) = t_wake(i,k) + d_t_adjwk(i,k)
-                   q_wake(i,k) = q_wake(i,k) + d_q_adjwk(i,k)
-                   wake_deltat(i,k) = wake_deltat(i,k) + d_t_adjwk(i,k)
-                   wake_deltaq(i,k) = wake_deltaq(i,k) + d_q_adjwk(i,k)
+                   t_w(i,k) = t_w(i,k) + d_t_adjwk(i,k)
+                   q_w(i,k) = q_w(i,k) + d_q_adjwk(i,k)
+                   d_deltat_ajs_cv(i,k) = d_t_adjwk(i,k)
+                   d_deltaq_ajs_cv(i,k) = d_q_adjwk(i,k)
+                ELSE
+                   d_deltat_ajs_cv(i,k) = 0.
+                   d_deltaq_ajs_cv(i,k) = 0.
                 ENDIF
              ENDDO
           ENDDO
+          CALL add_wake_tend &
+              (d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, wkoccur1, 'ajs_cv', abortphy) 
        ENDIF ! (iflag_wake>=1)
        !>jyg
@@ -2486,8 +2477,7 @@
           !jyg   iflag_con est dans clesphys
           !c          CALL concvl (iflag_con,iflag_clos,
-          clw=0.
           CALL concvl (iflag_clos, &
-               dtime, paprs, pplay, k_upper_cv, t_undi,q_undi, &
-               t_wake,q_wake,wake_s, &
+               dtime, paprs, pplay, k_upper_cv, t_x,q_x, &
+               t_w,q_w,wake_s, &
                u_seri,v_seri,tr_seri,nbtr_tmp, &
                ALE,ALP, &
@@ -2705,41 +2695,4 @@
           ENDDO
        ENDDO
-       !nrlmd+jyg<
-       DO k=1,klev
-          DO i=1,klon
-             wdt_PBL(i,k) =  0.
-             wdq_PBL(i,k) =  0.
-             udt_PBL(i,k) =  0.
-             udq_PBL(i,k) =  0.
-          ENDDO
-       ENDDO
-       !
-       IF (mod(iflag_pbl_split,2) .EQ. 1) THEN
-          DO k=1,klev
-             DO i=1,klon
-                wdt_PBL(i,k) = wdt_PBL(i,k) + d_t_vdf_w(i,k)/dtime
-                wdq_PBL(i,k) = wdq_PBL(i,k) + d_q_vdf_w(i,k)/dtime
-                udt_PBL(i,k) = udt_PBL(i,k) + d_t_vdf_x(i,k)/dtime
-                udq_PBL(i,k) = udq_PBL(i,k) + d_q_vdf_x(i,k)/dtime
-                !!        dt_dwn(i,k)  = dt_dwn(i,k) + d_t_vdf_w(i,k)/dtime
-                !!        dq_dwn(i,k)  = dq_dwn(i,k) + d_q_vdf_w(i,k)/dtime
-                !!        dt_a  (i,k)    = dt_a(i,k) + d_t_vdf_x(i,k)/dtime
-                !!        dq_a  (i,k)    = dq_a(i,k) + d_q_vdf_x(i,k)/dtime
-             ENDDO
-          ENDDO
-       ENDIF
-       IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
-          DO k=1,klev
-             DO i=1,klon
-                !!        dt_dwn(i,k)  = dt_dwn(i,k) + 0.
-                !!        dq_dwn(i,k)  = dq_dwn(i,k) + 0.
-                !!        dt_a(i,k)   = dt_a(i,k)   + d_t_ajs(i,k)/dtime
-                !!        dq_a(i,k)   = dq_a(i,k)   + d_q_ajs(i,k)/dtime
-                udt_PBL(i,k)   = udt_PBL(i,k)   + d_t_ajs(i,k)/dtime
-                udq_PBL(i,k)   = udq_PBL(i,k)   + d_q_ajs(i,k)/dtime
-             ENDDO
-          ENDDO
-       ENDIF
-       !>nrlmd+jyg
 
        IF (iflag_wake==2) THEN
@@ -2770,31 +2723,34 @@
        !
        !calcul caracteristiques de la poche froide
-       call calWAKE (paprs,pplay,dtime &
-            ,t_seri,q_seri,omega &
-            ,dt_dwn,dq_dwn,M_dwn,M_up &
-            ,dt_a,dq_a,sigd &
-            ,wdt_PBL,wdq_PBL &
-            ,udt_PBL,udq_PBL &
-            ,wake_deltat,wake_deltaq,wake_dth &
-            ,wake_h,wake_s,wake_dens &
-            ,wake_pe,wake_fip,wake_gfl &
-            ,dt_wake,dq_wake &
-            ,wake_k, t_undi,q_undi &
-            ,wake_omgbdth,wake_dp_omgb &
-            ,wake_dtKE,wake_dqKE &
-            ,wake_dtPBL,wake_dqPBL &
-            ,wake_omg,wake_dp_deltomg &
-            ,wake_spread,wake_Cstar,wake_d_deltat_gw &
-            ,wake_ddeltat,wake_ddeltaq)
+       call calWAKE (iflag_wake_tend, paprs, pplay, dtime, &
+            t_seri, q_seri, omega,  &
+            dt_dwn, dq_dwn, M_dwn, M_up,  &
+            dt_a, dq_a,  &
+            sigd,  &
+            wake_deltat, wake_deltaq, wake_s, wake_dens,  &
+            wake_dth, wake_h,  &
+            wake_pe, wake_fip, wake_gfl,  &
+            d_t_wake, d_q_wake,  &
+            wake_k, t_x, q_x,  &
+            wake_omgbdth, wake_dp_omgb,  &
+            wake_dtKE, wake_dqKE,  &
+            wake_omg, wake_dp_deltomg,  &
+            wake_spread, wake_Cstar, d_deltat_wk_gw,  &
+            d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk)
        !
        !-----------------------------------------------------------------------
        ! ajout des tendances des poches froides
-       ! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
-       ! coherent avec les autres d_t_...
-       d_t_wake(:,:)=dt_wake(:,:)*dtime
-       d_q_wake(:,:)=dq_wake(:,:)*dtime
        CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
             abortphy,flag_inhib_tend)
        !------------------------------------------------------------------------
+
+!       Increment Wake state variables
+       IF (iflag_wake_tend .GT. 0.) THEN
+
+         CALL add_wake_tend &
+            (d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk, wake_k, &
+             'wake', abortphy) 
+
+       ENDIF   ! (iflag_wake_tend .GT. 0.)
 
     endif  ! (iflag_wake>=1)
@@ -2919,9 +2875,6 @@
                 DO i=1,klon
                    !
-                   wake_deltat(i,k) = wake_deltat(i,k) - d_t_ajs(i,k)
-                   wake_deltaq(i,k) = wake_deltaq(i,k) - d_q_ajs(i,k)
-                   !
-                   !!!t_seri(i,k) = t_therm(i,k) + wake_s(i)*wake_deltat(i,k)
-                   !!!q_seri(i,k) = q_therm(i,k) + wake_s(i)*wake_deltaq(i,k)
+                   d_deltat_the(i,k) = - d_t_ajs(i,k) 
+                   d_deltaq_the(i,k) = - d_q_ajs(i,k) 
                    !
                    d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i)) 
@@ -2932,18 +2885,13 @@
                 ENDDO
              ENDDO
-!!!          ELSE
-!!!             DO k=1,klev
-!!!                DO i=1,klon
-!!!                   t_seri(i,k) = t_therm(i,k)
-!!!                   q_seri(i,k) = q_therm(i,k)
-!!!                ENDDO
-!!!             ENDDO
           ENDIF
+!
+          CALL add_wake_tend &
+              (d_deltat_the, d_deltaq_the, dsig0, ddens0, wkoccur1, 'the', abortphy) 
+!
           !
           CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs,  &
                              dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend)
           !
-          !>jyg
-!jyg<
 !
           CALL alpale_th( dtime, lmax_th, t_seri, cell_area,  &
@@ -2952,4 +2900,5 @@
                           alp_bl, alp_bl_stat, &
                           proba_notrig, random_notrig)
+          !>jyg
 
           ! ------------------------------------------------------------------
Index: /LMDZ5/trunk/libf/phylmd/wake.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/wake.F90	(revision 2634)
+++ /LMDZ5/trunk/libf/phylmd/wake.F90	(revision 2635)
@@ -2,12 +2,14 @@
 ! $Id$
 
-SUBROUTINE wake(p, ph, pi, dtime, sigd_con, &
+SUBROUTINE wake(p, ph, pi, dtime, &
                 te0, qe0, omgb, &
                 dtdwn, dqdwn, amdwn, amup, dta, dqa, &
-                wdtpbl, wdqpbl, udtpbl, udqpbl, &
-                deltatw, deltaqw, dth, hw, sigmaw, wape, fip, gfl, &
-                dtls, dqls, ktopw, omgbdth, dp_omgb, wdens, tu, qu, &
-                dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, spread, cstar, &
-                d_deltat_gw, d_deltatw2, d_deltaqw2)
+                sigd_con, &
+                deltatw, deltaqw, sigmaw, 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_wdens2)                ! tendencies
 
 
@@ -33,8 +35,9 @@
   ! le declenchement de nouvelles colonnes convectives.
 
-  ! Variables d'etat : deltatw    : ecart de temperature wake-undisturbed
-  ! area
-  ! deltaqw    : ecart d'humidite wake-undisturbed area
-  ! sigmaw     : fraction d'aire occupee par la poche.
+  ! 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 :
@@ -53,5 +56,5 @@
   ! omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
   ! dp_deltomg  : vertical gradient of omg (s-1)
-  ! spread  : spreading term in dt_wake and dq_wake
+  ! 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).
@@ -126,5 +129,4 @@
   REAL, DIMENSION (klon, klev),     INTENT(IN)          :: te0, qe0
   REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dtdwn, dqdwn
-  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: wdtpbl, wdqpbl, udtpbl, udqpbl ! UNUSED
   REAL, DIMENSION (klon, klev),     INTENT(IN)          :: amdwn, amup
   REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dta, dqa
@@ -133,6 +135,8 @@
   !
   ! Input/Output
+  ! State variables
   REAL, DIMENSION (klon, klev),     INTENT(INOUT)       :: deltatw, deltaqw
   REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw
+  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens
 
   ! Sorties
@@ -143,13 +147,13 @@
   REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtls, dqls
   REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke
-  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtpbl, dqpbl
   REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: spread
-  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
   REAL, DIMENSION (klon, klev+1),   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
-  REAL, DIMENSION (klon),           INTENT(OUT)         :: wdens
   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_wdens2
 
   ! Variables internes
@@ -157,94 +161,98 @@
 
   ! Variables à fixer
-  REAL alon
-  LOGICAL, SAVE :: first = .TRUE.
+  REAL                                                  :: alon
+  LOGICAL, SAVE                                         :: first = .TRUE.
   !$OMP THREADPRIVATE(first)
-  REAL, SAVE ::  stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol  
+  REAL, SAVE                                            :: stark, wdens_ref, coefgw, alpk
+  REAL, SAVE                                            :: crep_upper, crep_sol  
   !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
 
-  REAL delta_t_min
-  INTEGER nsub
-  REAL dtimesub
-  REAL sigmad, hwmin, wapecut
-  REAL :: sigmaw_max
-  REAL :: dens_rate
-  REAL wdens0
+  REAL                                                  :: delta_t_min
+  INTEGER                                               :: nsub
+  REAL                                                  :: dtimesub
+  REAL                                                  :: sigmad, hwmin, wapecut
+  REAL                                                  :: sigmaw_max
+  REAL                                                  :: dens_rate
+  REAL                                                  :: wdens0
   ! IM 080208
-  LOGICAL, DIMENSION (klon) :: gwake
+  LOGICAL, DIMENSION (klon)                             :: gwake
 
   ! Variables de sauvegarde
-  REAL, DIMENSION (klon, klev) :: deltatw0
-  REAL, DIMENSION (klon, klev) :: deltaqw0
-  REAL, DIMENSION (klon, klev) :: te, qe
-  REAL, DIMENSION (klon) :: sigmaw0, sigmaw1
+  REAL, DIMENSION (klon, klev)                          :: deltatw0
+  REAL, DIMENSION (klon, klev)                          :: deltaqw0
+  REAL, DIMENSION (klon, klev)                          :: te, qe
+  REAL, DIMENSION (klon)                                :: sigmaw0, sigmaw1
 
   ! Variables pour les GW
-  REAL, DIMENSION (klon) :: ll
-  REAL, DIMENSION (klon, klev) :: n2
-  REAL, DIMENSION (klon, klev) :: cgw
-  REAL, DIMENSION (klon, klev) :: tgw
+  REAL, DIMENSION (klon)                                :: ll
+  REAL, DIMENSION (klon, klev)                          :: n2
+  REAL, DIMENSION (klon, klev)                          :: cgw
+  REAL, DIMENSION (klon, klev)                          :: tgw
 
   ! Variables liées au calcul de hw
-  REAL, DIMENSION (klon) :: ptop_provis, ptop, ptop_new
-  REAL, DIMENSION (klon) :: sum_dth
-  REAL, DIMENSION (klon) :: dthmin
-  REAL, DIMENSION (klon) :: z, dz, hw0
-  INTEGER, DIMENSION (klon) :: ktop, kupper
+  REAL, DIMENSION (klon)                                :: ptop_provis, ptop, ptop_new
+  REAL, DIMENSION (klon)                                :: sum_dth
+  REAL, DIMENSION (klon)                                :: dthmin
+  REAL, DIMENSION (klon)                                :: z, dz, hw0
+  INTEGER, DIMENSION (klon)                             :: ktop, kupper
 
   ! Sub-timestep tendencies and related variables
-  REAL d_deltatw(klon, klev), d_deltaqw(klon, klev)
-  REAL d_te(klon, klev), d_qe(klon, klev)
-  REAL d_sigmaw(klon), alpha(klon)
-  REAL q0_min(klon), q1_min(klon)
-  LOGICAL wk_adv(klon), ok_qx_qw(klon)
-  REAL epsilon
+  REAL, DIMENSION (klon, klev)                          :: d_deltatw, d_deltaqw
+  REAL, DIMENSION (klon, klev)                          :: d_te, d_qe
+  REAL, DIMENSION (klon)                                :: d_sigmaw, alpha
+  REAL, DIMENSION (klon)                                :: q0_min, q1_min
+  LOGICAL, DIMENSION (klon)                             :: wk_adv, ok_qx_qw
+  REAL, SAVE                                            :: epsilon
+  !$OMP THREADPRIVATE(epsilon)
   DATA epsilon/1.E-15/
 
   ! Autres variables internes
-  INTEGER isubstep, k, i
-
-  REAL, DIMENSION (klon) :: sum_thu, sum_tu, sum_qu, sum_thvu
-  REAL, DIMENSION (klon) :: sum_dq, sum_rho
-  REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn
-  REAL, DIMENSION (klon) :: av_thu, av_tu, av_qu, av_thvu
-  REAL, DIMENSION (klon) :: av_dth, av_dq, av_rho
-  REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn
-
-  REAL, DIMENSION (klon, klev) :: rho, rhow
-  REAL, DIMENSION (klon, klev+1) :: rhoh
-  REAL, DIMENSION (klon, klev) :: rhow_moyen
-  REAL, DIMENSION (klon, klev) :: zh
-  REAL, DIMENSION (klon, klev+1) :: zhh
-  REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2
-
-  REAL, DIMENSION (klon, klev) :: the, thu
-
-  ! REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
-
-  REAL, DIMENSION (klon, klev+1) :: omgbw
-  REAL, DIMENSION (klon) :: pupper
-  REAL, DIMENSION (klon) :: omgtop
-  REAL, DIMENSION (klon, klev) :: dp_omgbw
-  REAL, DIMENSION (klon) :: ztop, dztop
-  REAL, DIMENSION (klon, klev) :: alpha_up
-
-  REAL, DIMENSION (klon) :: rre1, rre2
-  REAL :: rrd1, rrd2
-  REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2
-  REAL, DIMENSION (klon, klev) :: d_th1, d_th2, d_dth
-  REAL, DIMENSION (klon, klev) :: d_q1, d_q2, d_dq
-  REAL, DIMENSION (klon, klev) :: omgbdq
-
-  REAL, DIMENSION (klon) :: ff, gg
-  REAL, DIMENSION (klon) :: wape2, cstar2, heff
-
-  REAL, DIMENSION (klon, klev) :: crep
-
-  REAL, DIMENSION (klon, klev) :: ppi
+  INTEGER                                               ::isubstep, k, i
+
+  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
+
+  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+1)                        :: omgbw
+  REAL, DIMENSION (klon)                                :: pupper
+  REAL, DIMENSION (klon)                                :: omgtop
+  REAL, DIMENSION (klon, klev)                          :: dp_omgbw
+  REAL, DIMENSION (klon)                                :: ztop, dztop
+  REAL, DIMENSION (klon, klev)                          :: alpha_up
+
+  REAL, DIMENSION (klon)                                :: rre1, rre2
+  REAL                                                  :: rrd1, rrd2
+  REAL, DIMENSION (klon, klev)                          :: th1, th2, q1, q2
+  REAL, DIMENSION (klon, klev)                          :: d_th1, d_th2, d_dth
+  REAL, DIMENSION (klon, klev)                          :: d_q1, d_q2, d_dq
+  REAL, DIMENSION (klon, klev)                          :: omgbdq
+
+  REAL, DIMENSION (klon)                                :: ff, gg
+  REAL, DIMENSION (klon)                                :: wape2, cstar2, heff
+                                                        
+  REAL, DIMENSION (klon, klev)                          :: crep
+                                                        
+  REAL, DIMENSION (klon, klev)                          :: ppi
 
   ! cc nrlmd
-  REAL, DIMENSION (klon) :: death_rate, nat_rate
-  REAL, DIMENSION (klon, klev) :: entr
-  REAL, DIMENSION (klon, klev) :: detr
+  REAL, DIMENSION (klon)                                :: death_rate, nat_rate
+  REAL, DIMENSION (klon, klev)                          :: entr
+  REAL, DIMENSION (klon, klev)                          :: detr
+
+  REAL, DIMENSION(klon)                                 :: sigmaw_in   ! pour les prints
 
   ! -------------------------------------------------------------------------
@@ -312,5 +320,5 @@
   ! Les densites peuvent evoluer si les poches debordent
   ! (voir au tout debut de la boucle sur les substeps)
-  wdens = wdens_ref
+  wdens(:) = wdens_ref
 
   ! print*,'stark',stark
@@ -347,4 +355,7 @@
     END DO
   END DO
+  DO i = 1, klon
+   sigmaw_in(i) = sigmaw(i)
+  END DO
   ! sigmaw1=sigmaw
   ! IF (sigd_con.GT.sigmaw1) THEN
@@ -353,10 +364,16 @@
   DO i = 1, klon
     ! c      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
-    sigmaw(i) = amax1(sigmaw(i), sigmad)
-    sigmaw(i) = amin1(sigmaw(i), 0.99)
+!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
     sigmaw0(i) = sigmaw(i)
     wape(i) = 0.
     wape2(i) = 0.
     d_sigmaw(i) = 0.
+    d_wdens2(i) = 0.
     ktopw(i) = 0
   END DO
@@ -433,6 +450,6 @@
         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)))
+        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
@@ -504,6 +521,6 @@
       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))
+        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
@@ -602,4 +619,6 @@
         deltatw(i, k) = 0.
         deltaqw(i, k) = 0.
+        d_deltatw2(i,k) = -deltatw0(i,k)
+        d_deltaqw2(i,k) = -deltaqw0(i,k)
       END IF
     END DO
@@ -671,6 +690,7 @@
     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)
+    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
 
@@ -686,4 +706,6 @@
         deltaqw(i, k) = 0.
         dth(i, k) = 0.
+        d_deltatw2(i,k) = -deltatw0(i,k)
+        d_deltaqw2(i,k) = -deltaqw0(i,k)
       END IF
     END DO
@@ -695,5 +717,10 @@
       cstar(i) = 0.
       hw(i) = hwmin
-      sigmaw(i) = amax1(sigmad, sigd_con(i))
+!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.
@@ -708,11 +735,11 @@
   ! --------------------------
   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)))
+    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)))
+      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)
@@ -752,7 +779,7 @@
       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)
+                               rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
+        wdens0 = (sigmaw(i)/(4.*3.14))* &
+          ((1.-sigmaw(i))*omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))*cstar(i)))**(2)
         IF (wdens(i)<=wdens0*1.1) THEN
           wdens(i) = wdens0
@@ -770,7 +797,13 @@
       IF (wk_adv(i)) THEN
         gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
-        sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
-      END IF
-    END DO
+!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
+
     DO i = 1, klon
       IF (wk_adv(i)) THEN
@@ -783,4 +816,5 @@
           death_rate(i) = 0.
         END IF
+
         d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* &
           dtimesub
@@ -798,5 +832,4 @@
       END IF
     END DO
-
 
     ! calcul de la difference de vitesse verticale poche - zone non perturbee
@@ -1006,13 +1039,15 @@
 
           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)
+            (rrd1*omg(i,k)*sigmaw(i)*d_th1(i,k) - &
+             rrd2*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1)- &
+             (1.-alpha_up(i,k))*omgbdth(i,k)- &
+             alpha_up(i,k+1)*omgbdth(i,k+1))*ppi(i, k)
           ! print*,'d_deltatw=',d_deltatw(i,k)
 
           d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
-            (rrd1*omg(i,k)*sigmaw(i)*d_q1(i,k)-rrd2*omg(i,k+1)*(1.-sigmaw( &
-            i))*d_q2(i,k+1)-(1.-alpha_up(i,k))*omgbdq(i,k)-alpha_up(i,k+1)* &
-            omgbdq(i,k+1))
+            (rrd1*omg(i,k)*sigmaw(i)*d_q1(i,k)- &
+             rrd2*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1)- &
+             (1.-alpha_up(i,k))*omgbdq(i,k)- &
+             alpha_up(i,k+1)*omgbdq(i,k+1))
           ! print*,'d_deltaqw=',d_deltaqw(i,k)
 
@@ -1024,26 +1059,19 @@
           ! C
           ! -----------------------------------------------------------------
-          d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i, &
-            k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/(ph(i,k)-ph(i, &
-            k+1)) &                ! cc nrlmd     $
-                                   ! -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k)
-            -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/(ph(i, &
-            k)-ph(i,k+1)) &        ! cc
-            )*ppi(i, k)
-
-          d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i, &
-            k)-rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/(ph(i,k)-ph(i, &
-            k+1)) &                ! cc nrlmd     $
-                                   ! -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k)
-            -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i, &
-            k+1))/(ph(i,k)-ph(i,k+1)) & ! cc
-            )
-          ! cc nrlmd
+          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)) )
         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))))
+          d_te(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)/(ph(i,k)-ph(i,k+1)))*ppi(i, k)
+
+          d_qe(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)/(ph(i,k)-ph(i,k+1)))
 
         END IF
@@ -1067,6 +1095,6 @@
           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)))
+          crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/ &
+            (p(i,1)-ph(i,kupper(i)))
 
 
@@ -1091,21 +1119,4 @@
           ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
 
-!jyg<
-!!
-!!---------------------------------------------------------------
-!! The change of delta_T due to PBL (vertical diffusion plus thermal plumes)
-!! is accounted for by the PBL and the Thermals schemes. It is now set to zero 
-!! within the Wake scheme.
-!!---------------------------------------------------------------
-          dtPBL(i,k) = 0.
-          dqPBL(i,k) = 0.
-! 
-!!          dtPBL(i,k)=wdtPBL(i,k) - udtPBL(i,k)
-!!          dqPBL(i,k)=wdqPBL(i,k) - udqPBL(i,k)
-! 
-!!          dtpbl(i, k) = (wdtpbl(i,k)/sigmaw(i)-udtpbl(i,k)/(1.-sigmaw(i)))
-!!          dqpbl(i, k) = (wdqpbl(i,k)/sigmaw(i)-udqpbl(i,k)/(1.-sigmaw(i)))
-          ! print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k)
-!>jyg
 !
 
@@ -1145,17 +1156,14 @@
 
           IF (dtimesub*tgw(i,k)<1.E-10) THEN
-            d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k)+dtpbl(i,k) & ! cc
-                                                                     ! $
-                                                                     ! -spread(i,k)*deltatw(i,k)
-              -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( &
-              i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
-              -tgw(i,k)*deltatw(i,k))
+            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)+dtpbl(i,k) & ! cc     $
-                                                 ! -spread(i,k)*deltatw(i,k)
-              -entr(i,k)*deltatw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw( &
-              i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) & ! cc
-              -tgw(i,k)*deltatw(i,k))
+            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
 
@@ -1164,8 +1172,7 @@
           gg(i) = d_deltaqw(i, k)/dtimesub
 
-          d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k)+dqpbl(i,k) & ! cc     $
-                                                                   ! -spread(i,k)*deltaqw(i,k))
-            -entr(i,k)*deltaqw(i,k)/sigmaw(i)-(death_rate(i)*sigmaw(i)+detr( &
-            i,k))*deltaqw(i,k)/(1.-sigmaw(i)))
+          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)))
           ! cc
 
@@ -1239,4 +1246,7 @@
       IF (wk_adv(i)) THEN
         sigmaw(i) = sigmaw(i) + d_sigmaw(i)
+!jyg<
+        d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
+!>jyg
       END IF
     END DO
@@ -1258,6 +1268,6 @@
         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))
+          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
@@ -1331,6 +1341,6 @@
         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))
+          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
@@ -1359,4 +1369,6 @@
           deltatw(i, k) = 0.
           deltaqw(i, k) = 0.
+          d_deltatw2(i,k) = -deltatw0(i,k)
+          d_deltaqw2(i,k) = -deltaqw0(i,k)
         END IF
       END DO
@@ -1449,6 +1461,6 @@
         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)
+        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
@@ -1463,4 +1475,6 @@
             deltaqw(i, k) = 0.
             dth(i, k) = 0.
+            d_deltatw2(i,k) = -deltatw0(i,k)
+            d_deltaqw2(i,k) = -deltaqw0(i,k)
           END IF
         END IF
@@ -1474,5 +1488,10 @@
           cstar(i) = 0.
           hw(i) = hwmin
-          sigmaw(i) = max(sigmad, sigd_con(i))
+!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.
@@ -1486,23 +1505,4 @@
   END DO ! end sub-timestep loop
 
-  ! -----------------------------------------------------------------
-  ! Get back to tendencies per second
-
-  DO k = 1, klev
-    DO i = 1, klon
-
-      ! cc nrlmd        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
-      IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
-        ! cc
-        dtls(i, k) = dtls(i, k)/dtime
-        dqls(i, k) = dqls(i, k)/dtime
-        d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
-        d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
-        d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
-        ! c      print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
-        ! c     $         ,death_rate(i)*sigmaw(i)
-      END IF
-    END DO
-  END DO
 
 
@@ -1632,8 +1632,10 @@
       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)
+      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
@@ -1650,4 +1652,6 @@
         deltaqw(i, k) = 0.
         dth(i, k) = 0.
+        d_deltatw2(i,k) = -deltatw0(i,k)
+        d_deltaqw2(i,k) = -deltaqw0(i,k)
       END IF
     END DO
@@ -1663,5 +1667,10 @@
         cstar2(i) = 0.
         hw(i) = hwmin
-        sigmaw(i) = amax1(sigmad, sigd_con(i))
+!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.
@@ -1726,4 +1735,6 @@
         deltatw(i, k) = 0.
         deltaqw(i, k) = 0.
+        d_deltatw2(i,k) = -deltatw0(i,k)
+        d_deltaqw2(i,k) = -deltaqw0(i,k)
       END IF
     END DO
@@ -1734,4 +1745,5 @@
     IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=1.0)) .OR. (ktopw(i)<=2) .OR. &
         .NOT. ok_qx_qw(i)) THEN
+      ktopw(i) = 0
       wape(i) = 0.
       cstar(i) = 0.
@@ -1749,4 +1761,31 @@
     ! c     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
   END DO
+
+  ! -----------------------------------------------------------------
+  ! Get back to tendencies per second
+
+  DO k = 1, klev
+    DO i = 1, klon
+
+      ! cc nrlmd        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
+      IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
+        ! cc
+        dtls(i, k) = dtls(i, k)/dtime
+        dqls(i, k) = dqls(i, k)/dtime
+        d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
+        d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
+        d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
+        ! c      print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
+        ! c     $         ,death_rate(i)*sigmaw(i)
+      END IF
+    END DO
+  END DO
+!jyg<
+  DO i = 1, klon
+    d_sigmaw2(i) = d_sigmaw2(i)/dtime
+    d_wdens2(i) = d_wdens2(i)/dtime
+  ENDDO
+!>jyg
+
 
 
@@ -1792,6 +1831,6 @@
       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))
+          (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) - &
@@ -1807,12 +1846,12 @@
           ELSE
             IF (a>0.) THEN
-              alpha1(i) = 0.9*min((2.*c)/(-b+sqrt(discrim)), (-b+sqrt(discrim &
-                ))/(2.*a))
+              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))
+              alpha1(i) = 0.9*max( (2.*c)/(-b+sqrt(discrim)),  &
+                                   (-b+sqrt(discrim))/(2.*a))
             END IF
           END IF
@@ -1827,2 +1866,4 @@
 
 
+
+
