Index: LMDZ5/branches/testing/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/pbl_surface_mod.F90	(revision 2160)
+++ LMDZ5/branches/testing/libf/phylmd/pbl_surface_mod.F90	(revision 2187)
@@ -12,4 +12,5 @@
   USE dimphy
   USE mod_phys_lmdz_para,  ONLY : mpi_size
+  USE mod_grid_phy_lmdz,   ONLY : klon_glo
   USE ioipsl
   USE surface_data,        ONLY : type_ocean, ok_veget
@@ -174,4 +175,9 @@
        rain_f,    snow_f,    solsw_m,  sollw_m,       &
        t,         q,         u,        v,             &
+!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
+!!       t_x,       q_x,       t_w,      q_w,           &
+       wake_dlt,             wake_dlq,                &
+       wake_cstar,           wake_s,                  &
+!!!
        pplay,     paprs,     pctsrf,                  &
        ts,        alb1, alb2,ustar, u10m, v10m,wstar, &
@@ -181,6 +187,20 @@
        zxtsol,    zxfluxlat, zt2m,     qsat2m,        &
        d_t,       d_q,       d_u,      d_v, d_t_diss, & 
+!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
+       d_t_w,     d_q_w,                              &
+       d_t_x,     d_q_x,                              & 
+!!       d_wake_dlt,d_wake_dlq,                         &
+       zxsens_x,  zxfluxlat_x,zxsens_w,zxfluxlat_w,   &
+!!!
+!!! nrlmd le 13/06/2011
+       delta_tsurf,wake_dens,cdragh_x,cdragh_w,       &
+       cdragm_x,cdragm_w,kh,kh_x,kh_w,                &
+!!!
        zcoefh,    zcoefm,    slab_wfbils,             &
        qsol_d,    zq2m,      s_pblh,   s_plcl,        &
+!!!
+!!! jyg le 08/02/2012
+       s_pblh_x, s_plcl_x,   s_pblh_w, s_plcl_w,      &
+!!!
        s_capCL,   s_oliqCL,  s_cteiCL, s_pblT,        &
        s_therm,   s_trmb1,   s_trmb2,  s_trmb3,       &
@@ -191,5 +211,13 @@
        wfbils,    wfbilo,    flux_t,   flux_u, flux_v,&
        dflux_t,   dflux_q,   zxsnow,                  &
-       zxfluxt,   zxfluxq,   q2m,      flux_q, tke    )
+!jyg<
+!!       zxfluxt,   zxfluxq,   q2m,      flux_q, tke,   &
+       zxfluxt,   zxfluxq,   q2m,      flux_q, tke_x,   &
+!>jyg
+!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
+!!        tke_x,     tke_w                              &
+       wake_dltke                                     &
+!!!
+                        )
 !****************************************************************************************
 ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
@@ -221,4 +249,8 @@
 ! u--------input-R- vitesse u
 ! v--------input-R- vitesse v
+! wake_dlt-input-R- temperatre difference between (w) and (x) (K)
+! wake_dlq-input-R- humidity difference between (w) and (x) (kg/kg)
+!wake_cstar-input-R- wake gust front speed (m/s)
+! wake_s---input-R- wake fractionnal area
 ! ts-------input-R- temperature du sol (en Kelvin)
 ! paprs----input-R- pression a intercouche (Pa)
@@ -239,5 +271,6 @@
 ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
 !                    (orientation positive vers le bas)
-! tke---input/output-R- tke (kg/m**2/s)
+! tke_x---input/output-R- tke in the (x) region (kg/m**2/s)
+! wake_dltke-input/output-R- tke difference between (w) and (x) (kg/m**2/s)
 ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
 ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
@@ -299,14 +332,37 @@
 ! Martin
 
+!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
+!!    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: t_x       ! Température hors poche froide
+!!    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: t_w       ! Température dans la poches froide
+!!    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: q_x       ! 
+!!    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: q_w       ! Pareil pour l'humidité
+    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: wake_dlt  !temperature difference between (w) and (x) (K)
+    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: wake_dlq  !humidity difference between (w) and (x) (K)
+    REAL, DIMENSION(klon),        INTENT(IN)        :: wake_s    ! Fraction de poches froides
+    REAL, DIMENSION(klon),        INTENT(IN)        :: wake_cstar! Vitesse d'expansion des poches froides
+    REAL, DIMENSION(klon),        INTENT(IN)        :: wake_dens
+!!!
+
 ! Input/Output variables
 !****************************************************************************************
     REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ts      ! temperature at surface (K)
+    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: delta_tsurf !surface temperature difference between
+                                                                   !wake and off-wake regions
     REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: alb1    ! albedo in visible SW interval
     REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: alb2    ! albedo in near infra-red SW interval
+!jyg Pourquoi ustar et wstar sont-elles INOUT ?
     REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ustar   ! u* (m/s)
     REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT)   :: wstar   ! w* (m/s)
     REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: u10m    ! u speed at 10m
     REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: v10m    ! v speed at 10m
-    REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke 
+!jyg<
+!!    REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke
+    REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke_x
+!>jyg 
+
+!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
+    REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: wake_dltke ! TKE_w - TKE_x
+!!!
+
 ! Output variables
 !****************************************************************************************
@@ -325,4 +381,10 @@
     REAL, DIMENSION(klon),        INTENT(OUT)       :: zxevap     ! water vapour flux at surface, positiv upwards
     REAL, DIMENSION(klon),        INTENT(OUT)       :: zxtsol     ! temperature at surface, mean for each grid point
+!!! jyg le ???
+    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_t_w      !   !
+    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_q_w      !      !  Tendances dans les poches
+    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_t_x      !   !
+    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_q_x      !      !  Tendances hors des poches
+!!! jyg
     REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat  ! latent flux, mean for each grid point
     REAL, DIMENSION(klon),        INTENT(OUT)       :: zt2m       ! temperature at 2m, mean for each grid point
@@ -340,10 +402,34 @@
     ! coef for turbulent diffusion of U and V (?), mean for each grid point
 
+!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsens_x   ! Flux sensible hors poche
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsens_w   ! Flux sensible dans la poche
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat_x! Flux latent hors poche
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat_w! Flux latent dans la poche
+!!    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_wake_dlt
+!!    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_wake_dlq
+
 ! Output only for diagnostics
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh_x
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh_w
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragm_x
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragm_w
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: kh
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: kh_x
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: kh_w
+!!! 
     REAL, DIMENSION(klon),        INTENT(OUT)       :: slab_wfbils! heat balance at surface only for slab at ocean points
     REAL, DIMENSION(klon),        INTENT(OUT)       :: qsol_d     ! water height in the soil (mm)
     REAL, DIMENSION(klon),        INTENT(OUT)       :: zq2m       ! water vapour at 2m, mean for each grid point
     REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblh     ! height of the planetary boundary layer(HPBL)
+!!! jyg le 08/02/2012
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblh_x   ! height of the PBL in the off-wake region
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblh_w   ! height of the PBL in the wake region
+!!!
     REAL, DIMENSION(klon),        INTENT(OUT)       :: s_plcl     ! condensation level
+!!! jyg le 08/02/2012
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_plcl_x   ! condensation level in the off-wake region
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_plcl_w   ! condensation level in the wake region
+!!!
     REAL, DIMENSION(klon),        INTENT(OUT)       :: s_capCL    ! CAPE of PBL
     REAL, DIMENSION(klon),        INTENT(OUT)       :: s_oliqCL   ! liquid water intergral of PBL
@@ -409,8 +495,10 @@
 ! Other local variables
 !****************************************************************************************
+    INTEGER                            :: iflag_split 
     INTEGER                            :: i, k, nsrf 
     INTEGER                            :: knon, j
     INTEGER                            :: idayref
     INTEGER , DIMENSION(klon)          :: ni
+    REAL                               :: yt1_new
     REAL                               :: zx_alf1, zx_alf2 !valeur ambiante par extrapola
     REAL                               :: amn, amx
@@ -419,5 +507,5 @@
     REAL, DIMENSION(klon)              :: yts, yrugos, ypct, yz0_new
     REAL, DIMENSION(klon)              :: yalb, yalb1, yalb2
-    REAL, DIMENSION(klon)              :: yu1, yv1,ytoto
+    REAL, DIMENSION(klon)              :: yu1, yv1
     REAL, DIMENSION(klon)              :: ysnow, yqsurf, yagesno, yqsol
     REAL, DIMENSION(klon)              :: yrain_f, ysnow_f
@@ -474,6 +562,194 @@
     LOGICAL, PARAMETER                 :: zxli=.FALSE. ! utiliser un jeu de fonctions simples
     LOGICAL, PARAMETER                 :: check=.FALSE.
