Index: /LMDZ6/trunk/libf/phylmd/dyn1d/compar1d.h
===================================================================
--- /LMDZ6/trunk/libf/phylmd/dyn1d/compar1d.h	(revision 3887)
+++ /LMDZ6/trunk/libf/phylmd/dyn1d/compar1d.h	(revision 3888)
@@ -8,4 +8,5 @@
       real :: nat_surf
       real :: tsurf
+      real :: beta_surf
       real :: rugos
       real :: rugosh
@@ -45,5 +46,5 @@
       real    :: p_nudging_u, p_nudging_v, p_nudging_w, p_nudging_t, p_nudging_qv
       common/com_par1d/                                                 &
-     & nat_surf,tsurf,rugos,rugosh,                                     &
+     & nat_surf,tsurf,beta_surf,rugos,rugosh,                           &
      & xqsol,qsurf,psurf,zsurf,albedo,time,time_ini,xlat,xlon,airefi,   &
      & wtsurf,wqsurf,restart_runoff,xagesno,qsolinp,zpicinp,            &
Index: /LMDZ6/trunk/libf/phylmd/dyn1d/old_lmdz1d.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/dyn1d/old_lmdz1d.F90	(revision 3887)
+++ /LMDZ6/trunk/libf/phylmd/dyn1d/old_lmdz1d.F90	(revision 3888)
@@ -11,5 +11,5 @@
        du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
        falb_dir, falb_dif, &
-       ftsol, pbl_tke, pctsrf, radsol, rain_fall, snow_fall, ratqs, &
+       ftsol, beta_aridity, pbl_tke, pctsrf, radsol, rain_fall, snow_fall, ratqs, &
        rnebcon, rugoro, sig1, w01, solaire_etat0, sollw, sollwdown, &
        solsw, t_ancien, q_ancien, u_ancien, v_ancien, wake_cstar, &
@@ -656,4 +656,5 @@
       qsol = qsolinp
       qsurf = fq_sat(tsurf,psurf/100.)
+      beta_aridity(:,:) = beta_surf
       day1= day_ini
       time=daytime-day
@@ -795,6 +796,6 @@
 
         fder=0.
+        snsrf(1,:)=snowmass ! masse de neige des sous surface
         print *, 'snsrf', snsrf
-        snsrf(1,:)=snowmass ! masse de neige des sous surface
         qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
         fevap=0.
Index: /LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90	(revision 3887)
+++ /LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90	(revision 3888)
@@ -7,5 +7,5 @@
        du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
        falb_dir, falb_dif, &
-       ftsol, pbl_tke, pctsrf, radsol, rain_fall, snow_fall, ratqs, &
+       ftsol, beta_aridity, pbl_tke, pctsrf, radsol, rain_fall, snow_fall, ratqs, &
        rnebcon, rugoro, sig1, w01, solaire_etat0, sollw, sollwdown, &
        solsw, t_ancien, q_ancien, u_ancien, v_ancien, wake_cstar, &
@@ -429,4 +429,5 @@
       qsol = qsolinp
       qsurf = fq_sat(tsurf,psurf/100.)
+      beta_aridity(:,:) = beta_surf
       day1= day_ini
       time=daytime-day
Index: /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 3887)
+++ /LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 3888)
@@ -26,8 +26,8 @@
   USE cdrag_mod
   USE stdlevvar_mod
-  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
+  USE wx_pbl_var_mod,      ONLY : wx_pbl_init, wx_pbl_final, &
+                                  wx_pbl_prelim_0, wx_pbl_prelim_beta
+  USE wx_pbl_mod,          ONLY : wx_pbl0_merge, wx_pbl_split, wx_pbl_dts_merge, &
+                                  wx_pbl_check, wx_pbl_dts_check, wx_evappot
   use config_ocean_skin_m, only: activate_ocean_skin
 
@@ -37,10 +37,13 @@
   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: fder   ! flux drift
   !$OMP THREADPRIVATE(fder)
-  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: snow   ! snow at surface
+  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE    :: snow   ! snow at surface
   !$OMP THREADPRIVATE(snow)
   REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: qsurf  ! humidity at surface
   !$OMP THREADPRIVATE(qsurf)
-  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ftsoil ! soil temperature
+  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE          :: ftsoil ! soil temperature
   !$OMP THREADPRIVATE(ftsoil)
+  REAL, ALLOCATABLE, DIMENSION(:), SAVE              :: ydTs0, ydqs0  
+                                                     ! nul forced temperature and humidity differences
+  !$OMP THREADPRIVATE(ydTs0, ydqs0)
 
   INTEGER, SAVE :: iflag_pbl_surface_t2m_bug
@@ -99,8 +102,16 @@
     IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
 
+    ALLOCATE(ydTs0(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
+
+    ALLOCATE(ydqs0(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
+
     fder(:)       = fder_rst(:)
     snow(:,:)     = snow_rst(:,:)
     qsurf(:,:)    = qsurf_rst(:,:)
     ftsoil(:,:,:) = ftsoil_rst(:,:,:)
+    ydTs0(:) = 0.
+    ydqs0(:) = 0.
 
 !****************************************************************************************
@@ -183,4 +194,7 @@
        ts,SFRWL,   alb_dir, alb_dif,ustar, u10m, v10m,wstar, &
        cdragh,    cdragm,   zu1,    zv1,              &
+!jyg<   (26/09/2019)
+       beta, &
+!>jyg
        alb_dir_m,    alb_dif_m,  zxsens,   zxevap,    &
        alb3_lic,  runoff,    snowhgt,   qsnow,     to_ice,    sissnow,  &
@@ -206,5 +220,6 @@
        s_therm,   s_trmb1,   s_trmb2,  s_trmb3,       &
        zustar,zu10m,  zv10m,    fder_print,    &
-       zxqsurf,   rh2m,      zxfluxu,  zxfluxv,       &
+       zxqsurf, delta_qsurf,                       &
+       rh2m,      zxfluxu,  zxfluxv,               &
        z0m, z0h,   agesno,  sollw,    solsw,         &
        d_ts,      evap,    fluxlat,  t2m,           &
@@ -366,4 +381,7 @@
 ! Input/Output variables
 !****************************************************************************************
+!jyg<
+    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: beta    ! Aridity factor
+!>jyg
     REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ts      ! temperature at surface (K)
     REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: delta_tsurf !surface temperature difference between
@@ -468,4 +486,5 @@
     REAL, DIMENSION(klon),        INTENT(OUT)       :: fder_print ! fder for printing (=fder(i) + dflux_t(i) + dflux_q(i))
     REAL, DIMENSION(klon),        INTENT(OUT)       :: zxqsurf    ! humidity at surface, mean for each grid point
+    REAL, DIMENSION(klon),        INTENT(OUT)       :: delta_qsurf! humidity difference at surface, mean for each grid point
     REAL, DIMENSION(klon),        INTENT(OUT)       :: rh2m       ! relative humidity at 2m
     REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zxfluxu    ! u wind tension, mean for each grid point
@@ -540,8 +559,10 @@
     REAL, DIMENSION(klon)              :: r_co2_ppm     ! taux CO2 atmosphere
     REAL, DIMENSION(klon)              :: yts, yz0m, yz0h, ypct
+    REAL, DIMENSION(klon)              :: yz0h_old
 !albedo SB >>>
     REAL, DIMENSION(klon)              :: yalb,yalb_vis
 !albedo SB <<<
     REAL, DIMENSION(klon)              :: yt1, yq1, yu1, yv1
+    REAL, DIMENSION(klon)              :: yqa
     REAL, DIMENSION(klon)              :: ysnow, yqsurf, yagesno, yqsol
     REAL, DIMENSION(klon)              :: yrain_f, ysnow_f
@@ -577,8 +598,9 @@
     REAL, DIMENSION(klon)              :: yz0h_oupas
     REAL, DIMENSION(klon)              :: yfluxsens
+    REAL, DIMENSION(klon)              :: AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0
     REAL, DIMENSION(klon)              :: AcoefH, AcoefQ, BcoefH, BcoefQ
     REAL, DIMENSION(klon)              :: AcoefU, AcoefV, BcoefU, BcoefV
     REAL, DIMENSION(klon)              :: ypsref
-    REAL, DIMENSION(klon)              :: yevap, ytsurf_new, yalb3_new
+    REAL, DIMENSION(klon)              :: yevap, yevap_pot, ytsurf_new, yalb3_new
 !albedo SB >>>
     REAL, DIMENSION(klon,nsw)          :: yalb_dir_new, yalb_dif_new
@@ -592,5 +614,5 @@
     REAL, DIMENSION(klon,klev)         :: y_flux_u, y_flux_v
     REAL, DIMENSION(klon,klev)         :: ycoefh, ycoefm,ycoefq
-    REAL, DIMENSION(klon)              :: ycdragh, ycdragm
+    REAL, DIMENSION(klon)              :: ycdragh, ycdragq, ycdragm
     REAL, DIMENSION(klon,klev)         :: yu, yv
     REAL, DIMENSION(klon,klev)         :: yt, yq
@@ -624,5 +646,6 @@
     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)              :: ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w
+    REAL, DIMENSION(klon)              :: ycdragm_x, ycdragm_w
     REAL, DIMENSION(klon)              :: AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x
     REAL, DIMENSION(klon)              :: AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w
@@ -645,5 +668,8 @@
     REAL                               :: zx_qs_surf, zcor_surf, zdelta_surf
     REAL, DIMENSION(klon)              :: ytsurf_th, yqsatsurf
+!jyg<
     REAL, DIMENSION(klon)              :: ybeta
+    REAL, DIMENSION(klon)              :: ybeta_prev
+!>jyg
     REAL, DIMENSION(klon, klev)        :: d_u_x
     REAL, DIMENSION(klon, klev)        :: d_u_w
@@ -780,5 +806,14 @@
 !!! 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, DIMENSION(klon)              :: y_delta_tsurf, y_delta_tsurf_new
+    REAL, DIMENSION(klon)              :: delta_coef, tau_eq
+    REAL, DIMENSION(klon)              :: HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn
+    REAL, DIMENSION(klon)              :: phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0
+    REAL, DIMENSION(klon)              :: y_delta_qsurf
+    REAL, DIMENSION(klon)              :: y_delta_qsats
+    REAL, DIMENSION(klon)              :: yg_T, yg_Q
+    REAL, DIMENSION(klon)              :: yGamma_dTs_phiT, yGamma_dQs_phiQ
+    REAL, DIMENSION(klon)              :: ydTs_ins, ydqs_ins
+!
     REAL, PARAMETER                    :: facteur=2./sqrt(3.14)
     REAL, PARAMETER                    :: inertia=2000.
@@ -793,5 +828,7 @@
     REAL, DIMENSION(klon)              :: Kech_m
     REAL, DIMENSION(klon)              :: Kech_m_x, Kech_m_w 
-    REAL, DIMENSION(klon)              :: yts_x,yts_w
+    REAL, DIMENSION(klon)              :: yts_x, yts_w
+    REAL, DIMENSION(klon)              :: yqsatsrf0_x, yqsatsrf0_w
+    REAL, DIMENSION(klon)              :: yqsurf_x, yqsurf_w
 !jyg<
 !!    REAL, DIMENSION(klon)              :: Kech_Hp, Kech_H_xp, Kech_H_wp
@@ -800,7 +837,7 @@
 !!    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                               :: fact_cdrag
+    REAL                               :: z1lay
 
     REAL                               :: vent
@@ -950,4 +987,5 @@
  fder_print(:)=0.
  zxqsurf(:)=0.
+ delta_qsurf(:) = 0.
  zxfluxu(:,:)=0. ; zxfluxv(:,:)=0.
  solsw(:,:)=0. ; sollw(:,:)=0.
@@ -1032,6 +1070,11 @@
     y_delta_flux_t1=0.
     ydtsurf_th=0.
-    yts_x=0.      ; yts_w=0.
-    y_delta_tsurf=0.
+    yts_x(:)=0.      ; yts_w(:)=0.
+    y_delta_tsurf(:)=0. ; y_delta_qsurf(:)=0.
+    yqsurf_x(:)=0.      ; yqsurf_w(:)=0.
+    yg_T(:) = 0. ;        yg_Q(:) = 0.
+    yGamma_dTs_phiT(:) = 0. ; yGamma_dQs_phiQ(:) = 0.
+    ydTs_ins(:) = 0. ; ydqs_ins(:) = 0.
+
 !!!
     ytsoil = 999999. 
@@ -1302,4 +1345,21 @@
 !****************************************************************************************
 
+!
+!jyg<    (20190926)
+!   Provisional : set ybeta to standard values
+       IF (nsrf .NE. is_ter) THEN
+           ybeta(:) = 1.
+       ELSE
+           IF (iflag_split .EQ. 0) THEN
+              ybeta(:) = 1.
+           ELSE
+             DO j = 1, knon
+                i = ni(j)
+                ybeta(j)   = beta(i,nsrf)
+             ENDDO
+           ENDIF  ! (iflag_split .LE.1)
+       ENDIF !  (nsrf .NE. is_ter)
+!>jyg
+!
        DO j = 1, knon
           i = ni(j)
@@ -1521,5 +1581,5 @@
         CALL cdrag(knon, nsrf, &
             speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),&
-            yts_x, yqsurf, yz0m, yz0h, &
+            yts_x, yqsurf_x, yz0m, yz0h, &
             ycdragm_x, ycdragh_x, zri1_x, pref_x )
 
@@ -1548,5 +1608,5 @@
         CALL cdrag(knon, nsrf, &
             speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),&
-            yts_w, yqsurf, yz0m, yz0h, &
+            yts_w, yqsurf_w, yz0m, yz0h, &
             ycdragm_w, ycdragh_w, zri1_w, pref_w )
 !
@@ -1621,5 +1681,5 @@
        ENDIF
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
-            ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf, ycdragm_x, &
+            ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf_x, ycdragm_x, &
             ycoefm_x, ycoefh_x, ytke_x,y_treedrg)
 !            ycoefm_x, ycoefh_x, ytke_x)
@@ -1649,5 +1709,5 @@
        ENDIF
         CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
-            ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf, ycdragm_w, &
+            ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf_w, ycdragm_w, &
             ycoefm_w, ycoefh_w, ytke_w,y_treedrg)
 !            ycoefm_w, ycoefh_w, ytke_w)
@@ -1786,77 +1846,110 @@
          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, &
+!
+! Cdragq computation
+! ------------------
+    !******************************************************************************
+    ! Cdragq computed from cdrag
+    ! The difference comes only from a factor (f_z0qh_oce) on z0, so that
+    ! it can be computed inside wx_pbl0_merge
+    ! More complicated appraches may require the propagation through
+    ! pbl_surface of an independant cdragq variable.
+    !******************************************************************************
+!
+    IF ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce) THEN
+       ! Si on suit les formulations par exemple de Tessel, on 
+       ! a z0h=0.4*nu/u*, z0q=0.62*nu/u*, d'ou f_z0qh_oce=0.62/0.4=1.55
+!!       ycdragq_x(1:knon)=ycdragh_x(1:knon)*                                      &
+!!            log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
+!!       ycdragq_w(1:knon)=ycdragh_w(1:knon)*                                      &
+!!            log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
+!
+       DO j = 1,knon
+         z1lay = zgeo1(j)/RG
+         fact_cdrag = log(z1lay/yz0h(j))/log(z1lay/(f_z0qh_oce*yz0h(j)))
+         ycdragq_x(j)=ycdragh_x(j)*fact_cdrag
+         ycdragq_w(j)=ycdragh_w(j)*fact_cdrag
+!!     Print *,'YYYYpbl0: fact_cdrag ', fact_cdrag
+       ENDDO  ! j = 1,knon
+!
+!!  Print *,'YYYYpbl0: z1lay, yz0h, f_z0qh_oce, ycdragh_w, ycdragq_w ', &
+!!                z1lay, yz0h(1:knon), f_z0qh_oce, ycdragh_w(1:knon), ycdragq_w(1:knon)
+    ELSE
+       ycdragq_x(1:knon)=ycdragh_x(1:knon)
+       ycdragq_w(1:knon)=ycdragh_w(1:knon)
+    ENDIF  ! ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce)
+!
+         CALL wx_pbl_prelim_0(knon, nsrf, dtime, ypplay, ypaprs, ywake_s,  &
+                         yts, y_delta_tsurf, ygustiness, &
                          yt_x, yt_w, yq_x, yq_w, &
                          yu_x, yu_w, yv_x, yv_w, &
-                         ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w, &
+                         ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_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  &
+                         )
+         CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta,  &
+                         BcoefQ_x, BcoefQ_w  &
+                         )
+         CALL wx_pbl0_merge(knon, ypplay, ypaprs,  &
+                         ywake_s, ydTs0, ydqs0, &
+                         yt_x, yt_w, yq_x, yq_w, &
+                         yu_x, yu_w, yv_x, yv_w, &
+                         ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_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, &
+                         AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &
+                         BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
+                         ycdragh, ycdragq, 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 (iflag_split .eq. 2 .AND. nsrf .ne. is_oce) THEN
+           CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
+                           ywake_s, ybeta, ywake_cstar, ywake_dens, &
+                           AcoefH_x, AcoefH_w, &
+                           BcoefH_x, BcoefH_w, &
+                           AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
+                           AcoefH, AcoefQ, BcoefH, BcoefQ,  &
+                           HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
+                           phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
+                           yg_T, yg_Q, &
+                           yGamma_dTs_phiT, yGamma_dQs_phiQ, &
+                           ydTs_ins, ydqs_ins &
+                           )
+         ELSE !
+           AcoefH(:) = AcoefH_0(:)
+           AcoefQ(:) = AcoefQ_0(:)
+           BcoefH(:) = BcoefH_0(:)
+           BcoefQ(:) = BcoefQ_0(:)
+           yg_T(:) = 0.
+           yg_Q(:) = 0.
+           yGamma_dTs_phiT(:) = 0.
+           yGamma_dQs_phiQ(:) = 0.
+           ydTs_ins(:) = 0.
+           ydqs_ins(:) = 0.
+         ENDIF   ! (iflag_split .eq. 2)
+       ENDIF  ! (iflag_split .eq.0)
 !!!
        IF (prt_level >=10) THEN
-         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)
+         PRINT *,'pbl_surface (merge->): yt(1,:) ',yt(1,:)
+         PRINT *,'pbl_surface (merge->): yq(1,:) ',yq(1,:)
+         PRINT *,'pbl_surface (merge->): yu(1,:) ',yu(1,:)
+         PRINT *,'pbl_surface (merge->): yv(1,:) ',yv(1,:)
+         PRINT *,'pbl_surface (merge->): AcoefH(1), AcoefQ(1), AcoefU(1), AcoefV(1) ', &
+                                         AcoefH(1), AcoefQ(1), AcoefU(1), AcoefV(1)
+         PRINT *,'pbl_surface (merge->): BcoefH(1), BcoefQ(1), BcoefU(1), BcoefV(1) ', &
+                                         BcoefH(1), BcoefQ(1), BcoefU(1), BcoefV(1)
+
        ENDIF
 
+!  Save initial value of z0h for use in evappot (z0h wiil be computed again in the surface models)
+          yz0h_old(1:knon) = yz0h(1:knon)
+!
 !****************************************************************************************
 !
@@ -2117,35 +2210,102 @@
           y_flux_q1(j) = -yevap(j)
           ENDDO
-        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*,'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
-!!!!!
-!!!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
+        ENDIF ! (ok_flux_surf)
+!
+! ------------------------------------------------------------------------------
+! 12a)  Splitting
+! ------------------------------------------------------------------------------
+
        IF (iflag_split .GE. 1) THEN
