Index: LMDZ5/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/physiq_mod.F90	(revision 2630)
+++ 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
 
           ! ------------------------------------------------------------------