-    REAL, DIMENSION(klon)              :: Kech_h       ! Coefficient d'echange pour l'energie
+
+!!! nrlmd le 02/05/2011
+!!! jyg le 07/02/2012
+    REAL, DIMENSION(klon)              :: ywake_s, ywake_cstar, ywake_dens
+!!!
+    REAL, DIMENSION(klon,klev+1)       :: ytke_x, ytke_w
+    REAL, DIMENSION(klon,klev+1)       :: ywake_dltke
+    REAL, DIMENSION(klon,klev)         :: yu_x, yv_x, yu_w, yv_w
+    REAL, DIMENSION(klon,klev)         :: yt_x, yq_x, yt_w, yq_w
+    REAL, DIMENSION(klon,klev)         :: ycoefh_x, ycoefm_x, ycoefh_w, ycoefm_w
+    REAL, DIMENSION(klon,klev)         :: ycoefq_x, ycoefq_w
+    REAL, DIMENSION(klon)              :: ycdragh_x, ycdragm_x, ycdragh_w, ycdragm_w
+    REAL, DIMENSION(klon)              :: AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x
+    REAL, DIMENSION(klon)              :: AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w
+    REAL, DIMENSION(klon)              :: AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x
+    REAL, DIMENSION(klon)              :: AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w
+    REAL, DIMENSION(klon)              :: y_flux_t1_x, y_flux_q1_x, y_flux_t1_w, y_flux_q1_w
+    REAL, DIMENSION(klon)              :: y_flux_u1_x, y_flux_v1_x, y_flux_u1_w, y_flux_v1_w
+    REAL, DIMENSION(klon,klev)         :: y_flux_t_x, y_flux_q_x, y_flux_t_w, y_flux_q_w
+    REAL, DIMENSION(klon,klev)         :: y_flux_u_x, y_flux_v_x, y_flux_u_w, y_flux_v_w
+    REAL, DIMENSION(klon)              :: yfluxlat_x, yfluxlat_w
+    REAL, DIMENSION(klon,klev)         :: y_d_t_x, y_d_q_x, y_d_t_w, y_d_q_w
+    REAL, DIMENSION(klon,klev)         :: y_d_t_diss_x, y_d_t_diss_w
+    REAL, DIMENSION(klon,klev)         :: d_t_diss_x, d_t_diss_w
+    REAL, DIMENSION(klon,klev)         :: y_d_u_x, y_d_v_x, y_d_u_w, y_d_v_w
+    REAL, DIMENSION(klon, klev, nbsrf) :: flux_t_x, flux_q_x, flux_t_w, flux_q_w
+    REAL, DIMENSION(klon, klev, nbsrf) :: flux_u_x, flux_v_x, flux_u_w, flux_v_w
+    REAL, DIMENSION(klon, nbsrf)       :: fluxlat_x, fluxlat_w
+    REAL, DIMENSION(klon, klev)        :: zxfluxt_x, zxfluxq_x, zxfluxt_w, zxfluxq_w
+    REAL, DIMENSION(klon, klev)        :: zxfluxu_x, zxfluxv_x, zxfluxu_w, zxfluxv_w
+    REAL                               :: zx_qs_surf, zcor_surf, zdelta_surf
+    REAL, DIMENSION(klon)              :: ytsurf_th, yqsatsurf
+    REAL, DIMENSION(klon)              :: ybeta
+    REAL, DIMENSION(klon, klev)        :: d_u_x
+    REAL, DIMENSION(klon, klev)        :: d_u_w
+    REAL, DIMENSION(klon, klev)        :: d_v_x
+    REAL, DIMENSION(klon, klev)        :: d_v_w 
+
+    REAL, DIMENSION(klon,klev)         :: CcoefH, CcoefQ, DcoefH, DcoefQ
+    REAL, DIMENSION(klon,klev)         :: CcoefU, CcoefV, DcoefU, DcoefV
+    REAL, DIMENSION(klon,klev)         :: CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x
+    REAL, DIMENSION(klon,klev)         :: CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w
+    REAL, DIMENSION(klon,klev)         :: CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x
+    REAL, DIMENSION(klon,klev)         :: CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w
+    REAL, DIMENSION(klon,klev)         :: Kcoef_hq, Kcoef_m, gama_h, gama_q
+    REAL, DIMENSION(klon,klev)         :: Kcoef_hq_x, Kcoef_m_x, gama_h_x, gama_q_x
+    REAL, DIMENSION(klon,klev)         :: Kcoef_hq_w, Kcoef_m_w, gama_h_w, gama_q_w
+    REAL, DIMENSION(klon)              :: alf_1, alf_2, alf_1_x, alf_2_x, alf_1_w, alf_2_w
+!!!
+!!!jyg le 08/02/2012
+    REAL, DIMENSION(klon, nbsrf)       :: windsp
+!
+    REAL, DIMENSION(klon, nbsrf)       :: t2m_x
+    REAL, DIMENSION(klon, nbsrf)       :: q2m_x
+    REAL, DIMENSION(klon)              :: rh2m_x
+    REAL, DIMENSION(klon)              :: qsat2m_x
+    REAL, DIMENSION(klon, nbsrf)       :: u10m_x
+    REAL, DIMENSION(klon, nbsrf)       :: v10m_x
+    REAL, DIMENSION(klon, nbsrf)       :: ustar_x
+    REAL, DIMENSION(klon, nbsrf)       :: wstar_x
+!              
+    REAL, DIMENSION(klon, nbsrf)       :: pblh_x
+    REAL, DIMENSION(klon, nbsrf)       :: plcl_x
+    REAL, DIMENSION(klon, nbsrf)       :: capCL_x
+    REAL, DIMENSION(klon, nbsrf)       :: oliqCL_x
+    REAL, DIMENSION(klon, nbsrf)       :: cteiCL_x
+    REAL, DIMENSION(klon, nbsrf)       :: pblt_x
+    REAL, DIMENSION(klon, nbsrf)       :: therm_x
+    REAL, DIMENSION(klon, nbsrf)       :: trmb1_x
+    REAL, DIMENSION(klon, nbsrf)       :: trmb2_x
+    REAL, DIMENSION(klon, nbsrf)       :: trmb3_x
+!
+    REAL, DIMENSION(klon, nbsrf)       :: t2m_w
+    REAL, DIMENSION(klon, nbsrf)       :: q2m_w
+    REAL, DIMENSION(klon)              :: rh2m_w
+    REAL, DIMENSION(klon)              :: qsat2m_w
+    REAL, DIMENSION(klon, nbsrf)       :: u10m_w
+    REAL, DIMENSION(klon, nbsrf)       :: v10m_w
+    REAL, DIMENSION(klon, nbsrf)       :: ustar_w
+    REAL, DIMENSION(klon, nbsrf)       :: wstar_w
+!                           
+    REAL, DIMENSION(klon, nbsrf)       :: pblh_w
+    REAL, DIMENSION(klon, nbsrf)       :: plcl_w
+    REAL, DIMENSION(klon, nbsrf)       :: capCL_w
+    REAL, DIMENSION(klon, nbsrf)       :: oliqCL_w
+    REAL, DIMENSION(klon, nbsrf)       :: cteiCL_w
+    REAL, DIMENSION(klon, nbsrf)       :: pblt_w
+    REAL, DIMENSION(klon, nbsrf)       :: therm_w
+    REAL, DIMENSION(klon, nbsrf)       :: trmb1_w
+    REAL, DIMENSION(klon, nbsrf)       :: trmb2_w
+    REAL, DIMENSION(klon, nbsrf)       :: trmb3_w
+!
+    REAL, DIMENSION(klon)       :: yt2m_x
+    REAL, DIMENSION(klon)       :: yq2m_x
+    REAL, DIMENSION(klon)       :: yt10m_x
+    REAL, DIMENSION(klon)       :: yq10m_x
+    REAL, DIMENSION(klon)       :: yu10m_x
+    REAL, DIMENSION(klon)       :: yv10m_x
+    REAL, DIMENSION(klon)       :: yustar_x
+    REAL, DIMENSION(klon)       :: ywstar_x
+!              
+    REAL, DIMENSION(klon)       :: ypblh_x
+    REAL, DIMENSION(klon)       :: ylcl_x
+    REAL, DIMENSION(klon)       :: ycapCL_x
+    REAL, DIMENSION(klon)       :: yoliqCL_x
+    REAL, DIMENSION(klon)       :: ycteiCL_x
+    REAL, DIMENSION(klon)       :: ypblt_x
+    REAL, DIMENSION(klon)       :: ytherm_x
+    REAL, DIMENSION(klon)       :: ytrmb1_x
+    REAL, DIMENSION(klon)       :: ytrmb2_x
+    REAL, DIMENSION(klon)       :: ytrmb3_x
+!
+    REAL, DIMENSION(klon)       :: yt2m_w
+    REAL, DIMENSION(klon)       :: yq2m_w
+    REAL, DIMENSION(klon)       :: yt10m_w
+    REAL, DIMENSION(klon)       :: yq10m_w
+    REAL, DIMENSION(klon)       :: yu10m_w
+    REAL, DIMENSION(klon)       :: yv10m_w
+    REAL, DIMENSION(klon)       :: yustar_w
+    REAL, DIMENSION(klon)       :: ywstar_w
+!                       
+    REAL, DIMENSION(klon)       :: ypblh_w
+    REAL, DIMENSION(klon)       :: ylcl_w
+    REAL, DIMENSION(klon)       :: ycapCL_w
+    REAL, DIMENSION(klon)       :: yoliqCL_w
+    REAL, DIMENSION(klon)       :: ycteiCL_w
+    REAL, DIMENSION(klon)       :: ypblt_w
+    REAL, DIMENSION(klon)       :: ytherm_w
+    REAL, DIMENSION(klon)       :: ytrmb1_w
+    REAL, DIMENSION(klon)       :: ytrmb2_w
+    REAL, DIMENSION(klon)       :: ytrmb3_w
+!
+    REAL, DIMENSION(klon)              :: uzon_x, vmer_x
+    REAL, DIMENSION(klon)              :: zgeo1_x, tair1_x, qair1_x, tairsol_x
+!
+    REAL, DIMENSION(klon)              :: uzon_w, vmer_w
+    REAL, DIMENSION(klon)              :: zgeo1_w, tair1_w, qair1_w, tairsol_w
+
+!!! jyg le 25/03/2013
+!!    Variables intermediaires pour le raccord des deux colonnes à la surface
+    REAL   ::   dd_Ch
+    REAL   ::   dd_Cm
+    REAL   ::   dd_Kh
+    REAL   ::   dd_Km
+    REAL   ::   dd_u 
+    REAL   ::   dd_v 
+    REAL   ::   dd_t 
+    REAL   ::   dd_q 
+    REAL   ::   dd_AH
+    REAL   ::   dd_AQ
+    REAL   ::   dd_AU
+    REAL   ::   dd_AV
+    REAL   ::   dd_BH
+    REAL   ::   dd_BQ
+    REAL   ::   dd_BU
+    REAL   ::   dd_BV
+
+    REAL   ::   dd_KHp
+    REAL   ::   dd_KQp
+    REAL   ::   dd_KUp
+    REAL   ::   dd_KVp
+
+!!!
+!!! nrlmd le 13/06/2011
+    REAL, DIMENSION(klon)              :: y_delta_flux_t1, y_delta_flux_q1, y_delta_flux_u1, y_delta_flux_v1
+    REAL, DIMENSION(klon)              :: y_delta_tsurf,delta_coef,tau_eq
+    REAL, PARAMETER                    :: facteur=2./sqrt(3.14)
+    REAL, PARAMETER                    :: effusivity=2000.
+    REAL, DIMENSION(klon)              :: ytsurf_th_x,ytsurf_th_w,yqsatsurf_x,yqsatsurf_w
+    REAL, DIMENSION(klon)              :: ydtsurf_th
+    REAL                               :: zdelta_surf_x,zdelta_surf_w,zx_qs_surf_x,zx_qs_surf_w
+    REAL                               :: zcor_surf_x,zcor_surf_w
+    REAL                               :: mod_wind_x, mod_wind_w
+    REAL                               :: rho1
+    REAL, DIMENSION(klon)              :: Kech_h           ! Coefficient d'echange pour l'energie
+    REAL, DIMENSION(klon)              :: Kech_h_x, Kech_h_w 
+    REAL, DIMENSION(klon)              :: Kech_m
+    REAL, DIMENSION(klon)              :: Kech_m_x, Kech_m_w 
+    REAL, DIMENSION(klon)              :: yts_x,yts_w
+    REAL, DIMENSION(klon)              :: Kech_Hp, Kech_H_xp, Kech_H_wp
+    REAL, DIMENSION(klon)              :: Kech_Qp, Kech_Q_xp, Kech_Q_wp
+    REAL, DIMENSION(klon)              :: Kech_Up, Kech_U_xp, Kech_U_wp
+    REAL, DIMENSION(klon)              :: Kech_Vp, Kech_V_xp, Kech_V_wp
+
     REAL                               :: vent
+
+
+
+
+!!!
 
 ! For debugging with IOIPSL
@@ -514,8 +790,10 @@
 
 !****************************************************************************************
-
 ! End of declarations
 !****************************************************************************************
 
+      IF (prt_level >=10) print *,' -> pbl_surface, itap ',itap
+!
+      iflag_split = mod(iflag_pbl_split,2)
 
 !****************************************************************************************
@@ -529,5 +807,5 @@
       
        ! Initialize ok_flux_surf (for 1D model)
-       if (klon>1) ok_flux_surf=.FALSE.
+       if (klon_glo>1) ok_flux_surf=.FALSE.
        
        ! Initilize debug IO
@@ -573,27 +851,45 @@
 !****************************************************************************************
 ! 2) Initialization to zero 
-!    Done for all local variables that will be compressed later
-!    and argument with INTENT(OUT)
-!****************************************************************************************
-    cdragh = 0.0  ; cdragm = 0.0     ; dflux_t = 0.0   ; dflux_q = 0.0
-    ypct = 0.0    ; yts = 0.0        ; ysnow = 0.0
-    zv1 = 0.0     ; yqsurf = 0.0     ; yalb1 = 0.0     ; yalb2 = 0.0    
-    yrain_f = 0.0 ; ysnow_f = 0.0    ; yfder = 0.0     ; ysolsw = 0.0    
-    ysollw = 0.0  ; yrugos = 0.0     ; yu1 = 0.0    
-    yv1 = 0.0     ; ypaprs = 0.0     ; ypplay = 0.0
-    ydelp = 0.0   ; yu = 0.0         ; yv = 0.0        ; yt = 0.0         
-    yq = 0.0      ; y_dflux_t = 0.0  ; y_dflux_q = 0.0 
-    yrugoro = 0.0 ; ywindsp = 0.0   
-    d_ts = 0.0    ; yfluxlat=0.0     ; flux_t = 0.0    ; flux_q = 0.0     
-    flux_u = 0.0  ; flux_v = 0.0     ; d_t = 0.0       ; d_q = 0.0      
-    d_t_diss= 0.0 ;d_u = 0.0     ; d_v = 0.0        ; yqsol = 0.0    
-    ytherm = 0.0  ; ytke=0.
-    ! Martin
-    ysnowhgt = 0.0; yqsnow = 0.0     ; yrunoff = 0.0   ; ytoice =0.0
-    yalb3_new = 0.0  ; ysissnow = 0.0  ; ysollwd = 0.0
-    ypphi = 0.0   ; ycldt = 0.0      ; yrmu0 = 0.0
-    ! Martin
-    
-    tke(:,:,is_ave)=0.
+!****************************************************************************************
+!
+! 2a) Initialization of all argument variables with INTENT(OUT)
+!****************************************************************************************
+ lwdown_m(:)=0.
+ cdragh(:)=0. ; cdragm(:)=0.
+ zu1(:)=0. ; zv1(:)=0.
+ alb1_m(:)=0. ; alb2_m(:)=0. ; alb3_lic(:)=0.
+ zxsens(:)=0. ; zxevap(:)=0. ; zxtsol(:)=0.
+ d_t_w(:,:)=0. ; d_q_w(:,:)=0. ; d_t_x(:,:)=0. ; d_q_x(:,:)=0.
+ zxfluxlat(:)=0.
+ zt2m(:)=0. ; zq2m(:)=0. ; qsat2m(:)=0. ; rh2m(:)=0.
+ d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0.
+ zcoefh(:,:,:)=0. ; zcoefm(:,:,:)=0.
+ zxsens_x(:)=0. ; zxsens_w(:)=0. ; zxfluxlat_x(:)=0. ; zxfluxlat_w(:)=0.
+ cdragh_x(:)=0. ; cdragh_w(:)=0. ; cdragm_x(:)=0. ; cdragm_w(:)=0.
+ kh(:)=0. ; kh_x(:)=0. ; kh_w(:)=0.
+ slab_wfbils(:)=0.
+ qsol_d(:)=0.
+ s_pblh(:)=0. ; s_pblh_x(:)=0. ; s_pblh_w(:)=0.
+ s_plcl(:)=0. ; s_plcl_x(:)=0. ; s_plcl_w(:)=0.
+ s_capCL(:)=0. ; s_oliqCL(:)=0. ; s_cteiCL(:)=0. ; s_pblT(:)=0.
+ s_therm(:)=0.
+ s_trmb1(:)=0. ; s_trmb2(:)=0. ; s_trmb3(:)=0.
+ zxrugs(:)=0. ; zustar(:)=0.
+ zu10m(:)=0. ; zv10m(:)=0.
+ fder_print(:)=0.
+ zxqsurf(:)=0.
+ zxfluxu(:,:)=0. ; zxfluxv(:,:)=0.
+ rugos_d(:,:)=0. ; agesno_d(:,:)=0.
+ solsw(:,:)=0. ; sollw(:,:)=0.
+ d_ts(:,:)=0.
+ evap_d(:,:)=0.
+ fluxlat(:,:)=0.
+ wfbils(:,:)=0. ; wfbilo(:,:)=0.
+ flux_t(:,:,:)=0. ; flux_q(:,:,:)=0. ; flux_u(:,:,:)=0. ; flux_v(:,:,:)=0.
+ dflux_t(:)=0. ; dflux_q(:)=0.
+ zxsnow(:)=0.
+ zxfluxt(:,:)=0. ; zxfluxq(:,:)=0.
+ qsnow(:)=0. ; snowhgt(:)=0. ; to_ice(:)=0. ; sissnow(:)=0.
+ runoff(:)=0.
     IF (iflag_pbl<20.or.iflag_pbl>=30) THEN
        zcoefh(:,:,:) = 0.0
@@ -605,8 +901,139 @@
       zcoefh(:,:,is_ave)=0.
     ENDIF
