Index: LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 3102)
+++ LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 3179)
@@ -23,5 +23,8 @@
   USE climb_wind_mod,      ONLY : climb_wind_down, climb_wind_up
   USE coef_diff_turb_mod,  ONLY : coef_diff_turb
-
+  USE wx_pbl_mod,          ONLY : wx_pbl_init, wx_pbl_final, &
+!!                                  wx_pbl_fuse_no_dts, wx_pbl_split_no_dts, &
+!!                                  wx_pbl_fuse, wx_pbl_split
+                                  wx_pbl0_fuse, wx_pbl0_split
 
   IMPLICIT NONE
@@ -145,4 +148,13 @@
 !    CALL getin_p('iflag_frein',iflag_frein)
 !
+!jyg<
+!****************************************************************************************
+! Allocate variables for pbl splitting
+!
+!****************************************************************************************
+
+    CALL wx_pbl_init
+!>jyg
+
   END SUBROUTINE pbl_surface_init
 !  
@@ -201,6 +213,6 @@
 !!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
 !!        tke_x,     tke_w                              &
-       wake_dltke                                     &
-        , treedrg                                   &
+       wake_dltke,                                     &
+        treedrg                                   &
 !FC
 !!!
@@ -273,10 +285,12 @@
 ! pblh-----output-R- HCL
 ! pblT-----output-R- T au nveau HCL
-!
-    USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, co2_send
+! treedrg--output-R- tree drag (m)               
+!
+    USE carbon_cycle_mod,   ONLY : carbon_cycle_cpl, co2_send
     USE indice_sol_mod
-    USE time_phylmdz_mod, ONLY: day_ini,annee_ref,itau_phy
-    USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
-    USE print_control_mod, ONLY: prt_level,lunout
+    USE time_phylmdz_mod,   ONLY : day_ini,annee_ref,itau_phy
+    USE mod_grid_phy_lmdz,  ONLY : nbp_lon, nbp_lat
+    USE print_control_mod,  ONLY : prt_level,lunout
+    USE ioipsl_getin_p_mod, ONLY : getin_p
 
     IMPLICIT NONE
@@ -490,4 +504,6 @@
     CHARACTER(len=8), DIMENSION(nbsrf), SAVE :: cl_surf
 !$OMP THREADPRIVATE(cl_surf)
+    REAL, SAVE                               :: beta_land         ! beta for wx_dts
+!$OMP THREADPRIVATE(beta_land)
 
 ! Other local variables
@@ -507,5 +523,5 @@
     REAL, DIMENSION(klon)              :: yalb,yalb_vis
 !albedo SB <<<
-    REAL, DIMENSION(klon)              :: yu1, yv1
+    REAL, DIMENSION(klon)              :: yt1, yq1, yu1, yv1
     REAL, DIMENSION(klon)              :: ysnow, yqsurf, yagesno, yqsol
     REAL, DIMENSION(klon)              :: yrain_f, ysnow_f