-         CALL wx_pbl0_split(knon, dtime, ywake_s, &
+!
+         IF (nsrf .ne. is_oce) THEN
+!
+!         Compute potential evaporation and aridity factor  (jyg, 20200328)
+          ybeta_prev(:) = ybeta(:)
+             DO j = 1, knon
+               yqa(j) = AcoefQ(j) - BcoefQ(j)*yevap(j)*dtime
+             ENDDO
+!
+          CALL wx_evappot(knon, yqa, yTsurf_new, yevap_pot)
+!
+          ybeta(1:knon) = min(yevap(1:knon)/yevap_pot(1:knon), 1.)
+          
+          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_prev, beta, ytsurf_new', ybeta_prev(j), ybeta(j), ytsurf_new(j)
+            print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j)
+           ENDDO
+          ENDIF  ! (prt_level >=10)
+!
+! Second call to wx_pbl0_merge and wx_pbl_dts_merge in order to take into account 
+! the update of the aridity coeficient beta.
+!
+        CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta,  &
+                        BcoefQ_x, BcoefQ_w  &
+                        )
+        CALL wx_pbl0_merge(knon, ypplay, ypaprs,  &
+                          ywake_s, ydTs0, ydqs0, &
+                          yt_x, yt_w, yq_x, yq_w, &
+                          yu_x, yu_w, yv_x, yv_w, &
+                          ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_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_0, AcoefQ_0, AcoefU, AcoefV, &
+                          BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
+                          ycdragh, ycdragq, ycdragm, &
+                          yt1, yq1, yu1, yv1 &
+                          )
+          IF (iflag_split .eq. 2) THEN
+            CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
+                            ywake_s, ybeta, ywake_cstar, ywake_dens, &
+                            AcoefH_x, AcoefH_w, &
+                            BcoefH_x, BcoefH_w, &
+                            AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
+                            AcoefH, AcoefQ, BcoefH, BcoefQ,  &
+                            HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
+                            phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
+                            yg_T, yg_Q, &
+                            yGamma_dTs_phiT, yGamma_dQs_phiQ, &
+                            ydTs_ins, ydqs_ins &
+                            )
+          ELSE !
+            AcoefH(:) = AcoefH_0(:)
+            AcoefQ(:) = AcoefQ_0(:)
+            BcoefH(:) = BcoefH_0(:)
+            BcoefQ(:) = BcoefQ_0(:)
+            yg_T(:) = 0.
+            yg_Q(:) = 0.
+            yGamma_dTs_phiT(:) = 0.
+            yGamma_dQs_phiQ(:) = 0.
+            ydTs_ins(:) = 0.
+            ydqs_ins(:) = 0.
+          ENDIF   ! (iflag_split .eq. 2)
+!
+        ELSE    ! (nsrf .ne. is_oce)
+          ybeta(1:knon) = 1.
+          yevap_pot(1:knon) = yevap(1:knon)
+          AcoefH(:) = AcoefH_0(:)
+          AcoefQ(:) = AcoefQ_0(:)
+          BcoefH(:) = BcoefH_0(:)
+          BcoefQ(:) = BcoefQ_0(:)
+          yg_T(:) = 0.
+          yg_Q(:) = 0.
+          yGamma_dTs_phiT(:) = 0.
+          yGamma_dQs_phiQ(:) = 0.
+          ydTs_ins(:) = 0.
+          ydqs_ins(:) = 0.
+        ENDIF   ! (nsrf .ne. is_oce)
+! 
+        CALL wx_pbl_split(knon, nsrf, dtime, ywake_s, ybeta, iflag_split, &
+                       yg_T, yg_Q, &
+                       yGamma_dTs_phiT, yGamma_dQs_phiQ, &
+                       ydTs_ins, ydqs_ins, &
                        y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &
+!!!!                       HTRn_b, dd_HTRn, HTphiT_b, dd_HTphiT, &
+                       phiQ0_b, phiT0_b, &
                        y_flux_t1_x, y_flux_t1_w, &
                        y_flux_q1_x, y_flux_q1_w, &
@@ -2153,7 +2313,49 @@
                        y_flux_v1_x, y_flux_v1_w, &
                        yfluxlat_x, yfluxlat_w, &
-                       y_delta_tsurf &
+                       y_delta_qsats, &
+                       y_delta_tsurf_new, y_delta_qsurf &
                        )
+!
+         CALL wx_pbl_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
+                       yTs, y_delta_tsurf,  & 
+                       yqsurf, yTsurf_new,  &
+                       y_delta_tsurf_new, y_delta_qsats,  &
+                       AcoefH_x, AcoefH_w, &
+                       BcoefH_x, BcoefH_w, &
+                       AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
+                       AcoefH, AcoefQ, BcoefH, BcoefQ,  &
+                       y_flux_t1, y_flux_q1,  &
+                       y_flux_t1_x, y_flux_t1_w, &
+                       y_flux_q1_x, y_flux_q1_w)
+!
+         IF (nsrf .ne. is_oce) THEN
+           CALL wx_pbl_dts_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
+                         yTs, y_delta_tsurf,  & 
+                         yqsurf, yTsurf_new,  &
+                         y_delta_qsats, y_delta_tsurf_new, y_delta_qsurf,  &
+                         AcoefH_x, AcoefH_w, &
+                         BcoefH_x, BcoefH_w, &
+                         AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
+                         AcoefH, AcoefQ, BcoefH, BcoefQ,  &
+                         HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
+                         phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
+                         yg_T, yg_Q, &
+                         yGamma_dTs_phiT, yGamma_dQs_phiQ, &
+                         ydTs_ins, ydqs_ins, &
+                         y_flux_t1, y_flux_q1,  &
+                         y_flux_t1_x, y_flux_t1_w, &
+                         y_flux_q1_x, y_flux_q1_w )
+         ENDIF   ! (nsrf .ne. is_oce)
+!
+       ELSE  ! (iflag_split .ge. 1)
+         ybeta(1:knon) = 1.
+         yevap_pot(1:knon) = yevap(1:knon)
        ENDIF  ! (iflag_split .ge. 1)
+!
+       IF (prt_level >= 10) THEN
+         print *,'pbl_surface, ybeta , yevap, yevap_pot ', &
+                               ybeta , yevap, yevap_pot
+       ENDIF  ! (prt_level >= 10)
+!
 !>jyg
 !
@@ -2318,10 +2520,10 @@
        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
-
+!!
+!!        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
+!!
 !****************************************************************************************
 ! 13) Transform variables for output format : 
@@ -2438,4 +2640,5 @@
           i = ni(j)
           evap(i,nsrf) = - flux_q(i,1,nsrf)                  !jyg
+          beta(i,nsrf) = ybeta(j)                             !jyg
           d_ts(i,nsrf) = y_d_ts(j)
 !albedo SB >>>
@@ -2453,6 +2656,6 @@
           cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
           cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
-          dflux_t(i) = dflux_t(i) + y_dflux_t(j)
-          dflux_q(i) = dflux_q(i) + y_dflux_q(j)
+          dflux_t(i) = dflux_t(i) + y_dflux_t(j)*ypct(j)
+          dflux_q(i) = dflux_q(i) + y_dflux_q(j)*ypct(j)
        ENDDO
 
@@ -2470,5 +2673,8 @@
 !!! nrlmd le 13/06/2011
 !!jyg20170131          delta_tsurf(i,nsrf)=y_delta_tsurf(j)*ypct(j)
-          delta_tsurf(i,nsrf)=y_delta_tsurf(j)
+!!jyg20210118          delta_tsurf(i,nsrf)=y_delta_tsurf(j)
+          delta_tsurf(i,nsrf)=y_delta_tsurf_new(j)
+!
+          delta_qsurf(i) = delta_qsurf(i) + y_delta_qsurf(j)*ypct(j)
 !
           cdragh_x(i) = cdragh_x(i) + ycdragh_x(j)*ypct(j)
@@ -2671,5 +2877,6 @@
                * (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)
+!!          tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)
+          tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf_new(j)
           qairsol(j) = yqsurf(j)
         ENDDO
@@ -3268,4 +3475,6 @@
     IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
     IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
+    IF (ALLOCATED(ydTs0)) DEALLOCATE(ydTs0)
+    IF (ALLOCATED(ydqs0)) DEALLOCATE(ydqs0)
 
 !jyg<
Index: /LMDZ6/trunk/libf/phylmd/phyredem.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phyredem.F90	(revision 3887)
+++ /LMDZ6/trunk/libf/phylmd/phyredem.F90	(revision 3888)
@@ -12,5 +12,6 @@
   USE fonte_neige_mod,  ONLY : fonte_neige_final
   USE pbl_surface_mod,  ONLY : pbl_surface_final
-  USE phys_state_var_mod, ONLY: radpas, zmasq, pctsrf, ftsol, falb_dir,      &
+  USE phys_state_var_mod, ONLY: radpas, zmasq, pctsrf,                       &
+                                ftsol, beta_aridity, delta_tsurf, falb_dir,  &
                                 falb_dif, qsol, fevap, radsol, solsw, sollw, &
                                 sollwdown, rain_fall, snow_fall, z0m, z0h,   &
@@ -157,5 +158,11 @@
     END IF
 
+!    Surface variables 
     CALL put_field_srf1(pass,"TS","Temperature",ftsol(:,:))
+
+    CALL put_field_srf1(pass,"DELTA_TS","w-x surface temperature difference", delta_tsurf(:,:))
+
+    CALL put_field_srf1(pass,"BETA_S","Aridity factor", beta_aridity(:,:))
+!    End surface variables
 
 ! ================== Albedo =======================================
Index: /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 3887)
+++ /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 3888)
@@ -340,4 +340,6 @@
       REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxfluxlat_x, zxfluxlat_w
 !$OMP THREADPRIVATE(zxfluxlat_x, zxfluxlat_w)
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: delta_qsurf
+!$OMP THREADPRIVATE(delta_qsurf)
 !jyg<
 !!! Entrees supplementaires couche-limite
@@ -733,4 +735,5 @@
       ALLOCATE(sens_x(klon), sens_w(klon))
       ALLOCATE(zxfluxlat_x(klon), zxfluxlat_w(klon))
+      ALLOCATE(delta_qsurf(klon))
 !jyg<
 !!      ALLOCATE(t_x(klon,klev), t_w(klon,klev))
@@ -1032,4 +1035,5 @@
       DEALLOCATE(sens_x, sens_w)
       DEALLOCATE(zxfluxlat_x, zxfluxlat_w)
+      DEALLOCATE(delta_qsurf)
 !jyg<
 !!      DEALLOCATE(t_x, t_w)
Index: /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3887)
+++ /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3888)
@@ -814,6 +814,4 @@
 'flat_w', 'flat within_wake', 'W/m2', (/ ('', i=1, 10) /))
 !!
-  type(ctrl_out),save :: o_delta_tsurf    = ctrl_out((/ 1, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
-'delta_tsurf', 'Temperature difference (w-x)', 'K', (/ ('', i=1, 10) /))
   type(ctrl_out),save :: o_cdragh_x       = ctrl_out((/ 1, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
 'cdragh_x', 'cdragh off-wake', '', (/ ('', i=1, 10) /))
@@ -1094,4 +1092,16 @@
       ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'dltpbltke_sic',       &
       "TKE difference (w - x) "//clnsurf(4),"-", (/ ('', i=1, 10) /)) /)
+
+  TYPE(ctrl_out), SAVE :: o_delta_tsurf = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'delta_tsurf ', 'T_surf difference (w - x)', 'K', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_delta_tsurf_srf      = (/             &
+      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'delta_tsurf_ter',       &
+      "T_surf difference (w - x) "//clnsurf(1),"-", (/ ('', i=1, 10) /)), &
+      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'delta_tsurf_lic',       &
+      "T_surf difference (w - x) "//clnsurf(2),"-", (/ ('', i=1, 10) /)), &
+      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'delta_tsurf_oce',       &
+      "T_surf difference (w - x) "//clnsurf(3),"-", (/ ('', i=1, 10) /)), &
+      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'delta_tsurf_sic',       &
+      "T_surf difference (w - x) "//clnsurf(4),"-", (/ ('', i=1, 10) /)) /)
 
   TYPE(ctrl_out), SAVE :: o_kz = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
Index: /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 3887)
+++ /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 3888)
@@ -87,5 +87,5 @@
          o_dtvdf_x    , o_dtvdf_w    , o_dqvdf_x    , o_dqvdf_w    , &
          o_sens_x     , o_sens_w     , o_flat_x     , o_flat_w     , &
-         o_delta_tsurf, &
+         o_delta_tsurf, o_delta_tsurf_srf, &
          o_cdragh_x   , o_cdragh_w   , o_cdragm_x   , o_cdragm_w   , &
          o_kh         , o_kh_x       , o_kh_w       , &
@@ -1354,9 +1354,18 @@
 !
                CALL histwrite_phy(o_dqvdf_w    ,zx_tmp_fi3d) 
-               CALL histwrite_phy(o_sens_x     ,sens_x     ) 
-               CALL histwrite_phy(o_sens_w     ,sens_w     ) 
+       IF (vars_defined)  zx_tmp_fi2d(1:klon)=-1*sens_x(1:klon)
+               CALL histwrite_phy(o_sens_x     ,zx_tmp_fi2d) 
+       IF (vars_defined)  zx_tmp_fi2d(1:klon)=-1*sens_w(1:klon)
+               CALL histwrite_phy(o_sens_w     ,zx_tmp_fi2d) 
                CALL histwrite_phy(o_flat_x     ,zxfluxlat_x) 
                CALL histwrite_phy(o_flat_w     ,zxfluxlat_w) 
-               CALL histwrite_phy(o_delta_tsurf,delta_tsurf) 
+          zx_tmp_fi2d=0.
+          IF (vars_defined) THEN
+             DO nsrf=1,nbsrf
+                   zx_tmp_fi2d(:)=zx_tmp_fi2d(:) &
+                        +pctsrf(:,nsrf)*delta_tsurf(:,nsrf)
+             ENDDO
+          ENDIF
+               CALL histwrite_phy(o_delta_tsurf,zx_tmp_fi2d) 
                CALL histwrite_phy(o_cdragh_x   ,cdragh_x   ) 
                CALL histwrite_phy(o_cdragh_w   ,cdragh_w   ) 
Index: /LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 3887)
+++ /LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 3888)
@@ -32,4 +32,6 @@
       REAL, ALLOCATABLE, SAVE :: ftsol(:,:)
 !$OMP THREADPRIVATE(ftsol)
+      REAL, ALLOCATABLE, SAVE :: beta_aridity(:,:)
+!$OMP THREADPRIVATE(beta_aridity)
       REAL,ALLOCATABLE,SAVE :: qsol(:),fevap(:,:),z0m(:,:),z0h(:,:),agesno(:,:)
 !$OMP THREADPRIVATE(qsol,fevap,z0m,z0h,agesno)
@@ -96,8 +98,4 @@
       REAL, ALLOCATABLE, SAVE :: coefm(:,:,:) ! Kz momentum
 !$OMP THREADPRIVATE(pbl_tke, coefh,coefm)
-!nrlmd<
-      REAL, ALLOCATABLE, SAVE :: delta_tsurf(:,:) ! Surface temperature difference inside-outside cold pool
-!$OMP THREADPRIVATE(delta_tsurf)
-!>nrlmd
       REAL, ALLOCATABLE, SAVE :: zmax0(:), f0(:) ! 
 !$OMP THREADPRIVATE(zmax0,f0)
@@ -276,4 +274,8 @@
       REAL,ALLOCATABLE,SAVE :: wake_delta_pbl_TKE(:,:,:)
 !$OMP THREADPRIVATE(wake_delta_pbl_TKE)
+!nrlmd<
+      REAL, ALLOCATABLE, SAVE :: delta_tsurf(:,:) ! Surface temperature difference inside-outside cold pool
+!$OMP THREADPRIVATE(delta_tsurf)
+!>nrlmd
 !>jyg
 !
@@ -478,4 +480,5 @@
       ALLOCATE(pctsrf(klon,nbsrf))
       ALLOCATE(ftsol(klon,nbsrf))
+      ALLOCATE(beta_aridity(klon,nbsrf))
       ALLOCATE(qsol(klon),fevap(klon,nbsrf))
       ALLOCATE(z0m(klon,nbsrf+1),z0h(klon,nbsrf+1),agesno(klon,nbsrf))
@@ -674,4 +677,5 @@
 
       DEALLOCATE(pctsrf, ftsol, falb1, falb2)
+      DEALLOCATE(beta_aridity)
       DEALLOCATE(qsol,fevap,z0m,z0h,agesno)
 !FC
@@ -687,7 +691,4 @@
       DEALLOCATE(tr_ancien)                           !RomP
       DEALLOCATE(ratqs, pbl_tke,coefh,coefm)
-!nrlmd<
-      DEALLOCATE(delta_tsurf)
-!>nrlmd
       DEALLOCATE(zmax0, f0)
       DEALLOCATE(sig1, w01)
@@ -744,4 +745,7 @@
 !jyg<
       DEALLOCATE(wake_delta_pbl_TKE)
+!nrlmd<
+      DEALLOCATE(delta_tsurf)
+!>nrlmd
 !>jyg
       DEALLOCATE(pfrac_impa, pfrac_nucl)
Index: /LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 3887)
+++ /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 3888)
@@ -210,4 +210,5 @@
        zxrunofflic,                            &
        zxtsol, snow_lsc, zxfqfonte, zxqsurf,   &
+       delta_qsurf,                            &
        rain_lsc, rain_num,                     &
        !
@@ -2485,5 +2486,6 @@
     !   s_therm,   s_trmb1,   s_trmb2, s_trmb3,
     !   zu10m,     zv10m,   fder,
-    !   zxqsurf,   rh2m,      zxfluxu, zxfluxv,
+    !   zxqsurf,   delta_qsurf,
+    !   rh2m,      zxfluxu, zxfluxv,
     !   frugs,     agesno,    fsollw,  fsolsw,
     !   d_ts,      fevap,     fluxlat, t2m,
@@ -2546,4 +2548,5 @@
                                 !albedo SB <<<
             cdragh,    cdragm,  u1,    v1,            &
+            beta_aridity, &
                                 !albedo SB >>>
                                 ! albsol1,   albsol2,   sens,    evap,      &
@@ -2572,5 +2575,5 @@
             s_therm,   s_trmb1,   s_trmb2, s_trmb3, &
             zustar, zu10m,     zv10m,   fder, &
-            zxqsurf,   rh2m,      zxfluxu, zxfluxv, &
+            zxqsurf, delta_qsurf,   rh2m,      zxfluxu, zxfluxv, &
             z0m, z0h,     agesno,    fsollw,  fsolsw, &
             d_ts,      fevap,     fluxlat, t2m, &
@@ -4637,5 +4640,8 @@
 
     CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
-
+   !
+   ! Prevent pbl_tke_w from becoming negative
+    wake_delta_pbl_tke(:,:,:) = max(wake_delta_pbl_tke(:,:,:), -pbl_tke(:,:,:))
+   !
 
        ENDIF
Index: /LMDZ6/trunk/libf/phylmd/wx_pbl_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/wx_pbl_mod.F90	(revision 3887)
+++ /LMDZ6/trunk/libf/phylmd/wx_pbl_mod.F90	(revision 3888)
@@ -1,9 +1,8 @@
 MODULE wx_pbl_mod
 !
-! Planetary Boundary Layer and Surface module
-!
-! This module manage the calculation of turbulent diffusion in the boundary layer 
-! and all interactions towards the differents sub-surfaces.
-!
+! Split Planetary Boundary Layer
+!
+! This module manages the splitting of the boundary layer between two regions; the (w) 
+! region (inside cold pools) and the (x) region (outside cold pools)
 !
   USE dimphy
@@ -11,130 +10,14 @@
   IMPLICIT NONE
 
-  REAL, ALLOCATABLE, DIMENSION(:), SAVE        :: Kech_Tp, Kech_T_xp, Kech_T_wp
-  REAL, ALLOCATABLE, DIMENSION(:), SAVE        :: dd_KTp, KxKwTp, dd_AT, dd_BT
-!$OMP THREADPRIVATE(Kech_Tp, Kech_T_xp, Kech_T_wp, dd_KTp, KxKwTp, dd_AT, dd_BT)
-  REAL, ALLOCATABLE, DIMENSION(:), SAVE        :: Kech_Qp, Kech_Q_xp, Kech_Q_wp
-  REAL, ALLOCATABLE, DIMENSION(:), SAVE        :: dd_KQp, KxKwQp, dd_AQ, dd_BQ
-!$OMP THREADPRIVATE(Kech_Qp, Kech_Q_xp, Kech_Q_wp, dd_KQp, KxKwQp, dd_AQ, dd_BQ)
-  REAL, ALLOCATABLE, DIMENSION(:), SAVE        :: Kech_Up, Kech_U_xp, Kech_U_wp
-  REAL, ALLOCATABLE, DIMENSION(:), SAVE        :: dd_KUp, KxKwUp, dd_AU, dd_BU
-!$OMP THREADPRIVATE(Kech_Up, Kech_U_xp, Kech_U_wp, dd_KUp, KxKwUp, dd_AU, dd_BU)
-  REAL, ALLOCATABLE, DIMENSION(:), SAVE        :: Kech_Vp, Kech_V_xp, Kech_V_wp
-  REAL, ALLOCATABLE, DIMENSION(:), SAVE        :: dd_KVp, KxKwVp, dd_AV, dd_BV
-!$OMP THREADPRIVATE(Kech_Vp, Kech_V_xp, Kech_V_wp, dd_KVp, KxKwVp, dd_AV, dd_BV)
-
 CONTAINS
 !
 !****************************************************************************************
 !