+!!
+!  The components "is_ave" of tke_x and wake_deltke are "OUT" variables
+!jyg<
+!!    tke(:,:,is_ave)=0.
+    tke_x(:,:,is_ave)=0.
+    wake_dltke(:,:,is_ave)=0.
+!>jyg
+!!! jyg le 23/02/2013
+    t2m(:,:)       = 999999.     ! t2m and q2m are meaningfull only over sub-surfaces
+    q2m(:,:)       = 999999.     ! actually present in the grid cell.
+!!!
+    rh2m(:) = 0. ; qsat2m(:) = 0.
+!!!
+!!! jyg le 10/02/2012
+    rh2m_x(:) = 0. ; qsat2m_x(:) = 0. ; rh2m_w(:) = 0. ; qsat2m_w(:) = 0.
+!!!
+
+! 2b) Initialization of all local variables that will be compressed later
+!****************************************************************************************
+!!    cdragh = 0.0  ; cdragm = 0.0     ; dflux_t = 0.0   ; dflux_q = 0.0
+    ypct = 0.0    ; yts = 0.0        ; ysnow = 0.0
+!!    zv1 = 0.0     ; yqsurf = 0.0     ; yalb1 = 0.0     ; yalb2 = 0.0    
+    yqsurf = 0.0  ; yalb1 = 0.0      ; yalb2 = 0.0    
+    yrain_f = 0.0 ; ysnow_f = 0.0    ; yfder = 0.0     ; ysolsw = 0.0    
+    ysollw = 0.0  ; yrugos = 0.0     ; yu1 = 0.0    
+    yv1 = 0.0     ; ypaprs = 0.0     ; ypplay = 0.0
+    ydelp = 0.0   ; yu = 0.0         ; yv = 0.0        ; yt = 0.0         
+    yq = 0.0      ; y_dflux_t = 0.0  ; y_dflux_q = 0.0 
+    yrugoro = 0.0 ; ywindsp = 0.0   
+!!    d_ts = 0.0    ; yfluxlat=0.0     ; flux_t = 0.0    ; flux_q = 0.0     
+    yfluxlat=0.0 
+!!    flux_u = 0.0  ; flux_v = 0.0     ; d_t = 0.0       ; d_q = 0.0      
+!!    d_t_diss= 0.0 ;d_u = 0.0     ; d_v = 0.0 
+    yqsol = 0.0    
+    ytherm = 0.0  ; ytke=0.
+    ! Martin
+    ysnowhgt = 0.0; yqsnow = 0.0     ; yrunoff = 0.0   ; ytoice =0.0
+    yalb3_new = 0.0  ; ysissnow = 0.0  ; ysollwd = 0.0
+    ypphi = 0.0   ; ycldt = 0.0      ; yrmu0 = 0.0
+    ! Martin
+
+!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
+    ytke_x=0.     ; ytke_w=0.        ; ywake_dltke=0.
+    y_d_t_x=0.    ; y_d_t_w=0.       ; y_d_q_x=0.      ; y_d_q_w=0.
+!!    d_t_w=0.      ; d_q_w=0.         
+!!    d_t_x=0.      ; d_q_x=0.
+!!    d_wake_dlt=0.    ; d_wake_dlq=0.
+    yfluxlat_x=0. ; yfluxlat_w=0.
+    ywake_s=0.    ; ywake_cstar=0.   ;ywake_dens=0.
+!!!
+!!! nrlmd le 13/06/2011
+    tau_eq=0.     ; delta_coef=0.
+    y_delta_flux_t1=0.
+    ydtsurf_th=0.
+    yts_x=0.      ; yts_w=0.
+    y_delta_tsurf=0.
+!!!
     ytsoil = 999999. 
 
-    rh2m(:)        = 0.
-    qsat2m(:)      = 0.
+
+! 2c) Initialization of all local variables computed within the subsurface loop and used later on
+!****************************************************************************************
+    d_t_diss_x(:,:) = 0. ;        d_t_diss_w(:,:) = 0.
+    d_u_x(:,:)=0. ;               d_u_w(:,:)=0. 
+    d_v_x(:,:)=0. ;               d_v_w(:,:)=0.
+    flux_t_x(:,:,:)=0. ;          flux_t_w(:,:,:)=0. 
+    flux_q_x(:,:,:)=0. ;          flux_q_w(:,:,:)=0.
+!
+!jyg<
+    flux_u_x(:,:,:)=0. ;          flux_u_w(:,:,:)=0.
+    flux_v_x(:,:,:)=0. ;          flux_v_w(:,:,:)=0.
+    fluxlat_x(:,:)=0. ;           fluxlat_w(:,:)=0. 
+!>jyg
+!
+!jyg<
+! pblh,plcl,capCL,cteiCL ... are meaningfull only over sub-surfaces
+! actually present in the grid cell  ==> value set to 999999.
+!                           
+!jyg<
+       ustar(:,:)   = 999999.
+       wstar(:,:)   = 999999.
+       windsp(:,:)  = SQRT(u10m(:,:)**2 + v10m(:,:)**2 )
+       u10m(:,:)    = 999999. 
+       v10m(:,:)    = 999999. 
+!>jyg
+!
+       pblh(:,:)   = 999999.        ! Hauteur de couche limite
+       plcl(:,:)   = 999999.        ! Niveau de condensation de la CLA
+       capCL(:,:)  = 999999.        ! CAPE de couche limite
+       oliqCL(:,:) = 999999.        ! eau_liqu integree de couche limite
+       cteiCL(:,:) = 999999.        ! cloud top instab. crit. couche limite
+       pblt(:,:)   = 999999.        ! T a la Hauteur de couche limite
+       therm(:,:)  = 999999.
+       trmb1(:,:)  = 999999.        ! deep_cape
+       trmb2(:,:)  = 999999.        ! inhibition 
+       trmb3(:,:)  = 999999.        ! Point Omega
+!
+       t2m_x(:,:)    = 999999. 
+       q2m_x(:,:)    = 999999. 
+       ustar_x(:,:)   = 999999.
+       wstar_x(:,:)   = 999999.
+       u10m_x(:,:)   = 999999. 
+       v10m_x(:,:)   = 999999. 
+!                           
+       pblh_x(:,:)   = 999999.      ! Hauteur de couche limite
+       plcl_x(:,:)   = 999999.      ! Niveau de condensation de la CLA
+       capCL_x(:,:)  = 999999.      ! CAPE de couche limite
+       oliqCL_x(:,:) = 999999.      ! eau_liqu integree de couche limite
+       cteiCL_x(:,:) = 999999.      ! cloud top instab. crit. couche limite
+       pblt_x(:,:)   = 999999.      ! T a la Hauteur de couche limite
+       therm_x(:,:)  = 999999.      
+       trmb1_x(:,:)  = 999999.      ! deep_cape
+       trmb2_x(:,:)  = 999999.      ! inhibition 
+       trmb3_x(:,:)  = 999999.      ! Point Omega
+!
+       t2m_w(:,:)    = 999999. 
+       q2m_w(:,:)    = 999999. 
+       ustar_w(:,:)   = 999999.
+       wstar_w(:,:)   = 999999.
+       u10m_w(:,:)   = 999999. 
+       v10m_w(:,:)   = 999999. 
+                           
+       pblh_w(:,:)   = 999999.      ! Hauteur de couche limite
+       plcl_w(:,:)   = 999999.      ! Niveau de condensation de la CLA
+       capCL_w(:,:)  = 999999.      ! CAPE de couche limite
+       oliqCL_w(:,:) = 999999.      ! eau_liqu integree de couche limite
+       cteiCL_w(:,:) = 999999.      ! cloud top instab. crit. couche limite
+       pblt_w(:,:)   = 999999.      ! T a la Hauteur de couche limite
+       therm_w(:,:)  = 999999.      
+       trmb1_w(:,:)  = 999999.      ! deep_cape
+       trmb2_w(:,:)  = 999999.      ! inhibition 
+       trmb3_w(:,:)  = 999999.      ! Point Omega
+!!!      
+!
+!!!
 !****************************************************************************************
 ! 3) - Calculate pressure thickness of each layer
@@ -699,9 +1126,10 @@
 ! 4) Loop over different surfaces
 !
-! Only points containing a fraction of the sub surface will be threated.
+! Only points containing a fraction of the sub surface will be treated.
 ! 
 !****************************************************************************************
    
     loop_nbsrf: DO nsrf = 1, nbsrf
+       IF (prt_level >=10) print *,' Loop nsrf ',nsrf
 
 ! Search for index(ni) and size(knon) of domaine to treat
@@ -714,4 +1142,11 @@
           ENDIF
        ENDDO
+
+!!! jyg le 19/08/2012
+!       IF (knon <= 0) THEN
+!         IF (prt_level >= 10) print *,' no grid point for nsrf= ',nsrf
+!         cycle loop_nbsrf
+!       ENDIF
+!!!
 
        ! write index, with IOIPSL
@@ -752,5 +1187,8 @@
           yv1(j)     = v(i,1)
           ypaprs(j,klev+1) = paprs(i,klev+1)
-          ywindsp(j) = SQRT(u10m(i,nsrf)**2 + v10m(i,nsrf)**2 )
+!jyg<
+!!          ywindsp(j) = SQRT(u10m(i,nsrf)**2 + v10m(i,nsrf)**2 )
+          ywindsp(j) = windsp(i,nsrf)
+!>jyg
           ! Martin
           yzsig(j)   = zsig(i)
@@ -758,4 +1196,7 @@
           yrmu0(j)   = rmu0(i)
           ! Martin
+!!! nrlmd le 13/06/2011
+          y_delta_tsurf(j)=delta_tsurf(i,nsrf)
+!!!
        END DO
 
@@ -766,5 +1207,14 @@
              ypplay(j,k) = pplay(i,k)
              ydelp(j,k)  = delp(i,k)
-             ytke(j,k)   = tke(i,k,nsrf)
+          ENDDO
+       ENDDO
+!!! jyg le 07/02/2012 et le 10/04/2013
+        DO k = 1, klev
+          DO j = 1, knon
+             i = ni(j)
+!jyg<
+!!             ytke(j,k)   = tke(i,k,nsrf)
+             ytke(j,k)   = tke_x(i,k,nsrf)
+!>jyg
              yu(j,k) = u(i,k)
              yv(j,k) = v(i,k)
@@ -772,6 +1222,56 @@
              yq(j,k) = q(i,k)
           ENDDO
-       ENDDO
-
+        ENDDO
+!
+       IF (iflag_split .eq.1) THEN
+!!! nrlmd le 02/05/2011
+        DO k = 1, klev
+          DO j = 1, knon
+             i = ni(j)
+             yu_x(j,k) = u(i,k)
+             yv_x(j,k) = v(i,k)
+             yt_x(j,k) = t(i,k)-wake_s(i)*wake_dlt(i,k)
+             yq_x(j,k) = q(i,k)-wake_s(i)*wake_dlq(i,k)
+             yu_w(j,k) = u(i,k)
+             yv_w(j,k) = v(i,k)
+             yt_w(j,k) = t(i,k)+(1.-wake_s(i))*wake_dlt(i,k)
+             yq_w(j,k) = q(i,k)+(1.-wake_s(i))*wake_dlq(i,k)
+!!!
+          ENDDO
+        ENDDO
+!!! nrlmd le 02/05/2011
+        DO k = 1, klev+1
+          DO j = 1, knon
+             i = ni(j)
+!jyg<
+!!             ytke_x(j,k) = tke(i,k,nsrf)-wake_s(i)*wake_dltke(i,k,nsrf)
+!!             ytke_w(j,k) = tke(i,k,nsrf)+(1.-wake_s(i))*wake_dltke(i,k,nsrf)
+!!             ywake_dltke(j,k) = wake_dltke(i,k,nsrf)
+!!             ytke(j,k)     = tke(i,k,nsrf)
+!
+             ytke_x(j,k)      = tke_x(i,k,nsrf)
+             ytke(j,k)        = tke_x(i,k,nsrf)+wake_s(i)*wake_dltke(i,k,nsrf)
+             ytke_w(j,k)      = tke_x(i,k,nsrf)+wake_dltke(i,k,nsrf)
+             ywake_dltke(j,k) = wake_dltke(i,k,nsrf)
+!>jyg
+          ENDDO
+        ENDDO
+!!!
+!!! jyg le 07/02/2012
+        DO j = 1, knon
+          i = ni(j)
+          ywake_s(j)=wake_s(i)
+          ywake_cstar(j)=wake_cstar(i)
+          ywake_dens(j)=wake_dens(i)
+        ENDDO
+!!!
+!!! nrlmd le 13/06/2011
+        DO j=1,knon
+         yts_x(j)=yts(j)-ywake_s(j)*y_delta_tsurf(j)
+         yts_w(j)=yts(j)+(1.-ywake_s(j))*y_delta_tsurf(j)
+        ENDDO
+!!!
+       ENDIF  ! (iflag_split .eq.1)
+!!!
        DO k = 1, nsoilmx
           DO j = 1, knon
@@ -794,5 +1294,9 @@
 !****************************************************************************************
 
-       CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.0) THEN
+!!!
+!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
+        CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
             yu(:,1), yv(:,1), yt(:,1), yq(:,1), &
             yts, yqsurf, yrugos, &
@@ -810,15 +1314,64 @@
       ENDDO
      ENDIF
