Changeset 3888


Ignore:
Timestamp:
May 5, 2021, 12:50:37 PM (3 years ago)
Author:
jyg
Message:

New provisional version of the splitting of the
diffusive boundary layer into inwake and offwake
PBLs. The splitting of the diffuse BL should NOT
be activated yet for general purpose simulations.

The splitting is activated by:
mod(iflag_pbl_split,10)=1 for the option with
fixed surface temperature and
mod(iflag_pbl_split,10)=2 for the option with
coupled surface temperature.

iflag_pbl_split=0 ==> no splittingat all.
iflag_pbl_split=10 ==> splitting of thermals.
iflag_pbl_split=11 ==> splitting of thermals and
of vertical diffusion (fixed surf. temp.).
iflag_pbl_split=12 ==> splitting of thermals and
of vertical diffusion (coupled surf. temp.).

Location:
LMDZ6/trunk/libf/phylmd
Files:
1 added
11 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/dyn1d/compar1d.h

    r3594 r3888  
    88      real :: nat_surf
    99      real :: tsurf
     10      real :: beta_surf
    1011      real :: rugos
    1112      real :: rugosh
    … …  
    4546      real    :: p_nudging_u, p_nudging_v, p_nudging_w, p_nudging_t, p_nudging_qv
    4647      common/com_par1d/                                                 &
    47      & nat_surf,tsurf,rugos,rugosh,                                     &
     48     & nat_surf,tsurf,beta_surf,rugos,rugosh,                           &
    4849     & xqsol,qsurf,psurf,zsurf,albedo,time,time_ini,xlat,xlon,airefi,   &
    4950     & wtsurf,wqsurf,restart_runoff,xagesno,qsolinp,zpicinp,            &

  • LMDZ6/trunk/libf/phylmd/dyn1d/old_lmdz1d.F90

    r3781 r3888  
    1111       du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
    1212       falb_dir, falb_dif, &
    13        ftsol, pbl_tke, pctsrf, radsol, rain_fall, snow_fall, ratqs, &
     13       ftsol, beta_aridity, pbl_tke, pctsrf, radsol, rain_fall, snow_fall, ratqs, &
    1414       rnebcon, rugoro, sig1, w01, solaire_etat0, sollw, sollwdown, &
    1515       solsw, t_ancien, q_ancien, u_ancien, v_ancien, wake_cstar, &
    … …  
    656656      qsol = qsolinp
    657657      qsurf = fq_sat(tsurf,psurf/100.)
     658      beta_aridity(:,:) = beta_surf
    658659      day1= day_ini
    659660      time=daytime-day
    … …  
    795796
    796797        fder=0.
     798        snsrf(1,:)=snowmass ! masse de neige des sous surface
    797799        print *, 'snsrf', snsrf
    798         snsrf(1,:)=snowmass ! masse de neige des sous surface
    799800        qsurfsrf(1,:)=qsurf ! humidite de l'air des sous surface
    800801        fevap=0.

  • LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90

    r3781 r3888  
    77       du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
    88       falb_dir, falb_dif, &
    9        ftsol, pbl_tke, pctsrf, radsol, rain_fall, snow_fall, ratqs, &
     9       ftsol, beta_aridity, pbl_tke, pctsrf, radsol, rain_fall, snow_fall, ratqs, &
    1010       rnebcon, rugoro, sig1, w01, solaire_etat0, sollw, sollwdown, &
    1111       solsw, t_ancien, q_ancien, u_ancien, v_ancien, wake_cstar, &
    … …  
    429429      qsol = qsolinp
    430430      qsurf = fq_sat(tsurf,psurf/100.)
     431      beta_aridity(:,:) = beta_surf
    431432      day1= day_ini
    432433      time=daytime-day

  • LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90

    r3876 r3888  
    2626  USE cdrag_mod
    2727  USE stdlevvar_mod
    28   USE wx_pbl_mod,          ONLY : wx_pbl_init, wx_pbl_final, &
    29 !!                                  wx_pbl_fuse_no_dts, wx_pbl_split_no_dts, &
    30 !!                                  wx_pbl_fuse, wx_pbl_split
    31                                   wx_pbl0_fuse, wx_pbl0_split
     28  USE wx_pbl_var_mod,      ONLY : wx_pbl_init, wx_pbl_final, &
     29                                  wx_pbl_prelim_0, wx_pbl_prelim_beta
     30  USE wx_pbl_mod,          ONLY : wx_pbl0_merge, wx_pbl_split, wx_pbl_dts_merge, &
     31                                  wx_pbl_check, wx_pbl_dts_check, wx_evappot
    3232  use config_ocean_skin_m, only: activate_ocean_skin
    3333
    … …  
    3737  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: fder   ! flux drift
    3838  !$OMP THREADPRIVATE(fder)
    39   REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: snow   ! snow at surface
     39  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE    :: snow   ! snow at surface
    4040  !$OMP THREADPRIVATE(snow)
    4141  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: qsurf  ! humidity at surface
    4242  !$OMP THREADPRIVATE(qsurf)
    43   REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ftsoil ! soil temperature
     43  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE          :: ftsoil ! soil temperature
    4444  !$OMP THREADPRIVATE(ftsoil)
     45  REAL, ALLOCATABLE, DIMENSION(:), SAVE              :: ydTs0, ydqs0 
     46                                                     ! nul forced temperature and humidity differences
     47  !$OMP THREADPRIVATE(ydTs0, ydqs0)
    4548
    4649  INTEGER, SAVE :: iflag_pbl_surface_t2m_bug
    … …  
    99102    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
    100103
     104    ALLOCATE(ydTs0(klon), stat=ierr)
     105    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
     106
     107    ALLOCATE(ydqs0(klon), stat=ierr)
     108    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
     109
    101110    fder(:)       = fder_rst(:)
    102111    snow(:,:)     = snow_rst(:,:)
    103112    qsurf(:,:)    = qsurf_rst(:,:)
    104113    ftsoil(:,:,:) = ftsoil_rst(:,:,:)
     114    ydTs0(:) = 0.
     115    ydqs0(:) = 0.
    105116
    106117!****************************************************************************************
    … …  
    183194       ts,SFRWL,   alb_dir, alb_dif,ustar, u10m, v10m,wstar, &
    184195       cdragh,    cdragm,   zu1,    zv1,              &
     196!jyg<   (26/09/2019)
     197       beta, &
     198!>jyg
    185199       alb_dir_m,    alb_dif_m,  zxsens,   zxevap,    &
    186200       alb3_lic,  runoff,    snowhgt,   qsnow,     to_ice,    sissnow,  &
    … …  
    206220       s_therm,   s_trmb1,   s_trmb2,  s_trmb3,       &
    207221       zustar,zu10m,  zv10m,    fder_print,    &
    208        zxqsurf,   rh2m,      zxfluxu,  zxfluxv,       &
     222       zxqsurf, delta_qsurf,                       &
     223       rh2m,      zxfluxu,  zxfluxv,               &
    209224       z0m, z0h,   agesno,  sollw,    solsw,         &
    210225       d_ts,      evap,    fluxlat,  t2m,           &
    … …  
    366381! Input/Output variables
    367382!****************************************************************************************
     383!jyg<
     384    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: beta    ! Aridity factor
     385!>jyg
    368386    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ts      ! temperature at surface (K)
    369387    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: delta_tsurf !surface temperature difference between
    … …  
    468486    REAL, DIMENSION(klon),        INTENT(OUT)       :: fder_print ! fder for printing (=fder(i) + dflux_t(i) + dflux_q(i))
    469487    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxqsurf    ! humidity at surface, mean for each grid point
     488    REAL, DIMENSION(klon),        INTENT(OUT)       :: delta_qsurf! humidity difference at surface, mean for each grid point
    470489    REAL, DIMENSION(klon),        INTENT(OUT)       :: rh2m       ! relative humidity at 2m
    471490    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zxfluxu    ! u wind tension, mean for each grid point
    … …  
    540559    REAL, DIMENSION(klon)              :: r_co2_ppm     ! taux CO2 atmosphere
    541560    REAL, DIMENSION(klon)              :: yts, yz0m, yz0h, ypct
     561    REAL, DIMENSION(klon)              :: yz0h_old
    542562!albedo SB >>>
    543563    REAL, DIMENSION(klon)              :: yalb,yalb_vis
    544564!albedo SB <<<
    545565    REAL, DIMENSION(klon)              :: yt1, yq1, yu1, yv1
     566    REAL, DIMENSION(klon)              :: yqa
    546567    REAL, DIMENSION(klon)              :: ysnow, yqsurf, yagesno, yqsol
    547568    REAL, DIMENSION(klon)              :: yrain_f, ysnow_f
    … …  
    577598    REAL, DIMENSION(klon)              :: yz0h_oupas
    578599    REAL, DIMENSION(klon)              :: yfluxsens
     600    REAL, DIMENSION(klon)              :: AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0
    579601    REAL, DIMENSION(klon)              :: AcoefH, AcoefQ, BcoefH, BcoefQ
    580602    REAL, DIMENSION(klon)              :: AcoefU, AcoefV, BcoefU, BcoefV
    581603    REAL, DIMENSION(klon)              :: ypsref
    582     REAL, DIMENSION(klon)              :: yevap, ytsurf_new, yalb3_new
     604    REAL, DIMENSION(klon)              :: yevap, yevap_pot, ytsurf_new, yalb3_new
    583605!albedo SB >>>
    584606    REAL, DIMENSION(klon,nsw)          :: yalb_dir_new, yalb_dif_new
    … …  
    592614    REAL, DIMENSION(klon,klev)         :: y_flux_u, y_flux_v
    593615    REAL, DIMENSION(klon,klev)         :: ycoefh, ycoefm,ycoefq
    594     REAL, DIMENSION(klon)              :: ycdragh, ycdragm
     616    REAL, DIMENSION(klon)              :: ycdragh, ycdragq, ycdragm
    595617    REAL, DIMENSION(klon,klev)         :: yu, yv
    596618    REAL, DIMENSION(klon,klev)         :: yt, yq
    … …  
    624646    REAL, DIMENSION(klon,klev)         :: ycoefh_x, ycoefm_x, ycoefh_w, ycoefm_w
    625647    REAL, DIMENSION(klon,klev)         :: ycoefq_x, ycoefq_w
    626     REAL, DIMENSION(klon)              :: ycdragh_x, ycdragm_x, ycdragh_w, ycdragm_w
     648    REAL, DIMENSION(klon)              :: ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w
     649    REAL, DIMENSION(klon)              :: ycdragm_x, ycdragm_w
    627650    REAL, DIMENSION(klon)              :: AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x
    628651    REAL, DIMENSION(klon)              :: AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w
    … …  
    645668    REAL                               :: zx_qs_surf, zcor_surf, zdelta_surf
    646669    REAL, DIMENSION(klon)              :: ytsurf_th, yqsatsurf
     670!jyg<
    647671    REAL, DIMENSION(klon)              :: ybeta
     672    REAL, DIMENSION(klon)              :: ybeta_prev
     673!>jyg
    648674    REAL, DIMENSION(klon, klev)        :: d_u_x
    649675    REAL, DIMENSION(klon, klev)        :: d_u_w
    … …  
    780806!!! nrlmd le 13/06/2011
    781807    REAL, DIMENSION(klon)              :: y_delta_flux_t1, y_delta_flux_q1, y_delta_flux_u1, y_delta_flux_v1
    782     REAL, DIMENSION(klon)              :: y_delta_tsurf,delta_coef,tau_eq
     808    REAL, DIMENSION(klon)              :: y_delta_tsurf, y_delta_tsurf_new
     809    REAL, DIMENSION(klon)              :: delta_coef, tau_eq
     810    REAL, DIMENSION(klon)              :: HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn
     811    REAL, DIMENSION(klon)              :: phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0
     812    REAL, DIMENSION(klon)              :: y_delta_qsurf
     813    REAL, DIMENSION(klon)              :: y_delta_qsats
     814    REAL, DIMENSION(klon)              :: yg_T, yg_Q
     815    REAL, DIMENSION(klon)              :: yGamma_dTs_phiT, yGamma_dQs_phiQ
     816    REAL, DIMENSION(klon)              :: ydTs_ins, ydqs_ins
     817!
    783818    REAL, PARAMETER                    :: facteur=2./sqrt(3.14)
    784819    REAL, PARAMETER                    :: inertia=2000.
    … …  
    793828    REAL, DIMENSION(klon)              :: Kech_m
    794829    REAL, DIMENSION(klon)              :: Kech_m_x, Kech_m_w
    795     REAL, DIMENSION(klon)              :: yts_x,yts_w
     830    REAL, DIMENSION(klon)              :: yts_x, yts_w
     831    REAL, DIMENSION(klon)              :: yqsatsrf0_x, yqsatsrf0_w
     832    REAL, DIMENSION(klon)              :: yqsurf_x, yqsurf_w
    796833!jyg<
    797834!!    REAL, DIMENSION(klon)              :: Kech_Hp, Kech_H_xp, Kech_H_wp
    … …  
    800837!!    REAL, DIMENSION(klon)              :: Kech_Vp, Kech_V_xp, Kech_V_wp
    801838!>jyg
    802 !jyg<
    803     REAL, DIMENSION(klon)              :: ah, bh     ! coefficients of the delta_Tsurf equation
    804 !>jyg
     839
     840    REAL                               :: fact_cdrag
     841    REAL                               :: z1lay
    805842
    806843    REAL                               :: vent
    … …  
    950987 fder_print(:)=0.
    951988 zxqsurf(:)=0.
     989 delta_qsurf(:) = 0.
    952990 zxfluxu(:,:)=0. ; zxfluxv(:,:)=0.
    953991 solsw(:,:)=0. ; sollw(:,:)=0.
    … …  
    10321070    y_delta_flux_t1=0.
    10331071    ydtsurf_th=0.
    1034     yts_x=0.      ; yts_w=0.
    1035     y_delta_tsurf=0.
     1072    yts_x(:)=0.      ; yts_w(:)=0.
     1073    y_delta_tsurf(:)=0. ; y_delta_qsurf(:)=0.
     1074    yqsurf_x(:)=0.      ; yqsurf_w(:)=0.
     1075    yg_T(:) = 0. ;        yg_Q(:) = 0.
     1076    yGamma_dTs_phiT(:) = 0. ; yGamma_dQs_phiQ(:) = 0.
     1077    ydTs_ins(:) = 0. ; ydqs_ins(:) = 0.
     1078
    10361079!!!
    10371080    ytsoil = 999999.
    … …  
    13021345!****************************************************************************************
    13031346
     1347!
     1348!jyg<    (20190926)
     1349!   Provisional : set ybeta to standard values
     1350       IF (nsrf .NE. is_ter) THEN
     1351           ybeta(:) = 1.
     1352       ELSE
     1353           IF (iflag_split .EQ. 0) THEN
     1354              ybeta(:) = 1.
     1355           ELSE
     1356             DO j = 1, knon
     1357                i = ni(j)
     1358                ybeta(j)   = beta(i,nsrf)
     1359             ENDDO
     1360           ENDIF  ! (iflag_split .LE.1)
     1361       ENDIF !  (nsrf .NE. is_ter)
     1362!>jyg
     1363!
    13041364       DO j = 1, knon
    13051365          i = ni(j)
    … …  
    15211581        CALL cdrag(knon, nsrf, &
    15221582            speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),&
    1523             yts_x, yqsurf, yz0m, yz0h, &
     1583            yts_x, yqsurf_x, yz0m, yz0h, &
    15241584            ycdragm_x, ycdragh_x, zri1_x, pref_x )
    15251585
    … …  
    15481608        CALL cdrag(knon, nsrf, &
    15491609            speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),&
    1550             yts_w, yqsurf, yz0m, yz0h, &
     1610            yts_w, yqsurf_w, yz0m, yz0h, &
    15511611            ycdragm_w, ycdragh_w, zri1_w, pref_w )
    15521612!
    … …  
    16211681       ENDIF
    16221682        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
    1623             ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf, ycdragm_x, &
     1683            ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf_x, ycdragm_x, &
    16241684            ycoefm_x, ycoefh_x, ytke_x,y_treedrg)
    16251685!            ycoefm_x, ycoefh_x, ytke_x)
    … …  
    16491709       ENDIF
    16501710        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
    1651             ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf, ycdragm_w, &
     1711            ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf_w, ycdragm_w, &
    16521712            ycoefm_w, ycoefh_w, ytke_w,y_treedrg)
    16531713!            ycoefm_w, ycoefh_w, ytke_w)
    … …  
    17861846         yt1(:) = yt(:,1)
    17871847         yq1(:) = yq(:,1)
    1788 !!       ELSE IF (iflag_split .eq. 1) THEN
    1789 !!!
    1790 !jyg<
    1791 !!         CALL wx_pbl_fuse_no_dts(knon, dtime, ypplay, ywake_s, &
    1792 !!                                 yt_x, yt_w, yq_x, yq_w, &
    1793 !!                                 yu_x, yu_w, yv_x, yv_w, &
    1794 !!                                 ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w, &
    1795 !!                                 AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
    1796 !!                                 AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
    1797 !!                                 BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
    1798 !!                                 BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
    1799 !!                                 AcoefH, AcoefQ, AcoefU, AcoefV, &
    1800 !!                                 BcoefH, BcoefQ, BcoefU, BcoefV, &
    1801 !!                                 ycdragh, ycdragm, &
    1802 !!                                 yt1, yq1, yu1, yv1 &
    1803 !!                                 )
    18041848       ELSE IF (iflag_split .ge. 1) THEN
    1805          CALL wx_pbl0_fuse(knon, dtime, ypplay, ywake_s, &
     1849!
     1850! Cdragq computation
     1851! ------------------
     1852    !******************************************************************************
     1853    ! Cdragq computed from cdrag
     1854    ! The difference comes only from a factor (f_z0qh_oce) on z0, so that
     1855    ! it can be computed inside wx_pbl0_merge
     1856    ! More complicated appraches may require the propagation through
     1857    ! pbl_surface of an independant cdragq variable.
     1858    !******************************************************************************
     1859!
     1860    IF ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce) THEN
     1861       ! Si on suit les formulations par exemple de Tessel, on
     1862       ! a z0h=0.4*nu/u*, z0q=0.62*nu/u*, d'ou f_z0qh_oce=0.62/0.4=1.55
     1863!!       ycdragq_x(1:knon)=ycdragh_x(1:knon)*                                      &
     1864!!            log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
     1865!!       ycdragq_w(1:knon)=ycdragh_w(1:knon)*                                      &
     1866!!            log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
     1867!
     1868       DO j = 1,knon
     1869         z1lay = zgeo1(j)/RG
     1870         fact_cdrag = log(z1lay/yz0h(j))/log(z1lay/(f_z0qh_oce*yz0h(j)))
     1871         ycdragq_x(j)=ycdragh_x(j)*fact_cdrag
     1872         ycdragq_w(j)=ycdragh_w(j)*fact_cdrag
     1873!!     Print *,'YYYYpbl0: fact_cdrag ', fact_cdrag
     1874       ENDDO  ! j = 1,knon
     1875!
     1876!!  Print *,'YYYYpbl0: z1lay, yz0h, f_z0qh_oce, ycdragh_w, ycdragq_w ', &
     1877!!                z1lay, yz0h(1:knon), f_z0qh_oce, ycdragh_w(1:knon), ycdragq_w(1:knon)
     1878    ELSE
     1879       ycdragq_x(1:knon)=ycdragh_x(1:knon)
     1880       ycdragq_w(1:knon)=ycdragh_w(1:knon)
     1881    ENDIF  ! ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce)
     1882!
     1883         CALL wx_pbl_prelim_0(knon, nsrf, dtime, ypplay, ypaprs, ywake_s,  &
     1884                         yts, y_delta_tsurf, ygustiness, &
    18061885                         yt_x, yt_w, yq_x, yq_w, &
    18071886                         yu_x, yu_w, yv_x, yv_w, &
    1808                          ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w, &
     1887                         ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
     1888                         ycdragm_x, ycdragm_w, &
     1889                         AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
     1890                         AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
     1891                         BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
     1892                         BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w  &
     1893                         )
     1894         CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta,  &
     1895                         BcoefQ_x, BcoefQ_w  &
     1896                         )
     1897         CALL wx_pbl0_merge(knon, ypplay, ypaprs,  &
     1898                         ywake_s, ydTs0, ydqs0, &
     1899                         yt_x, yt_w, yq_x, yq_w, &
     1900                         yu_x, yu_w, yv_x, yv_w, &
     1901                         ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
     1902                         ycdragm_x, ycdragm_w, &
    18091903                         AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
    18101904                         AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
    18111905                         BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
    18121906                         BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
    1813                          AcoefH, AcoefQ, AcoefU, AcoefV, &
    1814                          BcoefH, BcoefQ, BcoefU, BcoefV, &
    1815                          ycdragh, ycdragm, &
     1907                         AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &
     1908                         BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
     1909                         ycdragh, ycdragq, ycdragm, &
    18161910                         yt1, yq1, yu1, yv1 &
    18171911                         )
    1818 !!       ELSE IF (iflag_split .ge.2) THEN
    1819 !!!    Provisoire
    1820 !!         ah(:) = 0.
    1821 !!         bh(:) = 0.
    1822 !!         IF (nsrf == is_oce) THEN
    1823 !!           ybeta(:) = 1.
    1824 !!         ELSE
    1825 !!           ybeta(:) = beta_land
    1826 !!         ENDIF
    1827 !!         ycdragh(:) = ywake_s(:)*ycdragh_w(:) + (1.-ywake_s(:))*ycdragh_x(:)
    1828 !!         CALL wx_dts(knon, nsrf, ywake_cstar, ywake_s, ywake_dens, &
    1829 !!                     yts, ypplay(:,1), ybeta, ycdragh , ypaprs(:,1), &
    1830 !!                     yq(:,1), yt(:,1), yu(:,1), yv(:,1), ygustiness, &
    1831 !!                     ah, bh &
    1832 !!                     )
    1833 !!!
    1834 !!         CALL wx_pbl_fuse(knon, dtime, ypplay, ywake_s, &
    1835 !!                         yt_x, yt_w, yq_x, yq_w, &
    1836 !!                         yu_x, yu_w, yv_x, yv_w, &
    1837 !!                         ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w, &
    1838 !!                         AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
    1839 !!                         AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
    1840 !!                         BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
    1841 !!                         BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
    1842 !!                         ah, bh, &
    1843 !!                         AcoefH, AcoefQ, AcoefU, AcoefV, &
    1844 !!                         BcoefH, BcoefQ, BcoefU, BcoefV, &
    1845 !!                         ycdragh, ycdragm, &
    1846 !!                         yt1, yq1, yu1, yv1 &
    1847 !!                         )
    1848 !>jyg
    1849 !!!
    1850          ENDIF  ! (iflag_split .eq.0)
     1912         IF (iflag_split .eq. 2 .AND. nsrf .ne. is_oce) THEN
     1913           CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
     1914                           ywake_s, ybeta, ywake_cstar, ywake_dens, &
     1915                           AcoefH_x, AcoefH_w, &
     1916                           BcoefH_x, BcoefH_w, &
     1917                           AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
     1918                           AcoefH, AcoefQ, BcoefH, BcoefQ,  &
     1919                           HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
     1920                           phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
     1921                           yg_T, yg_Q, &
     1922                           yGamma_dTs_phiT, yGamma_dQs_phiQ, &
     1923                           ydTs_ins, ydqs_ins &
     1924                           )
     1925         ELSE !
     1926           AcoefH(:) = AcoefH_0(:)
     1927           AcoefQ(:) = AcoefQ_0(:)
     1928           BcoefH(:) = BcoefH_0(:)
     1929           BcoefQ(:) = BcoefQ_0(:)
     1930           yg_T(:) = 0.
     1931           yg_Q(:) = 0.
     1932           yGamma_dTs_phiT(:) = 0.
     1933           yGamma_dQs_phiQ(:) = 0.
     1934           ydTs_ins(:) = 0.
     1935           ydqs_ins(:) = 0.
     1936         ENDIF   ! (iflag_split .eq. 2)
     1937       ENDIF  ! (iflag_split .eq.0)
    18511938!!!
    18521939       IF (prt_level >=10) THEN
    1853          PRINT *,'pbl_surface (fuse->): yt(1,:) ',yt(1,:)
    1854          PRINT *,'pbl_surface (fuse->): yq(1,:) ',yq(1,:)
    1855          PRINT *,'pbl_surface (fuse->): yu(1,:) ',yu(1,:)
    1856          PRINT *,'pbl_surface (fuse->): yv(1,:) ',yv(1,:)
    1857          PRINT *,'pbl_surface (fuse->): AcoefH(1) ',AcoefH(1)
    1858          PRINT *,'pbl_surface (fuse->): BcoefH(1) ',BcoefH(1)
     1940         PRINT *,'pbl_surface (merge->): yt(1,:) ',yt(1,:)
     1941         PRINT *,'pbl_surface (merge->): yq(1,:) ',yq(1,:)
     1942         PRINT *,'pbl_surface (merge->): yu(1,:) ',yu(1,:)
     1943         PRINT *,'pbl_surface (merge->): yv(1,:) ',yv(1,:)
     1944         PRINT *,'pbl_surface (merge->): AcoefH(1), AcoefQ(1), AcoefU(1), AcoefV(1) ', &
     1945                                         AcoefH(1), AcoefQ(1), AcoefU(1), AcoefV(1)
     1946         PRINT *,'pbl_surface (merge->): BcoefH(1), BcoefQ(1), BcoefU(1), BcoefV(1) ', &
     1947                                         BcoefH(1), BcoefQ(1), BcoefU(1), BcoefV(1)
     1948
    18591949       ENDIF
    18601950
     1951!  Save initial value of z0h for use in evappot (z0h wiil be computed again in the surface models)
     1952          yz0h_old(1:knon) = yz0h(1:knon)
     1953!
    18611954!****************************************************************************************
    18621955!
    … …  
    21172210          y_flux_q1(j) = -yevap(j)
    21182211          ENDDO
    2119         ENDIF
    2120 
    2121        IF (prt_level >=10) THEN
    2122         DO j=1,knon
    2123          print*,'y_flux_t1,yfluxlat,wakes' &
    2124  &             ,  y_flux_t1(j), yfluxlat(j), ywake_s(j)
    2125          print*,'beta,ytsurf_new', ybeta(j), ytsurf_new(j)
    2126          print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j)
    2127         ENDDO
    2128        ENDIF
    2129 
    2130 !!! jyg le 07/02/2012 puis le 10/04/2013
    2131 !!       IF (iflag_split .eq.1) THEN
    2132 !!!!!
    2133 !!!jyg<
    2134 !!         CALL wx_pbl_split_no_dts(knon, ywake_s, &
    2135 !!                                AcoefH_x, AcoefH_w, &
    2136 !!                                AcoefQ_x, AcoefQ_w, &
    2137 !!                                AcoefU_x, AcoefU_w, &
    2138 !!                                AcoefV_x, AcoefV_w, &
    2139 !!                                y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &
    2140 !!                                y_flux_t1_x, y_flux_t1_w, &
    2141 !!                                y_flux_q1_x, y_flux_q1_w, &
    2142 !!                                y_flux_u1_x, y_flux_u1_w, &
    2143 !!                                y_flux_v1_x, y_flux_v1_w, &
    2144 !!                                yfluxlat_x, yfluxlat_w &
    2145 !!                                )
    2146 !!       ELSE IF (iflag_split .ge. 2) THEN
     2212        ENDIF ! (ok_flux_surf)
     2213!
     2214! ------------------------------------------------------------------------------
     2215! 12a)  Splitting
     2216! ------------------------------------------------------------------------------
     2217
    21472218       IF (iflag_split .GE. 1) THEN
    2148          CALL wx_pbl0_split(knon, dtime, ywake_s, &
     2219!
     2220         IF (nsrf .ne. is_oce) THEN
     2221!
     2222!         Compute potential evaporation and aridity factor  (jyg, 20200328)
     2223          ybeta_prev(:) = ybeta(:)
     2224             DO j = 1, knon
     2225               yqa(j) = AcoefQ(j) - BcoefQ(j)*yevap(j)*dtime
     2226             ENDDO
     2227!
     2228          CALL wx_evappot(knon, yqa, yTsurf_new, yevap_pot)
     2229!
     2230          ybeta(1:knon) = min(yevap(1:knon)/yevap_pot(1:knon), 1.)
     2231         
     2232          IF (prt_level >=10) THEN
     2233           DO j=1,knon
     2234            print*,'y_flux_t1,yfluxlat,wakes' &
     2235 &                ,  y_flux_t1(j), yfluxlat(j), ywake_s(j)
     2236            print*,'beta_prev, beta, ytsurf_new', ybeta_prev(j), ybeta(j), ytsurf_new(j)
     2237            print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j)
     2238           ENDDO
     2239          ENDIF  ! (prt_level >=10)
     2240!
     2241! Second call to wx_pbl0_merge and wx_pbl_dts_merge in order to take into account
     2242! the update of the aridity coeficient beta.
     2243!
     2244        CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta,  &
     2245                        BcoefQ_x, BcoefQ_w  &
     2246                        )
     2247        CALL wx_pbl0_merge(knon, ypplay, ypaprs,  &
     2248                          ywake_s, ydTs0, ydqs0, &
     2249                          yt_x, yt_w, yq_x, yq_w, &
     2250                          yu_x, yu_w, yv_x, yv_w, &
     2251                          ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
     2252                          ycdragm_x, ycdragm_w, &
     2253                          AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
     2254                          AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
     2255                          BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
     2256                          BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
     2257                          AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &
     2258                          BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
     2259                          ycdragh, ycdragq, ycdragm, &
     2260                          yt1, yq1, yu1, yv1 &
     2261                          )
     2262          IF (iflag_split .eq. 2) THEN
     2263            CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
     2264                            ywake_s, ybeta, ywake_cstar, ywake_dens, &
     2265                            AcoefH_x, AcoefH_w, &
     2266                            BcoefH_x, BcoefH_w, &
     2267                            AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
     2268                            AcoefH, AcoefQ, BcoefH, BcoefQ,  &
     2269                            HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
     2270                            phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
     2271                            yg_T, yg_Q, &
     2272                            yGamma_dTs_phiT, yGamma_dQs_phiQ, &
     2273                            ydTs_ins, ydqs_ins &
     2274                            )
     2275          ELSE !
     2276            AcoefH(:) = AcoefH_0(:)
     2277            AcoefQ(:) = AcoefQ_0(:)
     2278            BcoefH(:) = BcoefH_0(:)
     2279            BcoefQ(:) = BcoefQ_0(:)
     2280            yg_T(:) = 0.
     2281            yg_Q(:) = 0.
     2282            yGamma_dTs_phiT(:) = 0.
     2283            yGamma_dQs_phiQ(:) = 0.
     2284            ydTs_ins(:) = 0.
     2285            ydqs_ins(:) = 0.
     2286          ENDIF   ! (iflag_split .eq. 2)
     2287!
     2288        ELSE    ! (nsrf .ne. is_oce)
     2289          ybeta(1:knon) = 1.
     2290          yevap_pot(1:knon) = yevap(1:knon)
     2291          AcoefH(:) = AcoefH_0(:)
     2292          AcoefQ(:) = AcoefQ_0(:)
     2293          BcoefH(:) = BcoefH_0(:)
     2294          BcoefQ(:) = BcoefQ_0(:)
     2295          yg_T(:) = 0.
     2296          yg_Q(:) = 0.
     2297          yGamma_dTs_phiT(:) = 0.
     2298          yGamma_dQs_phiQ(:) = 0.
     2299          ydTs_ins(:) = 0.
     2300          ydqs_ins(:) = 0.
     2301        ENDIF   ! (nsrf .ne. is_oce)
     2302!
     2303        CALL wx_pbl_split(knon, nsrf, dtime, ywake_s, ybeta, iflag_split, &
     2304                       yg_T, yg_Q, &
     2305                       yGamma_dTs_phiT, yGamma_dQs_phiQ, &
     2306                       ydTs_ins, ydqs_ins, &
    21492307                       y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &
     2308!!!!                       HTRn_b, dd_HTRn, HTphiT_b, dd_HTphiT, &
     2309                       phiQ0_b, phiT0_b, &
    21502310                       y_flux_t1_x, y_flux_t1_w, &
    21512311                       y_flux_q1_x, y_flux_q1_w, &
    … …  
    21532313                       y_flux_v1_x, y_flux_v1_w, &
    21542314                       yfluxlat_x, yfluxlat_w, &
    2155                        y_delta_tsurf &
     2315                       y_delta_qsats, &
     2316                       y_delta_tsurf_new, y_delta_qsurf &
    21562317                       )
     2318!
     2319         CALL wx_pbl_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
     2320                       yTs, y_delta_tsurf,  &
     2321                       yqsurf, yTsurf_new,  &
     2322                       y_delta_tsurf_new, y_delta_qsats,  &
     2323                       AcoefH_x, AcoefH_w, &
     2324                       BcoefH_x, BcoefH_w, &
     2325                       AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
     2326                       AcoefH, AcoefQ, BcoefH, BcoefQ,  &
     2327                       y_flux_t1, y_flux_q1,  &
     2328                       y_flux_t1_x, y_flux_t1_w, &
     2329                       y_flux_q1_x, y_flux_q1_w)
     2330!
     2331         IF (nsrf .ne. is_oce) THEN
     2332           CALL wx_pbl_dts_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
     2333                         yTs, y_delta_tsurf,  &
     2334                         yqsurf, yTsurf_new,  &
     2335                         y_delta_qsats, y_delta_tsurf_new, y_delta_qsurf,  &
     2336                         AcoefH_x, AcoefH_w, &
     2337                         BcoefH_x, BcoefH_w, &
     2338                         AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
     2339                         AcoefH, AcoefQ, BcoefH, BcoefQ,  &
     2340                         HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
     2341                         phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
     2342                         yg_T, yg_Q, &
     2343                         yGamma_dTs_phiT, yGamma_dQs_phiQ, &
     2344                         ydTs_ins, ydqs_ins, &
     2345                         y_flux_t1, y_flux_q1,  &
     2346                         y_flux_t1_x, y_flux_t1_w, &
     2347                         y_flux_q1_x, y_flux_q1_w )
     2348         ENDIF   ! (nsrf .ne. is_oce)
     2349!
     2350       ELSE  ! (iflag_split .ge. 1)
     2351         ybeta(1:knon) = 1.
     2352         yevap_pot(1:knon) = yevap(1:knon)
    21572353       ENDIF  ! (iflag_split .ge. 1)
     2354!
     2355       IF (prt_level >= 10) THEN
     2356         print *,'pbl_surface, ybeta , yevap, yevap_pot ', &
     2357                               ybeta , yevap, yevap_pot
     2358       ENDIF  ! (prt_level >= 10)
     2359!
    21582360!>jyg
    21592361!
    … …  
    23182520       ENDIF  ! (iflag_split .eq.0)
    23192521!!!
    2320 
    2321         DO j = 1, knon
    2322           y_dflux_t(j) = y_dflux_t(j) * ypct(j)
    2323           y_dflux_q(j) = y_dflux_q(j) * ypct(j)
    2324         ENDDO
    2325 
     2522!!
     2523!!        DO j = 1, knon
     2524!!          y_dflux_t(j) = y_dflux_t(j) * ypct(j)
     2525!!          y_dflux_q(j) = y_dflux_q(j) * ypct(j)
     2526!!        ENDDO
     2527!!
    23262528!****************************************************************************************
    23272529! 13) Transform variables for output format :
    … …  
    24382640          i = ni(j)
    24392641          evap(i,nsrf) = - flux_q(i,1,nsrf)                  !jyg
     2642          beta(i,nsrf) = ybeta(j)                             !jyg
    24402643          d_ts(i,nsrf) = y_d_ts(j)
    24412644!albedo SB >>>
    … …  
    24532656          cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
    24542657          cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
    2455           dflux_t(i) = dflux_t(i) + y_dflux_t(j)
    2456           dflux_q(i) = dflux_q(i) + y_dflux_q(j)
     2658          dflux_t(i) = dflux_t(i) + y_dflux_t(j)*ypct(j)
     2659          dflux_q(i) = dflux_q(i) + y_dflux_q(j)*ypct(j)
    24572660       ENDDO
    24582661
    … …  
    24702673!!! nrlmd le 13/06/2011
    24712674!!jyg20170131          delta_tsurf(i,nsrf)=y_delta_tsurf(j)*ypct(j)
    2472           delta_tsurf(i,nsrf)=y_delta_tsurf(j)
     2675!!jyg20210118          delta_tsurf(i,nsrf)=y_delta_tsurf(j)
     2676          delta_tsurf(i,nsrf)=y_delta_tsurf_new(j)
     2677!
     2678          delta_qsurf(i) = delta_qsurf(i) + y_delta_qsurf(j)*ypct(j)
    24732679!
    24742680          cdragh_x(i) = cdragh_x(i) + ycdragh_x(j)*ypct(j)
    … …  
    26712877               * (ypaprs(j,1)-ypplay(j,1))
    26722878          tairsol(j) = yts(j) + y_d_ts(j)
    2673           tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)
     2879!!          tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)
     2880          tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf_new(j)
    26742881          qairsol(j) = yqsurf(j)
    26752882        ENDDO
    … …  
    32683475    IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
    32693476    IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
     3477    IF (ALLOCATED(ydTs0)) DEALLOCATE(ydTs0)
     3478    IF (ALLOCATED(ydqs0)) DEALLOCATE(ydqs0)
    32703479
    32713480!jyg<

  • LMDZ6/trunk/libf/phylmd/phyredem.F90

    r3856 r3888  
    1212  USE fonte_neige_mod,  ONLY : fonte_neige_final
    1313  USE pbl_surface_mod,  ONLY : pbl_surface_final
    14   USE phys_state_var_mod, ONLY: radpas, zmasq, pctsrf, ftsol, falb_dir,      &
     14  USE phys_state_var_mod, ONLY: radpas, zmasq, pctsrf,                       &
     15                                ftsol, beta_aridity, delta_tsurf, falb_dir,  &
    1516                                falb_dif, qsol, fevap, radsol, solsw, sollw, &
    1617                                sollwdown, rain_fall, snow_fall, z0m, z0h,   &
    … …  
    157158    END IF
    158159
     160!    Surface variables
    159161    CALL put_field_srf1(pass,"TS","Temperature",ftsol(:,:))
     162
     163    CALL put_field_srf1(pass,"DELTA_TS","w-x surface temperature difference", delta_tsurf(:,:))
     164
     165    CALL put_field_srf1(pass,"BETA_S","Aridity factor", beta_aridity(:,:))
     166!    End surface variables
    160167
    161168! ================== Albedo =======================================

  • LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90

    r3817 r3888  
    340340      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxfluxlat_x, zxfluxlat_w
    341341!$OMP THREADPRIVATE(zxfluxlat_x, zxfluxlat_w)
     342      REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: delta_qsurf
     343!$OMP THREADPRIVATE(delta_qsurf)
    342344!jyg<
    343345!!! Entrees supplementaires couche-limite
    … …  
    733735      ALLOCATE(sens_x(klon), sens_w(klon))
    734736      ALLOCATE(zxfluxlat_x(klon), zxfluxlat_w(klon))
     737      ALLOCATE(delta_qsurf(klon))
    735738!jyg<
    736739!!      ALLOCATE(t_x(klon,klev), t_w(klon,klev))
    … …  
    10321035      DEALLOCATE(sens_x, sens_w)
    10331036      DEALLOCATE(zxfluxlat_x, zxfluxlat_w)
     1037      DEALLOCATE(delta_qsurf)
    10341038!jyg<
    10351039!!      DEALLOCATE(t_x, t_w)

  • LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90

    r3873 r3888  
    814814'flat_w', 'flat within_wake', 'W/m2', (/ ('', i=1, 10) /))
    815815!!
    816   type(ctrl_out),save :: o_delta_tsurf    = ctrl_out((/ 1, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
    817 'delta_tsurf', 'Temperature difference (w-x)', 'K', (/ ('', i=1, 10) /))
    818816  type(ctrl_out),save :: o_cdragh_x       = ctrl_out((/ 1, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
    819817'cdragh_x', 'cdragh off-wake', '', (/ ('', i=1, 10) /))
    … …  
    10941092      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'dltpbltke_sic',       &
    10951093      "TKE difference (w - x) "//clnsurf(4),"-", (/ ('', i=1, 10) /)) /)
     1094
     1095  TYPE(ctrl_out), SAVE :: o_delta_tsurf = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     1096    'delta_tsurf ', 'T_surf difference (w - x)', 'K', (/ ('', i=1, 10) /))
     1097  TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_delta_tsurf_srf      = (/             &
     1098      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'delta_tsurf_ter',       &
     1099      "T_surf difference (w - x) "//clnsurf(1),"-", (/ ('', i=1, 10) /)), &
     1100      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'delta_tsurf_lic',       &
     1101      "T_surf difference (w - x) "//clnsurf(2),"-", (/ ('', i=1, 10) /)), &
     1102      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'delta_tsurf_oce',       &
     1103      "T_surf difference (w - x) "//clnsurf(3),"-", (/ ('', i=1, 10) /)), &
     1104      ctrl_out((/ 10, 4, 10, 10, 10, 10, 11, 11, 11, 11/),'delta_tsurf_sic',       &
     1105      "T_surf difference (w - x) "//clnsurf(4),"-", (/ ('', i=1, 10) /)) /)
    10961106
    10971107  TYPE(ctrl_out), SAVE :: o_kz = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &

  • LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90

    r3865 r3888  
    8787         o_dtvdf_x    , o_dtvdf_w    , o_dqvdf_x    , o_dqvdf_w    , &
    8888         o_sens_x     , o_sens_w     , o_flat_x     , o_flat_w     , &
    89          o_delta_tsurf, &
     89         o_delta_tsurf, o_delta_tsurf_srf, &
    9090         o_cdragh_x   , o_cdragh_w   , o_cdragm_x   , o_cdragm_w   , &
    9191         o_kh         , o_kh_x       , o_kh_w       , &
    … …  
    13541354!
    13551355               CALL histwrite_phy(o_dqvdf_w    ,zx_tmp_fi3d)
    1356                CALL histwrite_phy(o_sens_x     ,sens_x     )
    1357                CALL histwrite_phy(o_sens_w     ,sens_w     )
     1356       IF (vars_defined)  zx_tmp_fi2d(1:klon)=-1*sens_x(1:klon)
     1357               CALL histwrite_phy(o_sens_x     ,zx_tmp_fi2d)
     1358       IF (vars_defined)  zx_tmp_fi2d(1:klon)=-1*sens_w(1:klon)
     1359               CALL histwrite_phy(o_sens_w     ,zx_tmp_fi2d)
    13581360               CALL histwrite_phy(o_flat_x     ,zxfluxlat_x)
    13591361               CALL histwrite_phy(o_flat_w     ,zxfluxlat_w)
    1360                CALL histwrite_phy(o_delta_tsurf,delta_tsurf)
     1362          zx_tmp_fi2d=0.
     1363          IF (vars_defined) THEN
     1364             DO nsrf=1,nbsrf
     1365                   zx_tmp_fi2d(:)=zx_tmp_fi2d(:) &
     1366                        +pctsrf(:,nsrf)*delta_tsurf(:,nsrf)
     1367             ENDDO
     1368          ENDIF
     1369               CALL histwrite_phy(o_delta_tsurf,zx_tmp_fi2d)
    13611370               CALL histwrite_phy(o_cdragh_x   ,cdragh_x   )
    13621371               CALL histwrite_phy(o_cdragh_w   ,cdragh_w   )

  • LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90

    r3856 r3888  
    3232      REAL, ALLOCATABLE, SAVE :: ftsol(:,:)
    3333!$OMP THREADPRIVATE(ftsol)
     34      REAL, ALLOCATABLE, SAVE :: beta_aridity(:,:)
     35!$OMP THREADPRIVATE(beta_aridity)
    3436      REAL,ALLOCATABLE,SAVE :: qsol(:),fevap(:,:),z0m(:,:),z0h(:,:),agesno(:,:)
    3537!$OMP THREADPRIVATE(qsol,fevap,z0m,z0h,agesno)
    … …  
    9698      REAL, ALLOCATABLE, SAVE :: coefm(:,:,:) ! Kz momentum
    9799!$OMP THREADPRIVATE(pbl_tke, coefh,coefm)
    98 !nrlmd<
    99       REAL, ALLOCATABLE, SAVE :: delta_tsurf(:,:) ! Surface temperature difference inside-outside cold pool
    100 !$OMP THREADPRIVATE(delta_tsurf)
    101 !>nrlmd
    102100      REAL, ALLOCATABLE, SAVE :: zmax0(:), f0(:) !
    103101!$OMP THREADPRIVATE(zmax0,f0)
    … …  
    276274      REAL,ALLOCATABLE,SAVE :: wake_delta_pbl_TKE(:,:,:)
    277275!$OMP THREADPRIVATE(wake_delta_pbl_TKE)
     276!nrlmd<
     277      REAL, ALLOCATABLE, SAVE :: delta_tsurf(:,:) ! Surface temperature difference inside-outside cold pool
     278!$OMP THREADPRIVATE(delta_tsurf)
     279!>nrlmd
    278280!>jyg
    279281!
    … …  
    478480      ALLOCATE(pctsrf(klon,nbsrf))
    479481      ALLOCATE(ftsol(klon,nbsrf))
     482      ALLOCATE(beta_aridity(klon,nbsrf))
    480483      ALLOCATE(qsol(klon),fevap(klon,nbsrf))
    481484      ALLOCATE(z0m(klon,nbsrf+1),z0h(klon,nbsrf+1),agesno(klon,nbsrf))
    … …  
    674677
    675678      DEALLOCATE(pctsrf, ftsol, falb1, falb2)
     679      DEALLOCATE(beta_aridity)
    676680      DEALLOCATE(qsol,fevap,z0m,z0h,agesno)
    677681!FC
    … …  
    687691      DEALLOCATE(tr_ancien)                           !RomP
    688692      DEALLOCATE(ratqs, pbl_tke,coefh,coefm)
    689 !nrlmd<
    690       DEALLOCATE(delta_tsurf)
    691 !>nrlmd
    692693      DEALLOCATE(zmax0, f0)
    693694      DEALLOCATE(sig1, w01)
    … …  
    744745!jyg<
    745746      DEALLOCATE(wake_delta_pbl_TKE)
     747!nrlmd<
     748      DEALLOCATE(delta_tsurf)
     749!>nrlmd
    746750!>jyg
    747751      DEALLOCATE(pfrac_impa, pfrac_nucl)

  • LMDZ6/trunk/libf/phylmd/physiq_mod.F90

    r3877 r3888  
    210210       zxrunofflic,                            &
    211211       zxtsol, snow_lsc, zxfqfonte, zxqsurf,   &
     212       delta_qsurf,                            &
    212213       rain_lsc, rain_num,                     &
    213214       !
    … …  
    24852486    !   s_therm,   s_trmb1,   s_trmb2, s_trmb3,
    24862487    !   zu10m,     zv10m,   fder,
    2487     !   zxqsurf,   rh2m,      zxfluxu, zxfluxv,
     2488    !   zxqsurf,   delta_qsurf,
     2489    !   rh2m,      zxfluxu, zxfluxv,
    24882490    !   frugs,     agesno,    fsollw,  fsolsw,
    24892491    !   d_ts,      fevap,     fluxlat, t2m,
    … …  
    25462548                                !albedo SB <<<
    25472549            cdragh,    cdragm,  u1,    v1,            &
     2550            beta_aridity, &
    25482551                                !albedo SB >>>
    25492552                                ! albsol1,   albsol2,   sens,    evap,      &
    … …  
    25722575            s_therm,   s_trmb1,   s_trmb2, s_trmb3, &
    25732576            zustar, zu10m,     zv10m,   fder, &
    2574             zxqsurf,   rh2m,      zxfluxu, zxfluxv, &
     2577            zxqsurf, delta_qsurf,   rh2m,      zxfluxu, zxfluxv, &
    25752578            z0m, z0h,     agesno,    fsollw,  fsolsw, &
    25762579            d_ts,      fevap,     fluxlat, t2m, &
    … …  
    46374640
    46384641    CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
    4639 
     4642   !
     4643   ! Prevent pbl_tke_w from becoming negative
     4644    wake_delta_pbl_tke(:,:,:) = max(wake_delta_pbl_tke(:,:,:), -pbl_tke(:,:,:))
     4645   !
    46404646
    46414647       ENDIF

  • LMDZ6/trunk/libf/phylmd/wx_pbl_mod.F90

    r3181 r3888  
    11MODULE wx_pbl_mod
    22!
    3 ! Planetary Boundary Layer and Surface module
    4 !
    5 ! This module manage the calculation of turbulent diffusion in the boundary layer
    6 ! and all interactions towards the differents sub-surfaces.
    7 !
     3! Split Planetary Boundary Layer
     4!
     5! This module manages the splitting of the boundary layer between two regions; the (w)
     6! region (inside cold pools) and the (x) region (outside cold pools)
    87!
    98  USE dimphy
    … …  
    1110  IMPLICIT NONE
    1211
    13   REAL, ALLOCATABLE, DIMENSION(:), SAVE        :: Kech_Tp, Kech_T_xp, Kech_T_wp
    14   REAL, ALLOCATABLE, DIMENSION(:), SAVE        :: dd_KTp, KxKwTp, dd_AT, dd_BT
    15 !$OMP THREADPRIVATE(Kech_Tp, Kech_T_xp, Kech_T_wp, dd_KTp, KxKwTp, dd_AT, dd_BT)
    16   REAL, ALLOCATABLE, DIMENSION(:), SAVE        :: Kech_Qp, Kech_Q_xp, Kech_Q_wp
    17   REAL, ALLOCATABLE, DIMENSION(:), SAVE        :: dd_KQp, KxKwQp, dd_AQ, dd_BQ
    18 !$OMP THREADPRIVATE(Kech_Qp, Kech_Q_xp, Kech_Q_wp, dd_KQp, KxKwQp, dd_AQ, dd_BQ)
    19   REAL, ALLOCATABLE, DIMENSION(:), SAVE        :: Kech_Up, Kech_U_xp, Kech_U_wp
    20   REAL, ALLOCATABLE, DIMENSION(:), SAVE        :: dd_KUp, KxKwUp, dd_AU, dd_BU
    21 !$OMP THREADPRIVATE(Kech_Up, Kech_U_xp, Kech_U_wp, dd_KUp, KxKwUp, dd_AU, dd_BU)
    22   REAL, ALLOCATABLE, DIMENSION(:), SAVE        :: Kech_Vp, Kech_V_xp, Kech_V_wp
    23   REAL, ALLOCATABLE, DIMENSION(:), SAVE        :: dd_KVp, KxKwVp, dd_AV, dd_BV
    24 !$OMP THREADPRIVATE(Kech_Vp, Kech_V_xp, Kech_V_wp, dd_KVp, KxKwVp, dd_AV, dd_BV)
    25 
    2612CONTAINS
    2713!
    2814!****************************************************************************************
    2915!
    30 SUBROUTINE wx_pbl_init
    31 
    32 ! Local variables
    33 !****************************************************************************************
    34     INTEGER                       :: ierr
    35  
    36 
    37 !****************************************************************************************
    38 ! Allocate module variables
    39 !
    40 !****************************************************************************************   
    41 
    42     ierr = 0
    43 
    44     ALLOCATE(Kech_Tp(klon), stat=ierr)
    45     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    46 
    47     ALLOCATE(Kech_T_xp(klon), stat=ierr)
    48     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    49 
    50     ALLOCATE(Kech_T_wp(klon), stat=ierr)
    51     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    52 
    53     ALLOCATE(dd_KTp(klon), stat=ierr)
    54     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    55 
    56     ALLOCATE(KxKwTp(klon), stat=ierr)
    57     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    58 
    59     ALLOCATE(dd_AT(klon), stat=ierr)
    60     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    61 
    62     ALLOCATE(dd_BT(klon), stat=ierr)
    63     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    64 
    65 !----------------------------------------------------------------------------
    66     ALLOCATE(Kech_Qp(klon), stat=ierr)
    67     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    68 
    69     ALLOCATE(Kech_Q_xp(klon), stat=ierr)
    70     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    71 
    72     ALLOCATE(Kech_Q_wp(klon), stat=ierr)
    73     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    74 
    75     ALLOCATE(dd_KQp(klon), stat=ierr)
    76     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    77 
    78     ALLOCATE(KxKwQp(klon), stat=ierr)
    79     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    80 
    81     ALLOCATE(dd_AQ(klon), stat=ierr)
    82     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    83 
    84     ALLOCATE(dd_BQ(klon), stat=ierr)
    85     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    86 
    87 !----------------------------------------------------------------------------
    88     ALLOCATE(Kech_Up(klon), stat=ierr)
    89     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    90 
    91     ALLOCATE(Kech_U_xp(klon), stat=ierr)
    92     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    93 
    94     ALLOCATE(Kech_U_wp(klon), stat=ierr)
    95     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    96 
    97     ALLOCATE(dd_KUp(klon), stat=ierr)
    98     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    99 
    100     ALLOCATE(KxKwUp(klon), stat=ierr)
    101     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    102 
    103     ALLOCATE(dd_AU(klon), stat=ierr)
    104     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    105 
    106     ALLOCATE(dd_BU(klon), stat=ierr)
    107     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    108 
    109 !----------------------------------------------------------------------------
    110     ALLOCATE(Kech_Vp(klon), stat=ierr)
    111     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    112 
    113     ALLOCATE(Kech_V_xp(klon), stat=ierr)
    114     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    115 
    116     ALLOCATE(Kech_V_wp(klon), stat=ierr)
    117     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    118 
    119     ALLOCATE(dd_KVp(klon), stat=ierr)
    120     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    121 
    122     ALLOCATE(KxKwVp(klon), stat=ierr)
    123     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    124 
    125     ALLOCATE(dd_AV(klon), stat=ierr)
    126     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    127 
    128     ALLOCATE(dd_BV(klon), stat=ierr)
    129     IF (ierr /= 0) CALL abort_physic('wx_pbl_init', 'pb in allocation',1)
    130 
    131 !----------------------------------------------------------------------------
    132 
    133 END SUBROUTINE wx_pbl_init
    134 
    135 SUBROUTINE wx_pbl0_fuse(knon, dtime, ypplay, ywake_s, &
     16SUBROUTINE wx_pbl0_merge(knon, ypplay, ypaprs,  &
     17                                 sigw, dTs_forcing, dqs_forcing,  &
    13618                                 yt_x, yt_w, yq_x, yq_w, &
    13719                                 yu_x, yu_w, yv_x, yv_w, &
    138                                  ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w, &
     20                                 ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
     21                                 ycdragm_x, ycdragm_w, &
    13922                                 AcoefT_x, AcoefT_w, AcoefQ_x, AcoefQ_w, &
    14023                                 AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
    … …  
    14326                                 AcoefT, AcoefQ, AcoefU, AcoefV, &
    14427                                 BcoefT, BcoefQ, BcoefU, BcoefV, &
    145                                  ycdragh, ycdragm, &
     28                                 ycdragh, ycdragq, ycdragm, &
    14629                                 yt1, yq1, yu1, yv1 &
    14730                                 )
    14831!
     32
     33    USE wx_pbl_var_mod
     34
    14935    USE print_control_mod, ONLY: prt_level,lunout
     36    USE indice_sol_mod, ONLY: is_oce
    15037!
    15138    INCLUDE "YOMCST.h"
     39    INCLUDE "FCTTRE.h"
     40    INCLUDE "YOETHF.h"
     41    INCLUDE "clesphys.h"
    15242!
    15343    INTEGER,                      INTENT(IN)        :: knon    ! number of grid cells
    154     REAL,                         INTENT(IN)        :: dtime   ! time step size (s)
    15544    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypplay  ! mid-layer pressure (Pa)
    156     REAL, DIMENSION(knon),        INTENT(IN)        :: ywake_s ! cold pools fractional area
     45    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypaprs  ! pressure at layer interfaces (pa)
     46    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw ! cold pools fractional area
     47    REAL, DIMENSION(knon),        INTENT(IN)        :: dTs_forcing ! forced temperature difference (w)-(x)
     48    REAL, DIMENSION(knon),        INTENT(IN)        :: dqs_forcing ! forced humidity difference (w)-(x)
    15749    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: yt_x, yt_w, yq_x, yq_w
    15850    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: yu_x, yu_w, yv_x, yv_w
    159     REAL, DIMENSION(knon),        INTENT(IN)        :: ycdragh_x, ycdragh_w, ycdragm_x, ycdragm_w
     51    REAL, DIMENSION(knon),        INTENT(IN)        :: ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w
     52    REAL, DIMENSION(knon),        INTENT(IN)        :: ycdragm_x, ycdragm_w
    16053    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT_x, AcoefT_w, AcoefQ_x, AcoefQ_w
    16154    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w
    … …  
    16457    REAL, DIMENSION(knon),        INTENT(OUT)       :: AcoefT, AcoefQ, AcoefU, AcoefV
    16558    REAL, DIMENSION(knon),        INTENT(OUT)       :: BcoefT, BcoefQ, BcoefU, BcoefV
    166     REAL, DIMENSION(knon),        INTENT(OUT)       :: ycdragh, ycdragm
     59    REAL, DIMENSION(knon),        INTENT(OUT)       :: ycdragh, ycdragq, ycdragm
    16760    REAL, DIMENSION(knon),        INTENT(OUT)       :: yt1, yq1, yu1, yv1  ! Apparent T, q, u, v at first level, as
    16861                                                                           !seen by surface modules
    … …  
    17063! Local variables
    17164    INTEGER                    :: j
    172     REAL                       :: rho1
    173     REAL                       :: mod_wind_x
    174     REAL                       :: mod_wind_w   
    175     REAL                       :: dd_Cdragh
    176     REAL                       :: dd_Cdragm
    17765    REAL                       :: dd_Kh
     66    REAL                       :: dd_Kq
    17867    REAL                       :: dd_Km
    17968    REAL                       :: dd_u
    … …  
    18271    REAL                       :: dd_q
    18372!
    184     REAL                       :: KCT, KCQ, KCU, KCV
    185 !
    186     REAL                       :: BBT, BBQ, BBU, BBV
    187     REAL                       :: DDT, DDQ, DDU, DDV
    188     REAL                       :: LambdaT, LambdaQ, LambdaU, LambdaV
    18973    REAL                       :: LambdaTs, LambdaQs, LambdaUs, LambdaVs
    19074!
    19175    REAL, DIMENSION(knon)      :: sigx       ! fractional area of (x) region
    192 
    193     REAL, DIMENSION(knon)      :: Kech_h    ! Energy exchange coefficient
    194     REAL, DIMENSION(knon)      :: Kech_h_x, Kech_h_w
    195     REAL, DIMENSION(knon)      :: Kech_m    ! Momentum exchange coefficient
    196     REAL, DIMENSION(knon)      :: Kech_m_x, Kech_m_w
    197 
    198 !!!
    199 !!! jyg le 09/04/2013 ; passage aux nouvelles expressions en differences
    200 
    201         sigx(:) = 1.-ywake_s(:)
    202 
     76!
     77!
     78   sigx(1:knon) = 1.-sigw(1:knon)
     79!                                           
     80!
    20381        DO j=1,knon
    20482!
    205 ! Calcul des coefficients d echange
    206          mod_wind_x = 1.0+SQRT(yu_x(j,1)**2+yv_x(j,1)**2)
    207          mod_wind_w = 1.0+SQRT(yu_w(j,1)**2+yv_w(j,1)**2)
    208 !!         rho1 = ypplay(j,1)/(RD*yt(j,1))
    209          rho1 = ypplay(j,1)/(RD*(yt_x(j,1) + ywake_s(j)*(yt_w(j,1)-yt_x(j,1))))
    210          Kech_h_x(j) = ycdragh_x(j) * mod_wind_x * rho1
    211          Kech_h_w(j) = ycdragh_w(j) * mod_wind_w * rho1
    212          Kech_m_x(j) = ycdragm_x(j) * mod_wind_x * rho1
    213          Kech_m_w(j) = ycdragm_w(j) * mod_wind_w * rho1
    214 !
    215          dd_Kh = Kech_h_w(j) - Kech_h_x(j)
    216          dd_Km = Kech_m_w(j) - Kech_m_x(j)
    217          IF (prt_level >=10) THEN
    218           print *,' mod_wind_x, mod_wind_w ', mod_wind_x, mod_wind_w
    219           print *,' rho1 ',rho1
    220           print *,' ycdragh_x(j),ycdragm_x(j) ',ycdragh_x(j),ycdragm_x(j)
    221           print *,' ycdragh_w(j),ycdragm_w(j) ',ycdragh_w(j),ycdragm_w(j)
    222           print *,' dd_Kh: ',dd_Kh
    223          ENDIF
    224 !
    225          Kech_h(j) = Kech_h_x(j) + ywake_s(j)*dd_Kh
    226          Kech_m(j) = Kech_m_x(j) + ywake_s(j)*dd_Km
    227 !
    228 ! Calcul des coefficients d echange corriges des retroactions
    229         Kech_T_xp(j) = Kech_h_x(j)/(1.-BcoefT_x(j)*Kech_h_x(j)*dtime)
    230         Kech_T_wp(j) = Kech_h_w(j)/(1.-BcoefT_w(j)*Kech_h_w(j)*dtime)
    231         Kech_Q_xp(j) = Kech_h_x(j)/(1.-BcoefQ_x(j)*Kech_h_x(j)*dtime)
    232         Kech_Q_wp(j) = Kech_h_w(j)/(1.-BcoefQ_w(j)*Kech_h_w(j)*dtime)
    233         Kech_U_xp(j) = Kech_m_x(j)/(1.-BcoefU_x(j)*Kech_m_x(j)*dtime)
    234         Kech_U_wp(j) = Kech_m_w(j)/(1.-BcoefU_w(j)*Kech_m_w(j)*dtime)
    235         Kech_V_xp(j) = Kech_m_x(j)/(1.-BcoefV_x(j)*Kech_m_x(j)*dtime)
    236         Kech_V_wp(j) = Kech_m_w(j)/(1.-BcoefV_w(j)*Kech_m_w(j)*dtime)
    237 !
    238          dd_KTp(j) = Kech_T_wp(j) - Kech_T_xp(j)
    239          dd_KQp(j) = Kech_Q_wp(j) - Kech_Q_xp(j)
    240          dd_KUp(j) = Kech_U_wp(j) - Kech_U_xp(j)
    241          dd_KVp(j) = Kech_V_wp(j) - Kech_V_xp(j)
    242 !
    243         Kech_Tp(j) = Kech_T_xp(j) + ywake_s(j)*dd_KTp(j)
    244         Kech_Qp(j) = Kech_Q_xp(j) + ywake_s(j)*dd_KQp(j)
    245         Kech_Up(j) = Kech_U_xp(j) + ywake_s(j)*dd_KUp(j)
    246         Kech_Vp(j) = Kech_V_xp(j) + ywake_s(j)*dd_KVp(j)
    247 !
    248 ! Calcul des differences w-x
    249        dd_Cdragm = ycdragm_w(j) - ycdragm_x(j)
    250        dd_Cdragh = ycdragh_w(j) - ycdragh_x(j)
     83!
     84! Compute w-x differences
     85       dd_t = yt_w(j,1) - yt_x(j,1)
     86       dd_q = yq_w(j,1) - yq_x(j,1)
    25187       dd_u = yu_w(j,1) - yu_x(j,1)
    25288       dd_v = yv_w(j,1) - yv_x(j,1)
    253        dd_t = yt_w(j,1) - yt_x(j,1)
    254        dd_q = yq_w(j,1) - yq_x(j,1)
    255        dd_AT(j) = AcoefT_w(j) - AcoefT_x(j)
    256        dd_AQ(j) = AcoefQ_w(j) - AcoefQ_x(j)
    257        dd_AU(j) = AcoefU_w(j) - AcoefU_x(j)
    258        dd_AV(j) = AcoefV_w(j) - AcoefV_x(j)
    259        dd_BT(j) = BcoefT_w(j) - BcoefT_x(j)
    260        dd_BQ(j) = BcoefQ_w(j) - BcoefQ_x(j)
    261        dd_BU(j) = BcoefU_w(j) - BcoefU_x(j)
    262        dd_BV(j) = BcoefV_w(j) - BcoefV_x(j)
    263 !
    264        KxKwTp(j) = Kech_T_xp(j)*Kech_T_wp(j)
    265        KxKwQp(j) = Kech_Q_xp(j)*Kech_Q_wp(j)
    266        KxKwUp(j) = Kech_U_xp(j)*Kech_U_wp(j)
    267        KxKwVp(j) = Kech_V_xp(j)*Kech_V_wp(j)
    268        BBT = (BcoefT_x(j) + sigx(j)*dd_BT(j))*dtime
    269        BBQ = (BcoefQ_x(j) + sigx(j)*dd_BQ(j))*dtime
    270        BBU = (BcoefU_x(j) + sigx(j)*dd_BU(j))*dtime
    271        BBV = (BcoefV_x(j) + sigx(j)*dd_BV(j))*dtime
    272        KCT = Kech_h(j)
    273        KCQ = Kech_h(j)
    274        KCU = Kech_m(j)
    275        KCV = Kech_m(j)
    276        DDT = Kech_Tp(j)
    277        DDQ = Kech_Qp(j)
    278        DDU = Kech_Up(j)
    279        DDV = Kech_Vp(j)
    280        LambdaT = dd_Kh/KCT
    281        LambdaQ = dd_Kh/KCQ
    282        LambdaU = dd_Km/KCU
    283        LambdaV = dd_Km/KCV
    284        LambdaTs = dd_KTp(j)/DDT
    285        LambdaQs = dd_KQp(j)/DDQ
    286        LambdaUs = dd_KUp(j)/DDU
    287        LambdaVs = dd_KVp(j)/DDV
    288 !
    289        IF (prt_level >=10) THEN
    290           print *,'Variables pour la fusion : Kech_T_xp(j)' ,Kech_T_xp(j)
    291           print *,'Variables pour la fusion : Kech_T_wp(j)' ,Kech_T_wp(j)
    292           print *,'Variables pour la fusion : Kech_Tp(j)' ,Kech_Tp(j)
    293           print *,'Variables pour la fusion : Kech_h(j)' ,Kech_h(j)
    294        ENDIF
     89!
     90! Merged exchange coefficients
     91         dd_Kh = Kech_h_w(j) - Kech_h_x(j)
     92         dd_Kq = Kech_q_w(j) - Kech_q_x(j)
     93         dd_Km = Kech_m_w(j) - Kech_m_x(j)
     94!
     95       LambdaTs = dd_KTp(j)/Kech_Tp(j)
     96       LambdaQs = dd_KQs(j)/Kech_Qs(j)
     97       LambdaUs = dd_KUp(j)/Kech_Up(j)
     98       LambdaVs = dd_KVp(j)/Kech_Vp(j)
    29599!
    296100! Calcul des coef A, B \'equivalents dans la couche 1
    297101!
    298        AcoefT(j) = AcoefT_x(j) + ywake_s(j)*dd_AT(j)*(1.+sigx(j)*LambdaTs)
    299        AcoefQ(j) = AcoefQ_x(j) + ywake_s(j)*dd_AQ(j)*(1.+sigx(j)*LambdaQs)
    300        AcoefU(j) = AcoefU_x(j) + ywake_s(j)*dd_AU(j)*(1.+sigx(j)*LambdaUs)
    301        AcoefV(j) = AcoefV_x(j) + ywake_s(j)*dd_AV(j)*(1.+sigx(j)*LambdaVs)
     102! The dTs_forcing and dqs_forcing terms are added for diagnostic purpose ; they should be zero in normal operation.
     103       AcoefT(j) = AcoefT_x(j) + sigw(j)*(1.+sigx(j)*LambdaTs)*(dd_AT(j) - C_p(j)*dTs_forcing(j))
     104       AcoefQ(j) = AcoefQ_x(j) + sigw(j)*(1.+sigx(j)*LambdaQs)*(dd_AQ(j) - dqs_forcing(j))
     105       AcoefU(j) = AcoefU_x(j) + sigw(j)*(1.+sigx(j)*LambdaUs)*dd_AU(j)
     106       AcoefV(j) = AcoefV_x(j) + sigw(j)*(1.+sigx(j)*LambdaVs)*dd_AV(j)
    302107!                                           
    303        BcoefT(j) = BcoefT_x(j) + ywake_s(j)*BcoefT_x(j)*sigx(j)*LambdaT*LambdaTs &
    304                                + ywake_s(j)*dd_BT(j)*(1.+sigx(j)*LambdaT)*(1.+sigx(j)*LambdaTs)
    305                                            
    306        BcoefQ(j) = BcoefQ_x(j) + ywake_s(j)*BcoefQ_x(j)*sigx(j)*LambdaQ*LambdaQs &
    307                                + ywake_s(j)*dd_BQ(j)*(1.+sigx(j)*LambdaQ)*(1.+sigx(j)*LambdaQs)
    308                                            
    309        BcoefU(j) = BcoefU_x(j) + ywake_s(j)*BcoefU_x(j)*sigx(j)*LambdaU*LambdaUs &
    310                                + ywake_s(j)*dd_BU(j)*(1.+sigx(j)*LambdaU)*(1.+sigx(j)*LambdaUs)
    311                                            
    312        BcoefV(j) = BcoefV_x(j) + ywake_s(j)*BcoefV_x(j)*sigx(j)*LambdaV*LambdaVs &
    313                                + ywake_s(j)*dd_BV(j)*(1.+sigx(j)*LambdaV)*(1.+sigx(j)*LambdaVs)
    314 
     108!
     109!!       BcoefT(j) = (sigw(j)*Kech_h_w(j)*Kech_T_pw(j)*BcoefT_w(j) + &
     110!!                sigx(j)*Kech_h_x(j)*Kech_T_px(j)*BcoefT_x(j) )/(Kech_h(j)*Kech_Tp(j))
     111!!       BcoefQ(j) = (sigw(j)*Kech_q_w(j)*Kech_Q_pw(j)*BcoefQ_w(j) + &
     112!!                sigx(j)*Kech_q_x(j)*Kech_Q_px(j)*BcoefQ_x(j) )/(Kech_q(j)*Kech_Qp(j))
     113!!       BcoefU(j) = (sigw(j)*Kech_m_w(j)*Kech_U_pw(j)*BcoefU_w(j) + &
     114!!                sigx(j)*Kech_m_x(j)*Kech_U_px(j)*BcoefU_x(j) )/(Kech_m(j)*Kech_Up(j))
     115!!       BcoefV(j) = (sigw(j)*Kech_m_w(j)*Kech_V_pw(j)*BcoefV_w(j) + &
     116!!                sigx(j)*Kech_m_x(j)*Kech_V_px(j)*BcoefV_x(j) )/(Kech_m(j)*Kech_Vp(j))
     117!
     118!!  Print *,'YYYYpbl0: BcoefT_x, sigw, sigx, dd_Kh, dd_KTp, Kech_h_w ', &
     119!!                     BcoefT_x, sigw, sigx, dd_Kh, dd_KTp, Kech_h_w
     120!!  Print *,'YYYYpbl0: Kech_T_pw, dd_BT, Kech_h, Kech_Tp ', &
     121!!                     Kech_T_pw, dd_BT, Kech_h, Kech_Tp
     122       BcoefT(j) = BcoefT_x(j) + sigw(j)*(sigx(j)*dd_Kh*dd_KTp(j)*BcoefT_x(j) + &
     123                                  Kech_h_w(j)*Kech_T_pw(j)*dd_BT(j))/(Kech_h(j)*Kech_Tp(j))
     124       BcoefQ(j) = BcoefQ_x(j) + sigw(j)*(sigx(j)*dd_Kh*dd_KQs(j)*BcoefQ_x(j) + &
     125                                  Kech_q_w(j)*Kech_Q_sw(j)*dd_BQ(j))/(Kech_q(j)*Kech_Qs(j))
     126       BcoefU(j) = BcoefU_x(j) + sigw(j)*(sigx(j)*dd_Km*dd_KUp(j)*BcoefU_x(j) + &
     127                                  Kech_m_w(j)*Kech_U_pw(j)*dd_BU(j))/(Kech_m(j)*Kech_Up(j))
     128       BcoefV(j) = BcoefV_x(j) + sigw(j)*(sigx(j)*dd_Km*dd_KVp(j)*BcoefV_x(j) + &
     129                                  Kech_m_w(j)*Kech_V_pw(j)*dd_BV(j))/(Kech_m(j)*Kech_Vp(j))
     130!>jyg
     131!
    315132!
    316133! Calcul des cdrag \'equivalents dans la couche
    317134!
    318        ycdragm(j) = ycdragm_x(j) + ywake_s(j)*dd_Cdragm
    319        ycdragh(j) = ycdragh_x(j) + ywake_s(j)*dd_Cdragh
     135       ycdragm(j) = ycdragm_x(j) + sigw(j)*dd_Cdragm(j)
     136       ycdragh(j) = ycdragh_x(j) + sigw(j)*dd_Cdragh(j)
     137       ycdragq(j) = ycdragq_x(j) + sigw(j)*dd_Cdragq(j)
    320138!
    321139! Calcul de T, q, u et v \'equivalents dans la couche 1
    322 !!       yt1(j) = yt_x(j,1) + ywake_s(j)*dd_t*(1.+sigx(j)*dd_Kh/KCT)
    323 !!       yq1(j) = yq_x(j,1) + ywake_s(j)*dd_q*(1.+sigx(j)*dd_Kh/KCQ)
    324 !!       yu1(j) = yu_x(j,1) + ywake_s(j)*dd_u*(1.+sigx(j)*dd_Km/KCU)
    325 !!       yv1(j) = yv_x(j,1) + ywake_s(j)*dd_v*(1.+sigx(j)*dd_Km/KCV)
    326        yt1(j) = yt_x(j,1) + ywake_s(j)*dd_t
    327        yq1(j) = yq_x(j,1) + ywake_s(j)*dd_q
    328        yu1(j) = yu_x(j,1) + ywake_s(j)*dd_u
    329        yv1(j) = yv_x(j,1) + ywake_s(j)*dd_v
     140!!       yt1(j) = yt_x(j,1) + sigw(j)*dd_t*(1.+sigx(j)*dd_Kh/KCT)
     141!!       yq1(j) = yq_x(j,1) + sigw(j)*dd_q*(1.+sigx(j)*dd_Kh/KCQ)
     142!!       yu1(j) = yu_x(j,1) + sigw(j)*dd_u*(1.+sigx(j)*dd_Km/KCU)
     143!!       yv1(j) = yv_x(j,1) + sigw(j)*dd_v*(1.+sigx(j)*dd_Km/KCV)
     144       yt1(j) = yt_x(j,1) + sigw(j)*dd_t
     145       yq1(j) = yq_x(j,1) + sigw(j)*dd_q
     146       yu1(j) = yu_x(j,1) + sigw(j)*dd_u
     147       yv1(j) = yv_x(j,1) + sigw(j)*dd_v
    330148
    331149
    … …  
    334152        RETURN
    335153
    336 END SUBROUTINE wx_pbl0_fuse
    337 
    338 SUBROUTINE wx_pbl0_split(knon, dtime, ywake_s, &
    339                        y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &
    340                        y_flux_t1_x, y_flux_t1_w, &
    341                        y_flux_q1_x, y_flux_q1_w, &
    342                        y_flux_u1_x, y_flux_u1_w, &
    343                        y_flux_v1_x, y_flux_v1_w, &
    344                        yfluxlat_x, yfluxlat_w, &
    345                        y_delta_tsurf &
    346                        )
    347 !
     154END SUBROUTINE wx_pbl0_merge
     155
     156SUBROUTINE wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
     157                                 sigw, beta, wcstar, wdens, &
     158                                 AT_x, AT_w, &
     159                                 BT_x, BT_w, &
     160                                 AcoefT0, AcoefQ0, BcoefT0, BcoefQ0, &
     161                                 AcoefT,  AcoefQ,  BcoefT,  BcoefQ, &
     162                                 HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
     163                                 phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
     164                                 g_T, g_Q, &
     165                                 Gamma_phiT, Gamma_phiQ, &
     166                                 dTs_ins, dqsatsrf_ins &
     167                                 )
     168!
     169
     170    USE wx_pbl_var_mod
     171
    348172    USE print_control_mod, ONLY: prt_level,lunout
    349173!
    350174    INCLUDE "YOMCST.h"
     175    INCLUDE "FCTTRE.h"
     176    INCLUDE "YOETHF.h"
    351177!
    352178    INTEGER,                      INTENT(IN)        :: knon    ! number of grid cells
    353179    REAL,                         INTENT(IN)        :: dtime   ! time step size (s)
    354     REAL, DIMENSION(knon),        INTENT(IN)        :: ywake_s ! cold pools fractional area
    355     REAL, DIMENSION(knon),        INTENT(IN)        :: y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1
    356 !
    357     REAL, DIMENSION(knon),        INTENT(OUT)       :: y_flux_t1_x, y_flux_t1_w
    358     REAL, DIMENSION(knon),        INTENT(OUT)       :: y_flux_q1_x, y_flux_q1_w
    359     REAL, DIMENSION(knon),        INTENT(OUT)       :: y_flux_u1_x, y_flux_u1_w
    360     REAL, DIMENSION(knon),        INTENT(OUT)       :: y_flux_v1_x, y_flux_v1_w
    361     REAL, DIMENSION(knon),        INTENT(OUT)       :: yfluxlat_x, yfluxlat_w
    362     REAL, DIMENSION(knon),        INTENT(OUT)       :: y_delta_tsurf
     180    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypplay  ! mid-layer pressure (Pa)
     181    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypaprs  ! pressure at layer interfaces (pa)
     182    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw    ! cold pool fractional area
     183    REAL, DIMENSION(knon),        INTENT(IN)        :: beta    ! evaporation by potential evaporation
     184    REAL, DIMENSION(knon),        INTENT(IN)        :: wcstar   ! cold pool gust front speed
     185    REAL, DIMENSION(knon),        INTENT(IN)        :: wdens    ! cold pool number density
     186    REAL, DIMENSION(knon),        INTENT(IN)        :: AT_x, AT_w
     187    REAL, DIMENSION(knon),        INTENT(IN)        :: BT_x, BT_w
     188    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT0, AcoefQ0, BcoefT0, BcoefQ0
     189!
     190    REAL, DIMENSION(knon),        INTENT(OUT)       :: AcoefT, AcoefQ, BcoefT, BcoefQ
     191    REAL, DIMENSION(knon),        INTENT(OUT)       :: HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn
     192    REAL, DIMENSION(knon),        INTENT(OUT)       :: phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0
     193    REAL, DIMENSION(knon),        INTENT(OUT)       :: g_T, g_Q
     194    REAL, DIMENSION(knon),        INTENT(OUT)       :: Gamma_phiT, Gamma_phiQ
     195    REAL, DIMENSION(knon),        INTENT(OUT)       :: dTs_ins, dqsatsrf_ins
     196!
     197! Local variables
     198    REAL, DIMENSION(knon)      :: qsat_x
     199    REAL, DIMENSION(knon)      :: qsat_w
     200    REAL, DIMENSION(knon)      :: dqsatdT_x
     201    REAL, DIMENSION(knon)      :: dqsatdT_w
     202!
     203    REAL, DIMENSION(knon)      :: T10_x
     204    REAL, DIMENSION(knon)      :: T10_w
     205    REAL, DIMENSION(knon)      :: phiT0_x
     206    REAL, DIMENSION(knon)      :: phiT0_w
     207    REAL, DIMENSION(knon)      :: phiQ0_x
     208    REAL, DIMENSION(knon)      :: phiQ0_w
     209    REAL, DIMENSION(knon)      :: Rn0_x
     210    REAL, DIMENSION(knon)      :: Rn0_w
     211    REAL, DIMENSION(knon)      :: Rp1_x
     212    REAL, DIMENSION(knon)      :: Rp1_w
     213    REAL, DIMENSION(knon)      :: Rps_x
     214    REAL, DIMENSION(knon)      :: Rps_w
     215!
     216    REAL, DIMENSION(knon)      :: HTphiT_x
     217    REAL, DIMENSION(knon)      :: HTphiT_w
     218    REAL, DIMENSION(knon)      :: HTphiQ_x
     219    REAL, DIMENSION(knon)      :: HTphiQ_w
     220    REAL, DIMENSION(knon)      :: HTRn_x
     221    REAL, DIMENSION(knon)      :: HTRn_w
     222!
     223    REAL, DIMENSION(knon)      :: HQphiT_x
     224    REAL, DIMENSION(knon)      :: HQphiT_w
     225    REAL, DIMENSION(knon)      :: HQphiQ_x
     226    REAL, DIMENSION(knon)      :: HQphiQ_w
     227    REAL, DIMENSION(knon)      :: HQRn_x
     228    REAL, DIMENSION(knon)      :: HQRn_w
     229!
     230    REAL, DIMENSION(knon)      :: HQphiT_b
     231    REAL, DIMENSION(knon)      :: dd_HQphiT
     232    REAL, DIMENSION(knon)      :: HQphiQ_b
     233    REAL, DIMENSION(knon)      :: dd_HQphiQ
     234    REAL, DIMENSION(knon)      :: HQRn_b
     235    REAL, DIMENSION(knon)      :: dd_HQRn
     236!
     237
     238    REAL, DIMENSION(knon)    :: sigx
     239!
     240    REAL, DIMENSION(knon)    :: Ts, T1
     241!!!    REAL, DIMENSION(knon)    :: qsat, dqsat_dT
     242!!!    REAL, DIMENSION(knon)    :: phiT0
     243!
     244!!!    REAL, DIMENSION(knon)    :: Cp, Lv
     245    REAL, DIMENSION(knon)    :: tau, Inert
     246!
     247    REAL                     :: dd_Kh
     248    REAL                     :: zdelta, zcvm5, zcor
     249    REAL                     :: qsat
     250!
     251    INTEGER                  :: j
     252
     253
     254!----------------------------------------------------------------------------
     255!  Reference state
     256!  ---------------
     257!   dqsat_dT_w = dqsat_dT(Ts0_w)                          dqsat_dT_x = dqsat_dT(Ts0_x)
     258!   T10_w = (AT_w/Cp - Kech_T_w BT_w dtime Ts0_w)/(1 - Kech_T_w BT_w dtime)
     259!                                                T10_x = (AT_x/Cp - Kech_T_x BT_x dtime Ts0_x)/(1 - Kech_T_x BT_x dtime)
     260!   phiT0_w = Kech_T_pw (AT_w - Cp Ts0_w)                 phiT0_x = Kech_T_px (AT_x - Cp Ts0_x)
     261!   phiQ0_w = Kech_Q_sw (beta AQ_w - qsatsrf0_w)          phiQ0_x = Kech_Q_sx (beta AQ_x - qsatsrf0_x)
     262!   Rn0_w = eps_1 Rsigma T10_w^4 - Rsigma Ts0_w^4         Rn0_x = eps_1 Rsigma T10_x^4 - Rsigma Ts0_x^4
     263!   Rp1_w = 4 eps_1 Rsigma T10_w^3                        Rp1_x = 4 eps_1 Rsigma T10_x^3
     264!   Rps_w = 4 Rsigma Ts0_w^3                              Rps_x = 4 Rsigma Ts0_x^3
     265!
     266!   phiT0_b = sigw phiT0_w + sigx phiT0_x
     267!   dphiT0 = phiT0_w - phiT0_x
     268!   phiQ0_b = sigw phiQ0_w + sigx phiQ0_x
     269!   dphiQ0 = phiQ0_w - phiQ0_x
     270!   Rn0_b = sigw Rn0_w + sigx Rn0_x
     271    dRn0 = Rn0_w - Rn0_x
     272!
     273!
     274!----------------------------------------------------------------------------
     275!  Elementary enthalpy equations
     276!  -----------------------------
     277!   phiT_w = phiT0_w - HTphiT_w (Ts_w-Ts0_w)            phiT_x = phiT0_x - HTphiT_x (Ts_x-Ts0_x)
     278!   phiQ_w = phiQ0_w - HTphiQ_w (Ts_w-Ts0_w)            phiQ_x = phiQ0_x - HTphiQ_x (Ts_x-Ts0_x)
     279!   Rn_w   = Rn0_w   - HTRn_w   (Ts_w-Ts0_w)            Rn_x   = Rn0_x   - HTRn_x   (Ts_x-Ts0_x)
     280!  DFlux_DT coefficients
     281!  ---------------------
     282!   Heat flux equation
     283!     HTphiT_w = Cp Kech_T_pw                            HTphiT_x = Cp Kech_T_px
     284!   Moisture flux equation
     285!     HTphiQ_w = beta Kech_Q_sw dqsat_dT_w               HTphiQ_x = beta Kech_Q_sx dqsat_dT_x
     286!   Radiation equation
     287!     HTRn_w = Rp1_w Kech_T_pw BcoefT_w dtime + Rps_w    HTRn_x = Rp1_x Kech_T_px BcoefT_x dtime + Rps_x
     288!
     289!----------------------------------------------------------------------------
     290!  Elementary moisture equations
     291!  -----------------------------
     292!   beta Ts_w   = beta Ts0_w    + QQ_w     (qsatsrf_w-qsatsrf0_w)    beta Ts_x   = beta Ts0_x    + QQ_x     (qsatsrf_x-qsatsrf0_x)
     293!   beta phiT_w = beta phiT0_w - HQphiT_w (qsatsrf_w-qsatsrf0_w)    beta phiQ_x = beta phiQ0_x - HTphiQ_x (qsatsrf_x-qsatsrf0_x)
     294!   beta phiQ_w = beta phiQ0_w - HQphiQ_w (qsatsrf_w-qsatsrf0_w)    beta phiQ_x = beta phiQ0_x - HTphiQ_x (qsatsrf_x-qsatsrf0_x)
     295!   beta Rn_w   = beta Rn0_w   - HQRn_w   (qsatsrf_w-qsatsrf0_w)    beta Rn_x   = beta Rn0_x   - HTRn_x   (qsatsrf_x-qsatsrf0_x)
     296!  DFluxDQ coefficients
     297!  ---------------------
     298!   dqsat_dT equation
     299!     QQ_w = 1. / dqsat_dT_w                             QQ_x = 1. / dqsat_dT_x
     300!   Heat flux equation
     301!     HQphiT_w = Cp Kech_T_pw QQ_w                       HQphiT_x = Cp Kech_T_px QQ_x
     302!   Moisture flux equation
     303!     HQphiQ_w = beta Kech_Q_sw                          HQphiQ_x = beta Kech_Q_sx
     304!   Radiation equation
     305!     HQRn_w = (Rp1_w Kech_T_pw BcoefT_w dtime + Rps_w) QQ_w
     306!                                         HQRn_x = (Rp1_x Kech_T_px BcoefT_x dtime + Rps_x) QQ_x
     307!
     308!----------------------------------------------------------------------------
     309! Mean values and w-x differences
     310! -------------------------------
     311!  HTphiT_b = sigw HTphiT_w + sigx HTphiT_x               dd_HTphiT = HTphiT_w - HTphiT_x
     312!  HTphiQ_b = sigw HTphiQ_w + sigx HTphiQ_x               dd_HTphiQ = HTphiQ_w - HTphiQ_x
     313!  HTRn_b   = sigw HTRn_w   + sigx HTRn_x                 dd_HTRn   = HTRn_w   - HTRn_x
     314!
     315!  QQ_b     = sigw QQ_w     + sigx QQ_x                   dd_QQ     = QQ_w     - QQ_x
     316!  HQphiT_b = sigw HQphiT_w + sigx HQphiT_x               dd_HQphiT = HQphiT_w - HQphiT_x
     317!  HQphiQ_b = sigw HQphiQ_w + sigx HQphiQ_x               dd_HQphiQ = HQphiQ_w - HQphiQ_x
     318!  HQRn_b   = sigw HQRn_w   + sigx HQRn_x                 dd_HQRn   = HQRn_w   - HQRn_x
     319!
     320!----------------------------------------------------------------------------
     321!  Equations
     322!  ---------
     323! (1 - g_T) dTs    = dTs_ins    + Gamma_phiT phiT
     324! (1 - g_Q) dqsatsrf = dqsatsrf_ins + Gamma_phiQ phiQ
     325!
     326! Feedback Gains
     327! --------------
     328! g_T = - (sqrt(tau)/I) [ HTphiT_b + Lv HTphiQ_b + HTRn_b +  &
     329!                        (dd_HTphiT + Lv dd_HTphiQ + dd_HTRn) (sigx - sigw - sigw sigx dd_HTphiT/HTphiT_b) ]
     330! g_Q = - (sqrt(tau)/(I QQ_b)) ( HQphiT_b + Lv HQphiQ_b + HQRn_b ) -  &
     331!         (sigx - sigw - sigw sigx dd_HQphiQ/HQphiQ_b)   &
     332!                          [ dd_QQ/QQ_b + (sqrt(tau)/(I QQ_b))(dd_HQphiT + Lv dd_HQphiQ + dd_HQRn) ]
     333!
     334!  Ts, qs Coupling coefficients                /
     335!  ----------------------------
     336! Gamma_phiT = (sqrt(tau)/(I HTphiT_b)) (dd_HTphiT + Lv dd_HTphiQ + dd_HTRn)
     337! Gamma_phiQ = (1/(HQphiQ_b QQ_b)) [ dd_QQ +  (sqrt(tau)/(I )) (dd_HQphiT + Lv dd_HQphiQ + dd_HQRn) ]
     338!
     339!  Insensitive changes
     340!  -------------------
     341! dTs_ins    = (1 - g_T) dTs0    - Gamma_phiT phiT0_b
     342! dqsatsrf_ins = (1 - g_Q) dqsatsrf0 - Gamma_phiQ phiQ0_b
     343!
     344!----------------------------------------------------------------------------
     345!  Effective coefficients Acoef and Bcoef
     346!  --------------------------------------
     347!  Equations
     348!  ---------
     349! Cp Ta = AcoefT + BcoefT phiT dtime
     350!    qa = AcoefQ + BcoefQ phiQ dtime
     351!  Coefficients
     352!  ------------
     353! AcoefT = AcoefT0 - sigw sigx (dd_KTp/Kech_Tp) Cp dTs_ins/(1 - g_T)
     354! BcoefT = BcoefT0 - sigw sigx (dd_KTp/Kech_Tp) Cp Gamma_phiT/(1 - g_T)/dtime
     355!
     356! AcoefQ = AcoefQ0 - sigw sigx (dd_KQp/Kech_Qp) dqs_ins/(1 - g_Q)
     357! BcoefQ = BcoefQ0 - sigw sigx (dd_KQp/Kech_Qp) Gamma_phiq/(1 - g_Q)/dtime
     358!
     359!==============================================================================
     360!
     361!
     362!  Parameters
     363!  ----------
     364   Inert(1:knon) = 2000.
     365   tau(1:knon) = sqrt(sigw(1:knon)/max(rpi*wdens(1:knon)*wcstar(1:knon)**2 , &
     366                                       sigw(1:knon)*1.e-12,smallestreal))
     367   sigx(1:knon) = 1.-sigw(1:knon)
     368!! Compute Cp, Lv, qsat, dqsat_dT.
     369!   C_p(1:knon) = RCpd
     370!   L_v(1:knon) = RLvtt
     371!
     372!      print *,' AAAA wx_pbl_dTs, C_p(j), qsat0(j), Ts0(j) : ', C_p(:), qsat0(:), Ts0(:)
     373!
     374!
     375   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))/  &
     376                   (1 - Kech_h_x(1:knon)*BT_x(1:knon)*dtime)
     377   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))/  &
     378                   (1 - Kech_h_w(1:knon)*BT_w(1:knon)*dtime)
     379!
     380   phiT0_x(1:knon) = Kech_T_px(1:knon)*(AT_x(1:knon) - C_p(1:knon)*Ts0_x(1:knon))
     381   phiT0_w(1:knon) = Kech_T_pw(1:knon)*(AT_w(1:knon) - C_p(1:knon)*Ts0_w(1:knon))
     382!
     383   phiQ0_x(1:knon) = Kech_Q_sx(1:knon)*(beta(1:knon)*AQ_x(1:knon) - qsatsrf0_x(1:knon))
     384   phiQ0_w(1:knon) = Kech_Q_sw(1:knon)*(beta(1:knon)*AQ_w(1:knon) - qsatsrf0_w(1:knon))
     385!
     386   Rn0_x(1:knon) = eps_1*Rsigma*T10_x(1:knon)**4 - Rsigma*Ts0_x(1:knon)**4
     387   Rn0_w(1:knon) = eps_1*Rsigma*T10_w(1:knon)**4 - Rsigma*Ts0_w(1:knon)**4
     388!
     389   Rp1_x(1:knon) = 4*eps_1*Rsigma*T10_x(1:knon)**3
     390   Rp1_w(1:knon) = 4*eps_1*Rsigma*T10_w(1:knon)**3
     391!
     392   Rps_x(1:knon) = 4*Rsigma*Ts0_x(1:knon)**3
     393   Rps_w(1:knon) = 4*Rsigma*Ts0_w(1:knon)**3
     394!
     395!  DFlux_DT coefficients
     396!  ---------------------
     397!   Heat flux equation
     398     HTphiT_x(1:knon) = C_p(1:knon)*Kech_T_px(1:knon)
     399     HTphiT_w(1:knon) = C_p(1:knon)*Kech_T_pw(1:knon)                       
     400!   Moisture flux equation
     401     HTphiQ_x(1:knon) = beta(1:knon)*Kech_Q_sx(1:knon)*dqsatdT0_x(1:knon)
     402     HTphiQ_w(1:knon) = beta(1:knon)*Kech_Q_sw(1:knon)*dqsatdT0_w(1:knon)         
     403!   Radiation equation
     404     HTRn_x(1:knon) = Rp1_x(1:knon)*Kech_T_px(1:knon)*BT_x(1:knon)*dtime + Rps_x(1:knon)
     405     HTRn_w(1:knon) = Rp1_w(1:knon)*Kech_T_pw(1:knon)*BT_w(1:knon)*dtime + Rps_w(1:knon) 
     406!
     407!  DFluxDQ coefficients
     408!  ---------------------
     409!   Heat flux equation
     410     HQphiT_x(1:knon) = C_p(1:knon)*Kech_T_px(1:knon)*QQ_x(1:knon)
     411     HQphiT_w(1:knon) = C_p(1:knon)*Kech_T_pw(1:knon)*QQ_w(1:knon)                 
     412!   Moisture flux equation
     413     HQphiQ_x(1:knon) = beta(1:knon)*Kech_Q_sx(1:knon)
     414     HQphiQ_w(1:knon) = beta(1:knon)*Kech_Q_sw(1:knon)                   
     415!   Radiation equation
     416     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)
     417     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)
     418!
     419! Mean values and w-x differences
     420! -------------------------------
     421  phiT0_b(1:knon) = sigw(1:knon)*phiT0_w(1:knon) + sigx(1:knon)*phiT0_x(1:knon)           
     422  phiQ0_b(1:knon) = sigw(1:knon)*phiQ0_w(1:knon) + sigx(1:knon)*phiQ0_x(1:knon)         
     423  Rn0_b(1:knon)   = sigw(1:knon)*Rn0_w(1:knon)   + sigx(1:knon)*Rn0_x(1:knon)         
     424!
     425  dphiT0(1:knon) = phiT0_w(1:knon) - phiT0_x(1:knon)           
     426  dphiQ0(1:knon) = phiQ0_w(1:knon) - phiQ0_x(1:knon)         
     427  dRn0(1:knon)   = Rn0_w(1:knon)   - Rn0_x(1:knon)         
     428!
     429  HTphiT_b(1:knon) = sigw(1:knon)*HTphiT_w(1:knon) + sigx(1:knon)*HTphiT_x(1:knon)           
     430  dd_HTphiT(1:knon) = HTphiT_w(1:knon) - HTphiT_x(1:knon)
     431!
     432  HTphiQ_b(1:knon) = sigw(1:knon)*HTphiQ_w(1:knon) + sigx(1:knon)*HTphiQ_x(1:knon)         
     433  dd_HTphiQ(1:knon) = HTphiQ_w(1:knon) - HTphiQ_x(1:knon)
     434!
     435  HTRn_b(1:knon)   = sigw(1:knon)*HTRn_w(1:knon)   + sigx(1:knon)*HTRn_x(1:knon)           
     436  dd_HTRn(1:knon)   = HTRn_w(1:knon)   - HTRn_x(1:knon)
     437!
     438  HQphiT_b(1:knon) = sigw(1:knon)*HQphiT_w(1:knon) + sigx(1:knon)*HQphiT_x(1:knon)         
     439  dd_HQphiT(1:knon) = HQphiT_w(1:knon) - HQphiT_x(1:knon)
     440!
     441  HQphiQ_b(1:knon) = sigw(1:knon)*HQphiQ_w(1:knon) + sigx(1:knon)*HQphiQ_x(1:knon)         
     442  dd_HQphiQ(1:knon) = HQphiQ_w - HQphiQ_x(1:knon)
     443!
     444  HQRn_b(1:knon)   = sigw(1:knon)*HQRn_w(1:knon)   + sigx(1:knon)*HQRn_x(1:knon)             
     445  dd_HQRn(1:knon)   = HQRn_w(1:knon)   - HQRn_x(1:knon)
     446!
     447! Feedback Gains
     448! --------------
     449 g_T(1:knon) = - (sqrt(tau(1:knon))/Inert(1:knon))  &
     450               * (HTphiT_b(1:knon) + L_v(1:knon)*HTphiQ_b(1:knon) + HTRn_b(1:knon)  &
     451                 + (dd_HTphiT(1:knon) + L_v(1:knon)*dd_HTphiQ(1:knon) + dd_HTRn(1:knon))  &
     452                 * (sigx(1:knon) - sigw(1:knon) - sigw(1:knon)*sigx(1:knon)*dd_HTphiT(1:knon)/HTphiT_b(1:knon)) )
     453!
     454!!!! DO j = 1,knon
     455!!!!  IF (mod(j,20) .eq.0) THEN
     456!!!!   print *, '   j     dd_QQ       QQ_b  dd_HQphiQ  dd_HQphiT   dd_HQRn   HQphiQ_b   HQphiT_b     HQRn_b '
     457!!!!  ENDIF
     458!!!!   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)
     459!!!! 1789 FORMAT( I4, 10(1X,E10.2))
     460!!!! ENDDO
     461   g_Q(1:knon) = - (dd_QQ(1:knon)/QQ_b(1:knon)) *  &
     462                    (sigx(1:knon)-sigw(1:knon)-sigw(1:knon)*sigx(1:knon)*dd_KQs(1:knon)/Kech_Qs(1:knon)) &
     463                 - sqrt(tau(1:knon))/(Inert(1:knon)*QQ_b(1:knon)) *  &
     464                   ( HQphiT_b(1:knon) + L_v(1:knon)*HQphiQ_b(1:knon) + HQRn_b(1:knon) +  &
     465                      (sigx(1:knon) - sigw(1:knon) - sigw(1:knon)*sigx(1:knon)*dd_KQs(1:knon)/Kech_Qs(1:knon)) *  &
     466                       (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
     467
     468!!   g_Q(1:knon) = - (dd_QQ(1:knon)/QQ_b(1:knon)) *  &
     469!!                    (sigx(1:knon)-sigw(1:knon)-sigw(1:knon)*sigx(1:knon)*dd_HQphiQ(1:knon)/HQphiQ_b(1:knon)) &
     470!!                 - sqrt(tau(1:knon))/(Inert(1:knon)*QQ_b(1:knon)) *  &
     471!!                   ( HQphiT_b(1:knon) + L_v(1:knon)*HQphiQ_b(1:knon) + HQRn_b(1:knon) +  &
     472!!                      (sigx(1:knon) - sigw(1:knon) - sigw(1:knon)*sigx(1:knon)*dd_HQphiQ(1:knon)/HQphiQ_b(1:knon)) *  &
     473!!                       (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
     474
     475!! g_Q(1:knon) = - (sqrt(tau(1:knon))/(Inert(1:knon)*QQ_b(1:knon))) *  &
     476!!                 ( HQphiT_b(1:knon) + L_v(1:knon)*HQphiQ_b(1:knon) + HQRn_b(1:knon) )  &
     477!!               - (sigx(1:knon) - sigw(1:knon) - sigw(1:knon)*sigx(1:knon)*dd_HQphiQ(1:knon)/HQphiQ_b(1:knon)) *   &
     478!!                 ( dd_QQ(1:knon)/QQ_b(1:knon)   &
     479!!                 + (sqrt(tau(1:knon))/(Inert(1:knon)*QQ_b(1:knon)))  &
     480!!                 * (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
     481
     482!  Ts, qs Coupling coefficients                /
     483!  ----------------------------
     484  Gamma_phiT(1:knon) = (sqrt(tau(1:knon))/(Inert(1:knon)*HTphiT_b(1:knon)))  &
     485                     * (dd_HTphiT(1:knon) + L_v(1:knon)*dd_HTphiQ(1:knon) + dd_HTRn(1:knon))
     486!
     487  Gamma_phiQ(1:knon) = (1./(Kech_Qs(1:knon)*QQ_b(1:knon))) * &
     488                        ( dd_QQ(1:knon)   &
     489                         + (sqrt(tau(1:knon))/(Inert(1:knon))) *  &
     490                          (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
     491
     492!!  Gamma_phiQ(1:knon) = (beta(1:knon)/(HQphiQ_b(1:knon)*QQ_b(1:knon))) * &
     493!!                        ( dd_QQ(1:knon)   &
     494!!                         + (sqrt(tau(1:knon))/(Inert(1:knon))) *  &
     495!!                          (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
     496
     497!!  Gamma_phiQ(1:knon) = (1/(HQphiQ_b(1:knon)*QQ_b(1:knon)))   &
     498!!                     * ( dd_QQ(1:knon)   &
     499!!                       + (sqrt(tau(1:knon))/(Inert(1:knon)))  &
     500!!                       * (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
     501!
     502!  Insensitive changes
     503!  -------------------
     504  dTs_ins(1:knon)    = (sqrt(tau(1:knon))/Inert(1:knon))*  &
     505                       (dphiT0(1:knon) + L_v(1:knon)*dphiQ0(1:knon) + dRn0(1:knon))
     506!
     507  dqsatsrf_ins(1:knon) = (beta(1:knon)/QQ_b(1:knon))*dTs_ins(1:knon)
     508!
     509   IF (prt_level .Ge. 10) THEN
     510      print *,'wx_pbl_merge, tau         ', tau
     511      print *,'wx_pbl_merge, AcoefT0     ', AcoefT0
     512      print *,'wx_pbl_merge, AcoefQ0     ', AcoefQ0
     513      print *,'wx_pbl_merge, BcoefT0     ', BcoefT0
     514      print *,'wx_pbl_merge, BcoefQ0     ', BcoefQ0
     515      print *,'wx_pbl_merge, qsat0_w, qsat0_x ', (qsat0_w(j), qsat0_x(j),j=1,knon)
     516      print *,'wx_pbl_merge, dqsatdT0_w, dqsatdT0_x ', (dqsatdT0_w(j), dqsatdT0_x(j),j=1,knon)
     517   ENDIF
     518!
     519!----------------------------------------------------------------------------
     520! 
     521!------------------------------------------------------------------------------
     522! 
     523!    Effective coefficients Acoef and Bcoef
     524!    --------------------------------------
     525   DO j = 1,knon
     526     AcoefT(j) = AcoefT0(j) - sigw(j)*sigx(j)*(dd_KTp(j)/Kech_Tp(j))*C_p(j)*   &
     527                 (dTs0(j) + (dTs_ins(j)-dTs0(j)-Gamma_phiT(j)*phiT0_b(j))/(1. - g_T(j)))
     528     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
     529     
     530     AcoefQ(j) = AcoefQ0(j) - sigw(j)*sigx(j)*(dd_KQs(j)/Kech_Qs(j))*    &
     531                 (dqsatsrf0(j) + (dqsatsrf_ins(j)-(beta(j)/QQ_b(j))*dTs0(j)-Gamma_phiQ(j)*phiQ0_b(j))/(1 - g_Q(j)))/ &
     532                 max(beta(j),1.e-4)
     533     BcoefQ(j) = BcoefQ0(j) - sigw(j)*sigx(j)*(dd_KQs(j)/Kech_Qs(j))*Gamma_phiQ(j)/(1 - g_Q(j))/ &
     534                 (max(beta(j),1.e-4)*dtime)
     535!!     AcoefQ(j) = AcoefQ0(j) - sigw(j)*sigx(j)*(dd_KQs(j)/Kech_Qs(j))*    &
     536!!                 (dqsatsrf0(j) + (dqsatsrf_ins(j)-(beta(j)/QQ_b(j))*dTs0(j)-Gamma_phiQ(j)*phiQ0_b(j))/(1 - g_Q(j)))/ &
     537!!                 beta(j)
     538!!     BcoefQ(j) = BcoefQ0(j) - sigw(j)*sigx(j)*(dd_KQs(j)/Kech_Qs(j))*Gamma_phiQ(j)/(1 - g_Q(j))/(beta(j)*dtime)
     539   ENDDO ! j = 1,knon
     540   
     541   IF (prt_level .Ge. 10) THEN
     542   print *,'wx_pbl_dts AAAA BcoefQ, BcoefQ0, sigw ', &
     543                            BcoefQ, BcoefQ0, sigw
     544      print *,'wx_pbl_dts_merge, dTs_ins      ', dTs_ins
     545      print *,'wx_pbl_dts_merge, dqs_ins      ', dqsatsrf_ins
     546   ENDIF
     547
     548     RETURN
     549
     550END SUBROUTINE wx_pbl_dts_merge
     551
     552SUBROUTINE wx_pbl_split(knon, nsrf, dtime, sigw, beta, iflag_split, &
     553                       g_T, g_Q, &
     554                       Gamma_phiT, Gamma_phiQ, &
     555                       dTs_ins, dqsatsrf_ins, &
     556                       phiT, phiQ, phiU, phiV, &
     557!!!!                       HTRn_b, dd_HTRn, HTphiT_b, dd_HTphiT, &
     558                       phiQ0_b, phiT0_b, &
     559                       phiT_x, phiT_w, &
     560                       phiQ_x, phiQ_w, &
     561                       phiU_x, phiU_w, &
     562                       phiV_x, phiV_w, &
     563                       philat_x, philat_w, &
     564!!!!                       Rn_b, dRn, &
     565                       dqsatsrf, &
     566                       dTs, delta_qsurf &
     567                       )
     568!
     569
     570    USE wx_pbl_var_mod
     571
     572    USE print_control_mod, ONLY: prt_level,lunout
     573    USE indice_sol_mod, ONLY: is_oce
     574!
     575    INCLUDE "YOMCST.h"
     576!
     577    INTEGER,                      INTENT(IN)        :: knon    ! number of grid cells
     578    INTEGER,                      INTENT(IN)        :: nsrf    ! surface type
     579    REAL,                         INTENT(IN)        :: dtime   ! time step size (s)
     580    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw ! cold pools fractional area
     581    REAL, DIMENSION(knon),        INTENT(IN)        :: beta ! aridity factor
     582    INTEGER,                      INTENT(IN)        :: iflag_split
     583    REAL, DIMENSION(knon),        INTENT(IN)        :: g_T, g_Q
     584    REAL, DIMENSION(knon),        INTENT(IN)        :: Gamma_phiT, Gamma_phiQ
     585    REAL, DIMENSION(knon),        INTENT(IN)        :: dTs_ins, dqsatsrf_ins
     586    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT, phiQ, phiU, phiV
     587    REAL, DIMENSION(knon),        INTENT(IN)        :: phiQ0_b, phiT0_b
     588!
     589    REAL, DIMENSION(knon),        INTENT(OUT)       :: phiT_x, phiT_w
     590    REAL, DIMENSION(knon),        INTENT(OUT)       :: phiQ_x, phiQ_w
     591    REAL, DIMENSION(knon),        INTENT(OUT)       :: phiU_x, phiU_w
     592    REAL, DIMENSION(knon),        INTENT(OUT)       :: phiV_x, phiV_w
     593    REAL, DIMENSION(knon),        INTENT(OUT)       :: philat_x, philat_w
     594    REAL, DIMENSION(knon),        INTENT(OUT)       :: dqsatsrf      ! beta delta(qsat(Ts))
     595    REAL, DIMENSION(knon),        INTENT(OUT)       :: dTs           ! Temperature difference at surface
     596    REAL, DIMENSION(knon),        INTENT(OUT)       :: delta_qsurf
    363597!
    364598!! Local variables
    365599    INTEGER                    :: j
    366     REAL, DIMENSION(knon)      :: y_delta_flux_t1, y_delta_flux_q1, y_delta_flux_u1, y_delta_flux_v1
    367 !
    368     REAL                       :: DDT, DDQ, DDU, DDV
    369     REAL                       :: LambdaTs, LambdaQs, LambdaUs, LambdaVs
     600    REAL, DIMENSION(knon)      :: dphiT, dphiQ, dphiU, dphiV
     601    REAL, DIMENSION(knon)      :: q1_x, q1_w
    370602!
    371603    REAL, DIMENSION(knon)      :: sigx       ! fractional area of (x) region
     604
     605!----------------------------------------------------------------------------
     606!  Equations
     607!  ---------
     608!!!!!! (1 - g_T) dTs    = dTs_ins    + Gamma_phiT phiT
     609!!!!!! (1 - g_Q) dqsatsrf = dqsatsrf_ins + Gamma_phiQ phiQ
     610!!!!!! dphiT = (dd_KTp/KTp) phiT + (     dd_AT - C_p dTs)*KxKwTp/KTp
     611!!!!!! dphiQ = (dd_KQs/KQs) phiQ + (beta dd_AQ - dqsatsrf )*KxKwQs/KQs
     612!!!!!! dphiU = (dd_KUp/KUp) phiU + (     dd_AU          )*KxKwUp/KUp
     613!!!!!! dphiV = (dd_KVp/KVp) phiV + (     dd_AV          )*KxKwVp/KVp
     614!
     615! (1 - g_T) (dTs-dTs0)    = dTs_ins-dTs0    + Gamma_phiT (phiT-phiT0)
     616! (1 - g_Q) dqsatsrf = dqsatsrf_ins + Gamma_phiQ phiQ
     617! dphiT = (dd_KTp/KTp) phiT + (     dd_AT - C_p dTs)*KxKwTp/KTp
     618! dphiQ = (dd_KQs/KQs) phiQ + (beta dd_AQ - dqsatsrf )*KxKwQs/KQs
     619! dphiU = (dd_KUp/KUp) phiU + (     dd_AU          )*KxKwUp/KUp
     620! dphiV = (dd_KVp/KVp) phiV + (     dd_AV          )*KxKwVp/KVp
     621!
    372622!!
    373         sigx(:) = 1.-ywake_s(:)
    374 
    375         DO j=1,knon
    376 !
    377        DDT = Kech_Tp(j)
    378        DDQ = Kech_Qp(j)
    379        DDU = Kech_Up(j)
    380        DDV = Kech_Vp(j)
    381 !
    382        LambdaTs =  dd_KTp(j)/DDT
    383        LambdaQs =  dd_KQp(j)/DDQ
    384        LambdaUs =  dd_KUp(j)/DDU
    385        LambdaVs =  dd_KVp(j)/DDV
    386 !
    387          y_delta_flux_t1(j) = y_flux_t1(j)*LambdaTs + dd_AT(j)*KxKwTp(j)/DDT
    388          y_delta_flux_q1(j) = y_flux_q1(j)*LambdaQs + dd_AQ(j)*KxKwQp(j)/DDQ
    389          y_delta_flux_u1(j) = y_flux_u1(j)*LambdaUs + dd_AU(j)*KxKwUp(j)/DDU
    390          y_delta_flux_v1(j) = y_flux_v1(j)*LambdaVs + dd_AV(j)*KxKwVp(j)/DDV
    391 !
    392          y_flux_t1_x(j)=y_flux_t1(j) - ywake_s(j)*y_delta_flux_t1(j)
    393          y_flux_t1_w(j)=y_flux_t1(j) + (1.-ywake_s(j))*y_delta_flux_t1(j)
    394          y_flux_q1_x(j)=y_flux_q1(j) - ywake_s(j)*y_delta_flux_q1(j)
    395          y_flux_q1_w(j)=y_flux_q1(j) + (1.-ywake_s(j))*y_delta_flux_q1(j)
    396          y_flux_u1_x(j)=y_flux_u1(j) - ywake_s(j)*y_delta_flux_u1(j)
    397          y_flux_u1_w(j)=y_flux_u1(j) + (1.-ywake_s(j))*y_delta_flux_u1(j)
    398          y_flux_v1_x(j)=y_flux_v1(j) - ywake_s(j)*y_delta_flux_v1(j)
    399          y_flux_v1_w(j)=y_flux_v1(j) + (1.-ywake_s(j))*y_delta_flux_v1(j)
    400 !
    401          yfluxlat_x(j)=y_flux_q1_x(j)*RLVTT
    402          yfluxlat_w(j)=y_flux_q1_w(j)*RLVTT
    403 !
    404 !       Delta_tsurf computation
    405 !!         y_delta_tsurf(j) = (1./RCPD)*(ah(j)*dd_AT(j) + &
    406 !!                                       ah(j)*y_flux_t1(j)*dd_BT(j)*dtime + &
    407 !!                                       y_delta_flux_t1(j)*(ah(j)*BBT+bh(j)) )
    408 !
    409            y_delta_tsurf(j) = 0.
    410 !
    411         ENDDO
     623        sigx(:) = 1.-sigw(:)
     624!
     625!      print *,' AAAA wx_pbl_split, C_p(j), qsat0(j), Ts0(j) : ', C_p(:), qsat0(:), Ts0(:)
     626!
     627   IF (iflag_split .EQ. 2 .AND. nsrf .NE. is_oce) THEN
     628!
     629!   Delta_tsurf and  Delta_qsurf computation
     630!   -----------------------------------------
     631      IF (prt_level >=10 ) THEN
     632        print *,' wx_pbl_split, dTs_ins, dTs0 , Gamma_phiT, g_T ', dTs_ins, dTs0, Gamma_phiT, g_T
     633        print *,' wx_pbl_split, dqsatsrf_ins, Gamma_phiQ, g_q ', dqsatsrf_ins, Gamma_phiQ, g_q
     634      ENDIF
     635!
     636      DO j = 1,knon
     637        dTs(j)    = dTs0(j) + (dTs_ins(j) - dTs0(j) + Gamma_phiT(j)*(phiT(j)-phiT0_b(j)) )/(1 - g_T(j))
     638        dqsatsrf(j) = dqsatsrf0(j) + (dqsatsrf_ins(j) - (beta(j)/QQ_b(j))*dTs0(j) + &
     639                       Gamma_phiQ(j)*(phiQ(j)-phiQ0_b(j)) )/(1 - g_Q(j))
     640      ENDDO ! j = 1,knon
     641!
     642        IF (prt_level >=10 ) THEN
     643          print *,' wx_pbl_split, dqsatsrf0, QQ_b ', dqsatsrf0, QQ_b
     644          print *,' wx_pbl_split, phiT0_b, phiT, dTs ', phiT0_b, phiT, dTs
     645          print *,' wx_pbl_split, phiQ0_b, phiQ, dqsatsrf ', phiQ0_b, phiQ, dqsatsrf
     646        ENDIF
     647   ELSE
     648        dTs(:) = 0.
     649        dqsatsrf(:) = 0.
     650   ENDIF ! (iflag_split .EQ. 2 .AND. nsrf .NE. is_oce)
     651!
     652     DO j = 1,knon
     653       dphiT(j) = (phiT(j)*dd_KTp(j) + (        dd_AT(j) - C_p(j)*dTs(j))*KxKwTp(j))/Kech_Tp(j)
     654       dphiQ(j) = (phiQ(j)*dd_KQs(j) + (beta(j)*dd_AQ(j) -        dqsatsrf(j))*KxKwQs(j))/Kech_Qs(j)
     655       dphiU(j) = (phiU(j)*dd_KUp(j) +          dd_AU(j)                 *KxKwUp(j))/Kech_Up(j)
     656       dphiV(j) = (phiV(j)*dd_KVp(j) +          dd_AV(j)                 *KxKwVp(j))/Kech_Vp(j)
     657!
     658       phiT_x(j)=phiT(j) - sigw(j)*dphiT(j)
     659       phiT_w(j)=phiT(j) + sigx(j)*dphiT(j)
     660       phiQ_x(j)=phiQ(j) - sigw(j)*dphiQ(j)
     661       phiQ_w(j)=phiQ(j) + sigx(j)*dphiQ(j)
     662       phiU_x(j)=phiU(j) - sigw(j)*dphiU(j)
     663       phiU_w(j)=phiU(j) + sigx(j)*dphiU(j)
     664       phiV_x(j)=phiV(j) - sigw(j)*dphiV(j)
     665       phiV_w(j)=phiV(j) + sigx(j)*dphiV(j)
     666!
     667       philat_x(j)=phiQ_x(j)*RLVTT
     668       philat_w(j)=phiQ_w(j)*RLVTT
     669     ENDDO ! j = 1,knon
     670!
     671     DO j = 1,knon
     672       q1_x(j) = AQ_x(j) + BQ_x(j)*phiQ_x(j)*dtime
     673       q1_w(j) = AQ_w(j) + BQ_w(j)*phiQ_w(j)*dtime
     674     ENDDO ! j = 1,knon
     675     DO j = 1,knon
     676       delta_qsurf(j) = (1.-beta(j))*(q1_w(j) - q1_x(j)) + dqsatsrf(j)
     677     ENDDO ! j = 1,knon
     678!
     679!!  Do j = 1,knon
     680!!     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)
     681!!     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)
     682!!  ENDDO
     683!
     684        IF (prt_level >=10 ) THEN
     685          print *,' wx_pbl_split, phiT, dphiT, dTs ', phiT, dphiT, dTs
     686          print *,' wx_pbl_split, phiQ, dphiQ, dqsatsrf ', phiQ, dphiQ, dqsatsrf
     687        ENDIF
     688!
     689        IF (prt_level >=10 ) THEN
     690!!          print *,' wx_pbl_split, verif dqsatsrf = beta dqsatdT0 dTs '
     691!!          print *,' wx_pbl_split, dqsatsrf, dqsatdT0*dTs ', dqsatsrf, dqsatdT0*dTs
     692        ENDIF
    412693!
    413694        RETURN
    414695
    415 END SUBROUTINE wx_pbl0_split
    416 
    417 SUBROUTINE wx_pbl_final
    418 !
    419 !****************************************************************************************
    420 ! Deallocate module variables
    421 !
    422 !****************************************************************************************   
    423 !
    424     IF (ALLOCATED(Kech_Tp))        DEALLOCATE(Kech_Tp)
    425     IF (ALLOCATED(Kech_T_xp))      DEALLOCATE(Kech_T_xp)
    426     IF (ALLOCATED(Kech_T_wp))      DEALLOCATE(Kech_T_wp)
    427     IF (ALLOCATED(dd_KTp))         DEALLOCATE(dd_KTp)
    428     IF (ALLOCATED(KxKwTp))         DEALLOCATE(KxKwTp)
    429     IF (ALLOCATED(dd_AT))          DEALLOCATE(dd_AT)
    430     IF (ALLOCATED(dd_BT))          DEALLOCATE(dd_BT)
    431     IF (ALLOCATED(Kech_Qp))        DEALLOCATE(Kech_Qp)
    432     IF (ALLOCATED(Kech_Q_xp))      DEALLOCATE(Kech_Q_xp)
    433     IF (ALLOCATED(Kech_Q_wp))      DEALLOCATE(Kech_Q_wp)
    434     IF (ALLOCATED(dd_KQp))         DEALLOCATE(dd_KQp)
    435     IF (ALLOCATED(KxKwQp))         DEALLOCATE(KxKwQp)
    436     IF (ALLOCATED(dd_AQ))          DEALLOCATE(dd_AQ)
    437     IF (ALLOCATED(dd_BQ))          DEALLOCATE(dd_BQ)
    438     IF (ALLOCATED(Kech_Up))        DEALLOCATE(Kech_Up)
    439     IF (ALLOCATED(Kech_U_xp))      DEALLOCATE(Kech_U_xp)
    440     IF (ALLOCATED(Kech_U_wp))      DEALLOCATE(Kech_U_wp)
    441     IF (ALLOCATED(dd_KUp))         DEALLOCATE(dd_KUp)
    442     IF (ALLOCATED(KxKwUp))         DEALLOCATE(KxKwUp)
    443     IF (ALLOCATED(dd_AU))          DEALLOCATE(dd_AU)
    444     IF (ALLOCATED(dd_BU))          DEALLOCATE(dd_BU)
    445     IF (ALLOCATED(Kech_Vp))        DEALLOCATE(Kech_Vp)
    446     IF (ALLOCATED(Kech_V_xp))      DEALLOCATE(Kech_V_xp)
    447     IF (ALLOCATED(Kech_V_wp))      DEALLOCATE(Kech_V_wp)
    448     IF (ALLOCATED(KxKwVp))         DEALLOCATE(KxKwVp)
    449     IF (ALLOCATED(dd_KVp))         DEALLOCATE(dd_KVp)
    450     IF (ALLOCATED(dd_AV))          DEALLOCATE(dd_AV)
    451     IF (ALLOCATED(dd_BV))          DEALLOCATE(dd_BV)
    452 
    453 END SUBROUTINE wx_pbl_final
     696END SUBROUTINE wx_pbl_split
     697
     698SUBROUTINE wx_pbl_check( knon, dtime, ypplay, ypaprs, &
     699                               sigw, beta, iflag_split,   &
     700                               Ts0_b9, dTs09,   &
     701                               qs_b9, Ts_b9,  &                         ! yqsurf, Tsurf_new
     702                               dTs9, dqsatsrf9,   &
     703                               AcoefT_x, AcoefT_w, &
     704                               BcoefT_x, BcoefT_w, &
     705                               AcoefT0, AcoefQ0, BcoefT0, BcoefQ0, &
     706                               AcoefT,  AcoefQ,  BcoefT,  BcoefQ, &
     707                               phiT_b9, phiQ_b9,  &
     708                               phiT_x9, phiT_w9, &
     709                               phiQ_x9, phiQ_w9 &
     710                               )
     711!
     712
     713    USE wx_pbl_var_mod
     714
     715    USE print_control_mod, ONLY: prt_level,lunout
     716!
     717    INCLUDE "YOMCST.h"
     718    INCLUDE "FCTTRE.h"
     719    INCLUDE "YOETHF.h"
     720!
     721    INTEGER,                      INTENT(IN)        :: knon         ! number of grid cells
     722    REAL,                         INTENT(IN)        :: dtime        ! time step size (s)
     723    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypplay       ! mid-layer pressure (Pa)
     724    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypaprs       ! pressure at layer interfaces (pa)
     725    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw         ! cold pools fractional area
     726    REAL, DIMENSION(knon),        INTENT(IN)        :: beta         ! aridity factor
     727    INTEGER,                      INTENT(IN)        :: iflag_split
     728    REAL, DIMENSION(knon),        INTENT(IN)        :: Ts0_b9, dTs09
     729    REAL, DIMENSION(knon),        INTENT(IN)        :: qs_b9, Ts_b9         ! yqsurf, Tsurf_new
     730    REAL, DIMENSION(knon),        INTENT(IN)        :: dTs9, dqsatsrf9
     731    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT_x, AcoefT_w
     732    REAL, DIMENSION(knon),        INTENT(IN)        :: BcoefT_x, BcoefT_w
     733    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT0, AcoefQ0, BcoefT0, BcoefQ0
     734!
     735    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT, AcoefQ, BcoefT, BcoefQ
     736    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT_b9, phiQ_b9
     737    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT_x9, phiT_w9
     738    REAL, DIMENSION(knon),        INTENT(IN)        :: phiQ_x9, phiQ_w9
     739!
     740!! Local variables
     741    INTEGER                    :: j
     742    REAL, DIMENSION(knon)      :: sigx                 ! fractional area of (x) region
     743    REAL, DIMENSION(knon)      :: AcoefT_b, AcoefQ_b   ! mean values of AcoefT and AcoefQ
     744    REAL                       :: zzt, zzq, zzqsat
     745    REAL                       :: zdelta, zcvm5, zcor, qsat
     746    REAL, DIMENSION(knon)      :: qsat_w, qsat_x
     747    REAL, DIMENSION(knon)      :: dqsatdT_w, dqsatdT_x
     748    REAL, DIMENSION(knon)      :: qsat_bs              ! qsat(Ts_b)
     749    REAL, DIMENSION(knon)      :: qsat01, dqsatdT01
     750    REAL, DIMENSION(knon)      :: Ts_x, Ts_w, qs_x, qs_w
     751    REAL, DIMENSION(knon)      :: T1_x, T1_w, q1_x, q1_w
     752    REAL, DIMENSION(knon)      :: Rn_x, Rn_w
     753    REAL, DIMENSION(knon)      :: phiQ0_x, phiQ0_w
     754    REAL, DIMENSION(knon)      :: Ta, qa
     755    REAL, DIMENSION(knon)      :: qsatsrf_w, qsatsrf_x, qsatsrf_b
     756    REAL, DIMENSION(knon)      :: qsurf_w, qsurf_x
     757    REAL                       :: dphiT, dphiQ
     758    REAL                       :: dqsatsrf1
     759    REAL                       :: phiT_w1, phiT_w2
     760    REAL                       :: phiT_x1, phiT_x2
     761    REAL                       :: phiQ_w1, phiQ_w2, phiQ_w3
     762    REAL                       :: phiQ_x1, phiQ_x2, phiQ_x3
     763    REAL                       :: phiT_b1, phiQ_b1
     764    REAL                       :: Kech_Q_sw1, Kech_Q_sx1
     765    REAL                       :: evap_pot
     766
     767!----------------------------------------------------------------------------
     768! Equations to be checked:
     769! -----------------------
     770!  Input : Ts0_b, dTs0, Ts_b, dTs, qsatsrf_b, dqsatsrf,
     771!          phiT_b, phiQ_b, phiT_w, phiT_x, phiQ_w, phiQ_x,
     772!         
     773!          AcoefT, AcoefQ, AcoefT_w, AcoefQ_w, AcoefT_x, AcoefQ_x,
     774!          BcoefT, BcoefQ, BcoefT_w, BcoefQ_w, BcoefT_x, BcoefQ_x
     775!
     776!  C_p T1_w = AcoefT_w + BcoefT_w phiT_w Delta t          C_p T1_x = AcoefT_x + BcoefT_x phiT_x Delta t
     777!  q1_w = AQ_w + BQ_w phiQ_w Delta t                      q1_x = AQ_x + BQ_x phiQ_x Delta t
     778!  qsatsrf_w = beta qsat(Ts_w)                            qsatsrf_x = beta qsat(Ts_x)
     779!  qsurf_w = (1-beta) q1_w + qsatsrf_w                    qsurf_x = (1-beta) q1_x + qsatsrf_x
     780!  phiT_w = Kech_h_w C_p ( T1_w - Ts_w)                   phiT_x = Kech_h_x C_p ( T1_x - Ts_x)             
     781!  phiT_w = Kech_T_pw ( AcoefT_w - C_p Ts_w)              phiT_x = Kech_T_px ( AcoefT_x - C_p Ts_x)
     782!  phiq_w = Kech_h_w ( beta q1_w - qsatsrf_w)             phiq_x = Kech_h_x ( beta q1_x - qsatsrf_x))
     783!  phiq_w = Kech_Q_sw (beta AQ_w -qsatsrf_w)              phiq_x = Kech_Q_sx (beta AQ_x -qsatsrf_x)
     784!  phiq_w = Kech_h_w (q1_w - qsurf_w)                     phiq_x = Kech_h_x (q1_x - qsurf_x)
     785!  phiT_b = sigw phiT_w + sigx phiT_x                     dphiT = phiT_w - phiT_x
     786!  phiQ_b = sigw phiQ_w + sigx phiQ_x                     dphiQ = phiQ_w - phiQ_x
     787!  Ts_b = sigw Ts_w + sigx Ts_x                           dTs = Ts_w - Ts_x
     788!  qsatsrf_b = sigw qsatsrf_w + sigx qsatsrf_x
     789!  C_p Ta = AcoefT + BcoefT phiT_b Delta t
     790!  qa = AcoefQ + BcoefQ phiQ_b Delta t
     791!  phiT_b = Kech_h C_p (Ta - Ts_b)
     792!  phiQ_b = beta Kech_h (qa - qsatsrf_b)
     793!  dTs = sqrt(tau)/I (dphit + L_v dphiq + dR)
     794
     795!----------------------------------------------------------------------------
     796!
     797!!
     798        sigx(:) = 1.-sigw(:)
     799        AcoefT_b(1:knon) = AcoefT_x(1:knon) + sigw(1:knon)*dd_AT(1:knon)
     800        AcoefQ_b(1:knon) = AQ_x(1:knon) + sigw(1:knon)*dd_AQ(1:knon)
     801
     802! Compute the three qsat and dqsatdTs
     803! ---------------------------------------------
     804!!   C_p(1:knon) = RCpd
     805!!   L_v(1:knon) = RLvtt
     806    IF (prt_level >=10 ) THEN
     807      print *,' AAAA wx_pbl_check, C_p(j), qsat0(j), Ts0(j) : ', C_p(:), qsat0(:), Ts0(:)
     808    ENDIF ! (prt_level >=10 )
     809!
     810   DO j = 1, knon
     811      zdelta = MAX(0.,SIGN(1.,RTT-Ts0_b9(j)))
     812      zcvm5 = R5LES*(1.-zdelta) + R5IES*zdelta
     813      qsat = R2ES*FOEEW(Ts0_b9(j),zdelta)/ypaprs(j,1)
     814      qsat = MIN(0.5,qsat)
     815      zcor = 1./(1.-RETV*qsat)
     816      qsat01(j) = fqsat*qsat*zcor
     817!!      dqsatdT0(j) = FOEDE(Ts0_b(j),zdelta,zcvm5,qsat0(j),zcor)/RLVTT    ! jyg 20210116
     818!!      dqsatdT0(j) = (RLvtt*(1.-zdelta)+RLSTT*zdelta)*qsat0(j)/(Rv*Ts0_b(j)*Ts0_b(j))
     819      dqsatdT01(j) = fqsat*FOEDE(Ts0_b9(j),zdelta,zcvm5,qsat01(j),zcor)
     820   ENDDO
     821!
     822!--------------------------------------------------------------------------------------------------
     823        IF (prt_level >=10 ) THEN
     824!
     825          DO j = 1, knon
     826!
     827   print *,'wx_pbl_check: Kech_h, Kech_q ', Kech_h(j), Kech_q(j)
     828!
     829  Ta(j) = (AcoefT(j) + BcoefT(j)*phiT_b9(j)*dtime)/C_p(j)
     830  qa(j) = AcoefQ(j) + BcoefQ(j)*phiQ_b9(j)*dtime
     831    print *, 'wx_pbl_check: j, Ta, qa ', Ta(j), qa(j)
     832!
     833  qsat_bs(j) = qsat01(j) + dqsatdT01(j)*(Ts_b9(j)-Ts0_b9(j))
     834!
     835   print *,'wx_pbl_check: qsat01, qsat_bs ', j,qsat01(j), qsat_bs(j)
     836!
     837  Ts_x(j) = Ts_b9(j) - sigw(j)*dTs9(j)
     838  Ts_w(j) = Ts_b9(j) + sigx(j)*dTs9(j)
     839    print *, 'wx_pbl_check: j, Ts_b9, Ts_w, Ts_x ', j, Ts_b9(j), Ts_w(j), Ts_x(j)
     840!
     841  qsat_x(j) = qsat0_x(j) + dqsatdT0_x(j)*(Ts_x(j)-Ts0_x(j))
     842  qsat_w(j) = qsat0_w(j) + dqsatdT0_w(j)*(Ts_w(j)-Ts0_w(j))
     843!
     844   print *,'wx_pbl_check: qsat0_w, qsat0_x, qsat_w, qsat_x ', qsat0_w(j), qsat0_x(j), qsat_w(j), qsat_x(j)
     845!
     846  T1_x(j) = (AcoefT_x(j) + BcoefT_x(j)*phiT_x9(j)*dtime) / C_p(j)
     847  T1_w(j) = (AcoefT_w(j) + BcoefT_w(j)*phiT_w9(j)*dtime) / C_p(j)
     848    print *, 'wx_pbl_check: j, T1_w, T1_x ', j, T1_w(j), T1_x(j)
     849!
     850  q1_x(j) = AQ_x(j) + BQ_x(j)*phiQ_x9(j)*dtime
     851  q1_w(j) = AQ_w(j) + BQ_w(j)*phiQ_w9(j)*dtime
     852    print *, 'wx_pbl_check: j, q1_w, q1_x ', j, q1_w(j), q1_x(j)
     853!
     854   qsatsrf_x(j) = beta(j)*qsat_x(j)
     855   qsatsrf_w(j) = beta(j)*qsat_w(j)
     856   qsatsrf_b(j) = sigw(j)*qsatsrf_w(j) + sigx(j)*qsatsrf_x(j)
     857!
     858   dqsatsrf1 = qsatsrf_w(j) - qsatsrf_x(j)
     859    print *, 'wx_pbl_check: j, qsatsrf_w, qsatsrf_x, dqsatsrf1, dqsatsrf9 ', &
     860                          qsatsrf_w(j), qsatsrf_x(j), dqsatsrf1, dqsatsrf9(j)
     861!
     862   qsurf_x(j) = (1-beta(j))*q1_x(j) + qsatsrf_x(j)
     863   qsurf_w(j) = (1-beta(j))*q1_w(j) + qsatsrf_w(j)
     864    print *, 'wx_pbl_check: j, qsurf_w, qsurf_x ', j, qsurf_w(j), qsurf_x(j)
     865!
     866!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     867!  Test qsat01 = qsat0    et   dqsatdT01 = dqsatdT0
     868!------------------------------------------------------------------------------------------------------
     869   print *, 'wx_pbl_check: j, qsat01(j), qsat0(j) ', j, qsat01(j), qsat0(j)
     870   print *, 'wx_pbl_check: j, dqsatdT01(j), dqsatdT0(j) ', j, dqsatdT01(j), dqsatdT0(j)
     871!
     872!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     873!  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)
     874!------------------------------------------------------------------------------------------------------
     875  Kech_Q_sx1 = Kech_q_x(j)/(1.-beta(j)*Kech_q_x(j)*BQ_x(j)*dtime)
     876  Kech_Q_sw1 = Kech_q_w(j)/(1.-beta(j)*Kech_q_w(j)*BQ_w(j)*dtime)
     877    print *, 'wx_pbl_check: j, Kech_Q_sx1, Kech_Q_sx(j)', j, Kech_Q_sx1, Kech_Q_sx(j)
     878    print *, 'wx_pbl_check: j, Kech_Q_sw1, Kech_Q_sw(j)', j, Kech_Q_sw1, Kech_Q_sw(j)
     879!
     880!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     881!  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))
     882!-----------------------------------------------------
     883    phiT_x1 = Kech_h_x(j)*C_p(j)*(T1_x(j)-Ts_x(j))
     884    phiT_w1 = Kech_h_w(j)*C_p(j)*(T1_w(j)-Ts_w(j))
     885!
     886!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     887!  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))
     888!-----------------------------------------------------
     889    phiT_x2 = Kech_T_px(j)*(AcoefT_x(j)-C_p(j)*Ts_x(j))
     890    phiT_w2 = Kech_T_pw(j)*(AcoefT_w(j)-C_p(j)*Ts_w(j))
     891    print *, 'wx_pbl_check: j, phiT_w1, phiT_w2, phiT_w9 ', j, phiT_w1, phiT_w2, phiT_w9(j)
     892    print *, 'wx_pbl_check: j, phiT_x1, phiT_x2, phiT_x9 ', j, phiT_x1, phiT_x2, phiT_x9(j)
     893!
     894!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     895!  Test phiq_w = Kech_q_w ( beta q1_w - qsatsrf_w)    phiq_x = Kech_q_x ( beta q1_x - qsatsrf_x))
     896!--------------------------------------------------------------
     897    phiq_x1 = Kech_q_x(j)*( beta(j)*q1_x(j) - qsatsrf_x(j))
     898    phiq_w1 = Kech_q_w(j)*( beta(j)*q1_w(j) - qsatsrf_w(j))
     899!
     900!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     901!  Test  phiq_w = Kech_Q_sw (beta AQ_w -qsatsrf_w)     phiq_x = Kech_Q_sx (beta AQ_x -qsatsrf_x)
     902!--------------------------------------------------------------
     903    phiq_x2 = Kech_Q_sx(j)*(beta(j)*AQ_x(j) -qsatsrf_x(j))
     904    phiq_w2 = Kech_Q_sw(j)*(beta(j)*AQ_w(j) -qsatsrf_w(j))
     905!
     906!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     907!  Test phiq_w = Kech_q_w ( q1_w - qsurf_w)    phiq_x = Kech_q_x ( q1_x - qsurf_x))
     908!--------------------------------------------------------------
     909    phiq_x3 = Kech_q_x(j)*( q1_x(j) - qsurf_x(j))
     910    phiq_w3 = Kech_q_w(j)*( q1_w(j) - qsurf_w(j))
     911    print *, 'wx_pbl_check: j, phiQ_w1, phiQ_w2, phiQ_w3, phiQ_w9 ', j, phiQ_w1, phiQ_w2, phiQ_w3, phiQ_w9(j)
     912    print *, 'wx_pbl_check: j, phiQ_x1, phiQ_x2, phiQ_x3, phiQ_x9 ', j, phiQ_x1, phiQ_x2, phiQ_x3, phiQ_x9(j)
     913!
     914!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     915!  Test phiT_b = Kech_h C_p (Ta - Ts_b)
     916!--------------------------------------------------------------
     917   phiT_b1 = Kech_h(j)*C_p(j)*(Ta(j) - Ts_b9(j))
     918   print *, 'wx_pbl_check: j, phiT_b1, PhiT_b9 ', j, phiT_b1, PhiT_b9(j)
     919!
     920!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     921!  Test phiQ_b = beta Kech_q (qa - qsat_bs)
     922!--------------------------------------------------------------
     923   evap_pot = Kech_q(j)*(qa(j) - qsat_bs(j))
     924   phiQ_b1 = beta(j)*Kech_q(j)*(qa(j) - qsat_bs(j))
     925   print *, 'wx_pbl_check: j, beta, evap_pot, phiQ_b1, PhiQ_b9 ', j, beta(j), evap_pot, phiQ_b1, PhiQ_b9(j)
     926!
     927!
     928          ENDDO  ! j = 1, knon
     929         
     930        ENDIF   ! (prt_level >=10 )
     931!--------------------------------------------------------------------------------------------------
     932
     933        RETURN
     934
     935END SUBROUTINE wx_pbl_check
     936
     937SUBROUTINE wx_pbl_dts_check( knon, dtime, ypplay, ypaprs, &
     938                               sigw, beta, iflag_split,   &
     939                               Ts0_b9, dTs09,   &
     940                               qs_b9, Ts_b9,  &                         ! yqsurf, Tsurf_new
     941                               dqsatsrf9, dTs9, delta_qsurf9,   &
     942                               AcoefT_x, AcoefT_w, &
     943                               BcoefT_x, BcoefT_w, &
     944                               AcoefT0, AcoefQ0, BcoefT0, BcoefQ0, &
     945                               AcoefT,  AcoefQ,  BcoefT,  BcoefQ, &
     946                               HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
     947                               phiT0_b9, dphiT09, phiQ0_b9, dphiQ09, Rn0_b9, dRn09, &
     948                               g_T, g_Q, &
     949                               Gamma_phiT, Gamma_phiQ, &
     950                               dTs_ins, dqsatsrf_ins, &
     951                               phiT_b9, phiQ_b9,  &
     952                               phiT_x9, phiT_w9, &
     953                               phiQ_x9, phiQ_w9  &
     954                               )
     955!
     956
     957    USE wx_pbl_var_mod
     958
     959    USE print_control_mod, ONLY: prt_level,lunout
     960!
     961    INCLUDE "YOMCST.h"
     962    INCLUDE "FCTTRE.h"
     963    INCLUDE "YOETHF.h"
     964!
     965    INTEGER,                      INTENT(IN)        :: knon         ! number of grid cells
     966    REAL,                         INTENT(IN)        :: dtime        ! time step size (s)
     967    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypplay       ! mid-layer pressure (Pa)
     968    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypaprs       ! pressure at layer interfaces (pa)
     969    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw         ! cold pools fractional area
     970    REAL, DIMENSION(knon),        INTENT(IN)        :: beta         ! aridity factor
     971    INTEGER,                      INTENT(IN)        :: iflag_split
     972    REAL, DIMENSION(knon),        INTENT(IN)        :: Ts0_b9, dTs09
     973    REAL, DIMENSION(knon),        INTENT(IN)        :: qs_b9, Ts_b9         ! yqsurf, Tsurf_new
     974    REAL, DIMENSION(knon),        INTENT(IN)        :: dTs9, dqsatsrf9
     975    REAL, DIMENSION(knon),        INTENT(IN)        :: delta_qsurf9
     976    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT_x, AcoefT_w
     977    REAL, DIMENSION(knon),        INTENT(IN)        :: BcoefT_x, BcoefT_w
     978    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT0, AcoefQ0, BcoefT0, BcoefQ0
     979!
     980    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT, AcoefQ, BcoefT, BcoefQ
     981    REAL, DIMENSION(knon),        INTENT(IN)        :: HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn
     982    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT0_b9, dphiT09, phiQ0_b9, dphiQ09, Rn0_b9, dRn09
     983    REAL, DIMENSION(knon),        INTENT(IN)        :: g_T, g_Q
     984    REAL, DIMENSION(knon),        INTENT(IN)        :: Gamma_phiT, Gamma_phiQ
     985    REAL, DIMENSION(knon),        INTENT(IN)        :: dTs_ins, dqsatsrf_ins
     986    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT_b9, phiQ_b9
     987    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT_x9, phiT_w9
     988    REAL, DIMENSION(knon),        INTENT(IN)        :: phiQ_x9, phiQ_w9
     989!
     990!! Local variables
     991    INTEGER                    :: j
     992    REAL, DIMENSION(knon)      :: sigx       ! fractional area of (x) region
     993    REAL, DIMENSION(knon)      :: AcoefT_b, AcoefQ_b   ! mean values of AcoefT and AcoefQ
     994    REAL                       :: zzt, zzq, zzqsat
     995    REAL                       :: zdelta, zcvm5, zcor, qsat
     996    REAL, DIMENSION(knon)      :: qsat_w, qsat_x
     997    REAL, DIMENSION(knon)      :: Ts_x, Ts_w, qs_x, qs_w
     998    REAL, DIMENSION(knon)      :: T1_x, T1_w, q1_x, q1_w
     999    REAL, DIMENSION(knon)      :: Rn_x, Rn_w
     1000    REAL, DIMENSION(knon)      :: Rn_b, dRn
     1001    REAL, DIMENSION(knon)      :: phiQ0_x, phiQ0_w
     1002    REAL, DIMENSION(knon)      :: Ta, qa
     1003    REAL, DIMENSION(knon)      :: err_phiT_w, err_phiT_x
     1004    REAL, DIMENSION(knon)      :: err_phiq_w, err_phiq_x
     1005    REAL, DIMENSION(knon)      :: err_phiT_b
     1006    REAL, DIMENSION(knon)      :: err_phiQ_b
     1007    REAL, DIMENSION(knon)      :: err2_phiT_b
     1008    REAL                       :: T1A_x, T1A_w, q1A_x, q1A_w
     1009    REAL                       :: qsatsrf_w, qsatsrf_x, qsatsrfb, qsbA
     1010    REAL                       :: dphiT, dphiQ
     1011    REAL                       :: dphiT_H, dphiQ_H
     1012    REAL                       :: phiQ_pot
     1013    REAL                       :: phiQ_w_m_phiQ0_w
     1014    REAL                       :: phiQ_x_m_phiQ0_x
     1015    REAL                       :: dphiQ_m_dphiQ0
     1016    REAL                       :: dphiT_m_dphiT0
     1017    REAL                       :: dRN_m_dRn0
     1018    REAL                       :: phiTb_m_phiT0b
     1019
     1020!----------------------------------------------------------------------------
     1021! Equations to be checked:
     1022! -----------------------
     1023!  Input : Ts0_b, dTs0, Ts_b, dTs, qsatsrf_b, dqsatsrf,
     1024!          phiT_b, phiQ_b, phiT_w, phiT_x, phiQ_w, phiQ_x,
     1025!         
     1026!          AcoefT, AcoefQ, AcoefT_w, AcoefQ_w, AcoefT_x, AcoefQ_x,
     1027!          BcoefT, BcoefQ, BcoefT_w, BcoefQ_w, BcoefT_x, BcoefQ_x
     1028!
     1029!  Ts_w = Ts_b + sigx dTs                                 Ts_x = Ts_b - sigw dTs
     1030!  T1_w = AcoefT_w + BcoefT_w phiT_w Delta t              T1_x = AcoefT_x + BcoefT_x phiT_x Delta t
     1031!  q1_w = AcoefQ_w + BcoefQ_w phiQ_w Delta t              q1_x = AcoefQ_x + BcoefQ_x phiQ_x Delta t
     1032!  phiT_w = Kech_h_w ( T1_w - Ts_w)                       phiT_x = Kech_h_x ( T1_x - Ts_x)             
     1033!  phiq_w = beta Kech_h_w ( q1_w - qsat(Ts_w))            phiq_x = beta Kech_h_x ( q1_x - qsat(Ts_x))
     1034!  phiT_b = sigw phiT_w + sigx phiT_x                     dphiT = phiT_w - phiT_x
     1035!  phiQ_b = sigw phiQ_w + sigx phiQ_x                     dphiQ = phiQ_w - phiQ_x
     1036!  Ts_b = sigw Ts_w + sigx Ts_x                           dTs = Ts_w - Ts_x
     1037!  Ta = AcoefT + BcoefT phiT_b Delta t
     1038!  qa = AcoefQ + BcoefQ phiQ_b Delta t
     1039!  phiT_b = Kech_h (Ta - Ts_b)
     1040!  phiQ_b = beta Kech_h (qa - qsat(Ts_b))
     1041!  dTs = sqrt(tau)/I (dphit + L_v dphiq + dR)
     1042
     1043!----------------------------------------------------------------------------
     1044!
     1045!!
     1046        sigx(:) = 1.-sigw(:)
     1047        AcoefT_b(1:knon) = AcoefT_x(1:knon) + sigw(1:knon)*dd_AT(1:knon)
     1048        AcoefQ_b(1:knon) = AQ_x(1:knon) + sigw(1:knon)*dd_AQ(1:knon)
     1049
     1050   IF (prt_level >=10 ) THEN
     1051    print *,'->wx_pbl_dts_check, HTphiT_b, HTphiQ_b, HTRn_b ', &
     1052                             HTphiT_b, HTphiQ_b, HTRn_b
     1053    print *,'->wx_pbl_dts_check, dd_HTphiT, dd_HTphiQ, dd_HTRn ', &
     1054                             dd_HTphiT, dd_HTphiQ, dd_HTRn
     1055   ENDIF ! (prt_level >=10 )
     1056!
     1057! Compute the three qsat and dqsatdTs
     1058! ---------------------------------------------
     1059!!      print *,' AAAA wx_pbl_dts_check, C_p(j), qsat0(j), Ts0(j) : ',  &
     1060!!                                      (C_p(j), qsat0(j), Ts0(j), j = 1,knon)
     1061!
     1062!
     1063!--------------------------------------------------------------------------------------------------
     1064        IF (prt_level >=10 ) THEN
     1065!
     1066          DO j = 1, knon
     1067  Ts_x(j) = Ts_b9(j) - sigw(j)*dTs9(j)
     1068  Ts_w(j) = Ts_b9(j) + sigx(j)*dTs9(j)
     1069    print *, 'wx_pbl_dts_check: j, Ts_b9, Ts_w, Ts_x ', j, Ts_b9(j), Ts_w(j), Ts_x(j)
     1070!
     1071  qsat_x(j) = qsat0_x(j) + dqsatdT0_x(j)*(Ts_x(j)-Ts0_x(j))
     1072  qsat_w(j) = qsat0_w(j) + dqsatdT0_w(j)*(Ts_w(j)-Ts0_w(j))
     1073!
     1074  T1_x(j) = (AcoefT_x(j) + BcoefT_x(j)*phiT_x9(j)*dtime) / C_p(j)
     1075  T1_w(j) = (AcoefT_w(j) + BcoefT_w(j)*phiT_w9(j)*dtime) / C_p(j)
     1076    print *, 'wx_pbl_dts_check: j, T1_w, T1_x ', j, T1_w(j), T1_x(j)
     1077!
     1078  q1_x(j) = AQ_x(j) + BQ_x(j)*phiQ_x9(j)*dtime
     1079  q1_w(j) = AQ_w(j) + BQ_w(j)*phiQ_w9(j)*dtime
     1080    print *, 'wx_pbl_dts_check: j, q1_w, q1_x ', j, q1_w(j), q1_x(j)
     1081!
     1082    Rn_x(j) = eps_1*Rsigma*T1_x(j)**4 - Rsigma*Ts_x(j)**4
     1083    Rn_w(j) = eps_1*Rsigma*T1_w(j)**4 - Rsigma*Ts_w(j)**4
     1084    Rn_b(j) = sigw(j)*Rn_w(j) + sigx(j)*Rn_x(j)
     1085    dRn(j) = dRn09(j) - ( HTRn_b(j) &
     1086                        +(sigx(j)-sigw(j))*dd_HTRn(j) &
     1087                        -sigw(j)*sigx(j)*dd_HTRn(j)*dd_HTphiT(j)/HTphiT_b(j) &
     1088                       )*(dTs9(j)-dTs09(j)) &
     1089                    + dd_HTRn(j)/HTphiT_b(j)*(phiT_b9(j)-phiT0_b9(j))
     1090!
     1091          print *,'wx_pbl_dts_check, dphiT, L_v*dphiQ, dRn, dTs ', &
     1092              phiT_w9(j)-phiT_x9(j), L_v(j)*(phiQ_w9(j)-phiQ_x9(j)), dRn(j), dTs9(j)
     1093!
     1094  phiQ0_x(j) = PhiQ0_b9(j) - sigw(j)*dphiQ09(j)
     1095  phiQ0_w(j) = PhiQ0_b9(j) + sigx(j)*dphiQ09(j)
     1096!
     1097!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1098!  Test phiQ_w-phiQ0_w = -beta*Kech_Q_sw*dqsatdT_w*(Ts_w-Ts0_w)
     1099!--------------------------------------------------------------
     1100  print *,'wx_pbl_dts_check: beta(j), Kech_Q_sw(j), dqsatdT0_w(j), Ts_w(j), Ts0_w(j) ', &
     1101                         beta(j), Kech_Q_sw(j), dqsatdT0_w(j), Ts_w(j), Ts0_w(j)
     1102  phiQ_w_m_phiQ0_w = -beta(j)*Kech_Q_sw(j)*dqsatdT0_w(j)*(Ts_w(j)-Ts0_w(j))
     1103    print *,'wx_pbl_dts_check: j, phiQ_w9-phiQ0_w, phiQ_w_m_phiQ0_w ', &
     1104                           j, phiQ_w9(j)-phiQ0_w(j), phiQ_w_m_phiQ0_w
     1105!
     1106!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1107!  Test phiQ_x-phiQ0_x = -beta*Kech_Q_sx*dqsatdT_x*(Ts_x-Ts0_x)
     1108!--------------------------------------------------------------
     1109  phiQ_x_m_phiQ0_x = -beta(j)*Kech_Q_sx(j)*dqsatdT0_x(j)*(Ts_x(j)-Ts0_x(j))
     1110    print *,'wx_pbl_dts_check: j, phiQ_x9-phiQ0_x, phiQ_x_m_phiQ0_x ', &
     1111                           j, phiQ_x9(j)-phiQ0_x(j), phiQ_x_m_phiQ0_x
     1112!
     1113!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1114!  Test dphiT-dphiT0 = -(HTphiT_b+(sigx-sigw)*dd_HTphiT)*(dTs-dTs0) - dd_HTphiT*(Ts_b-Ts0_b)
     1115!-------------------------------------------------------------------------------------------
     1116 dphiT = phiT_w9(j) - phiT_x9(j)
     1117 dphiT_m_dphiT0 = -(HTphiT_b(j)+(sigx(j)-sigw(j))*dd_HTphiT(j))*(dTs9(j)-dTs09(j)) &
     1118                  - dd_HTphiT(j)*(Ts_b9(j)-Ts0_b9(j))
     1119 print *,'wx_pbl_dts_check: j, dphiT-dphiT09, dphiT_m_dphiT0 ',j, dphiT-dphiT09(j), dphiT_m_dphiT0
     1120!
     1121!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1122!  Test dphiQ-dphiQ0 = -(HTphiQ_b+(sigx-sigw)*dd_HTphiQ)*(dTs-dTs0) - dd_HTphiQ*(Ts_b-Ts0_b)
     1123!-------------------------------------------------------------------------------------------
     1124 dphiQ = phiQ_w9(j) - phiQ_x9(j)
     1125 dphiQ_m_dphiQ0 = -(HTphiQ_b(j)+(sigx(j)-sigw(j))*dd_HTphiQ(j))*(dTs9(j)-dTs09(j)) &
     1126                  - dd_HTphiQ(j)*(Ts_b9(j)-Ts0_b9(j))
     1127 print *,'wx_pbl_dts_check: j, dphiQ-dphiQ09, dphiQ_m_dphiQ0 ',j, dphiQ-dphiQ09(j), dphiQ_m_dphiQ0
     1128!
     1129!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1130!  Test dRn-dRn0 = -(HTRn_b+(sigx-sigw)*dd_HTRn)*(dTs-dTs0) - dd_HTRn*(Ts_b-Ts0_b)
     1131!-------------------------------------------------------------------------------------------
     1132 dRn_m_dRn0 = -(HTRn_b(j)+(sigx(j)-sigw(j))*dd_HTRn(j))*(dTs9(j)-dTs09(j)) &
     1133                  - dd_HTRn(j)*(Ts_b9(j)-Ts0_b9(j))
     1134 print *,'wx_pbl_dts_check: j, dRn-dRn09, dRn_m_dRn0 ',j, dRn-dRn09(j), dRn_m_dRn0
     1135!
     1136!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1137!  Test phiT_b-phiT0_b = -sigx*sigw*dd_HTphiT*(dTs-dTs0) - HTphiT_b*(Ts_b-Ts0_b)
     1138!-------------------------------------------------------------------------------
     1139   phiTb_m_phiT0b = -sigx(j)*sigw(j)*dd_HTphiT(j)*(dTs9(j)-dTs09(j)) - HTphiT_b(j)*(Ts_b9(j)-Ts0_b9(j))
     1140   print *,'wx_pbl_dts_check: j, phiT_b9-phiT0_b9, phiTb_m_phiT0b ',j ,phiT_b9(j)-phiT0_b9(j), phiTb_m_phiT0b
     1141!
     1142!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1143!  Test phiT_w, phiT_x, dphiT from HTphiT
     1144!------------------------------------------
     1145!  phiT_w = Kech_h_w C_p ( T1_w - Ts_w)                   phiT_x = Kech_h_x C_p ( T1_x - Ts_x)             
     1146  err_phiT_x(j) = Kech_h_x(j)*C_p(j)*(T1_x(j) - Ts_x(j)) - phiT_x9(j)
     1147  err_phiT_w(j) = Kech_h_w(j)*C_p(j)*(T1_w(j) - Ts_w(j)) - phiT_w9(j)
     1148    print *, 'wx_pbl_dts_check: j, phiT_w9, phiT_x9, err_phiT_w, err_phiT_x ',   &
     1149                            j, phiT_w9(j), phiT_x9(j), err_phiT_w(j), err_phiT_x(j)
     1150  dphiT = phiT_w9(j) - phiT_x9(j)
     1151  dphiT_H = dphiT09(j) - ( HTphiT_b(j) &
     1152                            +(sigx(j)-sigw(j))*dd_HTphiT(j) &
     1153                            -sigw(j)*sigx(j)*dd_HTphiT(j)*dd_HTphiT(j)/HTphiT_b(j) &
     1154                           )*(dTs9(j)-dTs09(j)) &
     1155                         + dd_HTphiT(j)/HTphiT_b(j)*(phiT_b9(j)-phiT0_b9(j))
     1156    print *,'wx_pbl_dts_check: j, dphiT, dphiT_H ', j, dphiT, dphiT_H
     1157
     1158!
     1159!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1160!  Test phiq_w, phiq_x, dphiq from HTphiq
     1161!------------------------------------------
     1162!
     1163!  phiq_w = beta Kech_q_w ( q1_w - qsat(Ts_w))            phiq_x = beta Kech_q_x ( q1_x - qsat(Ts_x))
     1164  err_phiq_x(j) = beta(j)*Kech_q_x(j)*( q1_x(j) - qsat_x(j)) - phiq_x9(j)
     1165  err_phiq_w(j) = beta(j)*Kech_q_w(j)*( q1_w(j) - qsat_w(j)) - phiq_w9(j)
     1166  dphiQ = phiQ_w9(j) - phiQ_x9(j)
     1167  dphiQ_H = dphiQ09(j) - ( HTphiQ_b(j) &
     1168                            +(sigx(j)-sigw(j))*dd_HTphiQ(j) &
     1169                            -sigw(j)*sigx(j)*dd_HTphiQ(j)*dd_HTphiT(j)/HTphiT_b(j) &
     1170                           )*(dTs9(j)-dTs09(j)) &
     1171                         + dd_HTphiQ(j)/HTphiT_b(j)*(phiT_b9(j)-phiT0_b9(j))
     1172    print *,'wx_pbl_dts_check: j, dphiQ, dphiQ_H ', j, dphiQ, dphiQ_H
     1173!
     1174!  phiT_b = sigw phiT_w + sigx phiT_x                     dphiT = phiT_w - phiT_x
     1175  err_phiT_b(j) = sigw(j)*phiT_w9(j) + sigx(j)*phiT_x9(j) - phiT_b9(j)
     1176!
     1177!  phiQ_b = sigw phiQ_w + sigx phiQ_x                     dphiQ = phiQ_w - phiQ_x
     1178  err_phiQ_b(j) = sigw(j)*phiQ_w9(j) + sigx(j)*phiQ_x9(j) - phiQ_b9(j)
     1179!
     1180!  Ta = AcoefT + BcoefT phiT_b Delta t
     1181!  phiT_b = Kech_h C_p (Ta - Ts_b)
     1182  Ta(j) = (AcoefT(j) + BcoefT(j)*phiT_b9(j)*dtime) / C_p(j)
     1183  err2_phiT_b(j) = Kech_h(j)*C_p(j)*(Ta(j) - Ts_b9(j)) - phiT_b9(j)
     1184    print *, 'wx_pbl_dts_check: j, Ta, phiT_b9, err2_phiT_b ',   &
     1185                            j, Ta(j), phiT_b9(j), err2_phiT_b(j)
     1186!
     1187          ENDDO  ! j = 1, knon
     1188         
     1189        ENDIF   ! (prt_level >=10 )
     1190!--------------------------------------------------------------------------------------------------
     1191        RETURN
     1192
     1193END SUBROUTINE wx_pbl_dts_check
     1194
     1195SUBROUTINE wx_evappot(knon, q1, Ts, evap_pot)
     1196
     1197    USE wx_pbl_var_mod
     1198
     1199    INTEGER,                      INTENT(IN)        :: knon     ! number of grid cells
     1200    REAL, DIMENSION(knon),        INTENT(IN)        :: q1       ! specific humidity in layer 1
     1201    REAL, DIMENSION(knon),        INTENT(IN)        :: Ts       ! surface temperature
     1202!
     1203    REAL, DIMENSION(knon),        INTENT(OUT)       :: evap_pot ! potential evaporation
     1204!
     1205    INTEGER                   :: j
     1206    REAL                      :: qsat_bs
     1207!
     1208 DO j = 1,knon
     1209   evap_pot(j) = Kech_q(j)*(qsat0(j)+dqsatdT0(j)*(Ts(j)-Ts0(j))-q1(j))
     1210!
     1211  qsat_bs = qsat0(j)+dqsatdT0(j)*(Ts(j)-Ts0(j))
     1212!!  print *,'wx_evappot : Kech_q, qsat_bs, qa, evap_pot ', Kech_q(j), qsat_bs, q1(j), evap_pot(j)
     1213 ENDDO
     1214!
     1215 RETURN
     1216END SUBROUTINE wx_evappot
    4541217
    4551218END MODULE wx_pbl_mod
Note: See TracChangeset for help on using the changeset viewer.