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(qsat0, dqsatdT0)
  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, &
                                 Kech_h_x_out, Kech_h_w_out, Kech_h_out  &
                                 )

    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

    REAL, DIMENSION(knon),        INTENT(OUT)       :: Kech_h_x_out, Kech_h_w_out, Kech_h_out

! 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) = fqsat*qsat*zcor
      dqsatdT0_x(j) = fqsat*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) = fqsat*qsat*zcor
      dqsatdT0_w(j) = fqsat*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

         Kech_h_x_out(j) = Kech_h_x(j)
         Kech_h_w_out(j) = Kech_h_w(j)
         Kech_h_out(j)   = Kech_h(j)

! 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