-
-
-!****************************************************************************************
-! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefm et ycoefm.
-!
-!****************************************************************************************
-
-       CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
+        IF (prt_level >=10) print *,'clcdrag -> ycdragh ', ycdragh
+       ELSE  !(iflag_split .eq.0)
+        CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
+            yu_x(:,1), yv_x(:,1), yt_x(:,1), yq_x(:,1), &
+            yts_x, yqsurf, yrugos, &
+            ycdragm_x, ycdragh_x )
+! --- special Dice. JYG+MPL 25112013
+        IF (ok_prescr_ust) then
+         DO i = 1, knon
+          print *,'ycdragm_x avant=',ycdragm_x(i)
+          vent= sqrt(yu_x(i,1)*yu_x(i,1)+yv_x(i,1)*yv_x(i,1))
+          ycdragm_x(i) = ust*ust/(1.+vent)/vent
+          print *,'ycdragm_x ust yu yv apres=',ycdragm_x(i),ust,yu_x(i,1),yv_x(i,1)
+         ENDDO
+        ENDIF
+        IF (prt_level >=10) print *,'clcdrag -> ycdragh_x ', ycdragh_x
+!
+        CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
+            yu_w(:,1), yv_w(:,1), yt_w(:,1), yq_w(:,1), &
+            yts_w, yqsurf, yrugos, &
+            ycdragm_w, ycdragh_w )
+! --- special Dice. JYG+MPL 25112013
+        IF (ok_prescr_ust) then
+         DO i = 1, knon
+          print *,'ycdragm_w avant=',ycdragm_w(i)
+          vent= sqrt(yu_w(i,1)*yu_w(i,1)+yv_w(i,1)*yv_w(i,1))
+          ycdragm_w(i) = ust*ust/(1.+vent)/vent
+          print *,'ycdragm_w ust yu yv apres=',ycdragm_w(i),ust,yu_w(i,1),yv_w(i,1)
+         ENDDO
+        ENDIF
+        IF (prt_level >=10) print *,'clcdrag -> ycdragh_w ', ycdragh_w
+!!!
+       ENDIF  ! (iflag_split .eq.0)
+!!!
+       
+
+!****************************************************************************************
+! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefh et ycoefm.
+!
+!****************************************************************************************
+
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.0) THEN
+!!!
+!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
+      IF (prt_level >=10) THEN
+      print *,' args coef_diff_turb: yu ',  yu  
+      print *,' args coef_diff_turb: yv ',  yv  
+      print *,' args coef_diff_turb: yq ',  yq  
+      print *,' args coef_diff_turb: yt ',  yt  
+      print *,' args coef_diff_turb: yts ', yts  
+      print *,' args coef_diff_turb: yrugos ', yrugos  
+      print *,' args coef_diff_turb: yqsurf ', yqsurf  
+      print *,' args coef_diff_turb: ycdragm ', ycdragm
+      print *,' args coef_diff_turb: ycdragh ', ycdragh
+      print *,' args coef_diff_turb: ytke ', ytke
+       ENDIF
+        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
             ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, ycdragm, &
             ycoefm, ycoefh, ytke)
-
        IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
 ! In this case, coef_diff_turb is called for the Cd only
@@ -831,4 +1384,73 @@
        ENDDO
        ENDIF
+        IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh ',ycoefh
+!
+       ELSE  !(iflag_split .eq.0)
+      IF (prt_level >=10) THEN
+      print *,' args coef_diff_turb: yu_x ',  yu_x  
+      print *,' args coef_diff_turb: yv_x ',  yv_x  
+      print *,' args coef_diff_turb: yq_x ',  yq_x  
+      print *,' args coef_diff_turb: yt_x ',  yt_x  
+      print *,' args coef_diff_turb: yts_x ', yts_x  
+      print *,' args coef_diff_turb: yrugos ', yrugos  
+      print *,' args coef_diff_turb: yqsurf ', yqsurf  
+      print *,' args coef_diff_turb: ycdragm_x ', ycdragm_x
+      print *,' args coef_diff_turb: ycdragh_x ', ycdragh_x
+      print *,' args coef_diff_turb: ytke_x ', ytke_x
+       ENDIF
+        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
+            ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yrugos, yqsurf, ycdragm_x, &
+            ycoefm_x, ycoefh_x, ytke_x)
+       IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
+! In this case, coef_diff_turb is called for the Cd only
+       DO k = 2, klev
+          DO j = 1, knon
+             i = ni(j)
+             ycoefh_x(j,k)   = zcoefh(i,k,nsrf)
+             ycoefm_x(j,k)   = zcoefm(i,k,nsrf)
+          ENDDO
+       ENDDO
+       ENDIF
+        IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_x ',ycoefh_x
+!
+      IF (prt_level >=10) THEN
+      print *,' args coef_diff_turb: yu_w ',  yu_w  
+      print *,' args coef_diff_turb: yv_w ',  yv_w  
+      print *,' args coef_diff_turb: yq_w ',  yq_w  
+      print *,' args coef_diff_turb: yt_w ',  yt_w  
+      print *,' args coef_diff_turb: yts_w ', yts_w  
+      print *,' args coef_diff_turb: yrugos ', yrugos  
+      print *,' args coef_diff_turb: yqsurf ', yqsurf  
+      print *,' args coef_diff_turb: ycdragm_w ', ycdragm_w
+      print *,' args coef_diff_turb: ycdragh_w ', ycdragh_w
+      print *,' args coef_diff_turb: ytke_w ', ytke_w
+       ENDIF
+        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
+            ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yrugos, yqsurf, ycdragm_w, &
+            ycoefm_w, ycoefh_w, ytke_w)
+       IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
+! In this case, coef_diff_turb is called for the Cd only
+       DO k = 2, klev
+          DO j = 1, knon
+             i = ni(j)
+             ycoefh_w(j,k)   = zcoefh(i,k,nsrf)
+             ycoefm_w(j,k)   = zcoefm(i,k,nsrf)
+          ENDDO
+       ENDDO
+       ENDIF
+        IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_w ',ycoefh_w
+!
+!!!jyg le 10/04/2013
+!!   En attendant de traiter le transport des traceurs dans les poches froides, formule
+!!   arbitraire pour ycoefh et ycoefm
+      DO k = 2,klev
+        DO j = 1,knon
+         ycoefh(j,k) = ycoefh_x(j,k) + ywake_s(j)*(ycoefh_w(j,k) - ycoefh_x(j,k))
+         ycoefm(j,k) = ycoefm_x(j,k) + ywake_s(j)*(ycoefm_w(j,k) - ycoefm_x(j,k))
+        ENDDO
+      ENDDO
+!!!
+       ENDIF  ! (iflag_split .eq.0)
+!!!
        
 !****************************************************************************************
@@ -843,12 +1465,62 @@
 
 ! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q 
-       CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.0) THEN
+!!!
+!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
+        CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
             ydelp, yt, yq, dtime, &
+!!! jyg le 09/05/2011
+            CcoefH, CcoefQ, DcoefH, DcoefQ, &
+            Kcoef_hq, gama_q, gama_h, &
+!!!
             AcoefH, AcoefQ, BcoefH, BcoefQ)
+       ELSE  !(iflag_split .eq.0)
+        CALL climb_hq_down(knon, ycoefh_x, ypaprs, ypplay, &
+            ydelp, yt_x, yq_x, dtime, &
+!!! nrlmd le 02/05/2011
+            CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
+            Kcoef_hq_x, gama_q_x, gama_h_x, &
+!!!
+            AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x)
+!
+        CALL climb_hq_down(knon, ycoefh_w, ypaprs, ypplay, &
+            ydelp, yt_w, yq_w, dtime, &
+!!! nrlmd le 02/05/2011
+            CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
+            Kcoef_hq_w, gama_q_w, gama_h_w, &
+!!!
+            AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w)
+!!!
+       ENDIF  ! (iflag_split .eq.0)
+!!!
 
 ! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
-       CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.0) THEN
+!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
+        CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &
+!!! jyg le 09/05/2011
+            CcoefU, CcoefV, DcoefU, DcoefV, &
+            Kcoef_m, alf_1, alf_2, &
+!!!
             AcoefU, AcoefV, BcoefU, BcoefV)
-      
+       ELSE  ! (iflag_split .eq.0)
+        CALL climb_wind_down(knon, dtime, ycoefm_x, ypplay, ypaprs, yt_x, ydelp, yu_x, yv_x, &
+!!! nrlmd le 02/05/2011
+            CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
+            Kcoef_m_x, alf_1_x, alf_2_x, &
+!!!
+            AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x)
+!
+        CALL climb_wind_down(knon, dtime, ycoefm_w, ypplay, ypaprs, yt_w, ydelp, yu_w, yv_w, &
+!!! nrlmd le 02/05/2011
+            CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
+            Kcoef_m_w, alf_1_w, alf_2_w, &
+!!!
+            AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w)
+!!!      
+       ENDIF  ! (iflag_split .eq.0)
+!!!
 
 !****************************************************************************************
@@ -870,4 +1542,134 @@
        END IF
 
+!!! nrlmd le 13/06/2011
+!----- On finit le calcul des coefficients d'échange:on multiplie le cdrag par le module du vent et la densité dans la première couche 
+!          Kech_h_x(j) = ycdragh_x(j) * &
+!             (1.0+SQRT(yu_x(j,1)**2+yv_x(j,1)**2)) * &
+!             ypplay(j,1)/(RD*yt_x(j,1))
+!          Kech_h_w(j) = ycdragh_w(j) * &
+!             (1.0+SQRT(yu_w(j,1)**2+yv_w(j,1)**2)) * &
+!             ypplay(j,1)/(RD*yt_w(j,1))
+!          Kech_h(j) = (1.-ywake_s(j))*Kech_h_x(j)+ywake_s(j)*Kech_h_w(j)
+! 
+!          Kech_m_x(j) = ycdragm_x(j) * &
+!             (1.0+SQRT(yu_x(j,1)**2+yv_x(j,1)**2)) * &
+!             ypplay(j,1)/(RD*yt_x(j,1))
+!          Kech_m_w(j) = ycdragm_w(j) * &
+!             (1.0+SQRT(yu_w(j,1)**2+yv_w(j,1)**2)) * &
+!             ypplay(j,1)/(RD*yt_w(j,1))
+!          Kech_m(j) = (1.-ywake_s(j))*Kech_m_x(j)+ywake_s(j)*Kech_m_w(j)
+!!!
+
+!!! nrlmd le 02/05/2011  -----------------------On raccorde les 2 colonnes dans la couche 1 
+!----------------------------------------------------------------------------------------
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.1) THEN
+!!!
+!!! jyg le 09/04/2013 ; passage aux nouvelles expressions en differences
+
+        DO j=1,knon
+!
+! Calcul des coefficients d echange
+         mod_wind_x = 1.0+SQRT(yu_x(j,1)**2+yv_x(j,1)**2)
+         mod_wind_w = 1.0+SQRT(yu_w(j,1)**2+yv_w(j,1)**2)
+         rho1 = ypplay(j,1)/(RD*yt(j,1))
+         Kech_h_x(j) = ycdragh_x(j) * mod_wind_x * rho1
+         Kech_h_w(j) = ycdragh_w(j) * mod_wind_w * rho1
+         Kech_m_x(j) = ycdragm_x(j) * mod_wind_x * rho1
+         Kech_m_w(j) = ycdragm_w(j) * mod_wind_w * rho1
+!
+         dd_Kh = Kech_h_w(j) - Kech_h_x(j)
+         dd_Km = Kech_m_w(j) - Kech_m_x(j)
+         IF (prt_level >=10) THEN
+          print *,' mod_wind_x, mod_wind_w ', mod_wind_x, mod_wind_w
+          print *,' rho1 ',rho1
+          print *,' ycdragh_x(j),ycdragm_x(j) ',ycdragh_x(j),ycdragm_x(j)
+          print *,' ycdragh_w(j),ycdragm_w(j) ',ycdragh_w(j),ycdragm_w(j)
+          print *,' dd_Kh: ',dd_KH
+         ENDIF
+!
+         Kech_h(j) = Kech_h_x(j) + ywake_s(j)*dd_Kh
+         Kech_m(j) = Kech_m_x(j) + ywake_s(j)*dd_Km
+!
+! Calcul des coefficients d echange corriges des retroactions
+        Kech_H_xp(j) = Kech_h_x(j)/(1.-BcoefH_x(j)*Kech_h_x(j)*dtime)
+        Kech_H_wp(j) = Kech_h_w(j)/(1.-BcoefH_w(j)*Kech_h_w(j)*dtime)
+        Kech_Q_xp(j) = Kech_h_x(j)/(1.-BcoefQ_x(j)*Kech_h_x(j)*dtime)
+        Kech_Q_wp(j) = Kech_h_w(j)/(1.-BcoefQ_w(j)*Kech_h_w(j)*dtime)
+        Kech_U_xp(j) = Kech_m_x(j)/(1.-BcoefU_x(j)*Kech_m_x(j)*dtime)
+        Kech_U_wp(j) = Kech_m_w(j)/(1.-BcoefU_w(j)*Kech_m_w(j)*dtime)
+        Kech_V_xp(j) = Kech_m_x(j)/(1.-BcoefV_x(j)*Kech_m_x(j)*dtime)
+        Kech_V_wp(j) = Kech_m_w(j)/(1.-BcoefV_w(j)*Kech_m_w(j)*dtime)
+!
+         dd_KHp = Kech_H_wp(j) - Kech_H_xp(j)
+         dd_KQp = Kech_Q_wp(j) - Kech_Q_xp(j)
+         dd_KUp = Kech_U_wp(j) - Kech_U_xp(j)
+         dd_KVp = Kech_V_wp(j) - Kech_V_xp(j)
+!
+        Kech_Hp(j) = Kech_H_xp(j) + ywake_s(j)*dd_KHp
+        Kech_Qp(j) = Kech_Q_xp(j) + ywake_s(j)*dd_KQp
+        Kech_Up(j) = Kech_U_xp(j) + ywake_s(j)*dd_KUp
+        Kech_Vp(j) = Kech_V_xp(j) + ywake_s(j)*dd_KVp
+!
+! Calcul des differences w-x
+       dd_CM = ycdragm_w(j) - ycdragm_x(j)
+       dd_CH = ycdragh_w(j) - ycdragh_x(j)
+       dd_u = yu_w(j,1) - yu_x(j,1)
+       dd_v = yv_w(j,1) - yv_x(j,1)
+       dd_t = yt_w(j,1) - yt_x(j,1)
+       dd_q = yq_w(j,1) - yq_x(j,1)
+       dd_AH = AcoefH_w(j) - AcoefH_x(j)
+       dd_AQ = AcoefQ_w(j) - AcoefQ_x(j)
+       dd_AU = AcoefU_w(j) - AcoefU_x(j)
+       dd_AV = AcoefV_w(j) - AcoefV_x(j)
+       dd_BH = BcoefH_w(j) - BcoefH_x(j)
+       dd_BQ = BcoefQ_w(j) - BcoefQ_x(j)
+       dd_BU = BcoefU_w(j) - BcoefU_x(j)
+       dd_BV = BcoefV_w(j) - BcoefV_x(j)
+!
+       IF (prt_level >=10) THEN
+          print *,'Variables pour la fusion : Kech_H_xp(j)' ,Kech_H_xp(j)
+          print *,'Variables pour la fusion : Kech_H_wp(j)' ,Kech_H_wp(j)
+          print *,'Variables pour la fusion : Kech_Hp(j)' ,Kech_Hp(j)
+          print *,'Variables pour la fusion : Kech_h(j)' ,Kech_h(j)
+       ENDIF
+!
+! Calcul des coef A, B équivalents dans la couche 1
+!
+       AcoefH(j) = AcoefH_x(j) + ywake_s(j)*(Kech_H_wp(j)/Kech_Hp(j))*dd_AH
+       AcoefQ(j) = AcoefQ_x(j) + ywake_s(j)*(Kech_Q_wp(j)/Kech_Qp(j))*dd_AQ
+       AcoefU(j) = AcoefU_x(j) + ywake_s(j)*(Kech_U_wp(j)/Kech_Up(j))*dd_AU
+       AcoefV(j) = AcoefV_x(j) + ywake_s(j)*(Kech_V_wp(j)/Kech_Vp(j))*dd_AV
+!
+       BcoefH(j) = BcoefH_x(j) + ywake_s(j)*BcoefH_x(j)*(dd_Kh/Kech_h(j))*(1.+Kech_H_wp(j)/Kech_Hp(j)) &
+                               + ywake_s(j)*(Kech_H_wp(j)/Kech_Hp(j))*(Kech_h_w(j)/Kech_h(j))*dd_BH
+
+       BcoefQ(j) = BcoefQ_x(j) + ywake_s(j)*BcoefQ_x(j)*(dd_Kh/Kech_h(j))*(1.+Kech_Q_wp(j)/Kech_Qp(j)) &
+                               + ywake_s(j)*(Kech_Q_wp(j)/Kech_Qp(j))*(Kech_h_w(j)/Kech_h(j))*dd_BQ
+
+       BcoefU(j) = BcoefU_x(j) + ywake_s(j)*BcoefU_x(j)*(dd_Km/Kech_h(j))*(1.+Kech_U_wp(j)/Kech_Up(j)) &
+                               + ywake_s(j)*(Kech_U_wp(j)/Kech_Up(j))*(Kech_m_w(j)/Kech_m(j))*dd_BU
+
+       BcoefV(j) = BcoefV_x(j) + ywake_s(j)*BcoefV_x(j)*(dd_Km/Kech_h(j))*(1.+Kech_V_wp(j)/Kech_Vp(j)) &
+                               + ywake_s(j)*(Kech_V_wp(j)/Kech_Vp(j))*(Kech_m_w(j)/Kech_m(j))*dd_BV
+
+!
+! Calcul des cdrag équivalents dans la couche 
+!
+       ycdragm(j) = ycdragm_x(j) + ywake_s(j)*dd_CM
+       ycdragh(j) = ycdragh_x(j) + ywake_s(j)*dd_CH
+!
+! Calcul de T, q, u et v équivalents dans la couche 1
+       yt(j,1) = yt_x(j,1) + ywake_s(j)*(Kech_h_w(j)/Kech_h(j))*dd_t
+       yq(j,1) = yq_x(j,1) + ywake_s(j)*(Kech_h_w(j)/Kech_h(j))*dd_q
+       yu(j,1) = yu_x(j,1) + ywake_s(j)*(Kech_m_w(j)/Kech_m(j))*dd_u
+       yv(j,1) = yv_x(j,1) + ywake_s(j)*(Kech_m_w(j)/Kech_m(j))*dd_v
+
+
+        ENDDO
+!!!
+       ENDIF  ! (iflag_split .eq.1)
+!!!
+
 !****************************************************************************************
 !