-SUBROUTINE wx_pbl_init
-
-! Local variables
-!****************************************************************************************
-    INTEGER                       :: ierr
- 
-
-!****************************************************************************************
-! Allocate module variables
-!
-!****************************************************************************************    
-
-    ierr = 0
-
-    ALLOCATE(Kech_Tp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(Kech_T_xp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(Kech_T_wp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(dd_KTp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(KxKwTp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(dd_AT(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(dd_BT(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-!----------------------------------------------------------------------------
-    ALLOCATE(Kech_Qp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(Kech_Q_xp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(Kech_Q_wp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(dd_KQp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(KxKwQp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(dd_AQ(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(dd_BQ(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-!----------------------------------------------------------------------------
-    ALLOCATE(Kech_Up(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(Kech_U_xp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(Kech_U_wp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(dd_KUp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(KxKwUp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(dd_AU(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(dd_BU(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-!----------------------------------------------------------------------------
-    ALLOCATE(Kech_Vp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(Kech_V_xp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(Kech_V_wp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(dd_KVp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(KxKwVp(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(dd_AV(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-    ALLOCATE(dd_BV(klon), stat=ierr)
-    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
-
-!----------------------------------------------------------------------------
-
-END SUBROUTINE wx_pbl_init
-
-SUBROUTINE wx_pbl0_fuse(knon, dtime, ypplay, ywake_s, &
+SUBROUTINE wx_pbl0_merge(knon, ypplay, ypaprs,  &
+                                 sigw, dTs_forcing, dqs_forcing,  &
                                  yt_x, yt_w, yq_x, yq_w, &
                                  yu_x, yu_w, yv_x, yv_w, &
-                                 ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w, &
+                                 ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
+                                 ycdragm_x, ycdragm_w, &
                                  AcoefT_x, AcoefT_w, AcoefQ_x, AcoefQ_w, &
                                  AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
@@ -143,19 +26,29 @@
                                  AcoefT, AcoefQ, AcoefU, AcoefV, &
                                  BcoefT, BcoefQ, BcoefU, BcoefV, &
-                                 ycdragh, ycdragm, &
+                                 ycdragh, ycdragq, ycdragm, &
                                  yt1, yq1, yu1, yv1 &
                                  )
 !
+
+    USE wx_pbl_var_mod
+
     USE print_control_mod, ONLY: prt_level,lunout
+    USE indice_sol_mod, ONLY: is_oce
 !
     INCLUDE "YOMCST.h"
+    INCLUDE "FCTTRE.h"
+    INCLUDE "YOETHF.h"
+    INCLUDE "clesphys.h"
 !
     INTEGER,                      INTENT(IN)        :: knon    ! number of grid cells
-    REAL,                         INTENT(IN)        :: dtime   ! time step size (s)
     REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypplay  ! mid-layer pressure (Pa)
-    REAL, DIMENSION(knon),        INTENT(IN)        :: ywake_s ! cold pools fractional area
+    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypaprs  ! pressure at layer interfaces (pa)
+    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw ! cold pools fractional area
+    REAL, DIMENSION(knon),        INTENT(IN)        :: dTs_forcing ! forced temperature difference (w)-(x)
+    REAL, DIMENSION(knon),        INTENT(IN)        :: dqs_forcing ! forced humidity difference (w)-(x)
     REAL, DIMENSION(knon,klev),   INTENT(IN)        :: yt_x, yt_w, yq_x, yq_w
     REAL, DIMENSION(knon,klev),   INTENT(IN)        :: yu_x, yu_w, yv_x, yv_w
-    REAL, DIMENSION(knon),        INTENT(IN)        :: ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w
+    REAL, DIMENSION(knon),        INTENT(IN)        :: ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w
+    REAL, DIMENSION(knon),        INTENT(IN)        :: ycdragm_x, ycdragm_w
     REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT_x, AcoefT_w, AcoefQ_x, AcoefQ_w
     REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w
@@ -164,5 +57,5 @@
     REAL, DIMENSION(knon),        INTENT(OUT)       :: AcoefT, AcoefQ, AcoefU, AcoefV
     REAL, DIMENSION(knon),        INTENT(OUT)       :: BcoefT, BcoefQ, BcoefU, BcoefV
-    REAL, DIMENSION(knon),        INTENT(OUT)       :: ycdragh, ycdragm
+    REAL, DIMENSION(knon),        INTENT(OUT)       :: ycdragh, ycdragq, ycdragm
     REAL, DIMENSION(knon),        INTENT(OUT)       :: yt1, yq1, yu1, yv1  ! Apparent T, q, u, v at first level, as
                                                                            !seen by surface modules
@@ -170,10 +63,6 @@
 ! Local variables
     INTEGER                    :: j
-    REAL                       :: rho1
-    REAL                       :: mod_wind_x
-    REAL                       :: mod_wind_w   
-    REAL                       :: dd_Cdragh
-    REAL                       :: dd_Cdragm
     REAL                       :: dd_Kh
+    REAL                       :: dd_Kq
     REAL                       :: dd_Km
     REAL                       :: dd_u 
@@ -182,150 +71,79 @@
     REAL                       :: dd_q 
 !
-    REAL                       :: KCT, KCQ, KCU, KCV
-!
-    REAL                       :: BBT, BBQ, BBU, BBV
-    REAL                       :: DDT, DDQ, DDU, DDV
-    REAL                       :: LambdaT, LambdaQ, LambdaU, LambdaV 
     REAL                       :: LambdaTs, LambdaQs, LambdaUs, LambdaVs 
 !
     REAL, DIMENSION(knon)      :: sigx       ! fractional area of (x) region
-
-    REAL, DIMENSION(knon)      :: Kech_h    ! Energy exchange coefficient
-    REAL, DIMENSION(knon)      :: Kech_h_x, Kech_h_w
-    REAL, DIMENSION(knon)      :: Kech_m    ! Momentum exchange coefficient
-    REAL, DIMENSION(knon)      :: Kech_m_x, Kech_m_w
-
-!!!
-!!! jyg le 09/04/2013 ; passage aux nouvelles expressions en differences
-
-        sigx(:) = 1.-ywake_s(:)
-
+!
+!
+   sigx(1:knon) = 1.-sigw(1:knon)
+!                                           
+!
         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))
-         rho1 = ypplay(j,1)/(RD*(yt_x(j,1) + ywake_s(j)*(yt_w(j,1)-yt_x(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_T_xp(j) = Kech_h_x(j)/(1.-BcoefT_x(j)*Kech_h_x(j)*dtime)
-        Kech_T_wp(j) = Kech_h_w(j)/(1.-BcoefT_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_KTp(j) = Kech_T_wp(j) - Kech_T_xp(j)
-         dd_KQp(j) = Kech_Q_wp(j) - Kech_Q_xp(j)
-         dd_KUp(j) = Kech_U_wp(j) - Kech_U_xp(j)
-         dd_KVp(j) = Kech_V_wp(j) - Kech_V_xp(j)
-!
-        Kech_Tp(j) = Kech_T_xp(j) + ywake_s(j)*dd_KTp(j)
-        Kech_Qp(j) = Kech_Q_xp(j) + ywake_s(j)*dd_KQp(j)
-        Kech_Up(j) = Kech_U_xp(j) + ywake_s(j)*dd_KUp(j)
-        Kech_Vp(j) = Kech_V_xp(j) + ywake_s(j)*dd_KVp(j)
-!
-! Calcul des differences w-x
-       dd_Cdragm = ycdragm_w(j) - ycdragm_x(j)
-       dd_Cdragh = ycdragh_w(j) - ycdragh_x(j)
+!
+! Compute w-x differences
+       dd_t = yt_w(j,1) - yt_x(j,1)
+       dd_q = yq_w(j,1) - yq_x(j,1)
        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_AT(j) = AcoefT_w(j) - AcoefT_x(j)
-       dd_AQ(j) = AcoefQ_w(j) - AcoefQ_x(j)
-       dd_AU(j) = AcoefU_w(j) - AcoefU_x(j)
-       dd_AV(j) = AcoefV_w(j) - AcoefV_x(j)
-       dd_BT(j) = BcoefT_w(j) - BcoefT_x(j)
-       dd_BQ(j) = BcoefQ_w(j) - BcoefQ_x(j)
-       dd_BU(j) = BcoefU_w(j) - BcoefU_x(j)
-       dd_BV(j) = BcoefV_w(j) - BcoefV_x(j)
-!
-       KxKwTp(j) = Kech_T_xp(j)*Kech_T_wp(j)
-       KxKwQp(j) = Kech_Q_xp(j)*Kech_Q_wp(j)
-       KxKwUp(j) = Kech_U_xp(j)*Kech_U_wp(j)
-       KxKwVp(j) = Kech_V_xp(j)*Kech_V_wp(j)
-       BBT = (BcoefT_x(j) + sigx(j)*dd_BT(j))*dtime
-       BBQ = (BcoefQ_x(j) + sigx(j)*dd_BQ(j))*dtime
-       BBU = (BcoefU_x(j) + sigx(j)*dd_BU(j))*dtime
-       BBV = (BcoefV_x(j) + sigx(j)*dd_BV(j))*dtime
-       KCT = Kech_h(j)
-       KCQ = Kech_h(j)
-       KCU = Kech_m(j)
-       KCV = Kech_m(j)
-       DDT = Kech_Tp(j)
-       DDQ = Kech_Qp(j)
-       DDU = Kech_Up(j)
-       DDV = Kech_Vp(j)
-       LambdaT = dd_Kh/KCT
-       LambdaQ = dd_Kh/KCQ
-       LambdaU = dd_Km/KCU
-       LambdaV = dd_Km/KCV
-       LambdaTs = dd_KTp(j)/DDT
-       LambdaQs = dd_KQp(j)/DDQ
-       LambdaUs = dd_KUp(j)/DDU
-       LambdaVs = dd_KVp(j)/DDV
-!
-       IF (prt_level >=10) THEN
-          print *,'Variables pour la fusion : Kech_T_xp(j)' ,Kech_T_xp(j)
-          print *,'Variables pour la fusion : Kech_T_wp(j)' ,Kech_T_wp(j)
-          print *,'Variables pour la fusion : Kech_Tp(j)' ,Kech_Tp(j)
-          print *,'Variables pour la fusion : Kech_h(j)' ,Kech_h(j)
-       ENDIF
+!
+! Merged exchange coefficients
+         dd_Kh = Kech_h_w(j) - Kech_h_x(j)
+         dd_Kq = Kech_q_w(j) - Kech_q_x(j)
+         dd_Km = Kech_m_w(j) - Kech_m_x(j)
+!
+       LambdaTs = dd_KTp(j)/Kech_Tp(j)
+       LambdaQs = dd_KQs(j)/Kech_Qs(j)
+       LambdaUs = dd_KUp(j)/Kech_Up(j)
+       LambdaVs = dd_KVp(j)/Kech_Vp(j)
 !
 ! Calcul des coef A, B \'equivalents dans la couche 1
 !
-       AcoefT(j) = AcoefT_x(j) + ywake_s(j)*dd_AT(j)*(1.+sigx(j)*LambdaTs)
-       AcoefQ(j) = AcoefQ_x(j) + ywake_s(j)*dd_AQ(j)*(1.+sigx(j)*LambdaQs)
-       AcoefU(j) = AcoefU_x(j) + ywake_s(j)*dd_AU(j)*(1.+sigx(j)*LambdaUs)
-       AcoefV(j) = AcoefV_x(j) + ywake_s(j)*dd_AV(j)*(1.+sigx(j)*LambdaVs)
+! The dTs_forcing and dqs_forcing terms are added for diagnostic purpose ; they should be zero in normal operation.
+       AcoefT(j) = AcoefT_x(j) + sigw(j)*(1.+sigx(j)*LambdaTs)*(dd_AT(j) - C_p(j)*dTs_forcing(j))
+       AcoefQ(j) = AcoefQ_x(j) + sigw(j)*(1.+sigx(j)*LambdaQs)*(dd_AQ(j) - dqs_forcing(j))
+       AcoefU(j) = AcoefU_x(j) + sigw(j)*(1.+sigx(j)*LambdaUs)*dd_AU(j)
+       AcoefV(j) = AcoefV_x(j) + sigw(j)*(1.+sigx(j)*LambdaVs)*dd_AV(j)
 !                                           
-       BcoefT(j) = BcoefT_x(j) + ywake_s(j)*BcoefT_x(j)*sigx(j)*LambdaT*LambdaTs &
-                               + ywake_s(j)*dd_BT(j)*(1.+sigx(j)*LambdaT)*(1.+sigx(j)*LambdaTs)
-                                           
-       BcoefQ(j) = BcoefQ_x(j) + ywake_s(j)*BcoefQ_x(j)*sigx(j)*LambdaQ*LambdaQs &
-                               + ywake_s(j)*dd_BQ(j)*(1.+sigx(j)*LambdaQ)*(1.+sigx(j)*LambdaQs)
-                                           
-       BcoefU(j) = BcoefU_x(j) + ywake_s(j)*BcoefU_x(j)*sigx(j)*LambdaU*LambdaUs &
-                               + ywake_s(j)*dd_BU(j)*(1.+sigx(j)*LambdaU)*(1.+sigx(j)*LambdaUs)
-                                           
-       BcoefV(j) = BcoefV_x(j) + ywake_s(j)*BcoefV_x(j)*sigx(j)*LambdaV*LambdaVs &
-                               + ywake_s(j)*dd_BV(j)*(1.+sigx(j)*LambdaV)*(1.+sigx(j)*LambdaVs)
-
+!
+!!       BcoefT(j) = (sigw(j)*Kech_h_w(j)*Kech_T_pw(j)*BcoefT_w(j) + &
+!!                sigx(j)*Kech_h_x(j)*Kech_T_px(j)*BcoefT_x(j) )/(Kech_h(j)*Kech_Tp(j))
+!!       BcoefQ(j) = (sigw(j)*Kech_q_w(j)*Kech_Q_pw(j)*BcoefQ_w(j) + &
+!!                sigx(j)*Kech_q_x(j)*Kech_Q_px(j)*BcoefQ_x(j) )/(Kech_q(j)*Kech_Qp(j))
+!!       BcoefU(j) = (sigw(j)*Kech_m_w(j)*Kech_U_pw(j)*BcoefU_w(j) + &
+!!                sigx(j)*Kech_m_x(j)*Kech_U_px(j)*BcoefU_x(j) )/(Kech_m(j)*Kech_Up(j))
+!!       BcoefV(j) = (sigw(j)*Kech_m_w(j)*Kech_V_pw(j)*BcoefV_w(j) + &
+!!                sigx(j)*Kech_m_x(j)*Kech_V_px(j)*BcoefV_x(j) )/(Kech_m(j)*Kech_Vp(j))
+!
+!!  Print *,'YYYYpbl0: BcoefT_x, sigw, sigx, dd_Kh, dd_KTp, Kech_h_w ', &
+!!                     BcoefT_x, sigw, sigx, dd_Kh, dd_KTp, Kech_h_w
+!!  Print *,'YYYYpbl0: Kech_T_pw, dd_BT, Kech_h, Kech_Tp ', &
+!!                     Kech_T_pw, dd_BT, Kech_h, Kech_Tp
+       BcoefT(j) = BcoefT_x(j) + sigw(j)*(sigx(j)*dd_Kh*dd_KTp(j)*BcoefT_x(j) + &
+                                  Kech_h_w(j)*Kech_T_pw(j)*dd_BT(j))/(Kech_h(j)*Kech_Tp(j))
+       BcoefQ(j) = BcoefQ_x(j) + sigw(j)*(sigx(j)*dd_Kh*dd_KQs(j)*BcoefQ_x(j) + &
+                                  Kech_q_w(j)*Kech_Q_sw(j)*dd_BQ(j))/(Kech_q(j)*Kech_Qs(j))
+       BcoefU(j) = BcoefU_x(j) + sigw(j)*(sigx(j)*dd_Km*dd_KUp(j)*BcoefU_x(j) + &
+                                  Kech_m_w(j)*Kech_U_pw(j)*dd_BU(j))/(Kech_m(j)*Kech_Up(j))
+       BcoefV(j) = BcoefV_x(j) + sigw(j)*(sigx(j)*dd_Km*dd_KVp(j)*BcoefV_x(j) + &
+                                  Kech_m_w(j)*Kech_V_pw(j)*dd_BV(j))/(Kech_m(j)*Kech_Vp(j))
+!>jyg
+!
 !
 ! Calcul des cdrag \'equivalents dans la couche 
 !
-       ycdragm(j) = ycdragm_x(j) + ywake_s(j)*dd_Cdragm
-       ycdragh(j) = ycdragh_x(j) + ywake_s(j)*dd_Cdragh
+       ycdragm(j) = ycdragm_x(j) + sigw(j)*dd_Cdragm(j)
+       ycdragh(j) = ycdragh_x(j) + sigw(j)*dd_Cdragh(j)
+       ycdragq(j) = ycdragq_x(j) + sigw(j)*dd_Cdragq(j)
 !
 ! Calcul de T, q, u et v \'equivalents dans la couche 1
-!!       yt1(j) = yt_x(j,1) + ywake_s(j)*dd_t*(1.+sigx(j)*dd_Kh/KCT)
-!!       yq1(j) = yq_x(j,1) + ywake_s(j)*dd_q*(1.+sigx(j)*dd_Kh/KCQ)
-!!       yu1(j) = yu_x(j,1) + ywake_s(j)*dd_u*(1.+sigx(j)*dd_Km/KCU)
-!!       yv1(j) = yv_x(j,1) + ywake_s(j)*dd_v*(1.+sigx(j)*dd_Km/KCV)
-       yt1(j) = yt_x(j,1) + ywake_s(j)*dd_t
-       yq1(j) = yq_x(j,1) + ywake_s(j)*dd_q
-       yu1(j) = yu_x(j,1) + ywake_s(j)*dd_u
-       yv1(j) = yv_x(j,1) + ywake_s(j)*dd_v
+!!       yt1(j) = yt_x(j,1) + sigw(j)*dd_t*(1.+sigx(j)*dd_Kh/KCT)
+!!       yq1(j) = yq_x(j,1) + sigw(j)*dd_q*(1.+sigx(j)*dd_Kh/KCQ)
+!!       yu1(j) = yu_x(j,1) + sigw(j)*dd_u*(1.+sigx(j)*dd_Km/KCU)
+!!       yv1(j) = yv_x(j,1) + sigw(j)*dd_v*(1.+sigx(j)*dd_Km/KCV)
+       yt1(j) = yt_x(j,1) + sigw(j)*dd_t
+       yq1(j) = yq_x(j,1) + sigw(j)*dd_q
+       yu1(j) = yu_x(j,1) + sigw(j)*dd_u
+       yv1(j) = yv_x(j,1) + sigw(j)*dd_v
 
 
@@ -334,122 +152,1067 @@
         RETURN
 
-END SUBROUTINE wx_pbl0_fuse
-
-SUBROUTINE 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 &
-                       )
-!
+END SUBROUTINE wx_pbl0_merge
+
+SUBROUTINE wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
+                                 sigw, beta, wcstar, wdens, &
+                                 AT_x, AT_w, &
+                                 BT_x, BT_w, &
+                                 AcoefT0, AcoefQ0, BcoefT0, BcoefQ0, &
+                                 AcoefT,  AcoefQ,  BcoefT,  BcoefQ, &
+                                 HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
+                                 phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
+                                 g_T, g_Q, &
+                                 Gamma_phiT, Gamma_phiQ, &
+                                 dTs_ins, dqsatsrf_ins &
+                                 )
+!
+
+    USE wx_pbl_var_mod
+
     USE print_control_mod, ONLY: prt_level,lunout
 !
     INCLUDE "YOMCST.h"
+    INCLUDE "FCTTRE.h"
+    INCLUDE "YOETHF.h"
 !
     INTEGER,                      INTENT(IN)        :: knon    ! number of grid cells
     REAL,                         INTENT(IN)        :: dtime   ! time step size (s)
-    REAL, DIMENSION(knon),        INTENT(IN)        :: ywake_s ! cold pools fractional area
-    REAL, DIMENSION(knon),        INTENT(IN)        :: y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1
-!
-    REAL, DIMENSION(knon),        INTENT(OUT)       :: y_flux_t1_x, y_flux_t1_w
-    REAL, DIMENSION(knon),        INTENT(OUT)       :: y_flux_q1_x, y_flux_q1_w
-    REAL, DIMENSION(knon),        INTENT(OUT)       :: y_flux_u1_x, y_flux_u1_w
-    REAL, DIMENSION(knon),        INTENT(OUT)       :: y_flux_v1_x, y_flux_v1_w
-    REAL, DIMENSION(knon),        INTENT(OUT)       :: yfluxlat_x, yfluxlat_w
-    REAL, DIMENSION(knon),        INTENT(OUT)       :: y_delta_tsurf
+    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypplay  ! mid-layer pressure (Pa)
+    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypaprs  ! pressure at layer interfaces (pa)
+    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw    ! cold pool fractional area
+    REAL, DIMENSION(knon),        INTENT(IN)        :: beta    ! evaporation by potential evaporation
+    REAL, DIMENSION(knon),        INTENT(IN)        :: wcstar   ! cold pool gust front speed
+    REAL, DIMENSION(knon),        INTENT(IN)        :: wdens    ! cold pool number density
+    REAL, DIMENSION(knon),        INTENT(IN)        :: AT_x, AT_w
+    REAL, DIMENSION(knon),        INTENT(IN)        :: BT_x, BT_w
+    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT0, AcoefQ0, BcoefT0, BcoefQ0
+!
+    REAL, DIMENSION(knon),        INTENT(OUT)       :: AcoefT, AcoefQ, BcoefT, BcoefQ
+    REAL, DIMENSION(knon),        INTENT(OUT)       :: HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn
+    REAL, DIMENSION(knon),        INTENT(OUT)       :: phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0
+    REAL, DIMENSION(knon),        INTENT(OUT)       :: g_T, g_Q
+    REAL, DIMENSION(knon),        INTENT(OUT)       :: Gamma_phiT, Gamma_phiQ
+    REAL, DIMENSION(knon),        INTENT(OUT)       :: dTs_ins, dqsatsrf_ins
+!
+! Local variables
+    REAL, DIMENSION(knon)      :: qsat_x
+    REAL, DIMENSION(knon)      :: qsat_w
+    REAL, DIMENSION(knon)      :: dqsatdT_x
+    REAL, DIMENSION(knon)      :: dqsatdT_w
+!
+    REAL, DIMENSION(knon)      :: T10_x
+    REAL, DIMENSION(knon)      :: T10_w
+    REAL, DIMENSION(knon)      :: phiT0_x
+    REAL, DIMENSION(knon)      :: phiT0_w
+    REAL, DIMENSION(knon)      :: phiQ0_x
+    REAL, DIMENSION(knon)      :: phiQ0_w
+    REAL, DIMENSION(knon)      :: Rn0_x
+    REAL, DIMENSION(knon)      :: Rn0_w
+    REAL, DIMENSION(knon)      :: Rp1_x
+    REAL, DIMENSION(knon)      :: Rp1_w
+    REAL, DIMENSION(knon)      :: Rps_x
+    REAL, DIMENSION(knon)      :: Rps_w
+!
+    REAL, DIMENSION(knon)      :: HTphiT_x
+    REAL, DIMENSION(knon)      :: HTphiT_w
+    REAL, DIMENSION(knon)      :: HTphiQ_x
+    REAL, DIMENSION(knon)      :: HTphiQ_w
+    REAL, DIMENSION(knon)      :: HTRn_x
+    REAL, DIMENSION(knon)      :: HTRn_w
+!
+    REAL, DIMENSION(knon)      :: HQphiT_x
+    REAL, DIMENSION(knon)      :: HQphiT_w
+    REAL, DIMENSION(knon)      :: HQphiQ_x
+    REAL, DIMENSION(knon)      :: HQphiQ_w
+    REAL, DIMENSION(knon)      :: HQRn_x
+    REAL, DIMENSION(knon)      :: HQRn_w
+!
+    REAL, DIMENSION(knon)      :: HQphiT_b
+    REAL, DIMENSION(knon)      :: dd_HQphiT
+    REAL, DIMENSION(knon)      :: HQphiQ_b
+    REAL, DIMENSION(knon)      :: dd_HQphiQ
+    REAL, DIMENSION(knon)      :: HQRn_b
+    REAL, DIMENSION(knon)      :: dd_HQRn
+!
+
+    REAL, DIMENSION(knon)    :: sigx
+!
+    REAL, DIMENSION(knon)    :: Ts, T1
+!!!    REAL, DIMENSION(knon)    :: qsat, dqsat_dT
+!!!    REAL, DIMENSION(knon)    :: phiT0
+!
+!!!    REAL, DIMENSION(knon)    :: Cp, Lv
+    REAL, DIMENSION(knon)    :: tau, Inert 
+!
+    REAL                     :: dd_Kh
+    REAL                     :: zdelta, zcvm5, zcor
+    REAL                     :: qsat
+!
+    INTEGER                  :: j
+
+
+!----------------------------------------------------------------------------
+!  Reference state
+!  ---------------
+!   dqsat_dT_w = dqsat_dT(Ts0_w)                          dqsat_dT_x = dqsat_dT(Ts0_x)
+!   T10_w = (AT_w/Cp - Kech_T_w BT_w dtime Ts0_w)/(1 - Kech_T_w BT_w dtime)
+!                                                T10_x = (AT_x/Cp - Kech_T_x BT_x dtime Ts0_x)/(1 - Kech_T_x BT_x dtime)
+!   phiT0_w = Kech_T_pw (AT_w - Cp Ts0_w)                 phiT0_x = Kech_T_px (AT_x - Cp Ts0_x)
+!   phiQ0_w = Kech_Q_sw (beta AQ_w - qsatsrf0_w)          phiQ0_x = Kech_Q_sx (beta AQ_x - qsatsrf0_x)
+!   Rn0_w = eps_1 Rsigma T10_w^4 - Rsigma Ts0_w^4         Rn0_x = eps_1 Rsigma T10_x^4 - Rsigma Ts0_x^4
+!   Rp1_w = 4 eps_1 Rsigma T10_w^3                        Rp1_x = 4 eps_1 Rsigma T10_x^3
+!   Rps_w = 4 Rsigma Ts0_w^3                              Rps_x = 4 Rsigma Ts0_x^3
+!
+!   phiT0_b = sigw phiT0_w + sigx phiT0_x
+!   dphiT0 = phiT0_w - phiT0_x
+!   phiQ0_b = sigw phiQ0_w + sigx phiQ0_x
+!   dphiQ0 = phiQ0_w - phiQ0_x
+!   Rn0_b = sigw Rn0_w + sigx Rn0_x
+    dRn0 = Rn0_w - Rn0_x
+!
+!
+!----------------------------------------------------------------------------
+!  Elementary enthalpy equations
+!  -----------------------------
+!   phiT_w = phiT0_w - HTphiT_w (Ts_w-Ts0_w)            phiT_x = phiT0_x - HTphiT_x (Ts_x-Ts0_x)
+!   phiQ_w = phiQ0_w - HTphiQ_w (Ts_w-Ts0_w)            phiQ_x = phiQ0_x - HTphiQ_x (Ts_x-Ts0_x)
+!   Rn_w   = Rn0_w   - HTRn_w   (Ts_w-Ts0_w)            Rn_x   = Rn0_x   - HTRn_x   (Ts_x-Ts0_x)
+!  DFlux_DT coefficients
+!  ---------------------
+!   Heat flux equation
+!     HTphiT_w = Cp Kech_T_pw                            HTphiT_x = Cp Kech_T_px
+!   Moisture flux equation
+!     HTphiQ_w = beta Kech_Q_sw dqsat_dT_w               HTphiQ_x = beta Kech_Q_sx dqsat_dT_x
+!   Radiation equation
+!     HTRn_w = Rp1_w Kech_T_pw BcoefT_w dtime + Rps_w    HTRn_x = Rp1_x Kech_T_px BcoefT_x dtime + Rps_x
+!
+!----------------------------------------------------------------------------
+!  Elementary moisture equations
+!  -----------------------------
+!   beta Ts_w   = beta Ts0_w    + QQ_w     (qsatsrf_w-qsatsrf0_w)    beta Ts_x   = beta Ts0_x    + QQ_x     (qsatsrf_x-qsatsrf0_x)
+!   beta phiT_w = beta phiT0_w - HQphiT_w (qsatsrf_w-qsatsrf0_w)    beta phiQ_x = beta phiQ0_x - HTphiQ_x (qsatsrf_x-qsatsrf0_x)
+!   beta phiQ_w = beta phiQ0_w - HQphiQ_w (qsatsrf_w-qsatsrf0_w)    beta phiQ_x = beta phiQ0_x - HTphiQ_x (qsatsrf_x-qsatsrf0_x)
+!   beta Rn_w   = beta Rn0_w   - HQRn_w   (qsatsrf_w-qsatsrf0_w)    beta Rn_x   = beta Rn0_x   - HTRn_x   (qsatsrf_x-qsatsrf0_x)
+!  DFluxDQ coefficients
+!  ---------------------
+!   dqsat_dT equation
+!     QQ_w = 1. / dqsat_dT_w                             QQ_x = 1. / dqsat_dT_x
+!   Heat flux equation
+!     HQphiT_w = Cp Kech_T_pw QQ_w                       HQphiT_x = Cp Kech_T_px QQ_x
+!   Moisture flux equation
+!     HQphiQ_w = beta Kech_Q_sw                          HQphiQ_x = beta Kech_Q_sx
+!   Radiation equation
+!     HQRn_w = (Rp1_w Kech_T_pw BcoefT_w dtime + Rps_w) QQ_w
+!                                         HQRn_x = (Rp1_x Kech_T_px BcoefT_x dtime + Rps_x) QQ_x
+!
+!----------------------------------------------------------------------------
+! Mean values and w-x differences
+! -------------------------------
+!  HTphiT_b = sigw HTphiT_w + sigx HTphiT_x               dd_HTphiT = HTphiT_w - HTphiT_x
+!  HTphiQ_b = sigw HTphiQ_w + sigx HTphiQ_x               dd_HTphiQ = HTphiQ_w - HTphiQ_x
+!  HTRn_b   = sigw HTRn_w   + sigx HTRn_x                 dd_HTRn   = HTRn_w   - HTRn_x
+!
+!  QQ_b     = sigw QQ_w     + sigx QQ_x                   dd_QQ     = QQ_w     - QQ_x
+!  HQphiT_b = sigw HQphiT_w + sigx HQphiT_x               dd_HQphiT = HQphiT_w - HQphiT_x
+!  HQphiQ_b = sigw HQphiQ_w + sigx HQphiQ_x               dd_HQphiQ = HQphiQ_w - HQphiQ_x
+!  HQRn_b   = sigw HQRn_w   + sigx HQRn_x                 dd_HQRn   = HQRn_w   - HQRn_x
+!
+!----------------------------------------------------------------------------
+!  Equations
+!  ---------
+! (1 - g_T) dTs    = dTs_ins    + Gamma_phiT phiT
+! (1 - g_Q) dqsatsrf = dqsatsrf_ins + Gamma_phiQ phiQ
+!
+! Feedback Gains
+! --------------
+! g_T = - (sqrt(tau)/I) [ HTphiT_b + Lv HTphiQ_b + HTRn_b +  &
+!                        (dd_HTphiT + Lv dd_HTphiQ + dd_HTRn) (sigx - sigw - sigw sigx dd_HTphiT/HTphiT_b) ]
+! g_Q = - (sqrt(tau)/(I QQ_b)) ( HQphiT_b + Lv HQphiQ_b + HQRn_b ) -  &
+!         (sigx - sigw - sigw sigx dd_HQphiQ/HQphiQ_b)   &
+!                          [ dd_QQ/QQ_b + (sqrt(tau)/(I QQ_b))(dd_HQphiT + Lv dd_HQphiQ + dd_HQRn) ]
+!
+!  Ts, qs Coupling coefficients                /
+!  ----------------------------
+! Gamma_phiT = (sqrt(tau)/(I HTphiT_b)) (dd_HTphiT + Lv dd_HTphiQ + dd_HTRn)
+! Gamma_phiQ = (1/(HQphiQ_b QQ_b)) [ dd_QQ +  (sqrt(tau)/(I )) (dd_HQphiT + Lv dd_HQphiQ + dd_HQRn) ]
+!
+!  Insensitive changes
+!  -------------------
+! dTs_ins    = (1 - g_T) dTs0    - Gamma_phiT phiT0_b
+! dqsatsrf_ins = (1 - g_Q) dqsatsrf0 - Gamma_phiQ phiQ0_b
+!
+!----------------------------------------------------------------------------
+!  Effective coefficients Acoef and Bcoef
+!  --------------------------------------
+!  Equations
+!  ---------
+! Cp Ta = AcoefT + BcoefT phiT dtime
+!    qa = AcoefQ + BcoefQ phiQ dtime
+!  Coefficients
+!  ------------
+! AcoefT = AcoefT0 - sigw sigx (dd_KTp/Kech_Tp) Cp dTs_ins/(1 - g_T)
+! BcoefT = BcoefT0 - sigw sigx (dd_KTp/Kech_Tp) Cp Gamma_phiT/(1 - g_T)/dtime
+!
+! AcoefQ = AcoefQ0 - sigw sigx (dd_KQp/Kech_Qp) dqs_ins/(1 - g_Q)
+! BcoefQ = BcoefQ0 - sigw sigx (dd_KQp/Kech_Qp) Gamma_phiq/(1 - g_Q)/dtime
+!
+!==============================================================================
+!
+!
+!  Parameters
+!  ----------
+   Inert(1:knon) = 2000.
+   tau(1:knon) = sqrt(sigw(1:knon)/max(rpi*wdens(1:knon)*wcstar(1:knon)**2 , &
+                                       sigw(1:knon)*1.e-12,smallestreal))
+   sigx(1:knon) = 1.-sigw(1:knon)
+!! Compute Cp, Lv, qsat, dqsat_dT.
+!   C_p(1:knon) = RCpd
+!   L_v(1:knon) = RLvtt
+!
+!      print *,' AAAA wx_pbl_dTs, C_p(j), qsat0(j), Ts0(j) : ', C_p(:), qsat0(:), Ts0(:)
+!
+!
+   T10_x(1:knon) = (AT_x(1:knon)/C_p(1:knon) - Kech_h_x(1:knon)*BT_x(1:knon)*dtime*Ts0_x(1:knon))/  &
+                   (1 - Kech_h_x(1:knon)*BT_x(1:knon)*dtime)
+   T10_w(1:knon) = (AT_w(1:knon)/C_p(1:knon) - Kech_h_w(1:knon)*BT_w(1:knon)*dtime*Ts0_w(1:knon))/  &
+                   (1 - Kech_h_w(1:knon)*BT_w(1:knon)*dtime)
+!
+   phiT0_x(1:knon) = Kech_T_px(1:knon)*(AT_x(1:knon) - C_p(1:knon)*Ts0_x(1:knon))
+   phiT0_w(1:knon) = Kech_T_pw(1:knon)*(AT_w(1:knon) - C_p(1:knon)*Ts0_w(1:knon))
+!
+   phiQ0_x(1:knon) = Kech_Q_sx(1:knon)*(beta(1:knon)*AQ_x(1:knon) - qsatsrf0_x(1:knon))
+   phiQ0_w(1:knon) = Kech_Q_sw(1:knon)*(beta(1:knon)*AQ_w(1:knon) - qsatsrf0_w(1:knon))
+!
+   Rn0_x(1:knon) = eps_1*Rsigma*T10_x(1:knon)**4 - Rsigma*Ts0_x(1:knon)**4
+   Rn0_w(1:knon) = eps_1*Rsigma*T10_w(1:knon)**4 - Rsigma*Ts0_w(1:knon)**4
+!
+   Rp1_x(1:knon) = 4*eps_1*Rsigma*T10_x(1:knon)**3
+   Rp1_w(1:knon) = 4*eps_1*Rsigma*T10_w(1:knon)**3
+!
+   Rps_x(1:knon) = 4*Rsigma*Ts0_x(1:knon)**3
+   Rps_w(1:knon) = 4*Rsigma*Ts0_w(1:knon)**3
+!
+!  DFlux_DT coefficients
+!  ---------------------
+!   Heat flux equation
+     HTphiT_x(1:knon) = C_p(1:knon)*Kech_T_px(1:knon)
+     HTphiT_w(1:knon) = C_p(1:knon)*Kech_T_pw(1:knon)                       
+!   Moisture flux equation
+     HTphiQ_x(1:knon) = beta(1:knon)*Kech_Q_sx(1:knon)*dqsatdT0_x(1:knon)
+     HTphiQ_w(1:knon) = beta(1:knon)*Kech_Q_sw(1:knon)*dqsatdT0_w(1:knon)          
+!   Radiation equation
+     HTRn_x(1:knon) = Rp1_x(1:knon)*Kech_T_px(1:knon)*BT_x(1:knon)*dtime + Rps_x(1:knon)
+     HTRn_w(1:knon) = Rp1_w(1:knon)*Kech_T_pw(1:knon)*BT_w(1:knon)*dtime + Rps_w(1:knon)  
+!
+!  DFluxDQ coefficients
+!  ---------------------
+!   Heat flux equation
+     HQphiT_x(1:knon) = C_p(1:knon)*Kech_T_px(1:knon)*QQ_x(1:knon)
+     HQphiT_w(1:knon) = C_p(1:knon)*Kech_T_pw(1:knon)*QQ_w(1:knon)                 
+!   Moisture flux equation
+     HQphiQ_x(1:knon) = beta(1:knon)*Kech_Q_sx(1:knon)
+     HQphiQ_w(1:knon) = beta(1:knon)*Kech_Q_sw(1:knon)                   
+!   Radiation equation
+     HQRn_x(1:knon) = (Rp1_x(1:knon)*Kech_T_px(1:knon)*BT_x(1:knon)*dtime + Rps_x(1:knon))*QQ_x(1:knon)
+     HQRn_w(1:knon) = (Rp1_w(1:knon)*Kech_T_pw(1:knon)*BT_w(1:knon)*dtime + Rps_w(1:knon))*QQ_w(1:knon)
+!
+! Mean values and w-x differences
+! -------------------------------
+  phiT0_b(1:knon) = sigw(1:knon)*phiT0_w(1:knon) + sigx(1:knon)*phiT0_x(1:knon)           
+  phiQ0_b(1:knon) = sigw(1:knon)*phiQ0_w(1:knon) + sigx(1:knon)*phiQ0_x(1:knon)          
+  Rn0_b(1:knon)   = sigw(1:knon)*Rn0_w(1:knon)   + sigx(1:knon)*Rn0_x(1:knon)          
+!
+  dphiT0(1:knon) = phiT0_w(1:knon) - phiT0_x(1:knon)           
+  dphiQ0(1:knon) = phiQ0_w(1:knon) - phiQ0_x(1:knon)          
+  dRn0(1:knon)   = Rn0_w(1:knon)   - Rn0_x(1:knon)          
+!
+  HTphiT_b(1:knon) = sigw(1:knon)*HTphiT_w(1:knon) + sigx(1:knon)*HTphiT_x(1:knon)           
+  dd_HTphiT(1:knon) = HTphiT_w(1:knon) - HTphiT_x(1:knon)
+!
+  HTphiQ_b(1:knon) = sigw(1:knon)*HTphiQ_w(1:knon) + sigx(1:knon)*HTphiQ_x(1:knon)          
+  dd_HTphiQ(1:knon) = HTphiQ_w(1:knon) - HTphiQ_x(1:knon)
+!
+  HTRn_b(1:knon)   = sigw(1:knon)*HTRn_w(1:knon)   + sigx(1:knon)*HTRn_x(1:knon)           
+  dd_HTRn(1:knon)   = HTRn_w(1:knon)   - HTRn_x(1:knon)
+!
+  HQphiT_b(1:knon) = sigw(1:knon)*HQphiT_w(1:knon) + sigx(1:knon)*HQphiT_x(1:knon)          
+  dd_HQphiT(1:knon) = HQphiT_w(1:knon) - HQphiT_x(1:knon)
+!
+  HQphiQ_b(1:knon) = sigw(1:knon)*HQphiQ_w(1:knon) + sigx(1:knon)*HQphiQ_x(1:knon)          
+  dd_HQphiQ(1:knon) = HQphiQ_w - HQphiQ_x(1:knon)
+!
+  HQRn_b(1:knon)   = sigw(1:knon)*HQRn_w(1:knon)   + sigx(1:knon)*HQRn_x(1:knon)             
+  dd_HQRn(1:knon)   = HQRn_w(1:knon)   - HQRn_x(1:knon)
+!
+! Feedback Gains
+! --------------
+ g_T(1:knon) = - (sqrt(tau(1:knon))/Inert(1:knon))  &
+               * (HTphiT_b(1:knon) + L_v(1:knon)*HTphiQ_b(1:knon) + HTRn_b(1:knon)  &
+                 + (dd_HTphiT(1:knon) + L_v(1:knon)*dd_HTphiQ(1:knon) + dd_HTRn(1:knon))  &
+                 * (sigx(1:knon) - sigw(1:knon) - sigw(1:knon)*sigx(1:knon)*dd_HTphiT(1:knon)/HTphiT_b(1:knon)) )
+!
+!!!! DO j = 1,knon
+!!!!  IF (mod(j,20) .eq.0) THEN
+!!!!   print *, '   j     dd_QQ       QQ_b  dd_HQphiQ  dd_HQphiT   dd_HQRn   HQphiQ_b   HQphiT_b     HQRn_b '
+!!!!  ENDIF
+!!!!   print 1789, j, dd_QQ(j), QQ_b(j), dd_HQphiQ(j), dd_HQphiT(j), dd_HQRn(j), HQphiQ_b(j), HQphiT_b(j), HQRn_b(j)
+!!!! 1789 FORMAT( I4, 10(1X,E10.2))
+!!!! ENDDO
+   g_Q(1:knon) = - (dd_QQ(1:knon)/QQ_b(1:knon)) *  &
+                    (sigx(1:knon)-sigw(1:knon)-sigw(1:knon)*sigx(1:knon)*dd_KQs(1:knon)/Kech_Qs(1:knon)) &
+                 - sqrt(tau(1:knon))/(Inert(1:knon)*QQ_b(1:knon)) *  &
+                   ( HQphiT_b(1:knon) + L_v(1:knon)*HQphiQ_b(1:knon) + HQRn_b(1:knon) +  & 
+                      (sigx(1:knon) - sigw(1:knon) - sigw(1:knon)*sigx(1:knon)*dd_KQs(1:knon)/Kech_Qs(1:knon)) *  &
+                       (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
+
+!!   g_Q(1:knon) = - (dd_QQ(1:knon)/QQ_b(1:knon)) *  &
+!!                    (sigx(1:knon)-sigw(1:knon)-sigw(1:knon)*sigx(1:knon)*dd_HQphiQ(1:knon)/HQphiQ_b(1:knon)) &
+!!                 - sqrt(tau(1:knon))/(Inert(1:knon)*QQ_b(1:knon)) *  &
+!!                   ( HQphiT_b(1:knon) + L_v(1:knon)*HQphiQ_b(1:knon) + HQRn_b(1:knon) +  & 
+!!                      (sigx(1:knon) - sigw(1:knon) - sigw(1:knon)*sigx(1:knon)*dd_HQphiQ(1:knon)/HQphiQ_b(1:knon)) *  &
+!!                       (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
+
+!! g_Q(1:knon) = - (sqrt(tau(1:knon))/(Inert(1:knon)*QQ_b(1:knon))) *  &
+!!                 ( HQphiT_b(1:knon) + L_v(1:knon)*HQphiQ_b(1:knon) + HQRn_b(1:knon) )  & 
+!!               - (sigx(1:knon) - sigw(1:knon) - sigw(1:knon)*sigx(1:knon)*dd_HQphiQ(1:knon)/HQphiQ_b(1:knon)) *   &
+!!                 ( dd_QQ(1:knon)/QQ_b(1:knon)   &
+!!                 + (sqrt(tau(1:knon))/(Inert(1:knon)*QQ_b(1:knon)))  &
+!!                 * (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
+
+!  Ts, qs Coupling coefficients                /
+!  ----------------------------
+  Gamma_phiT(1:knon) = (sqrt(tau(1:knon))/(Inert(1:knon)*HTphiT_b(1:knon)))  &
+                     * (dd_HTphiT(1:knon) + L_v(1:knon)*dd_HTphiQ(1:knon) + dd_HTRn(1:knon))
+! 
+  Gamma_phiQ(1:knon) = (1./(Kech_Qs(1:knon)*QQ_b(1:knon))) * &
+                        ( dd_QQ(1:knon)   &
+                         + (sqrt(tau(1:knon))/(Inert(1:knon))) *  &
+                          (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
+
+!!  Gamma_phiQ(1:knon) = (beta(1:knon)/(HQphiQ_b(1:knon)*QQ_b(1:knon))) * &
+!!                        ( dd_QQ(1:knon)   &
+!!                         + (sqrt(tau(1:knon))/(Inert(1:knon))) *  &
+!!                          (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
+
+!!  Gamma_phiQ(1:knon) = (1/(HQphiQ_b(1:knon)*QQ_b(1:knon)))   &
+!!                     * ( dd_QQ(1:knon)   &
+!!                       + (sqrt(tau(1:knon))/(Inert(1:knon)))  &
+!!                       * (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
+!
+!  Insensitive changes
+!  -------------------
+  dTs_ins(1:knon)    = (sqrt(tau(1:knon))/Inert(1:knon))*  &
+                       (dphiT0(1:knon) + L_v(1:knon)*dphiQ0(1:knon) + dRn0(1:knon))
+!
+  dqsatsrf_ins(1:knon) = (beta(1:knon)/QQ_b(1:knon))*dTs_ins(1:knon)
+!
+   IF (prt_level .Ge. 10) THEN
+      print *,'wx_pbl_merge, tau         ', tau
+      print *,'wx_pbl_merge, AcoefT0     ', AcoefT0
+      print *,'wx_pbl_merge, AcoefQ0     ', AcoefQ0
+      print *,'wx_pbl_merge, BcoefT0     ', BcoefT0
+      print *,'wx_pbl_merge, BcoefQ0     ', BcoefQ0
+      print *,'wx_pbl_merge, qsat0_w, qsat0_x ', (qsat0_w(j), qsat0_x(j),j=1,knon)
+      print *,'wx_pbl_merge, dqsatdT0_w, dqsatdT0_x ', (dqsatdT0_w(j), dqsatdT0_x(j),j=1,knon)
+   ENDIF
+!
+!----------------------------------------------------------------------------
+!  
+!------------------------------------------------------------------------------
+!  
+!    Effective coefficients Acoef and Bcoef
+!    --------------------------------------
+   DO j = 1,knon
+     AcoefT(j) = AcoefT0(j) - sigw(j)*sigx(j)*(dd_KTp(j)/Kech_Tp(j))*C_p(j)*   &
+                 (dTs0(j) + (dTs_ins(j)-dTs0(j)-Gamma_phiT(j)*phiT0_b(j))/(1. - g_T(j)))
+     BcoefT(j) = BcoefT0(j) - sigw(j)*sigx(j)*(dd_KTp(j)/Kech_Tp(j))*C_p(j)*Gamma_phiT(j)/(1. - g_T(j))/dtime
+     
+     AcoefQ(j) = AcoefQ0(j) - sigw(j)*sigx(j)*(dd_KQs(j)/Kech_Qs(j))*    &
+                 (dqsatsrf0(j) + (dqsatsrf_ins(j)-(beta(j)/QQ_b(j))*dTs0(j)-Gamma_phiQ(j)*phiQ0_b(j))/(1 - g_Q(j)))/ &
+                 max(beta(j),1.e-4)
+     BcoefQ(j) = BcoefQ0(j) - sigw(j)*sigx(j)*(dd_KQs(j)/Kech_Qs(j))*Gamma_phiQ(j)/(1 - g_Q(j))/ &
+                 (max(beta(j),1.e-4)*dtime)
+!!     AcoefQ(j) = AcoefQ0(j) - sigw(j)*sigx(j)*(dd_KQs(j)/Kech_Qs(j))*    &
+!!                 (dqsatsrf0(j) + (dqsatsrf_ins(j)-(beta(j)/QQ_b(j))*dTs0(j)-Gamma_phiQ(j)*phiQ0_b(j))/(1 - g_Q(j)))/ &
+!!                 beta(j)
+!!     BcoefQ(j) = BcoefQ0(j) - sigw(j)*sigx(j)*(dd_KQs(j)/Kech_Qs(j))*Gamma_phiQ(j)/(1 - g_Q(j))/(beta(j)*dtime)
+   ENDDO ! j = 1,knon
+   
+   IF (prt_level .Ge. 10) THEN
+   print *,'wx_pbl_dts AAAA BcoefQ, BcoefQ0, sigw ', &
+                            BcoefQ, BcoefQ0, sigw
+      print *,'wx_pbl_dts_merge, dTs_ins      ', dTs_ins
+      print *,'wx_pbl_dts_merge, dqs_ins      ', dqsatsrf_ins
+   ENDIF
+
+     RETURN
+
+END SUBROUTINE wx_pbl_dts_merge
+
+SUBROUTINE wx_pbl_split(knon, nsrf, dtime, sigw, beta, iflag_split, &
+                       g_T, g_Q, &
+                       Gamma_phiT, Gamma_phiQ, &
+                       dTs_ins, dqsatsrf_ins, &
+                       phiT, phiQ, phiU, phiV, &
+!!!!                       HTRn_b, dd_HTRn, HTphiT_b, dd_HTphiT, &
+                       phiQ0_b, phiT0_b, &
+                       phiT_x, phiT_w, &
+                       phiQ_x, phiQ_w, &
+                       phiU_x, phiU_w, &
+                       phiV_x, phiV_w, &
+                       philat_x, philat_w, &
+!!!!                       Rn_b, dRn, &
+                       dqsatsrf, &
+                       dTs, delta_qsurf &
+                       )
+!
+
+    USE wx_pbl_var_mod
+
+    USE print_control_mod, ONLY: prt_level,lunout
+    USE indice_sol_mod, ONLY: is_oce
+!
+    INCLUDE "YOMCST.h"
+!
+    INTEGER,                      INTENT(IN)        :: knon    ! number of grid cells
+    INTEGER,                      INTENT(IN)        :: nsrf    ! surface type
+    REAL,                         INTENT(IN)        :: dtime   ! time step size (s)
+    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw ! cold pools fractional area
+    REAL, DIMENSION(knon),        INTENT(IN)        :: beta ! aridity factor
+    INTEGER,                      INTENT(IN)        :: iflag_split
+    REAL, DIMENSION(knon),        INTENT(IN)        :: g_T, g_Q
+    REAL, DIMENSION(knon),        INTENT(IN)        :: Gamma_phiT, Gamma_phiQ
+    REAL, DIMENSION(knon),        INTENT(IN)        :: dTs_ins, dqsatsrf_ins
+    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT, phiQ, phiU, phiV
+    REAL, DIMENSION(knon),        INTENT(IN)        :: phiQ0_b, phiT0_b
+!
+    REAL, DIMENSION(knon),        INTENT(OUT)       :: phiT_x, phiT_w
+    REAL, DIMENSION(knon),        INTENT(OUT)       :: phiQ_x, phiQ_w
+    REAL, DIMENSION(knon),        INTENT(OUT)       :: phiU_x, phiU_w
+    REAL, DIMENSION(knon),        INTENT(OUT)       :: phiV_x, phiV_w
+    REAL, DIMENSION(knon),        INTENT(OUT)       :: philat_x, philat_w
+    REAL, DIMENSION(knon),        INTENT(OUT)       :: dqsatsrf      ! beta delta(qsat(Ts))
+    REAL, DIMENSION(knon),        INTENT(OUT)       :: dTs           ! Temperature difference at surface
+    REAL, DIMENSION(knon),        INTENT(OUT)       :: delta_qsurf
 !
 !! Local variables
     INTEGER                    :: j
-    REAL, DIMENSION(knon)      :: y_delta_flux_t1, y_delta_flux_q1, y_delta_flux_u1, y_delta_flux_v1
-!
-    REAL                       :: DDT, DDQ, DDU, DDV
-    REAL                       :: LambdaTs, LambdaQs, LambdaUs, LambdaVs 
+    REAL, DIMENSION(knon)      :: dphiT, dphiQ, dphiU, dphiV
+    REAL, DIMENSION(knon)      :: q1_x, q1_w
 !
     REAL, DIMENSION(knon)      :: sigx       ! fractional area of (x) region
+
+!----------------------------------------------------------------------------
+!  Equations
+!  ---------
+!!!!!! (1 - g_T) dTs    = dTs_ins    + Gamma_phiT phiT
+!!!!!! (1 - g_Q) dqsatsrf = dqsatsrf_ins + Gamma_phiQ phiQ
+!!!!!! dphiT = (dd_KTp/KTp) phiT + (     dd_AT - C_p dTs)*KxKwTp/KTp
+!!!!!! dphiQ = (dd_KQs/KQs) phiQ + (beta dd_AQ - dqsatsrf )*KxKwQs/KQs
+!!!!!! dphiU = (dd_KUp/KUp) phiU + (     dd_AU          )*KxKwUp/KUp
+!!!!!! dphiV = (dd_KVp/KVp) phiV + (     dd_AV          )*KxKwVp/KVp
+! 
+! (1 - g_T) (dTs-dTs0)    = dTs_ins-dTs0    + Gamma_phiT (phiT-phiT0)
+! (1 - g_Q) dqsatsrf = dqsatsrf_ins + Gamma_phiQ phiQ
+! dphiT = (dd_KTp/KTp) phiT + (     dd_AT - C_p dTs)*KxKwTp/KTp
+! dphiQ = (dd_KQs/KQs) phiQ + (beta dd_AQ - dqsatsrf )*KxKwQs/KQs
+! dphiU = (dd_KUp/KUp) phiU + (     dd_AU          )*KxKwUp/KUp
+! dphiV = (dd_KVp/KVp) phiV + (     dd_AV          )*KxKwVp/KVp
+! 
 !!
-        sigx(:) = 1.-ywake_s(:)
-
-        DO j=1,knon
-!
-       DDT = Kech_Tp(j)
-       DDQ = Kech_Qp(j)
-       DDU = Kech_Up(j)
-       DDV = Kech_Vp(j)
-!
-       LambdaTs =  dd_KTp(j)/DDT
-       LambdaQs =  dd_KQp(j)/DDQ
-       LambdaUs =  dd_KUp(j)/DDU
-       LambdaVs =  dd_KVp(j)/DDV
-!
-         y_delta_flux_t1(j) = y_flux_t1(j)*LambdaTs + dd_AT(j)*KxKwTp(j)/DDT
-         y_delta_flux_q1(j) = y_flux_q1(j)*LambdaQs + dd_AQ(j)*KxKwQp(j)/DDQ
-         y_delta_flux_u1(j) = y_flux_u1(j)*LambdaUs + dd_AU(j)*KxKwUp(j)/DDU
-         y_delta_flux_v1(j) = y_flux_v1(j)*LambdaVs + dd_AV(j)*KxKwVp(j)/DDV
-!
-         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
-!
-!       Delta_tsurf computation
-!!         y_delta_tsurf(j) = (1./RCPD)*(ah(j)*dd_AT(j) + &
-!!                                       ah(j)*y_flux_t1(j)*dd_BT(j)*dtime + &
-!!                                       y_delta_flux_t1(j)*(ah(j)*BBT+bh(j)) )
-!
-           y_delta_tsurf(j) = 0.
-!
-        ENDDO
+        sigx(:) = 1.-sigw(:)
+!
+!      print *,' AAAA wx_pbl_split, C_p(j), qsat0(j), Ts0(j) : ', C_p(:), qsat0(:), Ts0(:)
+!
+   IF (iflag_split .EQ. 2 .AND. nsrf .NE. is_oce) THEN
+!
+!   Delta_tsurf and  Delta_qsurf computation
+!   -----------------------------------------
+      IF (prt_level >=10 ) THEN
+        print *,' wx_pbl_split, dTs_ins, dTs0 , Gamma_phiT, g_T ', dTs_ins, dTs0, Gamma_phiT, g_T
+        print *,' wx_pbl_split, dqsatsrf_ins, Gamma_phiQ, g_q ', dqsatsrf_ins, Gamma_phiQ, g_q
+      ENDIF
+!
+      DO j = 1,knon
+        dTs(j)    = dTs0(j) + (dTs_ins(j) - dTs0(j) + Gamma_phiT(j)*(phiT(j)-phiT0_b(j)) )/(1 - g_T(j)) 
+        dqsatsrf(j) = dqsatsrf0(j) + (dqsatsrf_ins(j) - (beta(j)/QQ_b(j))*dTs0(j) + &
+                       Gamma_phiQ(j)*(phiQ(j)-phiQ0_b(j)) )/(1 - g_Q(j)) 
+      ENDDO ! j = 1,knon
+!
+        IF (prt_level >=10 ) THEN
+          print *,' wx_pbl_split, dqsatsrf0, QQ_b ', dqsatsrf0, QQ_b
+          print *,' wx_pbl_split, phiT0_b, phiT, dTs ', phiT0_b, phiT, dTs
+          print *,' wx_pbl_split, phiQ0_b, phiQ, dqsatsrf ', phiQ0_b, phiQ, dqsatsrf
+        ENDIF
+   ELSE
+        dTs(:) = 0.
+        dqsatsrf(:) = 0.
+   ENDIF ! (iflag_split .EQ. 2 .AND. nsrf .NE. is_oce)
+!
+     DO j = 1,knon
+       dphiT(j) = (phiT(j)*dd_KTp(j) + (        dd_AT(j) - C_p(j)*dTs(j))*KxKwTp(j))/Kech_Tp(j)
+       dphiQ(j) = (phiQ(j)*dd_KQs(j) + (beta(j)*dd_AQ(j) -        dqsatsrf(j))*KxKwQs(j))/Kech_Qs(j)
+       dphiU(j) = (phiU(j)*dd_KUp(j) +          dd_AU(j)                 *KxKwUp(j))/Kech_Up(j)
+       dphiV(j) = (phiV(j)*dd_KVp(j) +          dd_AV(j)                 *KxKwVp(j))/Kech_Vp(j)
+!
+       phiT_x(j)=phiT(j) - sigw(j)*dphiT(j)
+       phiT_w(j)=phiT(j) + sigx(j)*dphiT(j)
+       phiQ_x(j)=phiQ(j) - sigw(j)*dphiQ(j)
+       phiQ_w(j)=phiQ(j) + sigx(j)*dphiQ(j)
+       phiU_x(j)=phiU(j) - sigw(j)*dphiU(j)
+       phiU_w(j)=phiU(j) + sigx(j)*dphiU(j)
+       phiV_x(j)=phiV(j) - sigw(j)*dphiV(j)
+       phiV_w(j)=phiV(j) + sigx(j)*dphiV(j)
+!
+       philat_x(j)=phiQ_x(j)*RLVTT
+       philat_w(j)=phiQ_w(j)*RLVTT
+     ENDDO ! j = 1,knon
+!
+     DO j = 1,knon
+       q1_x(j) = AQ_x(j) + BQ_x(j)*phiQ_x(j)*dtime
+       q1_w(j) = AQ_w(j) + BQ_w(j)*phiQ_w(j)*dtime
+     ENDDO ! j = 1,knon
+     DO j = 1,knon
+       delta_qsurf(j) = (1.-beta(j))*(q1_w(j) - q1_x(j)) + dqsatsrf(j)
+     ENDDO ! j = 1,knon
+!
+!!  Do j = 1,knon
+!!     print *,'XXXsplit : j, q1_x(j), AQ_x(j), BQ_x(j), phiQ_x(j) ', j, q1_x(j), AQ_x(j), BQ_x(j), phiQ_x(j)
+!!     print *,'XXXsplit : j, q1_w(j), AQ_w(j), BQ_w(j), phiQ_w(j) ', j, q1_w(j), AQ_w(j), BQ_w(j), phiQ_w(j)
+!!  ENDDO
+!
+        IF (prt_level >=10 ) THEN
+          print *,' wx_pbl_split, phiT, dphiT, dTs ', phiT, dphiT, dTs
+          print *,' wx_pbl_split, phiQ, dphiQ, dqsatsrf ', phiQ, dphiQ, dqsatsrf
+        ENDIF
+!
+        IF (prt_level >=10 ) THEN
+!!          print *,' wx_pbl_split, verif dqsatsrf = beta dqsatdT0 dTs '
+!!          print *,' wx_pbl_split, dqsatsrf, dqsatdT0*dTs ', dqsatsrf, dqsatdT0*dTs
+        ENDIF
 !
         RETURN
 
-END SUBROUTINE wx_pbl0_split
-
-SUBROUTINE wx_pbl_final
-!
-!****************************************************************************************
-! Deallocate module variables
-!
-!****************************************************************************************    
-!
-    IF (ALLOCATED(Kech_Tp))        DEALLOCATE(Kech_Tp)
-    IF (ALLOCATED(Kech_T_xp))      DEALLOCATE(Kech_T_xp)
-    IF (ALLOCATED(Kech_T_wp))      DEALLOCATE(Kech_T_wp)
-    IF (ALLOCATED(dd_KTp))         DEALLOCATE(dd_KTp)
-    IF (ALLOCATED(KxKwTp))         DEALLOCATE(KxKwTp)
-    IF (ALLOCATED(dd_AT))          DEALLOCATE(dd_AT)
-    IF (ALLOCATED(dd_BT))          DEALLOCATE(dd_BT)
-    IF (ALLOCATED(Kech_Qp))        DEALLOCATE(Kech_Qp)
-    IF (ALLOCATED(Kech_Q_xp))      DEALLOCATE(Kech_Q_xp)
-    IF (ALLOCATED(Kech_Q_wp))      DEALLOCATE(Kech_Q_wp)
-    IF (ALLOCATED(dd_KQp))         DEALLOCATE(dd_KQp)
-    IF (ALLOCATED(KxKwQp))         DEALLOCATE(KxKwQp)
-    IF (ALLOCATED(dd_AQ))          DEALLOCATE(dd_AQ)
-    IF (ALLOCATED(dd_BQ))          DEALLOCATE(dd_BQ)
-    IF (ALLOCATED(Kech_Up))        DEALLOCATE(Kech_Up)
-    IF (ALLOCATED(Kech_U_xp))      DEALLOCATE(Kech_U_xp)
-    IF (ALLOCATED(Kech_U_wp))      DEALLOCATE(Kech_U_wp)
-    IF (ALLOCATED(dd_KUp))         DEALLOCATE(dd_KUp)
-    IF (ALLOCATED(KxKwUp))         DEALLOCATE(KxKwUp)
-    IF (ALLOCATED(dd_AU))          DEALLOCATE(dd_AU)
-    IF (ALLOCATED(dd_BU))          DEALLOCATE(dd_BU)
-    IF (ALLOCATED(Kech_Vp))        DEALLOCATE(Kech_Vp)
-    IF (ALLOCATED(Kech_V_xp))      DEALLOCATE(Kech_V_xp)
-    IF (ALLOCATED(Kech_V_wp))      DEALLOCATE(Kech_V_wp)
-    IF (ALLOCATED(KxKwVp))         DEALLOCATE(KxKwVp)
-    IF (ALLOCATED(dd_KVp))         DEALLOCATE(dd_KVp)
-    IF (ALLOCATED(dd_AV))          DEALLOCATE(dd_AV)
-    IF (ALLOCATED(dd_BV))          DEALLOCATE(dd_BV)
-
-END SUBROUTINE wx_pbl_final
+END SUBROUTINE wx_pbl_split
+
+SUBROUTINE wx_pbl_check( knon, dtime, ypplay, ypaprs, &
+                               sigw, beta, iflag_split,   &
+                               Ts0_b9, dTs09,   &
+                               qs_b9, Ts_b9,  &                         ! yqsurf, Tsurf_new
+                               dTs9, dqsatsrf9,   &
+                               AcoefT_x, AcoefT_w, &
+                               BcoefT_x, BcoefT_w, &
+                               AcoefT0, AcoefQ0, BcoefT0, BcoefQ0, &
+                               AcoefT,  AcoefQ,  BcoefT,  BcoefQ, &
+                               phiT_b9, phiQ_b9,  &
+                               phiT_x9, phiT_w9, &
+                               phiQ_x9, phiQ_w9 &
+                               )
+!
+
+    USE wx_pbl_var_mod
+
+    USE print_control_mod, ONLY: prt_level,lunout
+!
+    INCLUDE "YOMCST.h"
+    INCLUDE "FCTTRE.h"
+    INCLUDE "YOETHF.h"
+!
+    INTEGER,                      INTENT(IN)        :: knon         ! number of grid cells
+    REAL,                         INTENT(IN)        :: dtime        ! time step size (s)
+    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypplay       ! mid-layer pressure (Pa)
+    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypaprs       ! pressure at layer interfaces (pa)
+    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw         ! cold pools fractional area
+    REAL, DIMENSION(knon),        INTENT(IN)        :: beta         ! aridity factor
+    INTEGER,                      INTENT(IN)        :: iflag_split
+    REAL, DIMENSION(knon),        INTENT(IN)        :: Ts0_b9, dTs09
+    REAL, DIMENSION(knon),        INTENT(IN)        :: qs_b9, Ts_b9         ! yqsurf, Tsurf_new
+    REAL, DIMENSION(knon),        INTENT(IN)        :: dTs9, dqsatsrf9
+    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT_x, AcoefT_w
+    REAL, DIMENSION(knon),        INTENT(IN)        :: BcoefT_x, BcoefT_w
+    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT0, AcoefQ0, BcoefT0, BcoefQ0
+!
+    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT, AcoefQ, BcoefT, BcoefQ
+    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT_b9, phiQ_b9
+    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT_x9, phiT_w9
+    REAL, DIMENSION(knon),        INTENT(IN)        :: phiQ_x9, phiQ_w9
+!
+!! Local variables
+    INTEGER                    :: j
+    REAL, DIMENSION(knon)      :: sigx                 ! fractional area of (x) region
+    REAL, DIMENSION(knon)      :: AcoefT_b, AcoefQ_b   ! mean values of AcoefT and AcoefQ 
+    REAL                       :: zzt, zzq, zzqsat
+    REAL                       :: zdelta, zcvm5, zcor, qsat
+    REAL, DIMENSION(knon)      :: qsat_w, qsat_x
+    REAL, DIMENSION(knon)      :: dqsatdT_w, dqsatdT_x
+    REAL, DIMENSION(knon)      :: qsat_bs              ! qsat(Ts_b)
+    REAL, DIMENSION(knon)      :: qsat01, dqsatdT01
+    REAL, DIMENSION(knon)      :: Ts_x, Ts_w, qs_x, qs_w
+    REAL, DIMENSION(knon)      :: T1_x, T1_w, q1_x, q1_w
+    REAL, DIMENSION(knon)      :: Rn_x, Rn_w
+    REAL, DIMENSION(knon)      :: phiQ0_x, phiQ0_w
+    REAL, DIMENSION(knon)      :: Ta, qa
+    REAL, DIMENSION(knon)      :: qsatsrf_w, qsatsrf_x, qsatsrf_b
+    REAL, DIMENSION(knon)      :: qsurf_w, qsurf_x
+    REAL                       :: dphiT, dphiQ
+    REAL                       :: dqsatsrf1
+    REAL                       :: phiT_w1, phiT_w2
+    REAL                       :: phiT_x1, phiT_x2
+    REAL                       :: phiQ_w1, phiQ_w2, phiQ_w3
+    REAL                       :: phiQ_x1, phiQ_x2, phiQ_x3
+    REAL                       :: phiT_b1, phiQ_b1
+    REAL                       :: Kech_Q_sw1, Kech_Q_sx1
+    REAL                       :: evap_pot
+
+!----------------------------------------------------------------------------
+! Equations to be checked:
+! -----------------------
+!  Input : Ts0_b, dTs0, Ts_b, dTs, qsatsrf_b, dqsatsrf, 
+!          phiT_b, phiQ_b, phiT_w, phiT_x, phiQ_w, phiQ_x, 
+!          
+!          AcoefT, AcoefQ, AcoefT_w, AcoefQ_w, AcoefT_x, AcoefQ_x,
+!          BcoefT, BcoefQ, BcoefT_w, BcoefQ_w, BcoefT_x, BcoefQ_x
+! 
+!  C_p T1_w = AcoefT_w + BcoefT_w phiT_w Delta t          C_p T1_x = AcoefT_x + BcoefT_x phiT_x Delta t
+!  q1_w = AQ_w + BQ_w phiQ_w Delta t                      q1_x = AQ_x + BQ_x phiQ_x Delta t
+!  qsatsrf_w = beta qsat(Ts_w)                            qsatsrf_x = beta qsat(Ts_x)
+!  qsurf_w = (1-beta) q1_w + qsatsrf_w                    qsurf_x = (1-beta) q1_x + qsatsrf_x
+!  phiT_w = Kech_h_w C_p ( T1_w - Ts_w)                   phiT_x = Kech_h_x C_p ( T1_x - Ts_x)             
+!  phiT_w = Kech_T_pw ( AcoefT_w - C_p Ts_w)              phiT_x = Kech_T_px ( AcoefT_x - C_p Ts_x)
+!  phiq_w = Kech_h_w ( beta q1_w - qsatsrf_w)             phiq_x = Kech_h_x ( beta q1_x - qsatsrf_x))
+!  phiq_w = Kech_Q_sw (beta AQ_w -qsatsrf_w)              phiq_x = Kech_Q_sx (beta AQ_x -qsatsrf_x)
+!  phiq_w = Kech_h_w (q1_w - qsurf_w)                     phiq_x = Kech_h_x (q1_x - qsurf_x)
+!  phiT_b = sigw phiT_w + sigx phiT_x                     dphiT = phiT_w - phiT_x
+!  phiQ_b = sigw phiQ_w + sigx phiQ_x                     dphiQ = phiQ_w - phiQ_x
+!  Ts_b = sigw Ts_w + sigx Ts_x                           dTs = Ts_w - Ts_x
+!  qsatsrf_b = sigw qsatsrf_w + sigx qsatsrf_x
+!  C_p Ta = AcoefT + BcoefT phiT_b Delta t
+!  qa = AcoefQ + BcoefQ phiQ_b Delta t
+!  phiT_b = Kech_h C_p (Ta - Ts_b)
+!  phiQ_b = beta Kech_h (qa - qsatsrf_b)
+!  dTs = sqrt(tau)/I (dphit + L_v dphiq + dR)
+
+!----------------------------------------------------------------------------
+!
+!!
+        sigx(:) = 1.-sigw(:)
+        AcoefT_b(1:knon) = AcoefT_x(1:knon) + sigw(1:knon)*dd_AT(1:knon)
+        AcoefQ_b(1:knon) = AQ_x(1:knon) + sigw(1:knon)*dd_AQ(1:knon)
+
+! Compute the three qsat and dqsatdTs
+! ---------------------------------------------
+!!   C_p(1:knon) = RCpd
+!!   L_v(1:knon) = RLvtt
+    IF (prt_level >=10 ) THEN
+      print *,' AAAA wx_pbl_check, C_p(j), qsat0(j), Ts0(j) : ', C_p(:), qsat0(:), Ts0(:)
+    ENDIF ! (prt_level >=10 )
+!
+   DO j = 1, knon
+      zdelta = MAX(0.,SIGN(1.,RTT-Ts0_b9(j)))
+      zcvm5 = R5LES*(1.-zdelta) + R5IES*zdelta
+      qsat = R2ES*FOEEW(Ts0_b9(j),zdelta)/ypaprs(j,1)
+      qsat = MIN(0.5,qsat)
+      zcor = 1./(1.-RETV*qsat)
+      qsat01(j) = fqsat*qsat*zcor
+!!      dqsatdT0(j) = FOEDE(Ts0_b(j),zdelta,zcvm5,qsat0(j),zcor)/RLVTT    ! jyg 20210116
+!!      dqsatdT0(j) = (RLvtt*(1.-zdelta)+RLSTT*zdelta)*qsat0(j)/(Rv*Ts0_b(j)*Ts0_b(j))
+      dqsatdT01(j) = fqsat*FOEDE(Ts0_b9(j),zdelta,zcvm5,qsat01(j),zcor)
+   ENDDO
+!
+!--------------------------------------------------------------------------------------------------
+        IF (prt_level >=10 ) THEN
+!
+          DO j = 1, knon
+!
+   print *,'wx_pbl_check: Kech_h, Kech_q ', Kech_h(j), Kech_q(j)
+!
+  Ta(j) = (AcoefT(j) + BcoefT(j)*phiT_b9(j)*dtime)/C_p(j)
+  qa(j) = AcoefQ(j) + BcoefQ(j)*phiQ_b9(j)*dtime
+    print *, 'wx_pbl_check: j, Ta, qa ', Ta(j), qa(j)
+!
+  qsat_bs(j) = qsat01(j) + dqsatdT01(j)*(Ts_b9(j)-Ts0_b9(j))
+!
+   print *,'wx_pbl_check: qsat01, qsat_bs ', j,qsat01(j), qsat_bs(j)
+!
+  Ts_x(j) = Ts_b9(j) - sigw(j)*dTs9(j)
+  Ts_w(j) = Ts_b9(j) + sigx(j)*dTs9(j)
+    print *, 'wx_pbl_check: j, Ts_b9, Ts_w, Ts_x ', j, Ts_b9(j), Ts_w(j), Ts_x(j)
+!
+  qsat_x(j) = qsat0_x(j) + dqsatdT0_x(j)*(Ts_x(j)-Ts0_x(j))
+  qsat_w(j) = qsat0_w(j) + dqsatdT0_w(j)*(Ts_w(j)-Ts0_w(j))
+!
+   print *,'wx_pbl_check: qsat0_w, qsat0_x, qsat_w, qsat_x ', qsat0_w(j), qsat0_x(j), qsat_w(j), qsat_x(j)
+!
+  T1_x(j) = (AcoefT_x(j) + BcoefT_x(j)*phiT_x9(j)*dtime) / C_p(j)
+  T1_w(j) = (AcoefT_w(j) + BcoefT_w(j)*phiT_w9(j)*dtime) / C_p(j)
+    print *, 'wx_pbl_check: j, T1_w, T1_x ', j, T1_w(j), T1_x(j)
+!
+  q1_x(j) = AQ_x(j) + BQ_x(j)*phiQ_x9(j)*dtime
+  q1_w(j) = AQ_w(j) + BQ_w(j)*phiQ_w9(j)*dtime
+    print *, 'wx_pbl_check: j, q1_w, q1_x ', j, q1_w(j), q1_x(j)
+!
+   qsatsrf_x(j) = beta(j)*qsat_x(j)
+   qsatsrf_w(j) = beta(j)*qsat_w(j) 
+   qsatsrf_b(j) = sigw(j)*qsatsrf_w(j) + sigx(j)*qsatsrf_x(j)
+!
+   dqsatsrf1 = qsatsrf_w(j) - qsatsrf_x(j)
+    print *, 'wx_pbl_check: j, qsatsrf_w, qsatsrf_x, dqsatsrf1, dqsatsrf9 ', &
+                          qsatsrf_w(j), qsatsrf_x(j), dqsatsrf1, dqsatsrf9(j)
+!
+   qsurf_x(j) = (1-beta(j))*q1_x(j) + qsatsrf_x(j)
+   qsurf_w(j) = (1-beta(j))*q1_w(j) + qsatsrf_w(j)
+    print *, 'wx_pbl_check: j, qsurf_w, qsurf_x ', j, qsurf_w(j), qsurf_x(j)
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test qsat01 = qsat0    et   dqsatdT01 = dqsatdT0
+!------------------------------------------------------------------------------------------------------
+   print *, 'wx_pbl_check: j, qsat01(j), qsat0(j) ', j, qsat01(j), qsat0(j)
+   print *, 'wx_pbl_check: j, dqsatdT01(j), dqsatdT0(j) ', j, dqsatdT01(j), dqsatdT0(j)
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test Kexh_Q_sw = Kech_q_w/(1.-beta*Kech_q_w*BcoefQ)   Kexh_Q_sx = Kech_q_x/(1.-beta*Kech_q_x*BcoefQ)
+!------------------------------------------------------------------------------------------------------
+  Kech_Q_sx1 = Kech_q_x(j)/(1.-beta(j)*Kech_q_x(j)*BQ_x(j)*dtime)
+  Kech_Q_sw1 = Kech_q_w(j)/(1.-beta(j)*Kech_q_w(j)*BQ_w(j)*dtime)
+    print *, 'wx_pbl_check: j, Kech_Q_sx1, Kech_Q_sx(j)', j, Kech_Q_sx1, Kech_Q_sx(j)
+    print *, 'wx_pbl_check: j, Kech_Q_sw1, Kech_Q_sw(j)', j, Kech_Q_sw1, Kech_Q_sw(j)
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test phiT_w = Kech_h_w*C_p(j)*(T1_w(j)-Ts_w(j))        phiT_x = Kech_h_x*C_p(j)*(T1_x(j)-Ts_x(j))
+!-----------------------------------------------------
+    phiT_x1 = Kech_h_x(j)*C_p(j)*(T1_x(j)-Ts_x(j))
+    phiT_w1 = Kech_h_w(j)*C_p(j)*(T1_w(j)-Ts_w(j))
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test phiT_w = Kech_T_pw*(AcoefT_w(j)-C_p(j)*Ts_w(j))   phiT_x = Kech_T_px*(AcoefT_x(j)-C_p(j)*Ts_x(j))
+!-----------------------------------------------------
+    phiT_x2 = Kech_T_px(j)*(AcoefT_x(j)-C_p(j)*Ts_x(j))
+    phiT_w2 = Kech_T_pw(j)*(AcoefT_w(j)-C_p(j)*Ts_w(j))
+    print *, 'wx_pbl_check: j, phiT_w1, phiT_w2, phiT_w9 ', j, phiT_w1, phiT_w2, phiT_w9(j)
+    print *, 'wx_pbl_check: j, phiT_x1, phiT_x2, phiT_x9 ', j, phiT_x1, phiT_x2, phiT_x9(j)
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test phiq_w = Kech_q_w ( beta q1_w - qsatsrf_w)    phiq_x = Kech_q_x ( beta q1_x - qsatsrf_x))
+!--------------------------------------------------------------
+    phiq_x1 = Kech_q_x(j)*( beta(j)*q1_x(j) - qsatsrf_x(j))
+    phiq_w1 = Kech_q_w(j)*( beta(j)*q1_w(j) - qsatsrf_w(j))
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test  phiq_w = Kech_Q_sw (beta AQ_w -qsatsrf_w)     phiq_x = Kech_Q_sx (beta AQ_x -qsatsrf_x)
+!--------------------------------------------------------------
+    phiq_x2 = Kech_Q_sx(j)*(beta(j)*AQ_x(j) -qsatsrf_x(j))
+    phiq_w2 = Kech_Q_sw(j)*(beta(j)*AQ_w(j) -qsatsrf_w(j))
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test phiq_w = Kech_q_w ( q1_w - qsurf_w)    phiq_x = Kech_q_x ( q1_x - qsurf_x))
+!--------------------------------------------------------------
+    phiq_x3 = Kech_q_x(j)*( q1_x(j) - qsurf_x(j))
+    phiq_w3 = Kech_q_w(j)*( q1_w(j) - qsurf_w(j))
+    print *, 'wx_pbl_check: j, phiQ_w1, phiQ_w2, phiQ_w3, phiQ_w9 ', j, phiQ_w1, phiQ_w2, phiQ_w3, phiQ_w9(j)
+    print *, 'wx_pbl_check: j, phiQ_x1, phiQ_x2, phiQ_x3, phiQ_x9 ', j, phiQ_x1, phiQ_x2, phiQ_x3, phiQ_x9(j)
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test phiT_b = Kech_h C_p (Ta - Ts_b)
+!--------------------------------------------------------------
+   phiT_b1 = Kech_h(j)*C_p(j)*(Ta(j) - Ts_b9(j))
+   print *, 'wx_pbl_check: j, phiT_b1, PhiT_b9 ', j, phiT_b1, PhiT_b9(j)
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test phiQ_b = beta Kech_q (qa - qsat_bs)
+!--------------------------------------------------------------
+   evap_pot = Kech_q(j)*(qa(j) - qsat_bs(j))
+   phiQ_b1 = beta(j)*Kech_q(j)*(qa(j) - qsat_bs(j))
+   print *, 'wx_pbl_check: j, beta, evap_pot, phiQ_b1, PhiQ_b9 ', j, beta(j), evap_pot, phiQ_b1, PhiQ_b9(j)
+!
+!
+          ENDDO  ! j = 1, knon
+          
+        ENDIF   ! (prt_level >=10 )
+!--------------------------------------------------------------------------------------------------
+
+        RETURN
+
+END SUBROUTINE wx_pbl_check
+
+SUBROUTINE wx_pbl_dts_check( knon, dtime, ypplay, ypaprs, &
+                               sigw, beta, iflag_split,   &
+                               Ts0_b9, dTs09,   &
+                               qs_b9, Ts_b9,  &                         ! yqsurf, Tsurf_new
+                               dqsatsrf9, dTs9, delta_qsurf9,   &
+                               AcoefT_x, AcoefT_w, &
+                               BcoefT_x, BcoefT_w, &
+                               AcoefT0, AcoefQ0, BcoefT0, BcoefQ0, &
+                               AcoefT,  AcoefQ,  BcoefT,  BcoefQ, &
+                               HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
+                               phiT0_b9, dphiT09, phiQ0_b9, dphiQ09, Rn0_b9, dRn09, &
+                               g_T, g_Q, &
+                               Gamma_phiT, Gamma_phiQ, &
+                               dTs_ins, dqsatsrf_ins, &
+                               phiT_b9, phiQ_b9,  &
+                               phiT_x9, phiT_w9, &
+                               phiQ_x9, phiQ_w9  &
+                               )
+!
+
+    USE wx_pbl_var_mod
+
+    USE print_control_mod, ONLY: prt_level,lunout
+!
+    INCLUDE "YOMCST.h"
+    INCLUDE "FCTTRE.h"
+    INCLUDE "YOETHF.h"
+!
+    INTEGER,                      INTENT(IN)        :: knon         ! number of grid cells
+    REAL,                         INTENT(IN)        :: dtime        ! time step size (s)
+    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypplay       ! mid-layer pressure (Pa)
+    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypaprs       ! pressure at layer interfaces (pa)
+    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw         ! cold pools fractional area
+    REAL, DIMENSION(knon),        INTENT(IN)        :: beta         ! aridity factor
+    INTEGER,                      INTENT(IN)        :: iflag_split
+    REAL, DIMENSION(knon),        INTENT(IN)        :: Ts0_b9, dTs09
+    REAL, DIMENSION(knon),        INTENT(IN)        :: qs_b9, Ts_b9         ! yqsurf, Tsurf_new
+    REAL, DIMENSION(knon),        INTENT(IN)        :: dTs9, dqsatsrf9
+    REAL, DIMENSION(knon),        INTENT(IN)        :: delta_qsurf9
+    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT_x, AcoefT_w
+    REAL, DIMENSION(knon),        INTENT(IN)        :: BcoefT_x, BcoefT_w
+    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT0, AcoefQ0, BcoefT0, BcoefQ0
+!
+    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT, AcoefQ, BcoefT, BcoefQ
+    REAL, DIMENSION(knon),        INTENT(IN)        :: HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn
+    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT0_b9, dphiT09, phiQ0_b9, dphiQ09, Rn0_b9, dRn09
+    REAL, DIMENSION(knon),        INTENT(IN)        :: g_T, g_Q
+    REAL, DIMENSION(knon),        INTENT(IN)        :: Gamma_phiT, Gamma_phiQ
+    REAL, DIMENSION(knon),        INTENT(IN)        :: dTs_ins, dqsatsrf_ins
+    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT_b9, phiQ_b9
+    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT_x9, phiT_w9
+    REAL, DIMENSION(knon),        INTENT(IN)        :: phiQ_x9, phiQ_w9
+!
+!! Local variables
+    INTEGER                    :: j
+    REAL, DIMENSION(knon)      :: sigx       ! fractional area of (x) region
+    REAL, DIMENSION(knon)      :: AcoefT_b, AcoefQ_b   ! mean values of AcoefT and AcoefQ 
+    REAL                       :: zzt, zzq, zzqsat
+    REAL                       :: zdelta, zcvm5, zcor, qsat
+    REAL, DIMENSION(knon)      :: qsat_w, qsat_x
+    REAL, DIMENSION(knon)      :: Ts_x, Ts_w, qs_x, qs_w
+    REAL, DIMENSION(knon)      :: T1_x, T1_w, q1_x, q1_w
+    REAL, DIMENSION(knon)      :: Rn_x, Rn_w
+    REAL, DIMENSION(knon)      :: Rn_b, dRn
+    REAL, DIMENSION(knon)      :: phiQ0_x, phiQ0_w
+    REAL, DIMENSION(knon)      :: Ta, qa
+    REAL, DIMENSION(knon)      :: err_phiT_w, err_phiT_x
+    REAL, DIMENSION(knon)      :: err_phiq_w, err_phiq_x
+    REAL, DIMENSION(knon)      :: err_phiT_b
+    REAL, DIMENSION(knon)      :: err_phiQ_b
+    REAL, DIMENSION(knon)      :: err2_phiT_b
+    REAL                       :: T1A_x, T1A_w, q1A_x, q1A_w
+    REAL                       :: qsatsrf_w, qsatsrf_x, qsatsrfb, qsbA
+    REAL                       :: dphiT, dphiQ
+    REAL                       :: dphiT_H, dphiQ_H
+    REAL                       :: phiQ_pot
+    REAL                       :: phiQ_w_m_phiQ0_w
+    REAL                       :: phiQ_x_m_phiQ0_x
+    REAL                       :: dphiQ_m_dphiQ0
+    REAL                       :: dphiT_m_dphiT0
+    REAL                       :: dRN_m_dRn0
+    REAL                       :: phiTb_m_phiT0b
+
+!----------------------------------------------------------------------------
+! Equations to be checked:
+! -----------------------
+!  Input : Ts0_b, dTs0, Ts_b, dTs, qsatsrf_b, dqsatsrf, 
+!          phiT_b, phiQ_b, phiT_w, phiT_x, phiQ_w, phiQ_x, 
+!          
+!          AcoefT, AcoefQ, AcoefT_w, AcoefQ_w, AcoefT_x, AcoefQ_x,
+!          BcoefT, BcoefQ, BcoefT_w, BcoefQ_w, BcoefT_x, BcoefQ_x
+! 
+!  Ts_w = Ts_b + sigx dTs                                 Ts_x = Ts_b - sigw dTs
+!  T1_w = AcoefT_w + BcoefT_w phiT_w Delta t              T1_x = AcoefT_x + BcoefT_x phiT_x Delta t
+!  q1_w = AcoefQ_w + BcoefQ_w phiQ_w Delta t              q1_x = AcoefQ_x + BcoefQ_x phiQ_x Delta t
+!  phiT_w = Kech_h_w ( T1_w - Ts_w)                       phiT_x = Kech_h_x ( T1_x - Ts_x)             
+!  phiq_w = beta Kech_h_w ( q1_w - qsat(Ts_w))            phiq_x = beta Kech_h_x ( q1_x - qsat(Ts_x))
+!  phiT_b = sigw phiT_w + sigx phiT_x                     dphiT = phiT_w - phiT_x
+!  phiQ_b = sigw phiQ_w + sigx phiQ_x                     dphiQ = phiQ_w - phiQ_x
+!  Ts_b = sigw Ts_w + sigx Ts_x                           dTs = Ts_w - Ts_x
+!  Ta = AcoefT + BcoefT phiT_b Delta t
+!  qa = AcoefQ + BcoefQ phiQ_b Delta t
+!  phiT_b = Kech_h (Ta - Ts_b)
+!  phiQ_b = beta Kech_h (qa - qsat(Ts_b))
+!  dTs = sqrt(tau)/I (dphit + L_v dphiq + dR)
+
+!----------------------------------------------------------------------------
+!
+!!
+        sigx(:) = 1.-sigw(:)
+        AcoefT_b(1:knon) = AcoefT_x(1:knon) + sigw(1:knon)*dd_AT(1:knon)
+        AcoefQ_b(1:knon) = AQ_x(1:knon) + sigw(1:knon)*dd_AQ(1:knon)
+
+   IF (prt_level >=10 ) THEN
+    print *,'->wx_pbl_dts_check, HTphiT_b, HTphiQ_b, HTRn_b ', &
+                             HTphiT_b, HTphiQ_b, HTRn_b 
+    print *,'->wx_pbl_dts_check, dd_HTphiT, dd_HTphiQ, dd_HTRn ', &
+                             dd_HTphiT, dd_HTphiQ, dd_HTRn 
+   ENDIF ! (prt_level >=10 )
+!
+! Compute the three qsat and dqsatdTs
+! ---------------------------------------------
+!!      print *,' AAAA wx_pbl_dts_check, C_p(j), qsat0(j), Ts0(j) : ',  &
+!!                                      (C_p(j), qsat0(j), Ts0(j), j = 1,knon)
+!
+!
+!--------------------------------------------------------------------------------------------------
+        IF (prt_level >=10 ) THEN
+!
+          DO j = 1, knon
+  Ts_x(j) = Ts_b9(j) - sigw(j)*dTs9(j)
+  Ts_w(j) = Ts_b9(j) + sigx(j)*dTs9(j)
+    print *, 'wx_pbl_dts_check: j, Ts_b9, Ts_w, Ts_x ', j, Ts_b9(j), Ts_w(j), Ts_x(j)
+!
+  qsat_x(j) = qsat0_x(j) + dqsatdT0_x(j)*(Ts_x(j)-Ts0_x(j))
+  qsat_w(j) = qsat0_w(j) + dqsatdT0_w(j)*(Ts_w(j)-Ts0_w(j))
+!
+  T1_x(j) = (AcoefT_x(j) + BcoefT_x(j)*phiT_x9(j)*dtime) / C_p(j)
+  T1_w(j) = (AcoefT_w(j) + BcoefT_w(j)*phiT_w9(j)*dtime) / C_p(j)
+    print *, 'wx_pbl_dts_check: j, T1_w, T1_x ', j, T1_w(j), T1_x(j)
+!
+  q1_x(j) = AQ_x(j) + BQ_x(j)*phiQ_x9(j)*dtime
+  q1_w(j) = AQ_w(j) + BQ_w(j)*phiQ_w9(j)*dtime
+    print *, 'wx_pbl_dts_check: j, q1_w, q1_x ', j, q1_w(j), q1_x(j)
+!
+    Rn_x(j) = eps_1*Rsigma*T1_x(j)**4 - Rsigma*Ts_x(j)**4
+    Rn_w(j) = eps_1*Rsigma*T1_w(j)**4 - Rsigma*Ts_w(j)**4
+    Rn_b(j) = sigw(j)*Rn_w(j) + sigx(j)*Rn_x(j)
+    dRn(j) = dRn09(j) - ( HTRn_b(j) &
+                        +(sigx(j)-sigw(j))*dd_HTRn(j) &
+                        -sigw(j)*sigx(j)*dd_HTRn(j)*dd_HTphiT(j)/HTphiT_b(j) &
+                       )*(dTs9(j)-dTs09(j)) &
+                    + dd_HTRn(j)/HTphiT_b(j)*(phiT_b9(j)-phiT0_b9(j))
+!
+          print *,'wx_pbl_dts_check, dphiT, L_v*dphiQ, dRn, dTs ', &
+              phiT_w9(j)-phiT_x9(j), L_v(j)*(phiQ_w9(j)-phiQ_x9(j)), dRn(j), dTs9(j)
+!
+  phiQ0_x(j) = PhiQ0_b9(j) - sigw(j)*dphiQ09(j)
+  phiQ0_w(j) = PhiQ0_b9(j) + sigx(j)*dphiQ09(j)
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test phiQ_w-phiQ0_w = -beta*Kech_Q_sw*dqsatdT_w*(Ts_w-Ts0_w)
+!--------------------------------------------------------------
+  print *,'wx_pbl_dts_check: beta(j), Kech_Q_sw(j), dqsatdT0_w(j), Ts_w(j), Ts0_w(j) ', &
+                         beta(j), Kech_Q_sw(j), dqsatdT0_w(j), Ts_w(j), Ts0_w(j)
+  phiQ_w_m_phiQ0_w = -beta(j)*Kech_Q_sw(j)*dqsatdT0_w(j)*(Ts_w(j)-Ts0_w(j))
+    print *,'wx_pbl_dts_check: j, phiQ_w9-phiQ0_w, phiQ_w_m_phiQ0_w ', &
+                           j, phiQ_w9(j)-phiQ0_w(j), phiQ_w_m_phiQ0_w
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test phiQ_x-phiQ0_x = -beta*Kech_Q_sx*dqsatdT_x*(Ts_x-Ts0_x)
+!--------------------------------------------------------------
+  phiQ_x_m_phiQ0_x = -beta(j)*Kech_Q_sx(j)*dqsatdT0_x(j)*(Ts_x(j)-Ts0_x(j))
+    print *,'wx_pbl_dts_check: j, phiQ_x9-phiQ0_x, phiQ_x_m_phiQ0_x ', &
+                           j, phiQ_x9(j)-phiQ0_x(j), phiQ_x_m_phiQ0_x
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test dphiT-dphiT0 = -(HTphiT_b+(sigx-sigw)*dd_HTphiT)*(dTs-dTs0) - dd_HTphiT*(Ts_b-Ts0_b)
+!-------------------------------------------------------------------------------------------
+ dphiT = phiT_w9(j) - phiT_x9(j)
+ dphiT_m_dphiT0 = -(HTphiT_b(j)+(sigx(j)-sigw(j))*dd_HTphiT(j))*(dTs9(j)-dTs09(j)) &
+                  - dd_HTphiT(j)*(Ts_b9(j)-Ts0_b9(j))
+ print *,'wx_pbl_dts_check: j, dphiT-dphiT09, dphiT_m_dphiT0 ',j, dphiT-dphiT09(j), dphiT_m_dphiT0
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test dphiQ-dphiQ0 = -(HTphiQ_b+(sigx-sigw)*dd_HTphiQ)*(dTs-dTs0) - dd_HTphiQ*(Ts_b-Ts0_b)
+!-------------------------------------------------------------------------------------------
+ dphiQ = phiQ_w9(j) - phiQ_x9(j)
+ dphiQ_m_dphiQ0 = -(HTphiQ_b(j)+(sigx(j)-sigw(j))*dd_HTphiQ(j))*(dTs9(j)-dTs09(j)) &
+                  - dd_HTphiQ(j)*(Ts_b9(j)-Ts0_b9(j))
+ print *,'wx_pbl_dts_check: j, dphiQ-dphiQ09, dphiQ_m_dphiQ0 ',j, dphiQ-dphiQ09(j), dphiQ_m_dphiQ0
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test dRn-dRn0 = -(HTRn_b+(sigx-sigw)*dd_HTRn)*(dTs-dTs0) - dd_HTRn*(Ts_b-Ts0_b)
+!-------------------------------------------------------------------------------------------
+ dRn_m_dRn0 = -(HTRn_b(j)+(sigx(j)-sigw(j))*dd_HTRn(j))*(dTs9(j)-dTs09(j)) &
+                  - dd_HTRn(j)*(Ts_b9(j)-Ts0_b9(j))
+ print *,'wx_pbl_dts_check: j, dRn-dRn09, dRn_m_dRn0 ',j, dRn-dRn09(j), dRn_m_dRn0
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test phiT_b-phiT0_b = -sigx*sigw*dd_HTphiT*(dTs-dTs0) - HTphiT_b*(Ts_b-Ts0_b)
+!-------------------------------------------------------------------------------
+   phiTb_m_phiT0b = -sigx(j)*sigw(j)*dd_HTphiT(j)*(dTs9(j)-dTs09(j)) - HTphiT_b(j)*(Ts_b9(j)-Ts0_b9(j))
+   print *,'wx_pbl_dts_check: j, phiT_b9-phiT0_b9, phiTb_m_phiT0b ',j ,phiT_b9(j)-phiT0_b9(j), phiTb_m_phiT0b
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test phiT_w, phiT_x, dphiT from HTphiT
+!------------------------------------------
+!  phiT_w = Kech_h_w C_p ( T1_w - Ts_w)                   phiT_x = Kech_h_x C_p ( T1_x - Ts_x)             
+  err_phiT_x(j) = Kech_h_x(j)*C_p(j)*(T1_x(j) - Ts_x(j)) - phiT_x9(j)
+  err_phiT_w(j) = Kech_h_w(j)*C_p(j)*(T1_w(j) - Ts_w(j)) - phiT_w9(j)
+    print *, 'wx_pbl_dts_check: j, phiT_w9, phiT_x9, err_phiT_w, err_phiT_x ',   &
+                            j, phiT_w9(j), phiT_x9(j), err_phiT_w(j), err_phiT_x(j)
+  dphiT = phiT_w9(j) - phiT_x9(j)
+  dphiT_H = dphiT09(j) - ( HTphiT_b(j) &
+                            +(sigx(j)-sigw(j))*dd_HTphiT(j) &
+                            -sigw(j)*sigx(j)*dd_HTphiT(j)*dd_HTphiT(j)/HTphiT_b(j) &
+                           )*(dTs9(j)-dTs09(j)) &
+                         + dd_HTphiT(j)/HTphiT_b(j)*(phiT_b9(j)-phiT0_b9(j))
+    print *,'wx_pbl_dts_check: j, dphiT, dphiT_H ', j, dphiT, dphiT_H
+
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Test phiq_w, phiq_x, dphiq from HTphiq
+!------------------------------------------
+!
+!  phiq_w = beta Kech_q_w ( q1_w - qsat(Ts_w))            phiq_x = beta Kech_q_x ( q1_x - qsat(Ts_x))
+  err_phiq_x(j) = beta(j)*Kech_q_x(j)*( q1_x(j) - qsat_x(j)) - phiq_x9(j)
+  err_phiq_w(j) = beta(j)*Kech_q_w(j)*( q1_w(j) - qsat_w(j)) - phiq_w9(j)
+  dphiQ = phiQ_w9(j) - phiQ_x9(j)
+  dphiQ_H = dphiQ09(j) - ( HTphiQ_b(j) &
+                            +(sigx(j)-sigw(j))*dd_HTphiQ(j) &
+                            -sigw(j)*sigx(j)*dd_HTphiQ(j)*dd_HTphiT(j)/HTphiT_b(j) &
+                           )*(dTs9(j)-dTs09(j)) &
+                         + dd_HTphiQ(j)/HTphiT_b(j)*(phiT_b9(j)-phiT0_b9(j))
+    print *,'wx_pbl_dts_check: j, dphiQ, dphiQ_H ', j, dphiQ, dphiQ_H
+!
+!  phiT_b = sigw phiT_w + sigx phiT_x                     dphiT = phiT_w - phiT_x
+  err_phiT_b(j) = sigw(j)*phiT_w9(j) + sigx(j)*phiT_x9(j) - phiT_b9(j)
+!
+!  phiQ_b = sigw phiQ_w + sigx phiQ_x                     dphiQ = phiQ_w - phiQ_x
+  err_phiQ_b(j) = sigw(j)*phiQ_w9(j) + sigx(j)*phiQ_x9(j) - phiQ_b9(j)
+!
+!  Ta = AcoefT + BcoefT phiT_b Delta t
+!  phiT_b = Kech_h C_p (Ta - Ts_b)
+  Ta(j) = (AcoefT(j) + BcoefT(j)*phiT_b9(j)*dtime) / C_p(j)
+  err2_phiT_b(j) = Kech_h(j)*C_p(j)*(Ta(j) - Ts_b9(j)) - phiT_b9(j)
+    print *, 'wx_pbl_dts_check: j, Ta, phiT_b9, err2_phiT_b ',   &
+                            j, Ta(j), phiT_b9(j), err2_phiT_b(j)
+!
+          ENDDO  ! j = 1, knon
+          
+        ENDIF   ! (prt_level >=10 )
+!--------------------------------------------------------------------------------------------------
+        RETURN
+
+END SUBROUTINE wx_pbl_dts_check
+
+SUBROUTINE wx_evappot(knon, q1, Ts, evap_pot)
+
+    USE wx_pbl_var_mod
+
+    INTEGER,                      INTENT(IN)        :: knon     ! number of grid cells
+    REAL, DIMENSION(knon),        INTENT(IN)        :: q1       ! specific humidity in layer 1
+    REAL, DIMENSION(knon),        INTENT(IN)        :: Ts       ! surface temperature
+!
+    REAL, DIMENSION(knon),        INTENT(OUT)       :: evap_pot ! potential evaporation
+!
+    INTEGER                   :: j
+    REAL                      :: qsat_bs
+!
+ DO j = 1,knon
+   evap_pot(j) = Kech_q(j)*(qsat0(j)+dqsatdT0(j)*(Ts(j)-Ts0(j))-q1(j))
+!
+  qsat_bs = qsat0(j)+dqsatdT0(j)*(Ts(j)-Ts0(j))
+!!  print *,'wx_evappot : Kech_q, qsat_bs, qa, evap_pot ', Kech_q(j), qsat_bs, q1(j), evap_pot(j)
+ ENDDO
+!
+ RETURN
+END SUBROUTINE wx_evappot
 
 END MODULE wx_pbl_mod
Index: /LMDZ6/trunk/libf/phylmd/wx_pbl_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/wx_pbl_var_mod.F90	(revision 3888)
+++ /LMDZ6/trunk/libf/phylmd/wx_pbl_var_mod.F90	(revision 3888)
@@ -0,0 +1,648 @@
+MODULE wx_pbl_var_mod
+!
+! Split Planetary Boundary Layer variables
+!
+! This module manages the variables necessary for the splitting of the boundary layer 
+! 
+!
+  USE dimphy
+
+  IMPLICIT NONE
+
+  REAL, PROTECTED, SAVE                             :: eps_1, fqsat, smallestreal
+!$OMP THREADPRIVATE(eps_1, fqsat, smallestreal)
+!
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: C_p, L_v
+!$OMP THREADPRIVATE(C_p, L_v)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: Ts0, dTs0
+!$OMP THREADPRIVATE(Ts0, dTs0)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: Ts0_x, Ts0_w
+!$OMP THREADPRIVATE(Ts0_x, Ts0_w)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: qsat0, dqsatdT0
+!$OMP THREADPRIVATE(qsat, dqsat_dT)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: qsat0_x, dqsatdT0_x
+!$OMP THREADPRIVATE(qsat0_x, dqsatdT0_x)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: qsat0_w, dqsatdT0_w
+!$OMP THREADPRIVATE(qsat0_w, dqsatdT0_w)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: QQ_b, dd_QQ
+!$OMP THREADPRIVATE(QQ_b, dd_QQ)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: QQ_x, QQ_w
+!$OMP THREADPRIVATE(QQ_x, QQ_w)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: qsatsrf0_x, qsatsrf0_w
+!$OMP THREADPRIVATE(qsatsrf0_x, qsatsrf0_w)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: dqsatsrf0
+!$OMP THREADPRIVATE(dqsatsrf0)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: q1_0b
+!$OMP THREADPRIVATE(q1_0b)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: dd_Cdragh, dd_Cdragm, dd_Cdragq 
+!$OMP THREADPRIVATE(dd_Cdragh, dd_Cdragm, dd_Cdragq )
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: Kech_h, Kech_h_x, Kech_h_w   ! Energy exchange coefficients
+!$OMP THREADPRIVATE(Kech_h, Kech_h_x, Kech_h_w)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: Kech_q, Kech_q_x, Kech_q_w   ! Moisture exchange coefficients
+!$OMP THREADPRIVATE(Kech_q, Kech_q_x, Kech_q_w)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: Kech_m, Kech_m_x, Kech_m_w   ! Momentum exchange coefficients
+!$OMP THREADPRIVATE(Kech_m, Kech_m_x, Kech_m_w)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: Kech_Tp, Kech_T_px, Kech_T_pw
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: dd_KTp, KxKwTp
+!$OMP THREADPRIVATE(Kech_Tp, Kech_T_px, Kech_T_pw, dd_KTp, KxKwTp)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: dd_AT, dd_BT
+!$OMP THREADPRIVATE(dd_AT, dd_BT)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: Kech_Qp, Kech_Q_px, Kech_Q_pw
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: dd_KQp, KxKwQp
+!$OMP THREADPRIVATE(Kech_Qp, Kech_Q_px, Kech_Q_pw, dd_KQp, KxKwQp)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: Kech_Qs, Kech_Q_sx, Kech_Q_sw
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: dd_KQs, KxKwQs
+!$OMP THREADPRIVATE(Kech_Qs, Kech_Q_sx, Kech_Q_sw, dd_KQs, KxKwQs)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: dd_AQ, dd_BQ
+!$OMP THREADPRIVATE(dd_AQ, dd_BQ)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: AQ_x, AQ_w, BQ_x, BQ_w
+!$OMP THREADPRIVATE(AQ_x, AQ_w, BQ_x, BQ_w)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: Kech_Up, Kech_U_px, Kech_U_pw
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: dd_KUp, KxKwUp
+!$OMP THREADPRIVATE(Kech_Up, Kech_U_px, Kech_U_pw, dd_KUp, KxKwUp)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: dd_AU, dd_BU
+!$OMP THREADPRIVATE(dd_AU, dd_BU)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: Kech_Vp, Kech_V_px, Kech_V_pw
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: dd_KVp, KxKwVp
+!$OMP THREADPRIVATE(Kech_Vp, Kech_V_px, Kech_V_pw, dd_KVp, KxKwVp)
+  REAL, ALLOCATABLE, PROTECTED, DIMENSION(:), SAVE  :: dd_AV, dd_BV
+!$OMP THREADPRIVATE(dd_AV, dd_BV)
+
+CONTAINS
+!
+!****************************************************************************************
+!
+SUBROUTINE wx_pbl_init
+
+! Local variables
+!****************************************************************************************
+    INTEGER                       :: ierr
+ 
+
+!****************************************************************************************
+! Allocate module variables
+!
+!****************************************************************************************    
+
+    ierr = 0
+
+    ALLOCATE(C_p(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(L_v(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Ts0(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dTs0(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Ts0_x(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Ts0_w(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(qsat0(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dqsatdT0(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(qsat0_x(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dqsatdT0_x(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(qsat0_w(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dqsatdT0_w(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(q1_0b(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(QQ_b(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_QQ(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(QQ_x(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(QQ_w(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(qsatsrf0_x(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(qsatsrf0_w(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dqsatsrf0(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_Cdragh(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_Cdragm(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_Cdragq(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_h(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_h_x(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_h_w(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_q(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_q_x(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_q_w(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_m(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_m_x(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_m_w(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_Tp(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_T_px(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_T_pw(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_KTp(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(KxKwTp(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_AT(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_BT(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+!----------------------------------------------------------------------------
+    ALLOCATE(Kech_Qp(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_Q_px(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_Q_pw(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_KQp(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(KxKwQp(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_Qs(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_Q_sx(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_Q_sw(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_KQs(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(KxKwQs(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+!!!!!!!!!!
+    ALLOCATE(AQ_x(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(AQ_w(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(BQ_x(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(BQ_w(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_AQ(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_BQ(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+!----------------------------------------------------------------------------
+    ALLOCATE(Kech_Up(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_U_px(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_U_pw(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_KUp(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(KxKwUp(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_AU(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_BU(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+!----------------------------------------------------------------------------
+    ALLOCATE(Kech_Vp(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_V_px(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(Kech_V_pw(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_KVp(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(KxKwVp(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_AV(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+    ALLOCATE(dd_BV(klon), stat=ierr)
+    IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
+
+!----------------------------------------------------------------------------
+
+END SUBROUTINE wx_pbl_init
+
+SUBROUTINE wx_pbl_prelim_0(knon, nsrf, dtime, ypplay, ypaprs, sigw,   &
+                                 yt_s, ydeltat_s, ygustiness, &
+                                 yt_x, yt_w, yq_x, yq_w, &
+                                 yu_x, yu_w, yv_x, yv_w, &
+                                 ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
+                                 ycdragm_x, ycdragm_w, &
+                                 AcoefT_x, AcoefT_w, AcoefQ_x, AcoefQ_w, &
+                                 AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
+                                 BcoefT_x, BcoefT_w, BcoefQ_x, BcoefQ_w, &
+                                 BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w  &
+                                 )
+!
+    USE print_control_mod, ONLY: prt_level,lunout
+    USE indice_sol_mod, ONLY: is_oce
+!
+    INCLUDE "YOMCST.h"
+    INCLUDE "FCTTRE.h"
+    INCLUDE "YOETHF.h"
+    INCLUDE "clesphys.h"
+!
+    INTEGER,                      INTENT(IN)        :: knon    ! number of grid cells
+    INTEGER,                      INTENT(IN)        :: nsrf    ! surface type
+    REAL,                         INTENT(IN)        :: dtime   ! time step size (s)
+    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypplay  ! mid-layer pressure (Pa)
+    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypaprs  ! pressure at layer interfaces (pa)
+    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw ! cold pools fractional area
+    REAL, DIMENSION(knon),        INTENT(IN)        :: yt_s
+    REAL, DIMENSION(knon),        INTENT(IN)        :: ydeltat_s
+    REAL, DIMENSION(knon),        INTENT(IN)        :: ygustiness
+    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: yt_x, yt_w, yq_x, yq_w
+    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: yu_x, yu_w, yv_x, yv_w
+    REAL, DIMENSION(knon),        INTENT(IN)        :: ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w
+    REAL, DIMENSION(knon),        INTENT(IN)        :: ycdragm_x, ycdragm_w
+    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT_x, AcoefT_w, AcoefQ_x, AcoefQ_w
+    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w
+    REAL, DIMENSION(knon),        INTENT(IN)        :: BcoefT_x, BcoefT_w, BcoefQ_x, BcoefQ_w
+    REAL, DIMENSION(knon),        INTENT(IN)        :: BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w
+!
+! Local variables
+    INTEGER                    :: j
+    REAL                       :: rho1
+    REAL                       :: mod_wind_x
+    REAL                       :: mod_wind_w   
+    REAL                       :: dd_Kh
+    REAL                       :: dd_Kq
+    REAL                       :: dd_Km
+!
+    REAL                       :: zdelta, zcvm5, zcor, qsat
+! 
+    REAL, DIMENSION(knon)      :: sigx       ! fractional area of (x) region
+!
+!!!
+!!! jyg le 09/04/2013 ; passage aux nouvelles expressions en differences
+
+!
+!  First computations
+!  ------------------
+   eps_1 = 0.5
+   smallestreal=tiny(smallestreal)
+!
+   sigx(1:knon) = 1.-sigw(1:knon)
+! Compute Cp, Lv, qsat, dqsat_dT.
+   L_v(1:knon) = RLvtt
+   Ts0(1:knon) = yt_s(1:knon)
+   dTs0(1:knon) = ydeltat_s(1:knon)
+   q1_0b(1:knon) = sigw(1:knon)*yq_w(1:knon,1)+sigx(1:knon)*yq_x(1:knon,1)
+!
+! fqsat determination 
+! -------------------
+   IF (nsrf == is_oce) THEN
+     fqsat = f_qsat_oce
+   ELSE
+     fqsat = 1.
+   ENDIF
+!
+!
+!  Reference state
+!  ---------------
+   DO j = 1, knon
+      zdelta = MAX(0.,SIGN(1.,RTT-Ts0(j)))
+      zcvm5 = R5LES*(1.-zdelta) + R5IES*zdelta
+      qsat = R2ES*FOEEW(Ts0(j),zdelta)/ypaprs(j,1)
+      qsat = MIN(0.5,qsat)
+      zcor = 1./(1.-RETV*qsat)
+      qsat0(j) = fqsat*qsat*zcor
+      dqsatdT0(j) = fqsat*FOEDE(Ts0(j),zdelta,zcvm5,qsat0(j),zcor)
+      C_p(j) = RCpd + qsat0(j)*(RCpv - RCpd)
+      C_p(j) = RCpd
+!
+!      print *,' AAAA wx_pbl0, C_p(j), qsat0(j), Ts0(j) : ', C_p(j), qsat0(j), Ts0(j)
+   ENDDO
+   DO j = 1, knon
+      Ts0_x(j) = Ts0(j) - sigw(j)*dTs0(j)
+      zdelta = MAX(0.,SIGN(1.,RTT-Ts0_x(j)))
+      zcvm5 = R5LES*(1.-zdelta) + R5IES*zdelta
+!!      zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
+      qsat = R2ES*FOEEW(Ts0_x(j),zdelta)/ypaprs(j,1)
+      qsat = MIN(0.5,qsat)
+      zcor = 1./(1.-RETV*qsat)
+      qsat0_x(j) = qsat*zcor
+      dqsatdT0_x(j) = FOEDE(Ts0_x(j),zdelta,zcvm5,qsat0_x(j),zcor)
+!!      dqsatdT0_x(j) = (RLvtt*(1.-zdelta)+RLSTT*zdelta)*qsat0_x(j)/(Rv*Ts0_x(j)*Ts0_x(j))
+   ENDDO
+   DO j = 1, knon
+      Ts0_w(j) = Ts0(j) + sigx(j)*dTs0(j)
+      zdelta = MAX(0.,SIGN(1.,RTT-Ts0_w(j)))
+      zcvm5 = R5LES*(1.-zdelta) + R5IES*zdelta
+!!      zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
+      qsat = R2ES*FOEEW(Ts0_w(j),zdelta)/ypaprs(j,1)
+      qsat = MIN(0.5,qsat)
+      zcor = 1./(1.-RETV*qsat)
+      qsat0_w(j) = qsat*zcor
+      dqsatdT0_w(j) = FOEDE(Ts0_w(j),zdelta,zcvm5,qsat0_w(j),zcor)
+!!      dqsatdT0_w(j) = (RLvtt*(1.-zdelta)+RLSTT*zdelta)*qsat0_w(j)/(Rv*Ts0_w(j)*Ts0_w(j))
+   ENDDO
+!
+   QQ_x(1:knon)  = 1./dqsatdT0_x(1:knon)
+   QQ_w(1:knon)  = 1./dqsatdT0_w(1:knon)
+   QQ_b(1:knon)  = sigw(1:knon)*QQ_w(1:knon) + sigx(1:knon)*QQ_x(1:knon)
+   dd_QQ(1:knon) = QQ_w(1:knon) - QQ_x(1:knon)
+!
+        DO j=1,knon
+!
+! Exchange coefficients computation
+! ---------------------------------
+!
+! Wind factor (Warning : this is not valid when using land_surf_orchidee)
+         mod_wind_x = min_wind_speed+SQRT(ygustiness(j)+yu_x(j,1)**2+yv_x(j,1)**2)
+         mod_wind_w = min_wind_speed+SQRT(ygustiness(j)+yu_w(j,1)**2+yv_w(j,1)**2)
+!
+!!         rho1 = ypplay(j,1)/(RD*yt(j,1))
+         rho1 = ypplay(j,1)/(RD*(yt_x(j,1) + sigw(j)*(yt_w(j,1)-yt_x(j,1))))
+! 
+! (w) and (x) Exchange coefficients
+         Kech_h_x(j) = ycdragh_x(j) * mod_wind_x * rho1
+         Kech_h_w(j) = ycdragh_w(j) * mod_wind_w * rho1
+         Kech_q_x(j) = ycdragq_x(j) * mod_wind_x * rho1
+         Kech_q_w(j) = ycdragq_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
+!!  Print *,'YYYYpbl0: ycdragh_x, ycdragq_x, mod_wind_x, rho1, Kech_q_x, Kech_h_x ', &
+!!                     ycdragh_x(j), ycdragq_x(j), mod_wind_x, rho1, Kech_q_x(j), Kech_h_x(j)
+!!  Print *,'YYYYpbl0: ycdragh_w, ycdragq_w, mod_wind_w, rho1, Kech_q_w, Kech_h_w ', &
+!!                     ycdragh_w(j), ycdragq_w(j), mod_wind_w, rho1, Kech_q_w(j), Kech_h_w(j)
+!
+! Merged exchange coefficients
+         dd_Kh = Kech_h_w(j) - Kech_h_x(j)
+         dd_Kq = Kech_q_w(j) - Kech_q_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) + sigw(j)*dd_Kh
+         Kech_q(j) = Kech_q_x(j) + sigw(j)*dd_Kq
+         Kech_m(j) = Kech_m_x(j) + sigw(j)*dd_Km
+!
+! Effective exchange coefficient computation
+! ------------------------------------------
+        Kech_T_px(j) = Kech_h_x(j)/(1.-BcoefT_x(j)*Kech_h_x(j)*dtime)
+        Kech_T_pw(j) = Kech_h_w(j)/(1.-BcoefT_w(j)*Kech_h_w(j)*dtime)
+!
+        Kech_Q_px(j) = Kech_q_x(j)/(1.-BcoefQ_x(j)*Kech_q_x(j)*dtime)
+        Kech_Q_pw(j) = Kech_q_w(j)/(1.-BcoefQ_w(j)*Kech_q_w(j)*dtime)
+!
+        Kech_U_px(j) = Kech_m_x(j)/(1.-BcoefU_x(j)*Kech_m_x(j)*dtime)
+        Kech_U_pw(j) = Kech_m_w(j)/(1.-BcoefU_w(j)*Kech_m_w(j)*dtime)
+!
+        Kech_V_px(j) = Kech_m_x(j)/(1.-BcoefV_x(j)*Kech_m_x(j)*dtime)
+        Kech_V_pw(j) = Kech_m_w(j)/(1.-BcoefV_w(j)*Kech_m_w(j)*dtime)
+!
+         dd_KTp(j) = Kech_T_pw(j) - Kech_T_px(j)
+         dd_KQp(j) = Kech_Q_pw(j) - Kech_Q_px(j)
+         dd_KUp(j) = Kech_U_pw(j) - Kech_U_px(j)
+         dd_KVp(j) = Kech_V_pw(j) - Kech_V_px(j)
+!
+        Kech_Tp(j) = Kech_T_px(j) + sigw(j)*dd_KTp(j)
+        Kech_Qp(j) = Kech_Q_px(j) + sigw(j)*dd_KQp(j)
+        Kech_Up(j) = Kech_U_px(j) + sigw(j)*dd_KUp(j)
+        Kech_Vp(j) = Kech_V_px(j) + sigw(j)*dd_KVp(j)
+!
+! Store AQ and BQ in the module variables
+       AQ_x(j) = AcoefQ_x(j)
+       AQ_w(j) = AcoefQ_w(j)
+       BQ_x(j) = BcoefQ_x(j)
+       BQ_w(j) = BcoefQ_w(j)
+!
+! Calcul des differences w-x
+       dd_Cdragm(j) = ycdragm_w(j) - ycdragm_x(j)
+       dd_Cdragh(j) = ycdragh_w(j) - ycdragh_x(j)
+       dd_Cdragq(j) = ycdragq_w(j) - ycdragq_x(j)
+!
+       dd_AT(j) = AcoefT_w(j) - AcoefT_x(j)
+       dd_AQ(j) = AcoefQ_w(j) - AcoefQ_x(j)
+       dd_AU(j) = AcoefU_w(j) - AcoefU_x(j)
+       dd_AV(j) = AcoefV_w(j) - AcoefV_x(j)
+       dd_BT(j) = BcoefT_w(j) - BcoefT_x(j)
+       dd_BQ(j) = BcoefQ_w(j) - BcoefQ_x(j)
+       dd_BU(j) = BcoefU_w(j) - BcoefU_x(j)
+       dd_BV(j) = BcoefV_w(j) - BcoefV_x(j)
+!
+       KxKwTp(j) = Kech_T_px(j)*Kech_T_pw(j)
+       KxKwQp(j) = Kech_Q_px(j)*Kech_Q_pw(j)
+       KxKwUp(j) = Kech_U_px(j)*Kech_U_pw(j)
+       KxKwVp(j) = Kech_V_px(j)*Kech_V_pw(j)
+!                                           
+!
+       IF (prt_level >=10) THEN
+          print *,'Variables pour la fusion : Kech_T_px(j)' ,Kech_T_px(j)
+          print *,'Variables pour la fusion : Kech_T_pw(j)' ,Kech_T_pw(j)
+          print *,'Variables pour la fusion : Kech_Tp(j)' ,Kech_Tp(j)
+          print *,'Variables pour la fusion : Kech_h(j)' ,Kech_h(j)
+       ENDIF
+
+     ENDDO  ! j = 1, knon
+
+     RETURN
+
+END SUBROUTINE wx_pbl_prelim_0
+
+SUBROUTINE wx_pbl_prelim_beta(knon, dtime,  &
+                                 sigw, beta,       &
+                                 BcoefQ_x, BcoefQ_w &
+                                 )
+!
+    USE print_control_mod, ONLY: prt_level,lunout
+    USE indice_sol_mod, ONLY: is_oce
+!
+    INTEGER,                      INTENT(IN)        :: knon    ! number of grid cells
+    REAL,                         INTENT(IN)        :: dtime   ! time step size (s)
+    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw ! cold pools fractional area
+    REAL, DIMENSION(knon),        INTENT(IN)        :: beta ! evaporation by potential evaporation
+    REAL, DIMENSION(knon),        INTENT(IN)        :: BcoefQ_x, BcoefQ_w
+!
+! Local variables
+    INTEGER                    :: j
+!
+   DO j = 1, knon
+!
+        qsatsrf0_x(j) = beta(j)*qsat0_x(j)
+        qsatsrf0_w(j) = beta(j)*qsat0_w(j)
+        dqsatsrf0(j)  = qsatsrf0_w(j) - qsatsrf0_x(j)
+!
+        Kech_Q_sx(j) = Kech_q_x(j)/(1.-beta(j)*BcoefQ_x(j)*Kech_q_x(j)*dtime)
+        Kech_Q_sw(j) = Kech_q_w(j)/(1.-beta(j)*BcoefQ_w(j)*Kech_q_w(j)*dtime)
+!
+        dd_KQs(j) = Kech_Q_sw(j) - Kech_Q_sx(j)
+!
+        Kech_Qs(j) = Kech_Q_sx(j) + sigw(j)*dd_KQs(j)
+!
+        KxKwQs(j) = Kech_Q_sx(j)*Kech_Q_sw(j)
+!
+!!  print *,'BBBBwx_prelim_beta : beta ', beta(j)
+!
+  ENDDO ! j = 1, knon
+
+  RETURN
+
+END SUBROUTINE wx_pbl_prelim_beta
+
+SUBROUTINE wx_pbl_final
+!
+!****************************************************************************************
+! Deallocate module variables
+!
+!****************************************************************************************    
+!
+    IF (ALLOCATED(C_p))           DEALLOCATE(C_p)
+    IF (ALLOCATED(L_v))           DEALLOCATE(L_v)
+    IF (ALLOCATED(Ts0))           DEALLOCATE(Ts0)
+    IF (ALLOCATED(dTs0))          DEALLOCATE(dTs0)
+    IF (ALLOCATED(Ts0_x))         DEALLOCATE(Ts0_x)
+    IF (ALLOCATED(Ts0_w))         DEALLOCATE(Ts0_w)
+    IF (ALLOCATED(qsat0))         DEALLOCATE(qsat0)
+    IF (ALLOCATED(dqsatdT0))      DEALLOCATE(dqsatdT0 )
+    IF (ALLOCATED(qsat0_x))       DEALLOCATE(qsat0_x)
+    IF (ALLOCATED(dqsatdT0_x))    DEALLOCATE(dqsatdT0_x )
+    IF (ALLOCATED(qsat0_w))       DEALLOCATE(qsat0_w)
+    IF (ALLOCATED(dqsatdT0_w))    DEALLOCATE(dqsatdT0_w )
+    IF (ALLOCATED(q1_0b))         DEALLOCATE(q1_0b)
+    IF (ALLOCATED(QQ_b))          DEALLOCATE(QQ_b)
+    IF (ALLOCATED(dd_QQ))         DEALLOCATE(dd_QQ)
+    IF (ALLOCATED(QQ_x))          DEALLOCATE(QQ_x)
+    IF (ALLOCATED(QQ_w))          DEALLOCATE(QQ_w)
+    IF (ALLOCATED(qsatsrf0_x))    DEALLOCATE(qsatsrf0_x)
+    IF (ALLOCATED(qsatsrf0_w))    DEALLOCATE(qsatsrf0_w)
+    IF (ALLOCATED(dqsatsrf0))     DEALLOCATE(dqsatsrf0)
+    IF (ALLOCATED(dd_Cdragh))     DEALLOCATE(dd_Cdragh)
+    IF (ALLOCATED(dd_Cdragm))     DEALLOCATE(dd_Cdragm)
+    IF (ALLOCATED(dd_Cdragq))     DEALLOCATE(dd_Cdragq)
+    IF (ALLOCATED(Kech_h))        DEALLOCATE(Kech_h)
+    IF (ALLOCATED(Kech_h_x))      DEALLOCATE(Kech_h_x)
+    IF (ALLOCATED(Kech_h_w))      DEALLOCATE(Kech_h_w)
+    IF (ALLOCATED(Kech_q))        DEALLOCATE(Kech_q)
+    IF (ALLOCATED(Kech_q_x))      DEALLOCATE(Kech_q_x)
+    IF (ALLOCATED(Kech_q_w))      DEALLOCATE(Kech_q_w)
+    IF (ALLOCATED(Kech_m))        DEALLOCATE(Kech_m)
+    IF (ALLOCATED(Kech_m_x))      DEALLOCATE(Kech_m_x)
+    IF (ALLOCATED(Kech_m_w))      DEALLOCATE(Kech_m_w)
+    IF (ALLOCATED(Kech_Tp))       DEALLOCATE(Kech_Tp)
+    IF (ALLOCATED(Kech_T_px))     DEALLOCATE(Kech_T_px)
+    IF (ALLOCATED(Kech_T_pw))     DEALLOCATE(Kech_T_pw)
+    IF (ALLOCATED(dd_KTp))        DEALLOCATE(dd_KTp)
+    IF (ALLOCATED(KxKwTp))        DEALLOCATE(KxKwTp)
+    IF (ALLOCATED(dd_AT))         DEALLOCATE(dd_AT)
+    IF (ALLOCATED(dd_BT))         DEALLOCATE(dd_BT)
+    IF (ALLOCATED(Kech_Qp))       DEALLOCATE(Kech_Qp)
+    IF (ALLOCATED(Kech_Q_px))     DEALLOCATE(Kech_Q_px)
+    IF (ALLOCATED(Kech_Q_pw))     DEALLOCATE(Kech_Q_pw)
+    IF (ALLOCATED(dd_KQp))        DEALLOCATE(dd_KQp)
+    IF (ALLOCATED(KxKwQp))        DEALLOCATE(KxKwQp)
+    IF (ALLOCATED(Kech_Qs))       DEALLOCATE(Kech_Qs)
+    IF (ALLOCATED(Kech_Q_sx))     DEALLOCATE(Kech_Q_sx)
+    IF (ALLOCATED(Kech_Q_sw))     DEALLOCATE(Kech_Q_sw)
+    IF (ALLOCATED(dd_KQs))        DEALLOCATE(dd_KQs)
+    IF (ALLOCATED(KxKwQs ))       DEALLOCATE(KxKwQs )
+    IF (ALLOCATED(AQ_x))          DEALLOCATE(AQ_x)
+    IF (ALLOCATED(AQ_w))          DEALLOCATE(AQ_w)
+    IF (ALLOCATED(BQ_x))          DEALLOCATE(BQ_x)
+    IF (ALLOCATED(BQ_w))          DEALLOCATE(BQ_w)
+    IF (ALLOCATED(dd_AQ))         DEALLOCATE(dd_AQ)
+    IF (ALLOCATED(dd_BQ ))        DEALLOCATE(dd_BQ )
+    IF (ALLOCATED(Kech_Up))       DEALLOCATE(Kech_Up)
+    IF (ALLOCATED(Kech_U_px))     DEALLOCATE(Kech_U_px)
+    IF (ALLOCATED(Kech_U_pw))     DEALLOCATE(Kech_U_pw)
+    IF (ALLOCATED(dd_KUp))        DEALLOCATE(dd_KUp)
+    IF (ALLOCATED(KxKwUp))        DEALLOCATE(KxKwUp)
+    IF (ALLOCATED(dd_AU))         DEALLOCATE(dd_AU)
+    IF (ALLOCATED(dd_BU))         DEALLOCATE(dd_BU)
+    IF (ALLOCATED(Kech_Vp))       DEALLOCATE(Kech_Vp)
+    IF (ALLOCATED(Kech_V_px))     DEALLOCATE(Kech_V_px)
+    IF (ALLOCATED(Kech_V_pw))     DEALLOCATE(Kech_V_pw)
+    IF (ALLOCATED(dd_KVp))        DEALLOCATE(dd_KVp)
+    IF (ALLOCATED(KxKwVp))        DEALLOCATE(KxKwVp)
+    IF (ALLOCATED(dd_AV))         DEALLOCATE(dd_AV)
+    IF (ALLOCATED(dd_BV))         DEALLOCATE(dd_BV)
+
+END SUBROUTINE wx_pbl_final
+
+END MODULE wx_pbl_var_mod
+