@@ -716,25 +732,27 @@
 !!! jyg le 25/03/2013
 !!    Variables intermediaires pour le raccord des deux colonnes \`a 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
+!jyg<
+!!    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
+!>jyg
 
 !!!
@@ -743,5 +761,5 @@
     REAL, DIMENSION(klon)              :: y_delta_tsurf,delta_coef,tau_eq
     REAL, PARAMETER                    :: facteur=2./sqrt(3.14)
-    REAL, PARAMETER                    :: effusivity=2000.
+    REAL, PARAMETER                    :: inertia=2000.
     REAL, DIMENSION(klon)              :: ytsurf_th_x,ytsurf_th_w,yqsatsurf_x,yqsatsurf_w
     REAL, DIMENSION(klon)              :: ydtsurf_th
@@ -755,8 +773,13 @@
     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
+!jyg<
+!!    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
+!>jyg
+!jyg<
+    REAL, DIMENSION(klon)              :: ah, bh     ! coefficients of the delta_Tsurf equation
+!>jyg
 
     REAL                               :: vent
@@ -824,4 +847,8 @@
        ! Initialize ok_flux_surf (for 1D model)
        if (klon_glo>1) ok_flux_surf=.FALSE.
+
+       ! intialize beta_land
+       beta_land = 0.5
+       call getin_p('beta_land', beta_land)
        
        ! Initilize debug IO
@@ -881,4 +908,5 @@
  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.
@@ -1301,4 +1329,8 @@
           ENDDO
         ENDDO
+        IF (prt_level .ge. 10) THEN
+          print *,'pbl_surface, wake_s(1), wake_dlt(1,:) ', wake_s(1), wake_dlt(1,:)
+          print *,'pbl_surface, wake_s(1), wake_dlq(1,:) ', wake_s(1), wake_dlq(1,:)
+        ENDIF
 !!! nrlmd le 02/05/2011
         DO k = 1, klev+1
@@ -1434,4 +1466,6 @@
             yts_w, yqsurf, yz0m, yz0h, &
             ycdragm_w, ycdragh_w, zri1_w, pref_w )
+!
+        zgeo1(:) = wake_s(:)*zgeo1_w(:) + (1.-wake_s(:))*zgeo1_x(:)
 
 ! --- special Dice. JYG+MPL 25112013 puis BOMEX
@@ -1475,5 +1509,5 @@
             ycoefm, ycoefh, ytke, y_treedrg)
 !            ycoefm, ycoefh, ytke)
-!FC y_treedrg ajouté
+!FC y_treedrg ajoute
        IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
 ! In this case, coef_diff_turb is called for the Cd only
@@ -1587,4 +1621,11 @@
 !!!
             AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x)
+!!!
+       IF (prt_level >=10) THEN
+         PRINT *,'pbl_surface (climb_hq_down.x->) AcoefH_x ',AcoefH_x
+         PRINT *,'pbl_surface (climb_hq_down.x->) AcoefQ_x ',AcoefQ_x
+         PRINT *,'pbl_surface (climb_hq_down.x->) BcoefH_x ',BcoefH_x
+         PRINT *,'pbl_surface (climb_hq_down.x->) BcoefQ_x ',BcoefQ_x
+       ENDIF
 !
         CALL climb_hq_down(knon, ycoefh_w, ypaprs, ypplay, &
@@ -1595,4 +1636,11 @@
 !!!
             AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w)
+!!!
+       IF (prt_level >=10) THEN
+         PRINT *,'pbl_surface (climb_hq_down.w->) AcoefH_w ',AcoefH_w
+         PRINT *,'pbl_surface (climb_hq_down.w->) AcoefQ_w ',AcoefQ_w
+         PRINT *,'pbl_surface (climb_hq_down.w->) BcoefH_w ',BcoefH_w
+         PRINT *,'pbl_surface (climb_hq_down.w->) BcoefQ_w ',BcoefQ_w
+       ENDIF
 !!!
        ENDIF  ! (iflag_split .eq.0)
@@ -1645,133 +1693,84 @@
        END IF
 
-!!! nrlmd le 13/06/2011
-!----- On finit le calcul des coefficients d'\'echange:on multiplie le cdrag par le module du vent et la densit\'e dans la premi\`ere 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)
-!
+!!! jyg le 01/02/2017
+       IF (iflag_split .eq. 0) THEN
+         yt1(:) = yt(:,1)
+         yq1(:) = yq(:,1)
+!!       ELSE IF (iflag_split .eq. 1) THEN
+!!!
+!jyg<
+!!         CALL wx_pbl_fuse_no_dts(knon, dtime, ypplay, ywake_s, &
+!!                                 yt_x, yt_w, yq_x, yq_w, &
+!!                                 yu_x, yu_w, yv_x, yv_w, &
+!!                                 ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w, &
+!!                                 AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
+!!                                 AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
+!!                                 BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
+!!                                 BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
+!!                                 AcoefH, AcoefQ, AcoefU, AcoefV, &
+!!                                 BcoefH, BcoefQ, BcoefU, BcoefV, &
+!!                                 ycdragh, ycdragm, &
+!!                                 yt1, yq1, yu1, yv1 &
+!!                                 )
+       ELSE IF (iflag_split .ge. 1) THEN
+         CALL wx_pbl0_fuse(knon, dtime, ypplay, ywake_s, &
+                         yt_x, yt_w, yq_x, yq_w, &
+                         yu_x, yu_w, yv_x, yv_w, &
+                         ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w, &
+                         AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
+                         AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
+                         BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
+                         BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
+                         AcoefH, AcoefQ, AcoefU, AcoefV, &
+                         BcoefH, BcoefQ, BcoefU, BcoefV, &
+                         ycdragh, ycdragm, &
+                         yt1, yq1, yu1, yv1 &
+                         )
+!!       ELSE IF (iflag_split .ge.2) THEN
+!!!    Provisoire
+!!         ah(:) = 0.
+!!         bh(:) = 0.
+!!         IF (nsrf == is_oce) THEN
+!!           ybeta(:) = 1.
+!!         ELSE
+!!           ybeta(:) = beta_land
+!!         ENDIF
+!!         ycdragh(:) = ywake_s(:)*ycdragh_w(:) + (1.-ywake_s(:))*ycdragh_x(:)
+!!         CALL wx_dts(knon, nsrf, ywake_cstar, ywake_s, ywake_dens, &
+!!                     yts, ypplay(:,1), ybeta, ycdragh , ypaprs(:,1), &
+!!                     yq(:,1), yt(:,1), yu(:,1), yv(:,1), ygustiness, &
+!!                     ah, bh &
+!!                     )
+!!!
+!!         CALL wx_pbl_fuse(knon, dtime, ypplay, ywake_s, &
+!!                         yt_x, yt_w, yq_x, yq_w, &
+!!                         yu_x, yu_w, yv_x, yv_w, &
+!!                         ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w, &
+!!                         AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
+!!                         AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
+!!                         BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
+!!                         BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
+!!                         ah, bh, &
+!!                         AcoefH, AcoefQ, AcoefU, AcoefV, &
+!!                         BcoefH, BcoefQ, BcoefU, BcoefV, &
+!!                         ycdragh, ycdragm, &
+!!                         yt1, yq1, yu1, yv1 &
+!!                         )
+!>jyg
+!!!
+         ENDIF  ! (iflag_split .eq.0)
+!!!
        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)
+         PRINT *,'pbl_surface (fuse->): yt(1,:) ',yt(1,:)
+         PRINT *,'pbl_surface (fuse->): yq(1,:) ',yq(1,:)
+         PRINT *,'pbl_surface (fuse->): yu(1,:) ',yu(1,:)
+         PRINT *,'pbl_surface (fuse->): yv(1,:) ',yv(1,:)
+         PRINT *,'pbl_surface (fuse->): AcoefH(1) ',AcoefH(1)
+         PRINT *,'pbl_surface (fuse->): BcoefH(1) ',BcoefH(1)
        ENDIF
-!
-! Calcul des coef A, B \'equivalents 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 \'equivalents 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 \'equivalents 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)
-!!!
 
 !****************************************************************************************
@@ -1810,5 +1809,6 @@
                rlon, rlat, yrmu0, &
                debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
-               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
+!!jyg               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
+               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
                AcoefH, AcoefQ, BcoefH, BcoefQ, & 
                AcoefU, AcoefV, BcoefU, BcoefV, & 
@@ -1856,5 +1856,6 @@
                yrmu0, ylwdown, yalb, ypphi(:,1), &
                ysolsw, ysollw, yts, ypplay(:,1), &
-               ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
+!!jyg               ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
+               ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
                AcoefH, AcoefQ, BcoefH, BcoefQ, &
                AcoefU, AcoefV, BcoefU, BcoefV, &
@@ -1893,5 +1894,6 @@
                ywindsp, rmu0, yfder, yts, &
                itap, dtime, jour, knon, ni, &
-               ypplay(:,1), zgeo1/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
+!!jyg               ypplay(:,1), zgeo1/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
+               ypplay(:,1), zgeo1/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
                AcoefH, AcoefQ, BcoefH, BcoefQ, &
                AcoefU, AcoefV, BcoefU, BcoefV, &
@@ -1931,5 +1933,6 @@
                itap, dtime, jour, knon, ni, &
                lafin, &
-               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
+!!jyg               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
+               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
                AcoefH, AcoefQ, BcoefH, BcoefQ, &
                AcoefU, AcoefV, BcoefU, BcoefV, &
@@ -1990,10 +1993,11 @@
           yfluxlat(:) =  flat
 !
-          IF (iflag_split .eq.0) THEN
+!!  Test sur iflag_split retire le 2/02/2018, sans vraiment comprendre la raison de ce test. (jyg)
+!!          IF (iflag_split .eq.0) THEN
              do j=1,knon
              Kech_h(j) = ycdragh(j) * (1.0+SQRT(yu(j,1)**2+yv(j,1)**2)) * &
                   ypplay(j,1)/(RD*yt(j,1))
              enddo
-          ENDIF ! (iflag_split .eq.0)
+!!          ENDIF ! (iflag_split .eq.0)
 
           DO j = 1, knon
@@ -2018,34 +2022,37 @@
  &             ,  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)
+         print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j)
         ENDDO
        ENDIF
 
 !!! 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
+!!       IF (iflag_split .eq.1) THEN
+!!!!!
+!!!jyg<
+!!         CALL wx_pbl_split_no_dts(knon, ywake_s, &
+!!                                AcoefH_x, AcoefH_w, &
+!!                                AcoefQ_x, AcoefQ_w, &
+!!                                AcoefU_x, AcoefU_w, &
+!!                                AcoefV_x, AcoefV_w, &
+!!                                y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &
+!!                                y_flux_t1_x, y_flux_t1_w, &
+!!                                y_flux_q1_x, y_flux_q1_w, &
+!!                                y_flux_u1_x, y_flux_u1_w, &
+!!                                y_flux_v1_x, y_flux_v1_w, &
+!!                                yfluxlat_x, yfluxlat_w &
+!!                                )
+!!       ELSE IF (iflag_split .ge. 2) THEN
+       IF (iflag_split .GE. 1) THEN
+         CALL wx_pbl0_split(knon, dtime, ywake_s, &
+                       y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &
+                       y_flux_t1_x, y_flux_t1_w, &
+                       y_flux_q1_x, y_flux_q1_w, &
+                       y_flux_u1_x, y_flux_u1_w, &
+                       y_flux_v1_x, y_flux_v1_w, &
+                       yfluxlat_x, yfluxlat_w, &
+                       y_delta_tsurf &
+                       )
+       ENDIF  ! (iflag_split .ge. 1)
+!>jyg
 !
  
@@ -2060,6 +2067,6 @@
 !!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)=dtime/(inertia*sqrt(tau_eq(j)))
+!!jyg!!          delta_coef(j)=facteur*sqrt(tau_eq(j))/inertia
 !!jyg!!!          delta_coef(j)=0.
 !!jyg!!       ELSE 
@@ -2090,4 +2097,5 @@
 !!jyg!!!!! fin nrlmd le 13/06/2011
 !!jyg!!
+       IF (iflag_split .ge. 1) THEN
        IF (prt_level >=10) THEN
         DO j = 1, knon
@@ -2105,10 +2113,10 @@
  &             , 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)
+         print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j)
         ENDDO
        ENDIF  ! (prt_level >=10)
 
 !!! jyg le 07/02/2012
-       ENDIF  ! (iflag_split .eq.1)
+       ENDIF  ! (iflag_split .ge.1)
 !!!
 
@@ -2363,5 +2371,7 @@
 !!!
 !!! nrlmd le 13/06/2011
-          delta_tsurf(i,nsrf)=y_delta_tsurf(j)*ypct(j)
+!!jyg20170131          delta_tsurf(i,nsrf)=y_delta_tsurf(j)*ypct(j)
+          delta_tsurf(i,nsrf)=y_delta_tsurf(j)
+!
           cdragh_x(i) = cdragh_x(i) + ycdragh_x(j)*ypct(j)
           cdragh_w(i) = cdragh_w(i) + ycdragh_w(j)*ypct(j)
@@ -2573,5 +2583,5 @@
 !      print*, tair1,yt(:,1),y_d_t(:,1)
 
-! Calculate the temperatureflag_pbl_surface_t2m_bugiflag_pbl_surface_t2m_bug et relative humidity at 2m and the wind at 10m 
+! Calculate the temperature and relative humidity at 2m and the wind at 10m 
 !!! jyg le 07/02/2012
        IF (iflag_split .eq.0) THEN
@@ -3089,4 +3099,13 @@
     IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
     IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
+
+!jyg<
+!****************************************************************************************
+! Deallocate variables for pbl splitting
+!
+!****************************************************************************************
+
+    CALL wx_pbl_final
+!>jyg
 
   END SUBROUTINE pbl_surface_final