@@ -893,5 +1695,5 @@
 !****************************************************************************************
 !
-! 10) Switch selon current surface
+! 10) Switch according to current surface
 !     It is necessary to start with the continental surfaces because the ocean
 !     needs their run-off.
@@ -968,5 +1770,7 @@
           !     y_flux_u1, y_flux_v1)
 
-          alb3_lic(:)=0.
+!jyg<
+!!          alb3_lic(:)=0.
+!>jyg
           DO j = 1, knon
              i = ni(j)
@@ -992,4 +1796,19 @@
                ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
                y_flux_u1, y_flux_v1)
+      IF (prt_level >=10) THEN
+          print *,'arg de surf_ocean: ycdragh ',ycdragh
+          print *,'arg de surf_ocean: ycdragm ',ycdragm
+          print *,'arg de surf_ocean: yt ', yt
+          print *,'arg de surf_ocean: yq ', yq
+          print *,'arg de surf_ocean: yts ', yts
+          print *,'arg de surf_ocean: AcoefH ',AcoefH
+          print *,'arg de surf_ocean: AcoefQ ',AcoefQ
+          print *,'arg de surf_ocean: BcoefH ',BcoefH
+          print *,'arg de surf_ocean: BcoefQ ',BcoefQ
+          print *,'arg de surf_ocean: yevap ',yevap
+          print *,'arg de surf_ocean: yfluxsens ',yfluxsens
+          print *,'arg de surf_ocean: yfluxlat ',yfluxlat
+          print *,'arg de surf_ocean: ytsurf_new ',ytsurf_new
+       ENDIF
           
        CASE(is_sic)
@@ -1036,31 +1855,178 @@
 !
 !****************************************************************************************
-! H and Q
-       IF (ok_flux_surf) THEN
-          PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
+
+!!!
+!!! jyg le 10/04/2013
+!!!
+        IF (ok_flux_surf) THEN
+          IF (prt_level >=10) THEN
+           PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
+          ENDIF
           y_flux_t1(:) =  fsens
           y_flux_q1(:) =  flat/RLVTT
           yfluxlat(:) =  flat
-
-          Kech_h(:) = ycdragh(:) * (1.0+SQRT(yu(:,1)**2+yv(:,1)**2)) * &
-               ypplay(:,1)/(RD*yt(:,1))
-          ytoto(:)=(1./RCPD)*(AcoefH(:)+BcoefH(:)*y_flux_t1(:)*dtime)
-          ytsurf_new(:)=ytoto(:)-y_flux_t1(:)/(Kech_h(:)*RCPD)
+!
+          IF (iflag_split .eq.0) THEN
+             Kech_h(:) = ycdragh(:) * (1.0+SQRT(yu(:,1)**2+yv(:,1)**2)) * &
+                  ypplay(:,1)/(RD*yt(:,1))
+          ENDIF ! (iflag_split .eq.0)
+
+          DO j = 1, knon
+            yt1_new=(1./RCPD)*(AcoefH(j)+BcoefH(j)*yfluxsens(j)*dtime)
+            ytsurf_new(j)=yt1_new-yfluxsens(j)/(Kech_h(j)*RCPD)
+          ENDDO
+
           y_d_ts(:) = ytsurf_new(:) - yts(:) 
 
-       ELSE
+        ELSE ! (ok_flux_surf)
           y_flux_t1(:) =  yfluxsens(:)
           y_flux_q1(:) = -yevap(:)
+        ENDIF
+
+       IF (prt_level >=10) THEN
+        DO j=1,knon
+         print*,'y_flux_t1,yfluxlat,wakes' &
+ &             ,  y_flux_t1(j), yfluxlat(j), ywake_s(j)
+         print*,'beta,ytsurf_new', ybeta(j), ytsurf_new(j)
+         print*,'effusivity,facteur,cstar', effusivity, facteur,wake_cstar(j)
+        ENDDO
        ENDIF
 
-       CALL climb_hq_up(knon, dtime, yt, yq, &
+!!! jyg le 07/02/2012 puis le 10/04/2013
+       IF (iflag_split .eq.1) THEN
+!!!
+        DO j=1,knon
+         y_delta_flux_t1(j) = ( Kech_H_wp(j)*Kech_H_xp(j)*(AcoefH_w(j)-AcoefH_x(j)) + &
+                                y_flux_t1(j)*(Kech_H_wp(j)-Kech_H_xp(j)) ) / Kech_Hp(j)
+         y_delta_flux_q1(j) = ( Kech_Q_wp(j)*Kech_Q_xp(j)*(AcoefQ_w(j)-AcoefQ_x(j)) + &
+                                y_flux_q1(j)*(Kech_Q_wp(j)-Kech_Q_xp(j)) ) / Kech_Qp(j)
+         y_delta_flux_u1(j) = ( Kech_U_wp(j)*Kech_U_xp(j)*(AcoefU_w(j)-AcoefU_x(j)) + &
+                                y_flux_u1(j)*(Kech_U_wp(j)-Kech_U_xp(j)) ) / Kech_Up(j)
+         y_delta_flux_v1(j) = ( Kech_V_wp(j)*Kech_V_xp(j)*(AcoefV_w(j)-AcoefV_x(j)) + &
+                                y_flux_v1(j)*(Kech_V_wp(j)-Kech_V_xp(j)) ) / Kech_Vp(j)
+!
+         y_flux_t1_x(j)=y_flux_t1(j) - ywake_s(j)*y_delta_flux_t1(j)
+         y_flux_t1_w(j)=y_flux_t1(j) + (1.-ywake_s(j))*y_delta_flux_t1(j)
+         y_flux_q1_x(j)=y_flux_q1(j) - ywake_s(j)*y_delta_flux_q1(j)
+         y_flux_q1_w(j)=y_flux_q1(j) + (1.-ywake_s(j))*y_delta_flux_q1(j)
+         y_flux_u1_x(j)=y_flux_u1(j) - ywake_s(j)*y_delta_flux_u1(j)
+         y_flux_u1_w(j)=y_flux_u1(j) + (1.-ywake_s(j))*y_delta_flux_u1(j)
+         y_flux_v1_x(j)=y_flux_v1(j) - ywake_s(j)*y_delta_flux_v1(j)
+         y_flux_v1_w(j)=y_flux_v1(j) + (1.-ywake_s(j))*y_delta_flux_v1(j)
+!
+         yfluxlat_x(j)=y_flux_q1_x(j)*RLVTT
+         yfluxlat_w(j)=y_flux_q1_w(j)*RLVTT
+
+        ENDDO
+!
+ 
+!!jyg!!   A reprendre apres reflexion   ===============================================
+!!jyg!!
+!!jyg!!        DO j=1,knon
+!!jyg!!!!! nrlmd le 13/06/2011
+!!jyg!!
+!!jyg!!!----Diffusion dans le sol dans le cas continental seulement
+!!jyg!!       IF (nsrf.eq.is_ter) THEN
+!!jyg!!!----Calcul du coefficient delta_coeff
+!!jyg!!          tau_eq(j)=(ywake_s(j)/2.)*(1./max(wake_cstar(j),0.01))*sqrt(0.4/(3.14*max(wake_dens(j),8e-12)))
+!!jyg!!
+!!jyg!!!          delta_coef(j)=dtime/(effusivity*sqrt(tau_eq(j)))
+!!jyg!!          delta_coef(j)=facteur*sqrt(tau_eq(j))/effusivity
+!!jyg!!!          delta_coef(j)=0.
+!!jyg!!       ELSE 
+!!jyg!!         delta_coef(j)=0.
+!!jyg!!       ENDIF
+!!jyg!!
+!!jyg!!!----Calcul de delta_tsurf
+!!jyg!!         y_delta_tsurf(j)=delta_coef(j)*y_delta_flux_t1(j)
+!!jyg!!
+!!jyg!!!----Si il n'y a pas des poches...
+!!jyg!!         IF (wake_cstar(j).le.0.01) THEN
+!!jyg!!           y_delta_tsurf(j)=0.
+!!jyg!!           y_delta_flux_t1(j)=0.
+!!jyg!!         ENDIF
+!!jyg!!
+!!jyg!!!-----Calcul de ybeta (evap_réelle/evap_potentielle)
+!!jyg!!!!!!! jyg le 23/02/2012
+!!jyg!!!!!!!
+!!jyg!!!!        ybeta(j)=y_flux_q1(j)   /    &
+!!jyg!!!! &        (Kech_h(j)*(yq(j,1)-yqsatsurf(j)))
+!!jyg!!!!!!        ybeta(j)=-1.*yevap(j)   /    &
+!!jyg!!!!!! &        (ywake_s(j)*Kech_h_w(j)*(yq_w(j,1)-yqsatsurf_w(j))+(1.-ywake_s(j))*Kech_h_x(j)*(yq_x(j,1)-yqsatsurf_x(j)))
+!!jyg!!!!!!! fin jyg
+!!jyg!!!!!
+!!jyg!!
+!!jyg!!       ENDDO
+!!jyg!!
+!!jyg!!!!! fin nrlmd le 13/06/2011
+!!jyg!!
+       IF (prt_level >=10) THEN
+        DO j = 1, knon
+         print*,'Chx,Chw,Ch', ycdragh_x(j), ycdragh_w(j), ycdragh(j)
+         print*,'Khx,Khw,Kh', Kech_h_x(j), Kech_h_w(j), Kech_h(j)
+!         print*,'tsurf_x,tsurf_w,tsurf,t1', ytsurf_th_x(j), ytsurf_th_w(j), ytsurf_th(j), yt(j,1)
+         print*,'tsurf_x,t1x,tsurf_w,t1w,tsurf,t1,t1_ancien', &
+ &               ytsurf_th_x(j), yt_x(j,1), ytsurf_th_w(j), yt_w(j,1), ytsurf_th(j), yt(j,1),t(j,1)
+         print*,'qsatsurf,qsatsurf_x,qsatsurf_w', yqsatsurf(j), yqsatsurf_x(j), yqsatsurf_w(j)
+         print*,'delta_coef,delta_flux,delta_tsurf,tau', delta_coef(j), y_delta_flux_t1(j), y_delta_tsurf(j), tau_eq(j) 
+        ENDDO
+
+        DO j=1,knon
+         print*,'fluxT_x, fluxT_w, y_flux_t1, fluxQ_x, fluxQ_w, yfluxlat, wakes' &
+ &             , y_flux_t1_x(j), y_flux_t1_w(j), y_flux_t1(j), y_flux_q1_x(j)*RLVTT, y_flux_q1_w(j)*RLVTT, yfluxlat(j), ywake_s(j)
+         print*,'beta,ytsurf_new,yqsatsurf', ybeta(j), ytsurf_new(j), yqsatsurf(j)
+         print*,'effusivity,facteur,cstar', effusivity, facteur,wake_cstar(j)
+        ENDDO
+       ENDIF
+
+!!! jyg le 07/02/2012
+       ENDIF  ! (iflag_split .eq.1)
+!!!
+
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.0) THEN
+!!!
+!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
+        CALL climb_hq_up(knon, dtime, yt, yq, &
             y_flux_q1, y_flux_t1, ypaprs, ypplay, &
+!!! jyg le 07/02/2012
+            AcoefH, AcoefQ, BcoefH, BcoefQ, &
+            CcoefH, CcoefQ, DcoefH, DcoefQ, &
+            Kcoef_hq, gama_q, gama_h, &
+!!!
             y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))    
