MODULE wx_pbl_mod

! 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

  IMPLICIT NONE

CONTAINS

!****************************************************************************************

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, 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, &
                                 AcoefT, AcoefQ, AcoefU, AcoefV, &
                                 BcoefT, BcoefQ, BcoefU, BcoefV, &
                                 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, 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)        :: 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, 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
    REAL, DIMENSION(knon),        INTENT(OUT)       :: AcoefT, AcoefQ, AcoefU, AcoefV
    REAL, DIMENSION(knon),        INTENT(OUT)       :: BcoefT, BcoefQ, BcoefU, BcoefV
    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

! Local variables
    INTEGER                    :: j
    REAL                       :: dd_Kh
    REAL                       :: dd_Kq
    REAL                       :: dd_Km
    REAL                       :: dd_u 
    REAL                       :: dd_v 
    REAL                       :: dd_t 
    REAL                       :: dd_q 

    REAL                       :: LambdaTs, LambdaQs, LambdaUs, LambdaVs 

    REAL, DIMENSION(knon)      :: sigx       ! fractional area of (x) region


   sigx(1:knon) = 1.-sigw(1:knon)


        DO j=1,knon


! 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)

! 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

! 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) = (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_Kq*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) + 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) + 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


        ENDDO

        RETURN

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,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 >= 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 >= 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)      :: 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.-sigw(:)

!      print *,' AAAA wx_pbl_split, C_p(j), qsat0(j), Ts0(j) : ', C_p(:), qsat0(:), Ts0(:)

   IF (iflag_split == 2 .AND. nsrf /= 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

!!    IF (knon .NE. 0) THEN
!!       call  iophys_ecrit('sigw', 1,'sigw', '.',sigw)
!!       call  iophys_ecrit('phit', 1,'phit', 'W/m2',phit)
!!       call  iophys_ecrit('phit_w', 1,'phit_w', 'W/m2',phit_w)
!!       call  iophys_ecrit('phit_x', 1,'phit_x', 'W/m2',phit_x)
!!       call  iophys_ecrit('phiq', 1,'phiq', 'kg/m2/s',phiq)
!!       call  iophys_ecrit('phiq_w', 1,'phiq_w', 'kg/m2/s',phiq_w)
!!       call  iophys_ecrit('phiq_x', 1,'phiq_x', 'kg/m2/s',phiq_x)
!!       call  iophys_ecrit('q1_w', 1,'q1_w', '.',q1_w)
!!       call  iophys_ecrit('q1_x', 1,'q1_x', '.',q1_x)
!!    ENDIF  ! (knon .NE. 0)

        RETURN

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