-       
-
-       CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
+       ELSE  !(iflag_split .eq.0)
+        CALL climb_hq_up(knon, dtime, yt_x, yq_x, &
+            y_flux_q1_x, y_flux_t1_x, ypaprs, ypplay, &
+!!! nrlmd le 02/05/2011
+            AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x, &
+            CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
+            Kcoef_hq_x, gama_q_x, gama_h_x, &
+!!!
+            y_flux_q_x(:,:), y_flux_t_x(:,:), y_d_q_x(:,:), y_d_t_x(:,:))    
+!
+       CALL climb_hq_up(knon, dtime, yt_w, yq_w, &
+            y_flux_q1_w, y_flux_t1_w, ypaprs, ypplay, &
+!!! nrlmd le 02/05/2011
+            AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w, &
+            CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
+            Kcoef_hq_w, gama_q_w, gama_h_w, &
+!!!
+            y_flux_q_w(:,:), y_flux_t_w(:,:), y_d_q_w(:,:), y_d_t_w(:,:))    
+!!!
+       ENDIF  ! (iflag_split .eq.0)
+!!!
+
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.0) THEN
+!!!
+!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
+        CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
+!!! jyg le 07/02/2012
+            AcoefU, AcoefV, BcoefU, BcoefV, &
+            CcoefU, CcoefV, DcoefU, DcoefV, &
+            Kcoef_m, &
+!!!
             y_flux_u, y_flux_v, y_d_u, y_d_v)
-
-
      y_d_t_diss(:,:)=0.
      IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
@@ -1071,8 +2037,51 @@
 !     print*,'yamada_c OK'
 
-       DO j = 1, knon
+       ELSE  !(iflag_split .eq.0)
+        CALL climb_wind_up(knon, dtime, yu_x, yv_x, y_flux_u1_x, y_flux_v1_x, &
+!!! nrlmd le 02/05/2011
+            AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x, &
+            CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
+            Kcoef_m_x, &
+!!!
+            y_flux_u_x, y_flux_v_x, y_d_u_x, y_d_v_x)
+!
+     y_d_t_diss_x(:,:)=0.
+     IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
+        CALL yamada_c(knon,dtime,ypaprs,ypplay &
+    &   ,yu_x,yv_x,yt_x,y_d_u_x,y_d_v_x,y_d_t_x,ycdragm_x,ytke_x,ycoefm_x,ycoefh_x &
+        ,ycoefq_x,y_d_t_diss_x,yustar_x &
+    &   ,iflag_pbl,nsrf)
+     ENDIF
+!     print*,'yamada_c OK'
+
+        CALL climb_wind_up(knon, dtime, yu_w, yv_w, y_flux_u1_w, y_flux_v1_w, &
+!!! nrlmd le 02/05/2011
+            AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w, &
+            CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
+            Kcoef_m_w, &
+!!!
+            y_flux_u_w, y_flux_v_w, y_d_u_w, y_d_v_w)
+!!!
+     y_d_t_diss_w(:,:)=0.
+     IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
+        CALL yamada_c(knon,dtime,ypaprs,ypplay &
+    &   ,yu_w,yv_w,yt_w,y_d_u_w,y_d_v_w,y_d_t_w,ycdragm_w,ytke_w,ycoefm_w,ycoefh_w &
+        ,ycoefq_w,y_d_t_diss_w,yustar_w &
+    &   ,iflag_pbl,nsrf)
+     ENDIF
+!     print*,'yamada_c OK'
+!
+        IF (prt_level >=10) THEN
+         print *, 'After climbing up, lfuxlat_x, fluxlat_w ', &
+               yfluxlat_x, yfluxlat_w
+        ENDIF
+!
+       ENDIF  ! (iflag_split .eq.0)
+!!!
+
+        DO j = 1, knon
           y_dflux_t(j) = y_dflux_t(j) * ypct(j)
           y_dflux_q(j) = y_dflux_q(j) * ypct(j)
-       ENDDO
+        ENDDO
 
 !****************************************************************************************
@@ -1084,6 +2093,10 @@
 !****************************************************************************************
 
-       DO k = 1, klev
-          DO j = 1, knon
+
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.0) THEN
+!!!
+        DO k = 1, klev
+           DO j = 1, knon
              i = ni(j)
              y_d_t_diss(j,k)  = y_d_t_diss(j,k) * ypct(j)
@@ -1099,19 +2112,83 @@
 
 
+           ENDDO
+        ENDDO
+
+
+       ELSE  !(iflag_split .eq.0)
+
+! Tendances hors poches
+        DO k = 1, klev
+          DO j = 1, knon
+            i = ni(j)
+            y_d_t_diss_x(j,k)  = y_d_t_diss_x(j,k) * ypct(j)
+            y_d_t_x(j,k)  = y_d_t_x(j,k) * ypct(j)
+            y_d_q_x(j,k)  = y_d_q_x(j,k) * ypct(j)
+            y_d_u_x(j,k)  = y_d_u_x(j,k) * ypct(j)
+            y_d_v_x(j,k)  = y_d_v_x(j,k) * ypct(j)
+
+            flux_t_x(i,k,nsrf) = y_flux_t_x(j,k)
+            flux_q_x(i,k,nsrf) = y_flux_q_x(j,k)
+            flux_u_x(i,k,nsrf) = y_flux_u_x(j,k)
+            flux_v_x(i,k,nsrf) = y_flux_v_x(j,k)
           ENDDO
-       ENDDO
+        ENDDO
+
+! Tendances dans les poches
+        DO k = 1, klev
+          DO j = 1, knon
+            i = ni(j)
+            y_d_t_diss_w(j,k)  = y_d_t_diss_w(j,k) * ypct(j)
+            y_d_t_w(j,k)  = y_d_t_w(j,k) * ypct(j)
+            y_d_q_w(j,k)  = y_d_q_w(j,k) * ypct(j)
+            y_d_u_w(j,k)  = y_d_u_w(j,k) * ypct(j)
+            y_d_v_w(j,k)  = y_d_v_w(j,k) * ypct(j)
+
+            flux_t_w(i,k,nsrf) = y_flux_t_w(j,k)
+            flux_q_w(i,k,nsrf) = y_flux_q_w(j,k)
+            flux_u_w(i,k,nsrf) = y_flux_u_w(j,k)
+            flux_v_w(i,k,nsrf) = y_flux_v_w(j,k)
+          ENDDO
+        ENDDO
+
+! Flux, tendances et Tke moyenne dans la maille
+        DO k = 1, klev
+          DO j = 1, knon
+            i = ni(j)
+            flux_t(i,k,nsrf) = flux_t_x(i,k,nsrf)+ywake_s(j)*(flux_t_w(i,k,nsrf)-flux_t_x(i,k,nsrf))
+            flux_q(i,k,nsrf) = flux_q_x(i,k,nsrf)+ywake_s(j)*(flux_q_w(i,k,nsrf)-flux_q_x(i,k,nsrf))
+            flux_u(i,k,nsrf) = flux_u_x(i,k,nsrf)+ywake_s(j)*(flux_u_w(i,k,nsrf)-flux_u_x(i,k,nsrf))
+            flux_v(i,k,nsrf) = flux_v_x(i,k,nsrf)+ywake_s(j)*(flux_v_w(i,k,nsrf)-flux_v_x(i,k,nsrf))
+          ENDDO
+        ENDDO
+        DO j=1,knon
+          yfluxlat(j)=yfluxlat_x(j)+ywake_s(j)*(yfluxlat_w(j)-yfluxlat_x(j))
+        ENDDO
+        IF (prt_level >=10) THEN
+          print *,' nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf) ', &
+                    nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf)
+        ENDIF
+
+        DO k = 1, klev
+          DO j = 1, knon
+            y_d_t_diss(j,k) = y_d_t_diss_x(j,k)+ywake_s(j)*(y_d_t_diss_w(j,k) -y_d_t_diss_x(j,k))
+            y_d_t(j,k) = y_d_t_x(j,k)+ywake_s(j)*(y_d_t_w(j,k) -y_d_t_x(j,k))
+            y_d_q(j,k) = y_d_q_x(j,k)+ywake_s(j)*(y_d_q_w(j,k) -y_d_q_x(j,k))
+            y_d_u(j,k) = y_d_u_x(j,k)+ywake_s(j)*(y_d_u_w(j,k) -y_d_u_x(j,k))
+            y_d_v(j,k) = y_d_v_x(j,k)+ywake_s(j)*(y_d_v_w(j,k) -y_d_v_x(j,k))
+          ENDDO
+        ENDDO
+
+       ENDIF  ! (iflag_split .eq.0)
+!!!
 
 !      print*,'Dans pbl OK1'
 
-       evap(:,nsrf) = - flux_q(:,1,nsrf)
-       
-       alb1(:, nsrf) = 0.
-       alb2(:, nsrf) = 0.
-       snow(:, nsrf) = 0.
-       qsurf(:, nsrf) = 0.
-       rugos(:, nsrf) = 0.
-       fluxlat(:,nsrf) = 0.
+!jyg<
+!!       evap(:,nsrf) = - flux_q(:,1,nsrf)
+!>jyg
        DO j = 1, knon
           i = ni(j)
+          evap(i,nsrf) = - flux_q(i,1,nsrf)                  !jyg
           d_ts(i,nsrf) = y_d_ts(j)
           alb1(i,nsrf) = yalb1_new(j)  
@@ -1130,11 +2207,82 @@
 !      print*,'Dans pbl OK2'
 
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.1) THEN
+!!!
+!!! nrlmd le 02/05/2011
+        DO j = 1, knon
+          i = ni(j)
+          fluxlat_x(i,nsrf) = yfluxlat_x(j)
+          fluxlat_w(i,nsrf) = yfluxlat_w(j)
+!!!
+!!! nrlmd le 13/06/2011
+          delta_tsurf(i,nsrf)=y_delta_tsurf(j)*ypct(j)
+          cdragh_x(i) = cdragh_x(i) + ycdragh_x(j)*ypct(j)
+          cdragh_w(i) = cdragh_w(i) + ycdragh_w(j)*ypct(j)
+          cdragm_x(i) = cdragm_x(i) + ycdragm_x(j)*ypct(j)
+          cdragm_w(i) = cdragm_w(i) + ycdragm_w(j)*ypct(j)
+          kh(i) = kh(i) + Kech_h(j)*ypct(j)
+          kh_x(i) = kh_x(i) + Kech_h_x(j)*ypct(j)
+          kh_w(i) = kh_w(i) + Kech_h_w(j)*ypct(j)
+!!!
+        END DO
+!!!      
+       ENDIF  ! (iflag_split .eq.1)
+!!!
+!!! nrlmd le 02/05/2011
+!!jyg le 20/02/2011
+!!        tke_x(:,:,nsrf)=0.
+!!        tke_w(:,:,nsrf)=0.
+!!jyg le 20/02/2011
+!!        DO k = 1, klev+1
+!!          DO j = 1, knon
+!!            i = ni(j)
+!!            wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
+!!            tke(i,k,nsrf)   = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
+!!          ENDDO
+!!        ENDDO
+!!jyg le 20/02/2011
+!!        DO k = 1, klev+1
+!!          DO j = 1, knon
+!!            i = ni(j)
+!!            tke(i,k,nsrf)=(1.-ywake_s(j))*tke_x(i,k,nsrf)+ywake_s(j)*tke_w(i,k,nsrf)
+!!          ENDDO
+!!        ENDDO
+!!!
+       IF (iflag_split .eq.0) THEN
+        DO k = 2, klev
+           DO j = 1, knon
+              i = ni(j)
+!jyg<
+!!              tke(i,k,nsrf)    = ytke(j,k)
+!!              tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
+              tke_x(i,k,nsrf)    = ytke(j,k)
+              tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)
+!>jyg
+           END DO
+        END DO
+
+       ELSE
+        DO k = 2, klev
+          DO j = 1, knon
+            i = ni(j)
+            wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
+!jyg<
+!!            tke(i,k,nsrf)   = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
+!!            tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j)
+            tke_x(i,k,nsrf)   = ytke_x(j,k)
+            tke_x(i,k,is_ave)   = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)
+            wake_dltke(i,k,is_ave)   = wake_dltke(i,k,is_ave) + wake_dltke(i,k,nsrf)*ypct(j)
+
+!>jyg
+          ENDDO
+        ENDDO
+       ENDIF  ! (iflag_split .eq.0)
+!!!
        DO k = 2, klev
           DO j = 1, knon
              i = ni(j)
-             tke(i,k,nsrf)    = ytke(j,k)
              zcoefh(i,k,nsrf) = ycoefh(j,k)
              zcoefm(i,k,nsrf) = ycoefm(j,k)
-             tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
              zcoefh(i,k,is_ave) = zcoefh(i,k,is_ave) + ycoefh(j,k)*ypct(j)
              zcoefm(i,k,is_ave) = zcoefm(i,k,is_ave) + ycoefm(j,k)*ypct(j)
@@ -1151,5 +2299,7 @@
        END IF
        
-       ftsoil(:,:,nsrf) = 0.
+!jyg<
+!!       ftsoil(:,:,nsrf) = 0.
+!>jyg
        DO k = 1, nsoilmx
           DO j = 1, knon
@@ -1159,4 +2309,30 @@
        END DO
        
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.1) THEN
+!!!
+!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
+        DO k = 1, klev
+          DO j = 1, knon
+           i = ni(j)
+           d_t_diss_x(i,k) = d_t_diss_x(i,k) + y_d_t_diss_x(j,k)
+           d_t_x(i,k) = d_t_x(i,k) + y_d_t_x(j,k)
+           d_q_x(i,k) = d_q_x(i,k) + y_d_q_x(j,k)
+           d_u_x(i,k) = d_u_x(i,k) + y_d_u_x(j,k)
+           d_v_x(i,k) = d_v_x(i,k) + y_d_v_x(j,k)
+!
+           d_t_diss_w(i,k) = d_t_diss_w(i,k) + y_d_t_diss_w(j,k)
+           d_t_w(i,k) = d_t_w(i,k) + y_d_t_w(j,k)
+           d_q_w(i,k) = d_q_w(i,k) + y_d_q_w(j,k)
+           d_u_w(i,k) = d_u_w(i,k) + y_d_u_w(j,k)
+           d_v_w(i,k) = d_v_w(i,k) + y_d_v_w(j,k)
+!
+!!           d_wake_dlt(i,k) = d_wake_dlt(i,k) + y_d_t_w(i,k)-y_d_t_x(i,k)
+!!           d_wake_dlq(i,k) = d_wake_dlq(i,k) + y_d_q_w(i,k)-y_d_q_x(i,k)
+          END DO
+        END DO
+!!!
+       ENDIF  ! (iflag_split .eq.1)
+!!!
        
        DO k = 1, klev
@@ -1173,26 +2349,16 @@
 !      print*,'Dans pbl OK4'
 
-!****************************************************************************************
-! 14) Calculate the temperature et relative humidity at 2m and the wind at 10m 
+       IF (prt_level >=10) THEN
+         print *, 'pbl_surface tendencies for w: d_t_w, d_t_x, d_t ', &
+          d_t_w(:,1), d_t_x(:,1), d_t(:,1)
+       ENDIF
+
+!****************************************************************************************
+! 14) Calculate the temperature and relative humidity at 2m and the wind at 10m 
 !     Call HBTM
 !
 !****************************************************************************************
-       t2m(:,nsrf)    = 0.
-       q2m(:,nsrf)    = 0.
-       ustar(:,nsrf)   = 0.
-       wstar(:,nsrf)   = 0.
-       u10m(:,nsrf)   = 0.
-       v10m(:,nsrf)   = 0.
-       pblh(:,nsrf)   = 0.        ! Hauteur de couche limite
-       plcl(:,nsrf)   = 0.        ! Niveau de condensation de la CLA
-       capCL(:,nsrf)  = 0.        ! CAPE de couche limite
-       oliqCL(:,nsrf) = 0.        ! eau_liqu integree de couche limite
-       cteiCL(:,nsrf) = 0.        ! cloud top instab. crit. couche limite
-       pblt(:,nsrf)   = 0.        ! T a la Hauteur de couche limite
-       therm(:,nsrf)  = 0.
-       trmb1(:,nsrf)  = 0.        ! deep_cape
-       trmb2(:,nsrf)  = 0.        ! inhibition 
-       trmb3(:,nsrf)  = 0.        ! Point Omega
-
+!!!
+!
 #undef T2m     
 #define T2m     
@@ -1203,6 +2369,7 @@
 !      print*,'tair1,yt(:,1),y_d_t(:,1)'
 !      print*, tair1,yt(:,1),y_d_t(:,1)
-       DO j=1, knon
-          i = ni(j)
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.0) THEN
+        DO j=1, knon
           uzon(j) = yu(j,1) + y_d_u(j,1)
           vmer(j) = yv(j,1) + y_d_v(j,1)
@@ -1212,4 +2379,33 @@
                * (ypaprs(j,1)-ypplay(j,1))
           tairsol(j) = yts(j) + y_d_ts(j)
+          qairsol(j) = yqsurf(j)
+        END DO
+       ELSE  ! (iflag_split .eq.0)
+        DO j=1, knon
+          uzon_x(j) = yu_x(j,1) + y_d_u_x(j,1)
+          vmer_x(j) = yv_x(j,1) + y_d_v_x(j,1)
+          tair1_x(j) = yt_x(j,1) + y_d_t_x(j,1) + y_d_t_diss_x(j,1)
+          qair1_x(j) = yq_x(j,1) + y_d_q_x(j,1)
+          zgeo1_x(j) = RD * tair1_x(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
+               * (ypaprs(j,1)-ypplay(j,1))
+          tairsol(j) = yts(j) + y_d_ts(j)
+          tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)
+          qairsol(j) = yqsurf(j)
+        END DO
+        DO j=1, knon
+          uzon_w(j) = yu_w(j,1) + y_d_u_w(j,1)
+          vmer_w(j) = yv_w(j,1) + y_d_v_w(j,1)
+          tair1_w(j) = yt_w(j,1) + y_d_t_w(j,1) + y_d_t_diss_w(j,1)
+          qair1_w(j) = yq_w(j,1) + y_d_q_w(j,1)
+          zgeo1_w(j) = RD * tair1_w(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
+               * (ypaprs(j,1)-ypplay(j,1))
+          tairsol_w(j) = tairsol(j) + (1.- ywake_s(j))*y_delta_tsurf(j)
+          qairsol(j) = yqsurf(j)
+        END DO
+!!!      
+       ENDIF  ! (iflag_split .eq.0)
+!!!
+       DO j=1, knon
+          i = ni(j)
           rugo1(j) = yrugos(j)
           IF(nsrf.EQ.is_oce) THEN
@@ -1218,5 +2414,4 @@
           psfce(j)=ypaprs(j,1)
           patm(j)=ypplay(j,1)
-          qairsol(j) = yqsurf(j)
        END DO
        
@@ -1226,21 +2421,59 @@
 
 ! Calculate the temperature et relative humidity at 2m and the wind at 10m 
-       CALL stdlevvar(klon, knon, nsrf, zxli, &
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.0) THEN
+        CALL stdlevvar(klon, knon, nsrf, zxli, &
             uzon, vmer, tair1, qair1, zgeo1, &
             tairsol, qairsol, rugo1, psfce, patm, &
             yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
-!      print*,'Dans pbl OK42B'
-
-       DO j=1, knon
+       ELSE  !(iflag_split .eq.0)
+        CALL stdlevvar(klon, knon, nsrf, zxli, &
+            uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
+            tairsol_x, qairsol, rugo1, psfce, patm, &
+            yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x)
+        CALL stdlevvar(klon, knon, nsrf, zxli, &
+            uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
+            tairsol_w, qairsol, rugo1, psfce, patm, &
+            yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w)
+!!!
+       ENDIF  ! (iflag_split .eq.0)
+!!!
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.0) THEN
+        DO j=1, knon
           i = ni(j)
           t2m(i,nsrf)=yt2m(j)
           q2m(i,nsrf)=yq2m(j)
-          
-          ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
+     ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
           ustar(i,nsrf)=yustar(j)
           u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)**2+vmer(j)**2)
           v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)**2+vmer(j)**2)
-
-       END DO
+        END DO
+       ELSE  !(iflag_split .eq.0)
+        DO j=1, knon
+          i = ni(j)
+          t2m_x(i,nsrf)=yt2m_x(j)
+          q2m_x(i,nsrf)=yq2m_x(j)
+     ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
+          ustar_x(i,nsrf)=yustar_x(j)
+          u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/SQRT(uzon_x(j)**2+vmer_x(j)**2)
+          v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/SQRT(uzon_x(j)**2+vmer_x(j)**2)
+        END DO
+        DO j=1, knon
+          i = ni(j)
+          t2m_w(i,nsrf)=yt2m_w(j)
+          q2m_w(i,nsrf)=yq2m_w(j)
+     ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
+          ustar_w(i,nsrf)=yustar_w(j)
+          u10m_w(i,nsrf)=(yu10m_w(j) * uzon_w(j))/SQRT(uzon_w(j)**2+vmer_w(j)**2)
+          v10m_w(i,nsrf)=(yu10m_w(j) * vmer_w(j))/SQRT(uzon_w(j)**2+vmer_w(j)**2)
+!
+          ustar(i,nsrf) = ustar_x(i,nsrf) + wake_s(i)*(ustar_w(i,nsrf)-ustar_x(i,nsrf))
+          u10m(i,nsrf) = u10m_x(i,nsrf) + wake_s(i)*(u10m_w(i,nsrf)-u10m_x(i,nsrf))
+          v10m(i,nsrf) = v10m_x(i,nsrf) + wake_s(i)*(v10m_w(i,nsrf)-v10m_x(i,nsrf))
+        END DO
+!!!
+       ENDIF  ! (iflag_split .eq.0)
+!!!
 
 !      print*,'Dans pbl OK43'
@@ -1248,4 +2481,6 @@
 !IM Ajoute dependance type surface
        IF (thermcep) THEN
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.0) THEN
           DO j = 1, knon
              i=ni(j)
@@ -1259,14 +2494,91 @@
              qsat2m(i) = qsat2m(i) + zx_qs1  * pctsrf(i,nsrf)
           END DO
+       ELSE  ! (iflag_split .eq.0)
+          DO j = 1, knon
+             i=ni(j)
+             zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_x(j) ))
+             zx_qs1  = r2es * FOEEW(yt2m_x(j),zdelta1)/paprs(i,1)
+             zx_qs1  = MIN(0.5,zx_qs1)
+             zcor1   = 1./(1.-RETV*zx_qs1)
+             zx_qs1  = zx_qs1*zcor1
+             
+             rh2m_x(i)   = rh2m_x(i)   + yq2m_x(j)/zx_qs1 * pctsrf(i,nsrf)
+             qsat2m_x(i) = qsat2m_x(i) + zx_qs1  * pctsrf(i,nsrf)
+          END DO
+          DO j = 1, knon
+             i=ni(j)
+             zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_w(j) ))
+             zx_qs1  = r2es * FOEEW(yt2m_w(j),zdelta1)/paprs(i,1)
+             zx_qs1  = MIN(0.5,zx_qs1)
+             zcor1   = 1./(1.-RETV*zx_qs1)
+             zx_qs1  = zx_qs1*zcor1
+             
+             rh2m_w(i)   = rh2m_w(i)   + yq2m_w(j)/zx_qs1 * pctsrf(i,nsrf)
+             qsat2m_w(i) = qsat2m_w(i) + zx_qs1  * pctsrf(i,nsrf)
+          END DO
+!!!      
+       ENDIF  ! (iflag_split .eq.0)
+!!!
        END IF
+!
+       IF (prt_level >=10) THEN
+         print *, 'T2m, q2m, RH2m ', &
+          t2m, q2m, rh2m
+       ENDIF
 
 !   print*,'OK pbl 5'
-       CALL hbtm(knon, ypaprs, ypplay, &
+!
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.0) THEN
+        CALL hbtm(knon, ypaprs, ypplay, &
             yt2m,yt10m,yq2m,yq10m,yustar,ywstar, &
             y_flux_t,y_flux_q,yu,yv,yt,yq, &
             ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
             ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
+          IF (prt_level >=10) THEN
+       print *,' Arg. de HBTM: yt2m ',yt2m
+       print *,' Arg. de HBTM: yt10m ',yt10m
+       print *,' Arg. de HBTM: yq2m ',yq2m
+       print *,' Arg. de HBTM: yq10m ',yq10m
+       print *,' Arg. de HBTM: yustar ',yustar
+       print *,' Arg. de HBTM: y_flux_t ',y_flux_t
+       print *,' Arg. de HBTM: y_flux_q ',y_flux_q
+       print *,' Arg. de HBTM: yu ',yu
+       print *,' Arg. de HBTM: yv ',yv
+       print *,' Arg. de HBTM: yt ',yt
+       print *,' Arg. de HBTM: yq ',yq
+          ENDIF
+       ELSE  ! (iflag_split .eq.0)
+        CALL HBTM(knon, ypaprs, ypplay, &
+            yt2m_x,yt10m_x,yq2m_x,yq10m_x,yustar_x,ywstar_x, &
+            y_flux_t_x,y_flux_q_x,yu_x,yv_x,yt_x,yq_x, &
+            ypblh_x,ycapCL_x,yoliqCL_x,ycteiCL_x,ypblT_x, &
+            ytherm_x,ytrmb1_x,ytrmb2_x,ytrmb3_x,ylcl_x)
+          IF (prt_level >=10) THEN
+       print *,' Arg. de HBTM: yt2m_x ',yt2m_x
+       print *,' Arg. de HBTM: yt10m_x ',yt10m_x
+       print *,' Arg. de HBTM: yq2m_x ',yq2m_x
+       print *,' Arg. de HBTM: yq10m_x ',yq10m_x
+       print *,' Arg. de HBTM: yustar_x ',yustar_x
+       print *,' Arg. de HBTM: y_flux_t_x ',y_flux_t_x
+       print *,' Arg. de HBTM: y_flux_q_x ',y_flux_q_x
+       print *,' Arg. de HBTM: yu_x ',yu_x
+       print *,' Arg. de HBTM: yv_x ',yv_x
+       print *,' Arg. de HBTM: yt_x ',yt_x
+       print *,' Arg. de HBTM: yq_x ',yq_x
+          ENDIF
+        CALL HBTM(knon, ypaprs, ypplay, &
+            yt2m_w,yt10m_w,yq2m_w,yq10m_w,yustar_w,ywstar_w, &
+            y_flux_t_w,y_flux_q_w,yu_w,yv_w,yt_w,yq_w, &
+            ypblh_w,ycapCL_w,yoliqCL_w,ycteiCL_w,ypblT_w, &
+            ytherm_w,ytrmb1_w,ytrmb2_w,ytrmb3_w,ylcl_w)
+!!!      
+       ENDIF  ! (iflag_split .eq.0)
+!!!
        
-       DO j=1, knon
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.0) THEN
+!!!
+        DO j=1, knon
           i = ni(j)
           pblh(i,nsrf)   = ypblh(j)
@@ -1281,6 +2593,53 @@
           trmb2(i,nsrf)  = ytrmb2(j)
           trmb3(i,nsrf)  = ytrmb3(j)
-       END DO
-       
+        END DO
+        IF (prt_level >=10) THEN
+          print *, 'After HBTM: pblh ', pblh
+          print *, 'After HBTM: plcl ', plcl
+          print *, 'After HBTM: cteiCL ', cteiCL
+        ENDIF
+       ELSE  !(iflag_split .eq.0)
+        DO j=1, knon
+          i = ni(j)
+          pblh_x(i,nsrf)   = ypblh_x(j)
+          wstar_x(i,nsrf)  = ywstar_x(j)
+          plcl_x(i,nsrf)   = ylcl_x(j)
+          capCL_x(i,nsrf)  = ycapCL_x(j)
+          oliqCL_x(i,nsrf) = yoliqCL_x(j)
+          cteiCL_x(i,nsrf) = ycteiCL_x(j)
+          pblT_x(i,nsrf)   = ypblT_x(j)
+          therm_x(i,nsrf)  = ytherm_x(j)
+          trmb1_x(i,nsrf)  = ytrmb1_x(j)
+          trmb2_x(i,nsrf)  = ytrmb2_x(j)
+          trmb3_x(i,nsrf)  = ytrmb3_x(j)
+        END DO
+        IF (prt_level >=10) THEN
+          print *, 'After HBTM: pblh_x ', pblh_x
+          print *, 'After HBTM: plcl_x ', plcl_x
+          print *, 'After HBTM: cteiCL_x ', cteiCL_x
+        ENDIF
+        DO j=1, knon
+          i = ni(j)
+          pblh_w(i,nsrf)   = ypblh_w(j)
+          wstar_w(i,nsrf)  = ywstar_w(j)
+          plcl_w(i,nsrf)   = ylcl_w(j)
+          capCL_w(i,nsrf)  = ycapCL_w(j)
+          oliqCL_w(i,nsrf) = yoliqCL_w(j)
+          cteiCL_w(i,nsrf) = ycteiCL_w(j)
+          pblT_w(i,nsrf)   = ypblT_w(j)
+          therm_w(i,nsrf)  = ytherm_w(j)
+          trmb1_w(i,nsrf)  = ytrmb1_w(j)
+          trmb2_w(i,nsrf)  = ytrmb2_w(j)
+          trmb3_w(i,nsrf)  = ytrmb3_w(j)
+        END DO
+        IF (prt_level >=10) THEN
+          print *, 'After HBTM: pblh_w ', pblh_w
+          print *, 'After HBTM: plcl_w ', plcl_w
+          print *, 'After HBTM: cteiCL_w ', cteiCL_w
+        ENDIF
+!!!
+       ENDIF  ! (iflag_split .eq.0)
+!!!
+
 !   print*,'OK pbl 6'
 #else 
@@ -1297,5 +2656,5 @@
 
 !****************************************************************************************
-! 16) Calculate the mean value over all sub-surfaces for som variables
+! 16) Calculate the mean value over all sub-surfaces for some variables
 !
 !****************************************************************************************
@@ -1304,4 +2663,38 @@
     zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
     zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
+    zxfluxt_x(:,:) = 0.0 ; zxfluxq_x(:,:) = 0.0
+    zxfluxu_x(:,:) = 0.0 ; zxfluxv_x(:,:) = 0.0
+    zxfluxt_w(:,:) = 0.0 ; zxfluxq_w(:,:) = 0.0
+    zxfluxu_w(:,:) = 0.0 ; zxfluxv_w(:,:) = 0.0
+
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.1) THEN
+!!!
+!!! nrlmd & jyg les 02/05/2011, 05/02/2012
+
+        DO nsrf = 1, nbsrf
+          DO k = 1, klev
+            DO i = 1, klon
+              zxfluxt_x(i,k) = zxfluxt_x(i,k) + flux_t_x(i,k,nsrf) * pctsrf(i,nsrf)
+              zxfluxq_x(i,k) = zxfluxq_x(i,k) + flux_q_x(i,k,nsrf) * pctsrf(i,nsrf)
+              zxfluxu_x(i,k) = zxfluxu_x(i,k) + flux_u_x(i,k,nsrf) * pctsrf(i,nsrf)
+              zxfluxv_x(i,k) = zxfluxv_x(i,k) + flux_v_x(i,k,nsrf) * pctsrf(i,nsrf)
+!
+              zxfluxt_w(i,k) = zxfluxt_w(i,k) + flux_t_w(i,k,nsrf) * pctsrf(i,nsrf)
+              zxfluxq_w(i,k) = zxfluxq_w(i,k) + flux_q_w(i,k,nsrf) * pctsrf(i,nsrf)
+              zxfluxu_w(i,k) = zxfluxu_w(i,k) + flux_u_w(i,k,nsrf) * pctsrf(i,nsrf)
+              zxfluxv_w(i,k) = zxfluxv_w(i,k) + flux_v_w(i,k,nsrf) * pctsrf(i,nsrf)
+            END DO
+          END DO
+        END DO
+
+    DO i = 1, klon
+      zxsens_x(i) = - zxfluxt_x(i,1)
+      zxsens_w(i) = - zxfluxt_w(i,1)
+    END DO
+!!!
+       ENDIF  ! (iflag_split .eq.1)
+!!!
+
     DO nsrf = 1, nbsrf
        DO k = 1, klev
@@ -1315,5 +2708,4 @@
     END DO
 
-!   print*,'OK pbl 8'
     DO i = 1, klon
        zxsens(i)     = - zxfluxt(i,1) ! flux de chaleur sensible au sol
@@ -1321,4 +2713,5 @@
        fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
     ENDDO
+!!!
    
 !
@@ -1329,4 +2722,8 @@
     zustar(:)=0.0 ; zu10m(:) = 0.0   ; zv10m(:) = 0.0
     s_pblh(:) = 0.0  ; s_plcl(:) = 0.0 
+!!! jyg le 07/02/2012
+     s_pblh_x(:) = 0.0  ; s_plcl_x(:) = 0.0 
+     s_pblh_w(:) = 0.0  ; s_plcl_w(:) = 0.0 
+!!!
     s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
     s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
@@ -1336,4 +2733,8 @@
     
 !   print*,'OK pbl 9'
+    
+!!! nrlmd le 02/05/2011
+    zxfluxlat_x(:) = 0.0  ;  zxfluxlat_w(:) = 0.0
+!!!
     
     DO nsrf = 1, nbsrf
@@ -1348,5 +2749,11 @@
           zxtsol(i)    = zxtsol(i)    + ts(i,nsrf)      * pctsrf(i,nsrf)
           zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
+       END DO
+    END DO
           
+!!! jyg le 07/02/2012
+       IF (iflag_split .eq.0) THEN
+        DO nsrf = 1, nbsrf
+         DO i = 1, klon          
           zt2m(i)  = zt2m(i)  + t2m(i,nsrf)  * pctsrf(i,nsrf)
           zq2m(i)  = zq2m(i)  + q2m(i,nsrf)  * pctsrf(i,nsrf)
@@ -1366,7 +2773,49 @@
           s_trmb2(i)  = s_trmb2(i)  + trmb2(i,nsrf) * pctsrf(i,nsrf)
           s_trmb3(i)  = s_trmb3(i)  + trmb3(i,nsrf) * pctsrf(i,nsrf)
-       END DO
-    END DO
-!   print*,'OK pbl 10'
+         END DO
+        END DO
+       ELSE  !(iflag_split .eq.0)
+        DO nsrf = 1, nbsrf
+         DO i = 1, klon          
+!!! nrlmd le 02/05/2011
+          zxfluxlat_x(i) = zxfluxlat_x(i) + fluxlat_x(i,nsrf) * pctsrf(i,nsrf)
+          zxfluxlat_w(i) = zxfluxlat_w(i) + fluxlat_w(i,nsrf) * pctsrf(i,nsrf)
+!!!
+!!! jyg le 08/02/2012
+!!  Pour le moment, on sort les valeurs dans (x) et (w) de pblh et de plcl ; 
+!!  pour zt2m, on fait la moyenne surfacique sur les sous-surfaces ;
+!!  pour qsat2m, on fait la moyenne surfacique sur (x) et (w) ;
+!!  pour les autres variables, on sort les valeurs de la region (x).
+          zt2m(i)  = zt2m(i)  + (t2m_x(i,nsrf)+wake_s(i)*(t2m_w(i,nsrf)-t2m_x(i,nsrf))) * pctsrf(i,nsrf)
+          zq2m(i)  = zq2m(i)  + q2m_x(i,nsrf)  * pctsrf(i,nsrf)
+          zustar(i) = zustar(i) + ustar_x(i,nsrf) * pctsrf(i,nsrf)
+          wstar(i,is_ave)=wstar(i,is_ave)+wstar_x(i,nsrf)*pctsrf(i,nsrf)
+          zu10m(i) = zu10m(i) + u10m_x(i,nsrf) * pctsrf(i,nsrf)
+          zv10m(i) = zv10m(i) + v10m_x(i,nsrf) * pctsrf(i,nsrf)
+!
+          s_pblh(i)     = s_pblh(i)     + pblh_x(i,nsrf)  * pctsrf(i,nsrf)
+          s_pblh_x(i)   = s_pblh_x(i)   + pblh_x(i,nsrf)  * pctsrf(i,nsrf)
+          s_pblh_w(i)   = s_pblh_w(i)   + pblh_w(i,nsrf)  * pctsrf(i,nsrf)
+!
+          s_plcl(i)     = s_plcl(i)     + plcl_x(i,nsrf)  * pctsrf(i,nsrf)
+          s_plcl_x(i)   = s_plcl_x(i)   + plcl_x(i,nsrf)  * pctsrf(i,nsrf)
+          s_plcl_w(i)   = s_plcl_w(i)   + plcl_w(i,nsrf)  * pctsrf(i,nsrf)
+!
+          s_capCL(i)  = s_capCL(i)  + capCL_x(i,nsrf) * pctsrf(i,nsrf)
+          s_oliqCL(i) = s_oliqCL(i) + oliqCL_x(i,nsrf)* pctsrf(i,nsrf)
+          s_cteiCL(i) = s_cteiCL(i) + cteiCL_x(i,nsrf)* pctsrf(i,nsrf)
+          s_pblT(i)   = s_pblT(i)   + pblT_x(i,nsrf)  * pctsrf(i,nsrf)
+          s_therm(i)  = s_therm(i)  + therm_x(i,nsrf) * pctsrf(i,nsrf)
+          s_trmb1(i)  = s_trmb1(i)  + trmb1_x(i,nsrf) * pctsrf(i,nsrf)
+          s_trmb2(i)  = s_trmb2(i)  + trmb2_x(i,nsrf) * pctsrf(i,nsrf)
+          s_trmb3(i)  = s_trmb3(i)  + trmb3_x(i,nsrf) * pctsrf(i,nsrf)
+         END DO
+        END DO
+        DO i = 1, klon          
+          qsat2m(i)= qsat2m_x(i)+ wake_s(i)*(qsat2m_x(i)-qsat2m_w(i))
+        END DO
+!!!
+       ENDIF  ! (iflag_split .eq.0)
+!!!
 
     IF (check) THEN
@@ -1508,5 +2957,5 @@
     REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: alb1, alb2
     REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: ustar,u10m, v10m
-    REAL, DIMENSION(klon,klev+1,nbsrf), INTENT(INOUT) :: tke
+    REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT) :: tke
 
 ! Local variables
@@ -1597,2 +3046,3 @@
 
 END MODULE pbl_surface_mod
+
