MODULE pbl_surface_subsrf_mod
!
! Planetary Boundary Layer and Surface module
!
! This module manages the calculation of turbulent diffusion in the boundary layer 
! and all interactions towards the differents sub-surfaces.
!
!
  USE pbl_surface_data
  USE dimphy
  USE mod_phys_lmdz_para,  ONLY : mpi_size
  USE mod_grid_phy_lmdz,   ONLY : klon_glo
  USE ioipsl
  USE surface_data,        ONLY : type_ocean, ok_veget, landice_opt, iflag_leads
  USE surf_land_mod,       ONLY : surf_land
  USE surf_landice_mod,    ONLY : surf_landice
  USE surf_ocean_mod,      ONLY : surf_ocean
  USE surf_seaice_mod,     ONLY : surf_seaice
  USE cpl_mod,             ONLY : gath2cpl
  USE climb_hq_mod,        ONLY : climb_hq_down, climb_hq_up
  USE climb_qbs_mod,       ONLY : climb_qbs_down, climb_qbs_up
  USE climb_wind_mod,      ONLY : climb_wind_down, climb_wind_up
  USE coef_diff_turb_mod,  ONLY : coef_diff_turb
  USE lmdz_call_atke,      ONLY : call_atke
  USE ioipsl_getin_p_mod,  ONLY : getin_p
  USE cdrag_mod
  USE stdlevvar_mod
  USE wx_pbl_var_mod,      ONLY : wx_pbl_init, wx_pbl_final, &
                                  wx_pbl_prelim_0, wx_pbl_prelim_beta
  USE wx_pbl_mod,          ONLY : wx_pbl0_merge, wx_pbl_split, wx_pbl_dts_merge, &
                                  wx_pbl_check, wx_pbl_dts_check, wx_evappot
  use config_ocean_skin_m, only: activate_ocean_skin
#ifdef ISO
  USE infotrac_phy, ONLY: niso,ntraciso=>ntiso    
#endif

  IMPLICIT NONE
  PRIVATE
  PUBLIC :: pbl_surface_subsrf
CONTAINS

  SUBROUTINE pbl_surface_subsrf( nsrf, knon, ni,      &
       dtime,     date0,     itap,     jour,          &
       debut,     lafin,                              &
       rlon,      rlat,      rugoro,   rmu0,          &
       lwdown_m,  pphi, cldt,          &
       rain_f,    snow_f,    bs_f,                    &
       gustiness,                                     &
       t,         q,        qbs,  u,        v,        &
       wake_dlt,             wake_dlq,                &
       wake_cstar,           wake_s,                  &
       pplay,     paprs,     pctsrf,                  &
       ts,SFRWL,   alb_dir, alb_dif,ustar, u10m, v10m,wstar, &
       cdragh,    cdragm,                             &
       beta, &
       icesub_ice, icemelt_ice, alb3_lic,  runoff,    snowhgt,   qsnow,     to_ice,    sissnow,  &
       qsat2m,                 &
       d_t,       d_q,    d_qbs,    d_u,      d_v, d_t_diss,            &
       d_t_w,     d_q_w,                             &
       d_t_x,     d_q_x,                             & 
       delta_tsurf,wake_dens,cdragh_x,cdragh_w,      &
       cdragm_x,cdragm_w,kh,kh_x,kh_w,               &
       zcoefh,    zcoefm,    slab_wfbils,            &
       qsol,    s_pblh,         &
       s_pblh_x, s_pblh_w,     &
       delta_qsurf,                         &
       rh2m,                       &
       z0m, z0h,   agesno,  sollw,    solsw,         &
       d_ts,      evap,    fluxlat,   t2m,           &
       flux_t,   flux_u, flux_v,                     &
       dflux_t,   dflux_q,                   &
       q2m, flux_q, flux_qbs, tke_x, eps_x, &
       wake_dltke,                                     &
       treedrg,hice ,tice, bilg_cumul,            &
       fcds, fcdi, dh_basal_growth, dh_basal_melt, &
       dh_top_melt, dh_snow2sic, &
       dtice_melt, dtice_snow2sic , &
       tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &
       cdragm_tersrf, cdragh_tersrf, &
       swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf &
#ifdef ISO
     &   ,xtrain_f, xtsnow_f,xt, &
     &   wake_dlxt,zxxtevap,xtevap, &
     &   d_xt,d_xt_w,d_xt_x, &
     &   xtsol,dflux_xt,zxxtsnow,zxfluxxt,flux_xt, &
     &   h1_diag,runoff_diag,xtrunoff_diag &
#endif      
     , n2mout, n2mout_x, n2mout_w, d_u_x, d_u_w, d_v_x, d_v_w, windsp, t2m_x,       &
       q2m_x, rh2m_x, qsat2m_x, u10m_x, v10m_x, ustar_x, wstar_x, pblh_x, plcl_x, capCL_x,     &
       oliqCL_x, cteiCL_x, pblt_x, therm_x, trmb1_x, trmb2_x, trmb3_x, t2m_w, q2m_w, rh2m_w,   &
       qsat2m_w, u10m_w, v10m_w, ustar_w, wstar_w, pblh_w, plcl_w, capCL_w, oliqCL_w, cteiCL_w,&
       pblt_w, therm_w, trmb1_w, trmb2_w, trmb3_w, pblh, plcl, capCL, oliqCL, cteiCL, pblT, &
       therm, trmb1, trmb2, trmb3, alb, snowerosion, iflag_split_ref, &
       delp, d_t_diss_x, d_t_diss_w, flux_t_x, flux_q_x, flux_t_w, flux_q_w,&
       flux_u_x, flux_v_x, flux_u_w, flux_v_w, fluxlat_x, fluxlat_w)
!$gpum horizontal knon klon


!****************************************************************************************
! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
! Objet: interface de "couche limite" (diffusion verticale)
!
!AA REM:
!AA-----
!AA Tout ce qui a trait au traceurs est dans phytrac maintenant
!AA pour l'instant le calcul de la couche limite pour les traceurs
!AA se fait avec cltrac et ne tient pas compte de la differentiation
!AA des sous-fraction de sol.
!AA REM bis :
!AA----------
!AA Pour pouvoir extraire les coefficient d'echanges et le vent 
!AA dans la premiere couche, 3 champs supplementaires ont ete crees
!AA zcoefh, zu1 et zv1. Pour l'instant nous avons moyenne les valeurs
!AA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir 
!AA si les informations des subsurfaces doivent etre prises en compte
!AA il faudra sortir ces memes champs en leur ajoutant une dimension, 
!AA c'est a dire nbsrf (nbre de subsurface).
!
! Arguments:
!
! dtime----input-R- interval du temps (secondes)
! itap-----input-I- numero du pas de temps
! date0----input-R- jour initial
! t--------input-R- temperature (K)
! q--------input-R- vapeur d'eau (kg/kg)
! u--------input-R- vitesse u
! v--------input-R- vitesse v
! wake_dlt-input-R- temperatre difference between (w) and (x) (K)
! wake_dlq-input-R- humidity difference between (w) and (x) (kg/kg)
!wake_cstar-input-R- wake gust front speed (m/s)
! wake_s---input-R- wake fractionnal area
! ts-------input-R- temperature du sol (en Kelvin)
! paprs----input-R- pression a intercouche (Pa)
! pplay----input-R- pression au milieu de couche (Pa)
! rlat-----input-R- latitude en degree
! z0m, z0h ----input-R- longeur de rugosite (en m)
! Martin
! cldt-----input-R- total cloud fraction
! Martin
!GG
! pphi-----input-R- geopotentiel de chaque couche (g z) (reference sol)
!GG
!
! d_t------output-R- le changement pour "t"
! d_q------output-R- le changement pour "q"
! d_u------output-R- le changement pour "u"
! d_v------output-R- le changement pour "v"
! d_ts-----output-R- le changement pour "ts"
! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
!                    (orientation positive vers le bas)
! tke_x---input/output-R- tke in the (x) region (kg/m**2/s)
! wake_dltke-input/output-R- tke difference between (w) and (x) (kg/m**2/s)
! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
! dflux_t--output-R- derive du flux sensible
! dflux_q--output-R- derive du flux latent
! zu1------output-R- le vent dans la premiere couche
! zv1------output-R- le vent dans la premiere couche
! trmb1----output-R- deep_cape
! trmb2----output-R- inhibition 
! trmb3----output-R- Point Omega
! cteiCL---output-R- Critere d'instab d'entrainmt des nuages de CL
! plcl-----output-R- Niveau de condensation
! pblh-----output-R- HCL
! pblT-----output-R- T au nveau HCL
! treedrg--output-R- tree drag (m)               
! qsurf_tersrf--output-R- surface specific humidity of continental sub-surfaces
! cdragm_tersrf--output-R- momentum drag coefficient of continental sub-surfaces
! cdragh_tersrf--output-R- heat drag coefficient of continental sub-surfaces
! tsurf_new_tersrf--output-R- surface temperature of continental sub-surfaces
! swnet_tersrf--output-R- net shortwave radiation of continental sub-surfaces
! lwnet_tersrf--output-R- net longwave radiation of continental sub-surfaces
! fluxsens_tersrf--output-R- sensible heat flux of continental sub-surfaces
! fluxlat_tersrf--output-R- latent heat flux of continental sub-surfaces
    USE carbon_cycle_mod,   ONLY : carbon_cycle_cpl  
    USE carbon_cycle_mod,   ONLY : co2_send, nbcf_out, fields_out, yfields_out, yfields_in
    use hbtm_mod, only: hbtm
    USE indice_sol_mod
    USE mod_grid_phy_lmdz,  ONLY : grid1dto2d_glo
    USE print_control_mod,  ONLY : prt_level,lunout
#ifdef ISO
  USE isotopes_mod, ONLY: Rdefault,iso_eau
#ifdef ISOVERIF
        USE isotopes_verif_mod
#endif
#ifdef ISOTRAC
        USE isotrac_mod, only: index_iso
#endif
#endif
USE dimpft_mod_h
    USE flux_arp_mod_h
    USE compbl_mod_h
    USE yoethf_mod_h
        USE clesphys_mod_h
    USE ioipsl_getin_p_mod, ONLY : getin_p
    use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal, dter, &
         dser, dt_ds, zsig, zmea, &
         frac_tersrf, z0m_tersrf, ratio_z0m_z0h_tersrf !AM 
    use phys_output_var_mod, only: tkt, tks, taur, sss
    use lmdz_blosno_ini, only : zeta_bs
    USE dimsoil_mod_h, ONLY: nsoilmx
    USE surf_param_mod, ONLY: eff_surf_param  !AM
    use wxios_mod, ONLY: missing_val_xios => missing_val, using_xios
    USE netcdf, only: missing_val_netcdf => nf90_fill_real
    USE yomcst_mod_h
    USE mod_phys_lmdz_para, ONLY : is_master
    USE cdrag_mod, ONLY : cdrag
    USE freinage_mod, ONLY : freinage
    USE yamada_c_mod, ONLY : yamada_c

IMPLICIT NONE

    INCLUDE "FCTTRE.h"
!****************************************************************************************
    INTEGER,                      INTENT(IN)        :: nsrf    ! indice current subsurface
    INTEGER,                      INTENT(IN)        :: knon    ! number of compressed points for current subsurface
    INTEGER,                      INTENT(IN)        :: ni(knon)! index for compressed current sub-surface
    REAL,                         INTENT(IN)        :: dtime   ! time interval (s)
    REAL,                         INTENT(IN)        :: date0   ! initial day
    INTEGER,                      INTENT(IN)        :: itap    ! time step
    INTEGER,                      INTENT(IN)        :: jour    ! current day of the year
    LOGICAL,                      INTENT(IN)        :: debut   ! true if first run step
    LOGICAL,                      INTENT(IN)        :: lafin   ! true if last run step
    REAL, DIMENSION(klon),        INTENT(IN)        :: rlon    ! longitudes in degrees
    REAL, DIMENSION(klon),        INTENT(IN)        :: rlat    ! latitudes in degrees
    REAL, DIMENSION(klon),        INTENT(IN)        :: rugoro  ! rugosity length
    REAL, DIMENSION(klon),        INTENT(IN)        :: rmu0    ! cosine of solar zenith angle
    REAL, DIMENSION(klon),        INTENT(IN)        :: rain_f  ! rain fall
    REAL, DIMENSION(klon),        INTENT(IN)        :: snow_f  ! snow fall
    REAL, DIMENSION(klon),        INTENT(IN)        :: bs_f  ! blowing snow fall
    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: t       ! temperature (K)
    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: q       ! water vapour (kg/kg)
    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: qbs       ! blowing snow specific content (kg/kg)
    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: u       ! u speed
    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: v       ! v speed
    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: pplay   ! mid-layer pression (Pa)
    REAL, DIMENSION(klon,klev+1), INTENT(IN)        :: paprs   ! pression between layers (Pa) 
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)        :: pctsrf  ! sub-surface fraction
    REAL, DIMENSION(klon),        INTENT(IN)        :: lwdown_m ! downward longwave radiation at mean s    
    REAL, DIMENSION(klon),        INTENT(IN)        :: gustiness ! gustiness
    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: pphi    ! geopotential (m2/s2)
    REAL, DIMENSION(klon),        INTENT(IN)        :: cldt    ! total cloud 

#ifdef ISO
    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(IN)        :: xt       ! water vapour (kg/kg)
    REAL, DIMENSION(ntraciso,klon),        INTENT(IN)        :: xtrain_f  ! rain fall
    REAL, DIMENSION(ntraciso,klon),        INTENT(IN)        :: xtsnow_f  ! snow fall
#endif
    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: wake_dlt  !temperature difference between (w) and (x) (K)
    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: wake_dlq  !humidity difference between (w) and (x) (K)
    REAL, DIMENSION(klon),        INTENT(IN)        :: wake_s    ! Fraction de poches froides
    REAL, DIMENSION(klon),        INTENT(IN)        :: wake_cstar! Vitesse d'expansion des poches froides
    REAL, DIMENSION(klon),        INTENT(IN)        :: wake_dens
#ifdef ISO
    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(IN)        :: wake_dlxt   
#endif
! Input/Output variables
!****************************************************************************************
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: beta    ! Aridity factor
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ts      ! temperature at surface (K)
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: delta_tsurf !surface temperature difference between
    REAL, DIMENSIOn(6),intent(in) :: SFRWL
    REAL, DIMENSION(klon, nsw, nbsrf), INTENT(INOUT)     :: alb_dir,alb_dif
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ustar   ! u* (m/s)
    REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT)   :: wstar   ! w* (m/s)
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: u10m    ! u speed at 10m
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: v10m    ! v speed at 10m
    REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke_x
    REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: wake_dltke ! TKE_w - TKE_x

! Output variables
!****************************************************************************************
    REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT)   :: eps_x      ! TKE dissipation rate

    REAL, DIMENSION(klon),        INTENT(INOUT)       :: cdragh     ! drag coefficient for T and Q
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: cdragm     ! drag coefficient for wind
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: alb3_lic
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: icesub_ice ! ice (no snow!) sublimation flux over ice sheet and sea ice 
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: icemelt_ice ! ice (no snow!) melting flux over ice sheet and sea ice 
    REAL, DIMENSION(klon,klev),   INTENT(INOUT)       :: d_t_w      !   !
    REAL, DIMENSION(klon,klev),   INTENT(INOUT)       :: d_q_w      !      !  Tendances dans les poches
    REAL, DIMENSION(klon,klev),   INTENT(INOUT)       :: d_t_x      !   !
    REAL, DIMENSION(klon,klev),   INTENT(INOUT)       :: d_q_x      !      !  Tendances hors des poches
    REAL, DIMENSION(klon),        INTENT(INOUT)     :: qsat2m
    REAL, DIMENSION(klon, klev),  INTENT(INOUT)       :: d_t        ! change in temperature 
    REAL, DIMENSION(klon, klev),  INTENT(INOUT)     :: d_t_diss       ! change in temperature 
    REAL, DIMENSION(klon, klev),  INTENT(INOUT)       :: d_q        ! change in water vapour
    REAL, DIMENSION(klon, klev),  INTENT(INOUT)       :: d_u        ! change in u speed
    REAL, DIMENSION(klon, klev),  INTENT(INOUT)       :: d_v        ! change in v speed
    REAL, DIMENSION(klon, klev),  INTENT(INOUT)       :: d_qbs        ! change in blowing snow specific content

    REAL, INTENT(INOUT):: zcoefh(klon, klev+1, nbsrf + 1) ! (klon, klev, nbsrf + 1) => only use (klon, klev, nbsrf+1)
    ! coef for turbulent diffusion of T and Q, mean for each grid point

    REAL, INTENT(INOUT):: zcoefm(klon, klev+1, nbsrf + 1) ! (klon, klev, nbsrf + 1) => only use (klon, klev, nbsrf+1)
    ! coef for turbulent diffusion of U and V (?), mean for each grid point
#ifdef ISO
    REAL, DIMENSION(ntraciso,klon),        INTENT(INOUT)       :: zxxtevap     ! water vapour flux at surface, positiv upwards
    REAL, DIMENSION(ntraciso,klon, klev),  INTENT(INOUT)       :: d_xt        ! change in water vapour
    REAL, DIMENSION(klon),                 INTENT(INOUT)       :: runoff_diag
    REAL, DIMENSION(niso,klon),            INTENT(INOUT)       :: xtrunoff_diag
    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(INOUT)       :: d_xt_w
    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(INOUT)       :: d_xt_x
#endif
! Output only for diagnostics
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: cdragh_x
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: cdragh_w
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: cdragm_x
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: cdragm_w
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: kh
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: kh_x
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: kh_w
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: slab_wfbils! heat balance at surface only for slab at ocean points
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: qsol     ! water height in the soil (mm)
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: s_pblh     ! height of the planetary boundary layer(HPBL)
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: s_pblh_x   ! height of the PBL in the off-wake region
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: s_pblh_w   ! height of the PBL in the wake region
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: delta_qsurf! humidity difference at surface, mean for each grid point
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: rh2m       ! relative humidity at 2m
    REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT)   :: z0m,z0h      ! rugosity length (m)
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: agesno   ! age of snow at surface
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: solsw      ! net shortwave radiation at surface 
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: sollw      ! net longwave radiation at surface
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: d_ts       ! change in temperature at surface
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: evap       ! evaporation at surface
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: fluxlat    ! latent flux
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: t2m        ! temperature at 2 meter height
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: flux_t     ! sensible heat flux (CpT) J/m**2/s (W/m**2)
                                                                  ! positve orientation downwards
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: flux_u     ! u wind tension (kg m/s)/(m**2 s) or Pascal
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: flux_v     ! v wind tension (kg m/s)/(m**2 s) or Pascal
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: treedrg  ! tree drag (m)     
!AM heterogeneous continental sub-surfaces
    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: tsurf_tersrf     ! surface temperature of continental sub-surfaces (K)               
    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: qsurf_tersrf     ! surface specific humidity of continental sub-surfaces (kg/kg)               
    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: tsurf_new_tersrf ! surface temperature of continental sub-surfaces (K)               
    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: cdragm_tersrf    ! momentum drag coefficient of continental sub-surfaces (-)               
    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: cdragh_tersrf    ! heat drag coefficient of continental sub-surfaces (-)               
    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: swnet_tersrf     ! net shortwave radiation of continental sub-surfaces (W/m2)               
    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: lwnet_tersrf     ! net longwave radiation of continental sub-surfaces (W/m2)               
    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: fluxsens_tersrf  ! sensible heat flux of continental sub-surfaces (W/m2)               
    REAL, DIMENSION(klon, nbtersrf), INTENT(INOUT) :: fluxlat_tersrf   ! latent heat flux of continental sub-surfaces (W/m2)               
    REAL, DIMENSION(klon, nsoilmx, nbtersrf), INTENT(INOUT) :: tsoil_tersrf ! soil temperature of continental sub-surfaces (K)               
#ifdef ISO        
    REAL, DIMENSION(niso,klon),   INTENT(INOUT)       :: xtsol      ! water height in the soil (mm)
    REAL, DIMENSION(ntraciso,klon, nbsrf)           :: xtevap     ! evaporation at surface
    REAL, DIMENSION(klon),        INTENT(INOUT)       :: h1_diag    ! just diagnostic, not useful
#endif

! Output not needed
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: dflux_t    ! change of sensible heat flux 
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: dflux_q    ! change of water vapour flux
    REAL, DIMENSION(klon, nbsrf),INTENT(INOUT)        :: q2m        ! water vapour at 2 meter height
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: flux_q     ! water vapour flux(latent flux) (kg/m**2/s)
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: flux_qbs   ! blowind snow vertical flux (kg/m**2

#ifdef ISO   
    REAL, DIMENSION(ntraciso,klon),              INTENT(INOUT) :: dflux_xt    ! change of water vapour flux
    REAL, DIMENSION(niso,klon),                  INTENT(INOUT) :: zxxtsnow    ! snow at surface, mean for each grid point
    REAL, DIMENSION(ntraciso,klon, klev),        INTENT(INOUT) :: zxfluxxt    ! water vapour flux, mean for each grid point 
    REAL, DIMENSION(ntraciso,klon, klev, nbsrf), INTENT(INOUT) :: flux_xt     ! water vapour flux(latent flux) (kg/m**2/s)  
#endif

    REAL, DIMENSION(klon),       INTENT(INOUT)        :: qsnow      ! snow water content
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: snowhgt    ! snow height
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: to_ice     ! snow passed to ice
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: sissnow    ! snow in snow model
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: runoff     ! runoff on land ice
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: hice      ! hice
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: tice      ! tice
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: bilg_cumul      ! flux cumulated
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: fcds
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: fcdi
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: dh_basal_growth
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: dh_basal_melt
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: dh_top_melt
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: dh_snow2sic
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: dtice_melt
    REAL, DIMENSION(klon),       INTENT(INOUT)        :: dtice_snow2sic

! variables temporaires en "klon" (nom compressée) passée en argument pour les sous-surface

    INTEGER, DIMENSION(klon, nbsrf, 6), INTENT(INOUT) :: n2mout
    INTEGER, DIMENSION(klon, nbsrf, 6), INTENT(INOUT) :: n2mout_x
    INTEGER, DIMENSION(klon, nbsrf, 6), INTENT(INOUT) :: n2mout_w
    REAL, DIMENSION(klon, klev), INTENT(INOUT)        :: d_u_x
    REAL, DIMENSION(klon, klev), INTENT(INOUT)        :: d_u_w
    REAL, DIMENSION(klon, klev), INTENT(INOUT)        :: d_v_x
    REAL, DIMENSION(klon, klev), INTENT(INOUT)        :: d_v_w 
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: windsp
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: t2m_x
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: q2m_x
    REAL, DIMENSION(klon), INTENT(INOUT)              :: rh2m_x
    REAL, DIMENSION(klon), INTENT(INOUT)              :: qsat2m_x
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: u10m_x
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: v10m_x
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: ustar_x
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: wstar_x
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: pblh_x
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: plcl_x
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: capCL_x
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: oliqCL_x
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: cteiCL_x
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: pblt_x
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: therm_x
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: trmb1_x
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: trmb2_x
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: trmb3_x
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: t2m_w
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: q2m_w
    REAL, DIMENSION(klon), INTENT(INOUT)              :: rh2m_w
    REAL, DIMENSION(klon), INTENT(INOUT)              :: qsat2m_w
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: u10m_w
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: v10m_w
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: ustar_w
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: wstar_w
!                           
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: pblh_w
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: plcl_w
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: capCL_w
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: oliqCL_w
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: cteiCL_w
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: pblt_w
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: therm_w
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: trmb1_w
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: trmb2_w
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: trmb3_w
!
    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: pblh         ! height of the planetary boundary layer
    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: plcl         ! condensation level
    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: capCL
    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: oliqCL
    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: cteiCL
    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: pblT
    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: therm
    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: trmb1        ! deep cape
    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: trmb2        ! inhibition
    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: trmb3        ! point Omega
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: alb          ! mean albedo for whole SW interval
    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: snowerosion    
    INTEGER,                     INTENT(INOUT)      :: iflag_split_ref 

!
! Other local variables
!****************************************************************************************
    INTEGER                            :: n 
    INTEGER                            :: iflag_split
    INTEGER                            :: i, k 
    INTEGER                            :: j
    REAL, DIMENSION(knon)              :: r_co2_ppm     ! taux CO2 atmosphere
    REAL                               :: yt1_new
    REAL, DIMENSION(knon)              :: yts, yz0m, yz0h, ypct
    REAL, DIMENSION(knon)              :: yz0h_old
    REAL, DIMENSION(knon)              :: yalb,yalb_vis
    REAL, DIMENSION(knon)              :: yt1, yq1, yu1, yv1, yqbs1
    REAL, DIMENSION(knon)              :: yqa
    REAL, DIMENSION(knon)              :: ysnow, yqsurf, yagesno, yqsol
    REAL, DIMENSION(knon)              :: yrain_f, ysnow_f, ybs_f
#ifdef ISO
    REAL, DIMENSION(ntraciso,knon)     :: yxt1
    REAL, DIMENSION(niso,knon)         :: yxtsnow, yxtsol   
    REAL, DIMENSION(ntraciso,knon)     :: yxtrain_f, yxtsnow_f 
    REAL, DIMENSION(knon)              :: yrunoff_diag
    REAL, DIMENSION(niso,knon)         :: yxtrunoff_diag
    REAL, DIMENSION(niso,knon)         :: yRland_ice    
#endif
    REAL, DIMENSION(knon)              :: ysolsw, ysollw
    REAL, DIMENSION(knon)              :: yfder
    REAL, DIMENSION(knon)              :: yrugoro
    REAL, DIMENSION(knon)              :: yfluxlat
    REAL, DIMENSION(knon)              :: yfluxbs
    REAL, DIMENSION(knon)              :: y_d_ts
    REAL, DIMENSION(knon)              :: y_flux_t1, y_flux_q1
    REAL, DIMENSION(knon)              :: y_dflux_t, y_dflux_q
#ifdef ISO
    REAL, DIMENSION(ntraciso,knon)     ::  y_flux_xt1
    REAL, DIMENSION(ntraciso,knon)     ::  y_dflux_xt
#endif
    REAL, DIMENSION(knon)              :: y_flux_u1, y_flux_v1
    REAL, DIMENSION(knon)              :: y_flux_bs, y_flux0
    REAL, DIMENSION(knon)              :: yt2m, yq2m, yu10m
    INTEGER, DIMENSION(knon, nbsrf, 6) :: yn2mout, yn2mout_x, yn2mout_w
    REAL, DIMENSION(knon)              :: yustar
    REAL, DIMENSION(knon)              :: ywstar
    REAL, DIMENSION(knon)              :: ywindsp
    REAL, DIMENSION(knon)              :: yt10m, yq10m
    REAL, DIMENSION(knon)              :: ypblh
    REAL, DIMENSION(knon)              :: ylcl
    REAL, DIMENSION(knon)              :: ycapCL
    REAL, DIMENSION(knon)              :: yoliqCL
    REAL, DIMENSION(knon)              :: ycteiCL
    REAL, DIMENSION(knon)              :: ypblT
    REAL, DIMENSION(knon)              :: ytherm
    REAL, DIMENSION(knon)              :: ytrmb1
    REAL, DIMENSION(knon)              :: ytrmb2
    REAL, DIMENSION(knon)              :: ytrmb3

    REAL, DIMENSION(knon)       :: yt2m_x
    REAL, DIMENSION(knon)       :: yq2m_x
    REAL, DIMENSION(knon)       :: yt10m_x
    REAL, DIMENSION(knon)       :: yq10m_x
    REAL, DIMENSION(knon)       :: yu10m_x
    REAL, DIMENSION(knon)       :: yustar_x
    REAL, DIMENSION(knon)       :: ywstar_x
!              
    REAL, DIMENSION(knon)       :: ypblh_x
    REAL, DIMENSION(knon)       :: ylcl_x
    REAL, DIMENSION(knon)       :: ycapCL_x
    REAL, DIMENSION(knon)       :: yoliqCL_x
    REAL, DIMENSION(knon)       :: ycteiCL_x
    REAL, DIMENSION(knon)       :: ypblt_x
    REAL, DIMENSION(knon)       :: ytherm_x
    REAL, DIMENSION(knon)       :: ytrmb1_x
    REAL, DIMENSION(knon)       :: ytrmb2_x
    REAL, DIMENSION(knon)       :: ytrmb3_x

    REAL, DIMENSION(knon)              :: uzon, vmer
    REAL, DIMENSION(knon)              :: tair1, qair1, tairsol
    REAL, DIMENSION(knon)              :: psfce, patm
    REAL, DIMENSION(knon)              :: qairsol, zgeo1, speed, zri1, pref !speed, zri1, pref, added by Fuxing WANG, 04/03/2015
    REAL, DIMENSION(knon)              :: yz0h_oupas
    REAL, DIMENSION(knon)              :: yfluxsens
    REAL, DIMENSION(knon)              :: AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0
    REAL, DIMENSION(knon)              :: AcoefH, AcoefQ, BcoefH, BcoefQ
#ifdef ISO
    REAL, DIMENSION(ntraciso,knon)     :: AcoefXT, BcoefXT
#endif
    REAL, DIMENSION(knon)              :: AcoefU, AcoefV, BcoefU, BcoefV
    REAL, DIMENSION(knon)              :: AcoefQBS, BcoefQBS
    REAL, DIMENSION(knon)              :: ypsref
    REAL, DIMENSION(knon)              :: yevap, yevap_pot, ytsurf_new, yalb3_new, yicesub, yicemelt
    REAL, DIMENSION(knon,nsw)          :: yalb_dir_new, yalb_dif_new
    REAL, DIMENSION(knon,klev)         :: y_d_t, y_d_q, y_d_t_diss, y_d_qbs
    REAL, DIMENSION(knon,klev)         :: y_d_u, y_d_v
    REAL, DIMENSION(knon,klev)         :: y_flux_t, y_flux_q, y_flux_qbs
    REAL, DIMENSION(knon,klev)         :: y_flux_u, y_flux_v
    REAL, DIMENSION(knon,klev)         :: ycoefh,ycoefm,ycoefq,ycoefqbs
    REAL, DIMENSION(knon)              :: ycdragh, ycdragq, ycdragm
    REAL, DIMENSION(knon,klev)         :: yu, yv
    REAL, DIMENSION(knon,klev)         :: yt, yq, yqbs
#ifdef ISO
    REAL, DIMENSION(ntraciso,knon)      :: yxtevap
    REAL, DIMENSION(ntraciso,knon,klev) :: y_d_xt
    REAL, DIMENSION(ntraciso,knon,klev) :: y_flux_xt
    REAL, DIMENSION(ntraciso,knon,klev) :: yxt   
#endif
    REAL, DIMENSION(knon,klev)         :: ypplay, ydelp
    REAL, DIMENSION(klon,klev),INTENT(IN)         :: delp
    REAL, DIMENSION(knon,klev+1)       :: ypaprs
    REAL, DIMENSION(knon,klev+1)       :: ytke, yeps
    REAL, DIMENSION(knon,nsoilmx)      :: ytsoil
    REAL, DIMENSION(knon,nvm_lmdz)          :: yveget
    REAL, DIMENSION(knon,nvm_lmdz)          :: ylai
    REAL, DIMENSION(knon,nvm_lmdz)          :: yheight
    REAL, DIMENSION(knon,klev)              :: y_d_u_frein
    REAL, DIMENSION(knon,klev)              :: y_d_v_frein
    REAL, DIMENSION(knon,klev)              :: y_treedrg

    CHARACTER(len=80)                  :: abort_message
    CHARACTER(len=20)                  :: modname = 'pbl_surface'
    LOGICAL, PARAMETER                 :: zxli=.FALSE. ! utiliser un jeu de fonctions simples
    LOGICAL, PARAMETER                 :: check=.FALSE.

    REAL, DIMENSION(knon)              :: ywake_s, ywake_cstar, ywake_dens
    REAL, DIMENSION(knon,klev+1)       :: ytke_x, ytke_w, yeps_x, yeps_w
    REAL, DIMENSION(knon,klev+1)       :: ywake_dltke
    REAL, DIMENSION(knon,klev)         :: yu_x, yv_x, yu_w, yv_w
    REAL, DIMENSION(knon,klev)         :: yt_x, yq_x, yt_w, yq_w
    REAL, DIMENSION(knon,klev)         :: ycoefh_x, ycoefm_x, ycoefh_w, ycoefm_w
    REAL, DIMENSION(knon,klev)         :: ycoefq_x, ycoefq_w
    REAL, DIMENSION(knon)              :: ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w
    REAL, DIMENSION(knon)              :: ycdragm_x, ycdragm_w
    REAL, DIMENSION(knon)              :: AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x
    REAL, DIMENSION(knon)              :: AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w
    REAL, DIMENSION(knon)              :: AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x
    REAL, DIMENSION(knon)              :: AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w
    REAL, DIMENSION(knon)              :: y_flux_t1_x, y_flux_q1_x, y_flux_t1_w, y_flux_q1_w
    REAL, DIMENSION(knon)              :: y_flux_u1_x, y_flux_v1_x, y_flux_u1_w, y_flux_v1_w
    REAL, DIMENSION(knon,klev)         :: y_flux_t_x, y_flux_q_x, y_flux_t_w, y_flux_q_w
    REAL, DIMENSION(knon,klev)         :: y_flux_u_x, y_flux_v_x, y_flux_u_w, y_flux_v_w
    REAL, DIMENSION(knon)              :: yfluxlat_x, yfluxlat_w
    REAL, DIMENSION(knon,klev)         :: y_d_t_x, y_d_q_x, y_d_t_w, y_d_q_w
    REAL, DIMENSION(knon,klev)         :: y_d_t_diss_x, y_d_t_diss_w
    REAL, DIMENSION(klon,klev), INTENT(INOUT)         :: d_t_diss_x, d_t_diss_w
    REAL, DIMENSION(knon,klev)         :: y_d_u_x, y_d_v_x, y_d_u_w, y_d_v_w
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: flux_t_x, flux_q_x, flux_t_w, flux_q_w
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: flux_u_x, flux_v_x, flux_u_w, flux_v_w
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: fluxlat_x, fluxlat_w
    REAL, DIMENSION(knon)              :: ybeta
    REAL, DIMENSION(knon)              :: ybeta_prev
    REAL, DIMENSION(knon,klev)         :: CcoefH, CcoefQ, DcoefH, DcoefQ
    REAL, DIMENSION(knon,klev)         :: CcoefU, CcoefV, DcoefU, DcoefV
    REAL, DIMENSION(knon,klev)         :: CcoefQBS, DcoefQBS
    REAL, DIMENSION(knon,klev)         :: CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x
    REAL, DIMENSION(knon,klev)         :: CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w
    REAL, DIMENSION(knon,klev)         :: CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x
    REAL, DIMENSION(knon,klev)         :: CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w
    REAL, DIMENSION(knon,klev)         :: Kcoef_hq, Kcoef_m, gama_h, gama_q
    REAL, DIMENSION(knon,klev)         :: gama_qbs, Kcoef_qbs 
    REAL, DIMENSION(knon,klev)         :: Kcoef_hq_x, Kcoef_m_x, gama_h_x, gama_q_x
    REAL, DIMENSION(knon,klev)         :: Kcoef_hq_w, Kcoef_m_w, gama_h_w, gama_q_w
    REAL, DIMENSION(knon)              :: alf_1, alf_2, alf_1_x, alf_2_x, alf_1_w, alf_2_w
#ifdef ISO
    REAL, DIMENSION(ntraciso,knon,klev)         :: yxt_x, yxt_w
    REAL, DIMENSION(ntraciso,knon)              :: y_flux_xt1_x , y_flux_xt1_w   
    REAL, DIMENSION(ntraciso,knon,klev)         :: y_flux_xt_x,y_d_xt_x
    REAL, DIMENSION(ntraciso,knon,klev)         :: y_flux_xt_w,y_d_xt_w
    REAL, DIMENSION(ntraciso,klon,klev),INTENT(INOUT)    :: zxfluxxt_w, zxfluxxt_x
    REAL, DIMENSION(ntraciso,klon,klev,nbsrf), INTENT(INOUT)   :: flux_xt_x, flux_xt_w
    REAL, DIMENSION(ntraciso,knon)              :: AcoefXT_x, BcoefXT_x
    REAL, DIMENSION(ntraciso,knon)              :: AcoefXT_w, BcoefXT_w
    REAL, DIMENSION(ntraciso,knon,klev)         :: CcoefXT, DcoefXT
    REAL, DIMENSION(ntraciso,knon,klev)         :: CcoefXT_x, DcoefXT_x
    REAL, DIMENSION(ntraciso,knon,klev)         :: CcoefXT_w, DcoefXT_w
    REAL, DIMENSION(ntraciso,knon,klev)         :: gama_xt,gama_xt_x,gama_xt_w
#endif

    REAL, DIMENSION(knon)       :: yt2m_w
    REAL, DIMENSION(knon)       :: yq2m_w
    REAL, DIMENSION(knon)       :: yt10m_w
    REAL, DIMENSION(knon)       :: yq10m_w
    REAL, DIMENSION(knon)       :: yu10m_w
    REAL, DIMENSION(knon)       :: yustar_w
    REAL, DIMENSION(knon)       :: ywstar_w
!                       
    REAL, DIMENSION(knon)       :: ypblh_w
    REAL, DIMENSION(knon)       :: ylcl_w
    REAL, DIMENSION(knon)       :: ycapCL_w
    REAL, DIMENSION(knon)       :: yoliqCL_w
    REAL, DIMENSION(knon)       :: ycteiCL_w
    REAL, DIMENSION(knon)       :: ypblt_w
    REAL, DIMENSION(knon)       :: ytherm_w
    REAL, DIMENSION(knon)       :: ytrmb1_w
    REAL, DIMENSION(knon)       :: ytrmb2_w
    REAL, DIMENSION(knon)       :: ytrmb3_w
!
    REAL, DIMENSION(knon)       :: uzon_x, vmer_x, speed_x, zri1_x, pref_x !speed_x, zri1_x, pref_x, added by Fuxing WANG, 04/03/2015
    REAL, DIMENSION(knon)       :: zgeo1_x, tair1_x, qair1_x, tairsol_x
!
    REAL, DIMENSION(knon)       :: uzon_w, vmer_w, speed_w, zri1_w, pref_w !speed_w, zri1_w, pref_w, added by Fuxing WANG, 04/03/2015
    REAL, DIMENSION(knon)       :: zgeo1_w, tair1_w, qair1_w, tairsol_w
    REAL, DIMENSION(knon)       :: yus0, yvs0
!
    REAL, DIMENSION(knon)              :: y_delta_flux_t1
    REAL, DIMENSION(knon)              :: y_delta_tsurf, y_delta_tsurf_new
    REAL, DIMENSION(knon)              :: delta_coef, tau_eq
    REAL, DIMENSION(knon)              :: HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn
    REAL, DIMENSION(knon)              :: phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0
    REAL, DIMENSION(knon)              :: y_delta_qsurf
    REAL, DIMENSION(knon)              :: y_delta_qsats
    REAL, DIMENSION(knon)              :: yg_T, yg_Q
    REAL, DIMENSION(knon)              :: yGamma_dTs_phiT, yGamma_dQs_phiQ
    REAL, DIMENSION(knon)              :: ydTs_ins, ydqs_ins
!
    REAL, PARAMETER                    :: facteur = 2. / 1.772  ! ( == 2. / SQRT(3.14))
    REAL, PARAMETER                    :: inertia=2000.
    REAL, DIMENSION(knon)              :: ydtsurf_th
    REAL, DIMENSION(knon)              :: Kech_h           ! Coefficient d'echange pour l'energie
    REAL, DIMENSION(knon)              :: Kech_h_x, Kech_h_w 
    REAL, DIMENSION(knon)              :: yts_x, yts_w
    REAL, DIMENSION(knon)              :: yqsurf_x, yqsurf_w
    REAL                               :: fact_cdrag
    REAL                               :: z1lay
    REAL                               :: vent
    REAL, DIMENSION(knon)              :: ylwdown      ! jg : temporary (ysollwdown)
    REAL, DIMENSION(knon)              :: ygustiness      ! jg : temporary (ysollwdown)
    REAL                               :: zx_qs1, zcor1, zdelta1 
    REAL, DIMENSION(knon)              :: ytoice
    REAL, DIMENSION(knon)              :: ysnowhgt, yqsnow, ysissnow, yrunoff
    REAL, DIMENSION(knon)              :: yzmea
    REAL, DIMENSION(knon)              :: yzsig
    REAL, DIMENSION(knon)              :: ycldt
    REAL, DIMENSION(knon)              :: yrmu0
    REAL, DIMENSION(knon)              :: yri0

    REAL, DIMENSION(knon):: ydelta_sst, ydelta_sal, yds_ns, ydt_ns, ydter, &
         ydser, ydt_ds, ytkt, ytks, ytaur, ysss
    ! compression of delta_sst, delta_sal, ds_ns, dt_ns, dter, dser,
    ! dt_ds, tkt, tks, taur, sss on ocean points
    REAL :: missing_val

    REAL, DIMENSION(knon,klev)         :: ytheta
    REAL, DIMENSION(knon,klev)         :: ypphii
    REAL, DIMENSION(knon,klev)         :: ypphi
    REAL, DIMENSION(knon,klev)         :: ydthetadz
    REAL, DIMENSION(knon)              :: ydthetadz300
    REAL, DIMENSION(knon)              :: Ampl
    REAL, DIMENSION(knon, nbtersrf) :: yfrac_tersrf, yz0m_tersrf, yz0h_tersrf 
    REAL, DIMENSION(knon) :: yzxtsol     ! temperature at surface
    REAL, DIMENSION(knon)                   :: ypblh_tmp ! temporaire pblh compressed 
    REAL, DIMENSION(knon) :: zgeo1_over_RG
#ifdef ISO
    REAL, DIMENSION(knon)       :: h1
    INTEGER                     :: ixt
#endif

      IF (using_xios) THEN
        missing_val=missing_val_xios
      ELSE
        missing_val=missing_val_netcdf
      ENDIF

     yus0(:)=0. ; yvs0(:)=0.

!    loop_nbsrf: DO nsrf = 1, nbsrf                                        !<<<<<<<<<<<<<
                                                                          !<<<<<<<<<<<<<
       IF (prt_level >=10) print *,' Loop nsrf ',nsrf
!
       IF (iflag_split_ref == 3) THEN
         IF (nsrf == is_oce) THEN
            iflag_split = 1
         ELSE
            iflag_split=0
         ENDIF   !! (nsrf == is_oce)
       ELSE                     
         iflag_split = iflag_split_ref
       ENDIF   !! (iflag_split_ref == 3)

! Search for index(ni) and size(knon) of domaine to treat
!       ni(:) = 0
!       knon  = 0
!       DO i = 1, klon
!          IF (pctsrf(i,nsrf) > 0.) THEN
!             knon = knon + 1
!             ni(knon) = i
!          ENDIF
!       ENDDO

!!! jyg le 19/08/2012
!       IF (knon <= 0) THEN
!         IF (prt_level >= 10) print *,' no grid point for nsrf= ',nsrf
!         cycle loop_nbsrf
!       ENDIF

!!!
! 2b) Initialization of all local variables that will be compressed later
!****************************************************************************************

    ypct = 0.0    ; yts = 0.0        ; ysnow = 0.0
    yqsurf = 0.0  ; yalb = 0.0 ; yalb_vis = 0.0
    yrain_f = 0.0 ; ysnow_f = 0.0  ; ybs_f=0.0  ; yfder = 0.0     ; ysolsw = 0.0
    ysollw = 0.0  ; yz0m = 0.0 ; yz0h = 0.0    ; yz0h_oupas = 0.0 ; yu1 = 0.0    
    yv1 = 0.0     ; ypaprs = 0.0     ; ypplay = 0.0     ; yqbs1 = 0.0
    ydelp = 0.0   ; yu = 0.0         ; yv = 0.0        ; yt = 0.0         
    yq = 0.0      ; y_dflux_t = 0.0  ; y_dflux_q = 0.0 
    yqbs(:,:)=0.0  
    yrugoro = 0.0 ; ywindsp = 0.0   
    yfluxlat=0.0 ; y_flux0(:)=0.0
    yqsol = 0.0  ; yzxtsol = 0.0  

    ytke=0.
    yeps=0.
    yri0(:)=0.
    y_treedrg=0.

    ysnowhgt = 0.0; yqsnow = 0.0     ; yrunoff = 0.0   ; ytoice =0.0
    yalb3_new = 0.0  ; ysissnow = 0.0 
    ycldt = 0.0      ; yrmu0 = 0.0
    y_d_qbs(:,:)=0.0

    ytke_x=0.     ; ytke_w=0.        ; ywake_dltke=0.
    yeps_x=0.     ; yeps_w=0.
    y_d_t_x=0.    ; y_d_t_w=0.       ; y_d_q_x=0.      ; y_d_q_w=0.
    yfluxlat_x=0. ; yfluxlat_w=0.
    ywake_s=0.    ; ywake_cstar=0.   ;ywake_dens=0.

    tau_eq=0.     ; delta_coef=0.
    y_delta_flux_t1=0.
    ydtsurf_th=0.
    yts_x(:)=0.      ; yts_w(:)=0.
    y_delta_tsurf(:)=0. ; y_delta_qsurf(:)=0.
    yqsurf_x(:)=0.      ; yqsurf_w(:)=0.
    yg_T(:) = 0. ;        yg_Q(:) = 0.
    yGamma_dTs_phiT(:) = 0. ; yGamma_dQs_phiQ(:) = 0.
    ydTs_ins(:) = 0. ; ydqs_ins(:) = 0.

    ytsoil = 999999. 
    y_d_u_frein(:,:)=0.
    y_d_v_frein(:,:)=0.

#ifdef ISO
   yxtrain_f = 0.0 ; yxtsnow_f = 0.0
   yxtsnow  = 0.0
   yxt = 0.0
   yxtsol = 0.0
   flux_xt = 0.0
   yRland_ice = 0.0

   y_dflux_xt = 0.0  
   y_d_xt_x=0.      ; y_d_xt_w=0.       
#endif 

! >> PC
!the yfields_out variable is defined in (klon,nbcf_out) even if it is used on
!the ORCHIDEE grid and as such should be defined in yfields_out(knon,nbcf_out) but
!the knon variable is not known at that level of pbl_surface_mod

!the yfields_in variable is defined in (klon,nbcf_in) even if it is used on the
!ORCHIDEE grid and as such should be defined in yfields_in(knon,nbcf_in) but the
!knon variable is not known at that level of pbl_surface_mod
  IF (ALLOCATED(yfields_out)) yfields_out(:,:) = 0.
  IF (ALLOCATED(yfields_in))  yfields_in(:,:) = 0.
! >> PC

  ypphi = 0.0  



      
!****************************************************************************************
! 5) Compress variables 
!
!****************************************************************************************

!   Provisional : set ybeta to standard values
       IF (nsrf .NE. is_ter) THEN
           ybeta(1:knon) = 1.
       ELSE
           IF (iflag_split .EQ. 0) THEN
              ybeta(1:knon) = 1.
           ELSE
             DO j = 1, knon
                i = ni(j)
                ybeta(j)   = beta(i,nsrf)
             ENDDO
           ENDIF  ! (iflag_split .LE.1)
       ENDIF !  (nsrf .NE. is_ter)
!
       DO j = 1, knon
          i = ni(j)
          ypct(j)    = pctsrf(i,nsrf)
          yts(j)     = ts(i,nsrf)
          ysnow(j)   = snow(i,nsrf)
          yqsurf(j)  = qsurf(i,nsrf)
          yalb(j)    = alb(i,nsrf)
          yalb_vis(j) = alb_dir(i,1,nsrf)
          IF (nsw==6) THEN
            yalb_vis(j)=(alb_dir(i,1,nsrf)*SFRWL(1)+alb_dir(i,2,nsrf)*SFRWL(2) &
              +alb_dir(i,3,nsrf)*SFRWL(3))/(SFRWL(1)+SFRWL(2)+SFRWL(3))
          ENDIF
          yrain_f(j) = rain_f(i)
          ysnow_f(j) = snow_f(i)
          ybs_f(j)   = bs_f(i)
          yagesno(j) = agesno(i,nsrf)
          yfder(j)   = fder(i)
          ylwdown(j) = lwdown_m(i)
          ygustiness(j) = gustiness(i)
          ysolsw(j)  = solsw(i,nsrf)
          ysollw(j)  = sollw(i,nsrf)
          yz0m(j)  = z0m(i,nsrf)
          yz0h(j)  = z0h(i,nsrf)
          yrugoro(j) = rugoro(i)
          yu1(j)     = u(i,1)
          yv1(j)     = v(i,1)
          yqbs1(j)   = qbs(i,1)
          ypaprs(j,klev+1) = paprs(i,klev+1)
          ywindsp(j) = windsp(i,nsrf)
          yzmea(j)   = zmea(i)
          yzsig(j)   = zsig(i)
          ycldt(j)   = cldt(i)
          yrmu0(j)   = rmu0(i)
          y_delta_tsurf(j)=delta_tsurf(i,nsrf)
          yfluxbs(j)=0.0
          y_flux_bs(j) = 0.0
!!!
#ifdef ISO
          DO ixt=1,ntraciso
            yxtrain_f(ixt,j) = xtrain_f(ixt,i)
            yxtsnow_f(ixt,j) = xtsnow_f(ixt,i)  
          ENDDO
          DO ixt=1,niso
            yxtsnow(ixt,j)   = xtsnow(ixt,i,nsrf)
          ENDDO    
          DO ixt=1,niso
            yRland_ice(ixt,j)= Rland_ice(ixt,i)  
          ENDDO    
#ifdef ISOVERIF
          IF (iso_eau >= 0) THEN
              call iso_verif_egalite_choix(ysnow_f(j), &
     &          yxtsnow_f(iso_eau,j),'pbl_surf_mod 862', &
     &          errmax,errmaxrel)
              call iso_verif_egalite_choix(ysnow(j), &
     &          yxtsnow(iso_eau,j),'pbl_surf_mod 872', &
     &          errmax,errmaxrel)
          ENDIF
#endif
#ifdef ISOVERIF
         DO ixt=1,ntraciso
           call iso_verif_noNaN(yxtsnow_f(ixt,j),'pbl_surf_mod 921')
         ENDDO
#endif
#endif
       ENDDO
! >> PC
!--compressing fields_out onto ORCHIDEE grid
!--these fields are shared and used directly in surf_land_orchidee_mod
        IF (nbcf_out > 0 .AND. ALLOCATED(yfields_out) .AND. ALLOCATED(fields_out)) THEN
           DO n = 1, nbcf_out
              DO j = 1, knon
                 i = ni(j)
                 yfields_out(j,n) = fields_out(i,n)
              ENDDO
           ENDDO
        END IF
! << PC


       DO k = 1, klev
          DO j = 1, knon
             i = ni(j)
             ypaprs(j,k) = paprs(i,k)
             ypplay(j,k) = pplay(i,k)
             ydelp(j,k)  = delp(i,k)
          ENDDO
       ENDDO

        DO k = 1, klev+1
          DO j = 1, knon
             i = ni(j)
             ytke(j,k)   = tke_x(i,k,nsrf)
          ENDDO
        ENDDO

        DO k = 1, klev
          DO j = 1, knon
             i = ni(j)
             y_treedrg(j,k) =  treedrg(i,k,nsrf)
             yu(j,k) = u(i,k)
             yv(j,k) = v(i,k)
             yt(j,k) = t(i,k)
             yq(j,k) = q(i,k)
             yqbs(j,k)=qbs(i,k)
             ypphi(j,k) = pphi(i,k)

#ifdef ISO
             DO ixt=1,ntraciso   
               yxt(ixt,j,k) = xt(ixt,i,k)
             ENDDO !DO ixt=1,ntraciso
#endif
          ENDDO
        ENDDO
!
       IF (iflag_split.GE.1) THEN

        DO k = 1, klev
          DO j = 1, knon
             i = ni(j)
             yu_x(j,k) = u(i,k)
             yv_x(j,k) = v(i,k)
             yt_x(j,k) = t(i,k)-wake_s(i)*wake_dlt(i,k)
             yq_x(j,k) = q(i,k)-wake_s(i)*wake_dlq(i,k)
             yu_w(j,k) = u(i,k)
             yv_w(j,k) = v(i,k)
             yt_w(j,k) = t(i,k)+(1.-wake_s(i))*wake_dlt(i,k)
             yq_w(j,k) = q(i,k)+(1.-wake_s(i))*wake_dlq(i,k)
#ifdef ISO
             DO ixt=1,ntraciso
               yxt_x(ixt,j,k) = xt(ixt,i,k)-wake_s(i)*wake_dlxt(ixt,i,k)
               yxt_w(ixt,j,k) = xt(ixt,i,k)+(1.-wake_s(i))*wake_dlxt(ixt,i,k)
             ENDDO
#endif
          ENDDO
        ENDDO

        IF (prt_level .ge. 10) THEN
          print *,'pbl_surface, wake_s(1), wake_dlt(1,:) ', wake_s(1), wake_dlt(1,:)
          print *,'pbl_surface, wake_s(1), wake_dlq(1,:) ', wake_s(1), wake_dlq(1,:)
        ENDIF

        DO k = 1, klev+1
          DO j = 1, knon
             i = ni(j)
             ytke_x(j,k)      = tke_x(i,k,nsrf)
             ytke(j,k)        = tke_x(i,k,nsrf)+wake_s(i)*wake_dltke(i,k,nsrf)
             ytke_w(j,k)      = tke_x(i,k,nsrf)+wake_dltke(i,k,nsrf)
             ywake_dltke(j,k) = wake_dltke(i,k,nsrf)
          ENDDO
        ENDDO

        DO j = 1, knon
          i = ni(j)
          ywake_s(j)=wake_s(i)
          ywake_cstar(j)=wake_cstar(i)
          ywake_dens(j)=wake_dens(i)
        ENDDO

        DO j=1,knon
         yts_x(j)=yts(j)-ywake_s(j)*y_delta_tsurf(j)
         yts_w(j)=yts(j)+(1.-ywake_s(j))*y_delta_tsurf(j)
        ENDDO

       ENDIF  ! (iflag_split .ge.1)

       DO k = 1, nsoilmx
          DO j = 1, knon
             i = ni(j)
             ytsoil(j,k) = ftsoil(i,k,nsrf)
          ENDDO
       ENDDO
       
       ! qsol(water height in soil) only for bucket continental model
       IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN 
          DO j = 1, knon
             i = ni(j)
             yqsol(j) = qsol(i)
#ifdef ISO
             DO ixt=1,niso
               yxtsol(ixt,j) = xtsol(ixt,i)
             ENDDO
#endif
          ENDDO
       ENDIF

       if (nsrf == is_oce .and. activate_ocean_skin >= 1) then
          if (activate_ocean_skin == 2 .and. type_ocean == "couple") then
             ydelta_sal(:knon) = delta_sal(ni(:knon))
             ydelta_sst(:knon) = delta_sst(ni(:knon))
             ydter(:knon) = dter(ni(:knon))
             ydser(:knon) = dser(ni(:knon))
             ydt_ds(:knon) = dt_ds(ni(:knon))
          end if
          
          yds_ns(:knon) = ds_ns(ni(:knon))
          ydt_ns(:knon) = dt_ns(ni(:knon))
       end if
       
!****************************************************************************************
! 6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm.
!
!****************************************************************************************


       IF (iflag_split .eq.0) THEN

        DO i = 1, knon
           zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
                * (ypaprs(i,1)-ypplay(i,1))
           speed(i) = SQRT(yu(i,1)**2+yv(i,1)**2)
        ENDDO

        !!! AM heterogeneous continental subsurfaces
        IF (nsrf .EQ. is_ter) THEN
          ! compute time-dependent effective surface parameters (function of zgeo1) !! AM
          IF (iflag_hetero_surf .GT. 0) THEN
            DO n=1,nbtersrf
              DO j=1,knon
                i = ni(j)
                yfrac_tersrf(j,n) = frac_tersrf(i,n)
                yz0m_tersrf(j,n) = z0m_tersrf(i,n)
                IF (ratio_z0m_z0h_tersrf(i,n) .NE. 0.) THEN
                  yz0h_tersrf(j,n) = z0m_tersrf(i,n) / ratio_z0m_z0h_tersrf(i,n)
                ELSE
                  yz0h_tersrf(j,n) = 0.
                ENDIF
              ENDDO
            ENDDO
            !
            zgeo1_over_RG(1:knon) = zgeo1(1:knon)/RG
            CALL eff_surf_param(knon, nbtersrf, yz0m_tersrf, yfrac_tersrf, 'CDN', yz0m, zgeo1_over_RG)
            CALL eff_surf_param(knon, nbtersrf, yz0h_tersrf, yfrac_tersrf, 'CDN', yz0h, zgeo1_over_RG)
            !
          ENDIF
        ENDIF

!
        ypblh_tmp(:)=s_pblh(ni(:))

        CALL cdrag(knon, nsrf, &
            speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1), ypblh_tmp, &
            yts, yqsurf, yz0m, yz0h, yri0, 0, &
            ycdragm, ycdragh, zri1, pref, rain_f, yzxtsol, ypplay(:,1))
        s_pblh(ni(:)) = ypblh_tmp(:)
! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
     IF (ok_prescr_ust) THEN
      DO i = 1, knon
       print *,'ycdragm avant=',ycdragm(i)
       vent= sqrt(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1))
       ycdragm(i) = ust*ust/(1.+vent)/vent
      ENDDO
     ENDIF

        IF (prt_level >=10) print *,'cdrag -> ycdragh ', ycdragh(1:knon)
       ELSE  !(iflag_split .eq.0)

        DO i = 1, knon
           zgeo1_x(i) = RD * yt_x(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
                * (ypaprs(i,1)-ypplay(i,1))
           speed_x(i) = SQRT(yu_x(i,1)**2+yv_x(i,1)**2)
        ENDDO

            ypblh_tmp(:)=s_pblh_x(ni(:))

            CALL cdrag(knon, nsrf, &
            speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),ypblh,&
            yts_x, yqsurf_x, yz0m, yz0h, yri0, 0, &
            ycdragm_x, ycdragh_x, zri1_x, pref_x, rain_f, yzxtsol, ypplay(:,1) )
    
            s_pblh_x(ni(:)) = ypblh_tmp(:)
! --- special Dice. JYG+MPL 25112013
        IF (ok_prescr_ust) THEN
         DO i = 1, knon
          vent= sqrt(yu_x(i,1)*yu_x(i,1)+yv_x(i,1)*yv_x(i,1))
          ycdragm_x(i) = ust*ust/(1.+vent)/vent
         ENDDO
        ENDIF
        IF (prt_level >=10) print *,'cdrag -> ycdragh_x ', ycdragh_x(1:knon)

        DO i = 1, knon
           zgeo1_w(i) = RD * yt_w(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
                * (ypaprs(i,1)-ypplay(i,1))
           speed_w(i) = SQRT(yu_w(i,1)**2+yv_w(i,1)**2)
        ENDDO

        ypblh_tmp(:)=s_pblh_w(ni(:))

        CALL cdrag(knon, nsrf, &
            speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),s_pblh_w,&
            yts_w, yqsurf_w, yz0m, yz0h, yri0, 0, &
            ycdragm_w, ycdragh_w, zri1_w, pref_w, rain_f, yzxtsol, ypplay(:,1) )
        
        s_pblh_w(ni(:)) = ypblh_tmp(:)
!
        IF(ok_bug_zg_wk_pbl) THEN
         zgeo1(1:knon) = wake_s(1:knon)*zgeo1_w(1:knon) + (1.-wake_s(1:knon))*zgeo1_x(1:knon)
        ELSE
         zgeo1(1:knon) = ywake_s(1:knon)*zgeo1_w(1:knon) + (1.-ywake_s(1:knon))*zgeo1_x(1:knon)
        ENDIF

! --- special Dice. JYG+MPL 25112013 puis BOMEX
        IF (ok_prescr_ust) THEN
         DO i = 1, knon
          vent= sqrt(yu_w(i,1)*yu_w(i,1)+yv_w(i,1)*yv_w(i,1))
          ycdragm_w(i) = ust*ust/(1.+vent)/vent
         ENDDO
        ENDIF
        IF (prt_level >=10) print *,'cdrag -> ycdragh_w ', ycdragh_w(1:knon)
       ENDIF  ! (iflag_split .eq.0)


!****************************************************************************************
! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefh et ycoefm.
!
!****************************************************************************************

       IF (iflag_split .eq.0) THEN

      IF (prt_level >=10) THEN
      print *,' args coef_diff_turb: yu ',  yu(1:knon,:)  
      print *,' args coef_diff_turb: yv ',  yv(1:knon,:)    
      print *,' args coef_diff_turb: yq ',  yq(1:knon,:)    
      print *,' args coef_diff_turb: yt ',  yt(1:knon,:)    
      print *,' args coef_diff_turb: yts ', yts(1:knon)
      print *,' args coef_diff_turb: yz0m ', yz0m(1:knon)
      print *,' args coef_diff_turb: yqsurf ', yqsurf(1:knon)  
      print *,' args coef_diff_turb: ycdragm ', ycdragm(1:knon)
      print *,' args coef_diff_turb: ycdragh ', ycdragh(1:knon)
      print *,' args coef_diff_turb: ytke ', ytke(1:knon,:)    
       ENDIF

        IF (iflag_pbl>=50) THEN

        CALL call_atke(dtime,knon,klev,nsrf,ni,ycdragm, ycdragh,yus0,yvs0,yts, &
                  yu, yv,yt,yq,ypplay,ypaprs,       &
                  ytke,yeps, ycoefm, ycoefh)

        ELSE


        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
            ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, &
            ycoefm, ycoefh, ytke, yeps, y_treedrg)

       IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
! In this case, coef_diff_turb is called for the Cd only
       DO k = 2, klev
          DO j = 1, knon
             i = ni(j)
             ycoefh(j,k)   = zcoefh(i,k,nsrf)
             ycoefm(j,k)   = zcoefm(i,k,nsrf)
          ENDDO
       ENDDO
       ENDIF

       ENDIF ! iflag_pbl >= 50

        IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh ',ycoefh(1:knon,:)


       ELSE  !(iflag_split .eq.0)

      
      IF (prt_level >=10) THEN
      print *,' args coef_diff_turb: yu_x ',  yu_x(1:knon,:)      
      print *,' args coef_diff_turb: yv_x ',  yv_x(1:knon,:)      
      print *,' args coef_diff_turb: yq_x ',  yq_x(1:knon,:)      
      print *,' args coef_diff_turb: yt_x ',  yt_x(1:knon,:)      
      print *,' args coef_diff_turb: yts_x ', yts_x(1:knon)
      print *,' args coef_diff_turb: yqsurf ', yqsurf(1:knon)  
      print *,' args coef_diff_turb: ycdragm_x ', ycdragm_x(1:knon)
      print *,' args coef_diff_turb: ycdragh_x ', ycdragh_x(1:knon)
      print *,' args coef_diff_turb: ytke_x ', ytke_x(1:knon,:)    
      ENDIF


        IF (iflag_pbl>=50) THEN
     
        CALL call_atke(dtime,knon,klev,nsrf,ni,ycdragm_x(1:knon),ycdragh_x(1:knon),yus0(1:knon),yvs0(1:knon),yts_x(1:knon),    &
                       yu_x(1:knon,:),yv_x(1:knon,:),yt_x(1:knon,:),yq_x(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:),  &
                       ytke_x(1:knon,:),yeps_x(1:knon,:),ycoefm_x(1:knon,:), ycoefh_x(1:knon,:))

        ELSE


        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
            ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf_x, ycdragm_x, &
            ycoefm_x, ycoefh_x, ytke_x,yeps_x,y_treedrg)

!FC doit on le mettre ( on ne l utilise pas si il y a du spliting)
       IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
! In this case, coef_diff_turb is called for the Cd only
       DO k = 2, klev
          DO j = 1, knon
             i = ni(j)
             ycoefh_x(j,k)   = zcoefh(i,k,nsrf)
             ycoefm_x(j,k)   = zcoefm(i,k,nsrf)
          ENDDO
       ENDDO
       ENDIF

        ENDIF ! iflag_pbl >= 50

        IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_x ',ycoefh_x(1:knon,:)
!
      IF (prt_level >=10) THEN
      print *,' args coef_diff_turb: yu_w ',  yu_w(1:knon,:)
      print *,' args coef_diff_turb: yv_w ',  yv_w(1:knon,:)  
      print *,' args coef_diff_turb: yq_w ',  yq_w(1:knon,:)  
      print *,' args coef_diff_turb: yt_w ',  yt_w(1:knon,:)  
      print *,' args coef_diff_turb: yts_w ', yts_w(1:knon)
      print *,' args coef_diff_turb: yqsurf ', yqsurf(1:knon)  
      print *,' args coef_diff_turb: ycdragm_w ', ycdragm_w(1:knon)
      print *,' args coef_diff_turb: ycdragh_w ', ycdragh_w(1:knon)
      print *,' args coef_diff_turb: ytke_w ', ytke_w(1:knon,:)
      ENDIF
      
        IF (iflag_pbl>=50) THEN
        

        CALL call_atke(dtime,knon,klev,nsrf,ni,ycdragm_w,ycdragh_w,yus0,yvs0,yts_w, &
                yu_w, yv_w, yt_w, yq_w, ypplay, ypaprs,      &
                ytke_w, yeps_w, ycoefm_w, ycoefh_w)

        ELSE

        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
            ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf_w, ycdragm_w, &
            ycoefm_w, ycoefh_w, ytke_w,yeps_w,y_treedrg)

       IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
! In this case, coef_diff_turb is called for the Cd only
       DO k = 2, klev
          DO j = 1, knon
             i = ni(j)
             ycoefh_w(j,k)   = zcoefh(i,k,nsrf)
             ycoefm_w(j,k)   = zcoefm(i,k,nsrf)
          ENDDO
       ENDDO
       ENDIF

       ENDIF ! iflag_pbl >= 50


        IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_w ',ycoefh_w(1:knon,:)

!!!jyg le 10/04/2013
!!   En attendant de traiter le transport des traceurs dans les poches froides, formule
!!   arbitraire pour ycoefh et ycoefm
      DO k = 2,klev
        DO j = 1,knon
         ycoefh(j,k) = ycoefh_x(j,k) + ywake_s(j)*(ycoefh_w(j,k) - ycoefh_x(j,k))
         ycoefm(j,k) = ycoefm_x(j,k) + ywake_s(j)*(ycoefm_w(j,k) - ycoefm_x(j,k))
        ENDDO
      ENDDO


       ENDIF  ! (iflag_split .eq.0)

       
!****************************************************************************************
! 
! 8) "La descente" - "The downhill"
!  
!  climb_hq_down and climb_wind_down calculate the coefficients
!  Ccoef_X et Dcoef_X for X=[H, Q, U, V].
!  Only the coefficients at surface for H and Q are returned.
!
!****************************************************************************************

! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q 
       IF (iflag_split .eq.0) THEN

        CALL climb_hq_down(knon, ni, ycoefh, ypaprs, ypplay, &
            ydelp, yt, yq, dtime, &
            CcoefH, CcoefQ, DcoefH, DcoefQ, &
            Kcoef_hq, gama_q, gama_h, &
            AcoefH, AcoefQ, BcoefH, BcoefQ &
#ifdef ISO
         &   ,yxt, CcoefXT, DcoefXT, gama_xt, AcoefXT, BcoefXT & 
#endif               
         &   )
       ELSE  !(iflag_split .eq.0)

        CALL climb_hq_down(knon, ni, ycoefh_x, ypaprs, ypplay, &
            ydelp, yt_x, yq_x, dtime, &
            CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
            Kcoef_hq_x, gama_q_x, gama_h_x, &
            AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x &
#ifdef ISO
         &   ,yxt_x, CcoefXT_x, DcoefXT_x, gama_xt_x, AcoefXT_x, BcoefXT_x & 
#endif               
         &   )

       IF (prt_level >=10) THEN
         PRINT *,'pbl_surface (climb_hq_down.x->) AcoefH_x ',AcoefH_x
         PRINT *,'pbl_surface (climb_hq_down.x->) AcoefQ_x ',AcoefQ_x
         PRINT *,'pbl_surface (climb_hq_down.x->) BcoefH_x ',BcoefH_x
         PRINT *,'pbl_surface (climb_hq_down.x->) BcoefQ_x ',BcoefQ_x
       ENDIF


        CALL climb_hq_down(knon, ni, ycoefh_w, ypaprs, ypplay, &
            ydelp, yt_w, yq_w, dtime, &
            CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
            Kcoef_hq_w, gama_q_w, gama_h_w, &
            AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w &
#ifdef ISO
         &   ,yxt_w, CcoefXT_w, DcoefXT_w, gama_xt_w, AcoefXT_w, BcoefXT_w & 
#endif               
         &   )

       IF (prt_level >=10) THEN
         PRINT *,'pbl_surface (climb_hq_down.w->) AcoefH_w ',AcoefH_w
         PRINT *,'pbl_surface (climb_hq_down.w->) AcoefQ_w ',AcoefQ_w
         PRINT *,'pbl_surface (climb_hq_down.w->) BcoefH_w ',BcoefH_w
         PRINT *,'pbl_surface (climb_hq_down.w->) BcoefQ_w ',BcoefQ_w
       ENDIF

       ENDIF  ! (iflag_split .eq.0)


! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
       IF (iflag_split .eq.0) THEN

        CALL climb_wind_down(knon, ni, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &
            CcoefU, CcoefV, DcoefU, DcoefV, &
            Kcoef_m, alf_1, alf_2, &
            AcoefU, AcoefV, BcoefU, BcoefV)
       ELSE  ! (iflag_split .eq.0)

        CALL climb_wind_down(knon, ni, dtime, ycoefm_x, ypplay, ypaprs, yt_x, ydelp, yu_x, yv_x, &
            CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
            Kcoef_m_x, alf_1_x, alf_2_x, &
            AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x)

        CALL climb_wind_down(knon, ni, dtime, ycoefm_w, ypplay, ypaprs, yt_w, ydelp, yu_w, yv_w, &
            CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
            Kcoef_m_w, alf_1_w, alf_2_w, &
            AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w)
       ENDIF  ! (iflag_split .eq.0)

! For blowing snow:
    IF (ok_bs) THEN
     ! following Bintanja et al 2000, part II and Vionnet V PhD thesis
     ! we assume that the eddy diffsivity coefficient for
     ! suspended particles is a fraction of Kh 
     do k=1,klev
        do j=1,knon
           ycoefqbs(j,k)=ycoefh(j,k)*zeta_bs 
        enddo
     enddo

     CALL climb_qbs_down(knon, ni, ycoefqbs, ypaprs, ypplay, &
     ydelp, yt, yqbs, dtime, & 
     CcoefQBS, DcoefQBS, &
     Kcoef_qbs, gama_qbs, &
     AcoefQBS, BcoefQBS)
    ENDIF

!****************************************************************************************
! 9) Small calculations
!
!****************************************************************************************

! - Reference pressure is given the values at surface level          
       ypsref(:) = ypaprs(:,1)  

! - CO2 field on 2D grid to be sent to ORCHIDEE
!   Transform to compressed field
       IF (carbon_cycle_cpl) THEN
          DO i=1,knon
             r_co2_ppm(i) = co2_send(ni(i))
          ENDDO
       ELSE
          r_co2_ppm(:) = co2_ppm     ! Constant field
       ENDIF

!!! nrlmd le 02/05/2011  -----------------------On raccorde les 2 colonnes dans la couche 1 

       IF (iflag_split .eq. 0) THEN
         yt1(:) = yt(:,1)
         yq1(:) = yq(:,1)
#ifdef ISO
         yxt1(:,:) = yxt(:,:,1)
#endif

       ELSE IF (iflag_split .ge. 1) THEN
#ifdef ISO
        call abort_physic('pbl_surface_mod 2149','isos pas encore dans iflag_split=1',1)
#endif

!
! Cdragq computation
! ------------------
    !******************************************************************************
    ! Cdragq computed from cdrag
    ! The difference comes only from a factor (f_z0qh_oce) on z0, so that
    ! it can be computed inside wx_pbl0_merge
    ! More complicated appraches may require the propagation through
    ! pbl_surface of an independant cdragq variable.
    !******************************************************************************
!
    IF ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce) THEN
       ! Si on suit les formulations par exemple de Tessel, on 
       ! a z0h=0.4*nu/u*, z0q=0.62*nu/u*, d'ou f_z0qh_oce=0.62/0.4=1.55
!!       ycdragq_x(1:knon)=ycdragh_x(1:knon)*                                      &
!!            log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
!!       ycdragq_w(1:knon)=ycdragh_w(1:knon)*                                      &
!!            log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
!
       DO j = 1,knon
         z1lay = zgeo1(j)/RG
         fact_cdrag = log(z1lay/yz0h(j))/log(z1lay/(f_z0qh_oce*yz0h(j)))
         ycdragq_x(j)=ycdragh_x(j)*fact_cdrag
         ycdragq_w(j)=ycdragh_w(j)*fact_cdrag
       ENDDO  ! j = 1,knon

    ELSE
       ycdragq_x(1:knon)=ycdragh_x(1:knon)
       ycdragq_w(1:knon)=ycdragh_w(1:knon)
    ENDIF  ! ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce)
!

         CALL wx_pbl_prelim_0(knon, nsrf, dtime, ypplay, ypaprs, ywake_s,  &
                         yts, y_delta_tsurf, ygustiness, &
                         yt_x, yt_w, yq_x, yq_w, &
                         yu_x, yu_w, yv_x, yv_w, &
                         ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
                         ycdragm_x, ycdragm_w, &
                         AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
                         AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
                         BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
                         BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
                         Kech_h_x, Kech_h_w, Kech_h  &
                         )
         CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta,  &
                         BcoefQ_x, BcoefQ_w  &
                         )

         CALL wx_pbl0_merge(knon, ypplay, ypaprs,  &
                         ywake_s, ydTs0, ydqs0, &
                         yt_x, yt_w, yq_x, yq_w, &
                         yu_x, yu_w, yv_x, yv_w, &
                         ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
                         ycdragm_x, ycdragm_w, &
                         AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
                         AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
                         BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
                         BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
                         AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &
                         BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
                         ycdragh, ycdragq, ycdragm, &
                         yt1, yq1, yu1, yv1 &
                         )
         IF (iflag_split .eq. 2 .AND. nsrf .ne. is_oce) THEN

           CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
                           ywake_s, ybeta, ywake_cstar, ywake_dens, &
                           AcoefH_x, AcoefH_w, &
                           BcoefH_x, BcoefH_w, &
                           AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
                           AcoefH, AcoefQ, BcoefH, BcoefQ,  &
                           HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
                           phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
                           yg_T, yg_Q, &
                           yGamma_dTs_phiT, yGamma_dQs_phiQ, &
                           ydTs_ins, ydqs_ins &
                           )
         ELSE !
           AcoefH(:) = AcoefH_0(:)
           AcoefQ(:) = AcoefQ_0(:)
           BcoefH(:) = BcoefH_0(:)
           BcoefQ(:) = BcoefQ_0(:)
           yg_T(:) = 0.
           yg_Q(:) = 0.
           yGamma_dTs_phiT(:) = 0.
           yGamma_dQs_phiQ(:) = 0.
           ydTs_ins(:) = 0.
           ydqs_ins(:) = 0.
         ENDIF   ! (iflag_split .eq. 2)
       ENDIF  ! (iflag_split .eq.0)

       IF (prt_level >=10) THEN
         DO i = 1, min(1,knon)
           PRINT *,'pbl_surface (merge->): yt(1,:) ',yt(i,:)
           PRINT *,'pbl_surface (merge->): yq(1,:) ',yq(i,:)
           PRINT *,'pbl_surface (merge->): yu(1,:) ',yu(i,:)
           PRINT *,'pbl_surface (merge->): yv(1,:) ',yv(i,:)
           PRINT *,'pbl_surface (merge->): AcoefH(1), AcoefQ(1), AcoefU(1), AcoefV(1) ', &
                                           AcoefH(i), AcoefQ(i), AcoefU(i), AcoefV(i)
           PRINT *,'pbl_surface (merge->): BcoefH(1), BcoefQ(1), BcoefU(1), BcoefV(1) ', &
                                           BcoefH(i), BcoefQ(i), BcoefU(i), BcoefV(i)
         ENDDO

       ENDIF

!  Save initial value of z0h for use in evappot (z0h wiil be computed again in the surface models)
          yz0h_old(1:knon) = yz0h(1:knon)
!
!****************************************************************************************
!
! Calulate t2m and q2m for the case of calculation at land grid points 
! t2m and q2m are needed as input to ORCHIDEE
!
!****************************************************************************************
       IF (nsrf == is_ter) THEN

          DO i = 1, knon
             zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
                  * (ypaprs(i,1)-ypplay(i,1))
          ENDDO

          ! Calculate the temperature et relative humidity at 2m and the wind at 10m 
          IF (iflag_new_t2mq2m==1) THEN

         CALL stdlevvarn(knon, knon, is_ter, zxli, &
               yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
               yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
               yt2m, yq2m, yt10m, yq10m, yu10m, yustar, &
               yn2mout)
          ELSE 

          CALL stdlevvar(knon, knon, is_ter, zxli, &
               yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
               yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
               yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, yzxtsol)
          ENDIF
          
       ENDIF

!****************************************************************************************
!
! 10) Switch according to current surface
!     It is necessary to start with the continental surfaces because the ocean
!     needs their run-off.
!
!****************************************************************************************
       SELECT CASE(nsrf)
     
       CASE(is_ter)

          CALL surf_land(itap, dtime, date0, jour, knon, ni,&
               rlon, rlat, yrmu0, &
               debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt1, yq1,&
               AcoefH, AcoefQ, BcoefH, BcoefQ, & 
               AcoefU, AcoefV, BcoefU, BcoefV, & 
               ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
               ylwdown, yq2m, yt2m, &
               ysnow, yqsol, yagesno, ytsoil, &
               yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,yfluxbs,&
               yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
               y_flux_u1, y_flux_v1, &
               yveget,ylai,yheight, tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &
               cdragm_tersrf, cdragh_tersrf, &
               swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf  &
#ifdef ISO
         &      ,yxtrain_f, yxtsnow_f,yxt1, &
         &      yxtsnow,yxtsol,yxtevap,h1, &
         &      yrunoff_diag,yxtrunoff_diag,yRland_ice &
#endif               
         &      )

          tsurf_tersrf(:,:) =  tsurf_new_tersrf(:,:) ! for next time step

            IF (ifl_pbltree .ge. 1) THEN

              CALL   freinage(knon, knon, yu, yv, yt, &
                yveget,ylai, yheight,ypaprs,ypplay,y_treedrg, y_d_u_frein,y_d_v_frein)
            ENDIF

               
! Special DICE MPL 05082013 puis BOMEX
       IF (ok_prescr_ust) THEN
          DO j=1,knon
            y_flux_u1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yu(j,1)*ypplay(j,1)/RD/yt(j,1)
            y_flux_v1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yv(j,1)*ypplay(j,1)/RD/yt(j,1)
          ENDDO
      ENDIF

#ifdef ISOVERIF
        DO j=1,knon
          DO ixt=1,ntraciso
            CALL iso_verif_noNaN(yxtevap(ixt,j), &
         &      'pbl_surface 1056a: apres surf_land')
          ENDDO
          DO ixt=1,niso
            CALL iso_verif_noNaN(yxtsol(ixt,j), &
         &      'pbl_surface 1056b: apres surf_land')
          ENDDO
        ENDDO
#endif
#ifdef ISOVERIF

        DO j=1,knon
          IF (iso_eau >= 0) THEN     
                 CALL iso_verif_egalite(yxtsnow(iso_eau,j), &
     &                                  ysnow(j),'pbl_surf_mod 1043')
          ENDIF !if (iso_eau.gt.0) then
        ENDDO !DO i=1,klon
#endif
    
       CASE(is_lic)
          IF (landice_opt .LT. 2) THEN
             ! Land ice is treated by LMDZ and not by ORCHIDEE

             CALL surf_landice(itap, dtime, knon, ni, &
                  rlon, rlat, debut, lafin, &
                  yrmu0, ylwdown, yalb, zgeo1, &
                  ysolsw, ysollw, yts, ypplay(:,1), &
                  ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt1, yq1,&
                  AcoefH, AcoefQ, BcoefH, BcoefQ, &
                  AcoefU, AcoefV, BcoefU, BcoefV, &
                  AcoefQBS, BcoefQBS, &
                  ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
                  ysnow, yqsurf, yqsol,yqbs1, yagesno, &
                  ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yicesub, yicemelt, yfluxsens,yfluxlat, &
                  yfluxbs, ytsurf_new, y_dflux_t, y_dflux_q, &
                  yzmea, yzsig, ycldt, &
                  ysnowhgt, yqsnow, ytoice, ysissnow, &
                  yalb3_new, yrunoff, &
                  y_flux_u1, y_flux_v1 &
#ifdef ISO
                  &    ,yxtrain_f, yxtsnow_f,yxt1,yRland_ice &
                  &    ,yxtsnow,yxtsol,yxtevap &
#endif              
                  &    )
             
             DO j = 1, knon
                i = ni(j)
                alb3_lic(i) = yalb3_new(j)
                snowhgt(i)   = ysnowhgt(j)
                qsnow(i)     = yqsnow(j)
                to_ice(i)    = ytoice(j)
                sissnow(i)   = ysissnow(j)
                runoff(i)    = yrunoff(j)
                icesub_ice(i) = icesub_ice(i) + yicesub(j)*ypct(j)
                icemelt_ice(i) = icemelt_ice(i) + yicemelt(j)*ypct(j)
             ENDDO
             ! Martin
             ! Special DICE MPL 05082013 puis BOMEX MPL 20150410
             IF (ok_prescr_ust) THEN
                DO j=1,knon
                   y_flux_u1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yu(j,1)*ypplay(j,1)/RD/yt(j,1)
                   y_flux_v1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yv(j,1)*ypplay(j,1)/RD/yt(j,1)
                ENDDO
             ENDIF

#ifdef ISOVERIF
             DO j=1,knon
               DO ixt=1,ntraciso
                 CALL iso_verif_noNaN(yxtevap(ixt,j), &
                        &             'pbl_surface 1095a: apres surf_landice')
               ENDDO
                do ixt=1,niso
                   call iso_verif_noNaN(yxtsol(ixt,j), &
                        &      'pbl_surface 1095b: apres surf_landice')
                enddo
             enddo
#endif
#ifdef ISOVERIF

             do j=1,knon
               IF (iso_eau >= 0) THEN     
                 CALL iso_verif_egalite(yxtsnow(iso_eau,j), &
                        &               ysnow(j),'pbl_surf_mod 1064')
               ENDIF !if (iso_eau >= 0) THEN
             ENDDO !DO i=1,klon
#endif
            
          END IF
          
       CASE(is_oce)
! calculate length scale PBL

        if (iflag_leads == 1) then
        ydthetadz = 999999.
        ypphii = 999999.
        ytheta = 999999.

        DO k = 1, klev
          DO j = 1, knon
             ytheta(j,k) = yt(j,k)*(ypplay(j,k)/1.e5)**(RD/RCPD)
          ENDDO
        ENDDO

        DO k = 2, klev
          DO j = 1, knon
             ydthetadz(j,k) = RG*( ytheta(j,k) - ytheta(j,k-1) ) / ( ypphi(j,k) - ypphi(j,k-1) )
             ypphii(j,k) = (ypphi(j,k)+ypphi(j,k-1))/(RG*2.)
          ENDDO
        ENDDO

!ym  minloc does't work on GPU (nvfortran + openacc) 
!        DO j = 1, knon
!             k= minloc(abs(ypphii(j,:)-300),1)
!             ydthetadz300(j)=ydthetadz(j,k)
!        ENDDO

        ydthetadz300(1:knon) = ydthetadz(1:knon,1)
        DO k=2, klev
          DO j=1,knon
            IF (abs(ypphii(j,k)-300) < ydthetadz300(j)) ydthetadz300(j) = abs(ypphii(j,k)-300)
          ENDDO
        ENDDO

        end if
           zgeo1_over_RG(1:knon) = zgeo1(1:knon)/RG
           CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb_vis, &
               ywindsp, yrmu0, yfder, yts, &
               itap, dtime, jour, knon, ni, &
               ypplay(:,1), zgeo1_over_RG, ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt(:,1), yq(:,1),&    ! ym missing init
               AcoefH, AcoefQ, BcoefH, BcoefQ, &
               AcoefU, AcoefV, BcoefU, BcoefV, &
               ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
               ysnow, yqsurf, yagesno, &
               yz0m, yz0h, SFRWL,yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
               ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
               y_flux_u1, y_flux_v1, ydelta_sst(:knon), ydelta_sal(:knon), &
               yds_ns(:knon), ydt_ns(:knon), ydter(:knon), ydser(:knon), &
               ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss, &
               ydthetadz300,Ampl                 &

#ifdef ISO
         &      ,yxtrain_f, yxtsnow_f,yxt1,Roce, &
         &      yxtsnow,yxtevap,h1 &
#endif               
         &      )

      IF (prt_level >=10) THEN
          print *,'arg de surf_ocean: ycdragh ',ycdragh(1:knon)
          print *,'arg de surf_ocean: ycdragm ',ycdragm(1:knon)
          print *,'arg de surf_ocean: yt ', yt(1:knon,:)
          print *,'arg de surf_ocean: yq ', yq(1:knon,:)
          print *,'arg de surf_ocean: yts ', yts(1:knon)
          print *,'arg de surf_ocean: AcoefH ',AcoefH(1:knon)
          print *,'arg de surf_ocean: AcoefQ ',AcoefQ(1:knon)
          print *,'arg de surf_ocean: BcoefH ',BcoefH(1:knon)
          print *,'arg de surf_ocean: BcoefQ ',BcoefQ(1:knon)
          print *,'arg de surf_ocean: yevap ',yevap(1:knon)
          print *,'arg de surf_ocean: yfluxsens ',yfluxsens(1:knon)
          print *,'arg de surf_ocean: yfluxlat ',yfluxlat(1:knon)
          print *,'arg de surf_ocean: ytsurf_new ',ytsurf_new(1:knon)
       ENDIF
! Special DICE MPL 05082013 puis BOMEX MPL 20150410
       IF (ok_prescr_ust) THEN
          DO j=1,knon
          y_flux_u1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yu(j,1)*ypplay(j,1)/RD/yt(j,1)
          y_flux_v1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yv(j,1)*ypplay(j,1)/RD/yt(j,1)
          ENDDO
      ENDIF
          
       CASE(is_sic)

          CALL surf_seaice( &
               rlon, rlat, ysolsw, ysollw, yalb_vis, yfder, &
               itap, dtime, jour, knon, ni, &
               lafin, &
               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
               AcoefH, AcoefQ, BcoefH, BcoefQ, &
               AcoefU, AcoefV, BcoefU, BcoefV, &
               ypsref, yu1, yv1, ygustiness, pctsrf, &
               ysnow, yqsurf, yqsol, yagesno, ytsoil, &
               yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yicesub, yicemelt, yfluxsens,yfluxlat,&
               ytsurf_new, y_dflux_t, y_dflux_q, &
               y_flux_u1, y_flux_v1, &
               hice,tice,bilg_cumul, &
               fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, &
               dtice_melt, dtice_snow2sic     &
#ifdef ISO
         &      ,yxtrain_f, yxtsnow_f,yxt1,Roce, &
         &      yxtsnow,yxtsol,yxtevap,Rland_ice &
#endif               
         &      )

             DO j = 1, knon
                i = ni(j)
                icesub_ice(i) = icesub_ice(i) + yicesub(j)*ypct(j)
                icemelt_ice(i) = icemelt_ice(i) + yicemelt(j)*ypct(j)
             ENDDO
! Special DICE MPL 05082013 puis BOMEX MPL 20150410
       IF (ok_prescr_ust) THEN
          DO j=1,knon
          y_flux_u1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yu(j,1)*ypplay(j,1)/RD/yt(j,1)
          y_flux_v1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yv(j,1)*ypplay(j,1)/RD/yt(j,1)
          ENDDO
       ENDIF

#ifdef ISOVERIF
        DO j=1,knon
          DO ixt=1,ntraciso
            CALL iso_verif_noNaN(yxtevap(ixt,j), &
         &                       'pbl_surface 1165a: apres surf_seaice')
          ENDDO
          DO ixt=1,niso
            CALL iso_verif_noNaN(yxtsol(ixt,j), &
         &      'pbl_surface 1165b: apres surf_seaice')
          ENDDO
        ENDDO
#endif
#ifdef ISOVERIF
        DO j=1,knon
          IF (iso_eau >= 0) THEN     
                 CALL iso_verif_egalite(yxtsnow(iso_eau,j), &
     &                                  ysnow(j),'pbl_surf_mod 1106')
          ENDIF !IF (iso_eau >= 0) THEN
        ENDDO !DO i=1,klon
#endif

       CASE DEFAULT
          WRITE(lunout,*) 'Surface index = ', nsrf
          abort_message = 'Surface index not valid'
!ym          CALL abort_physic(modname,abort_message,1)
       END SELECT


!****************************************************************************************
! 11) - Calcul the increment of surface temperature
!
!****************************************************************************************

       IF (evap0>=0.) THEN
          yevap(1:knon)=evap0
          yevap(1:knon)=RLVTT*evap0
       ENDIF

       y_d_ts(1:knon)   = ytsurf_new(1:knon) - yts(1:knon)
 
!****************************************************************************************
!
! 12) "La remontee" - "The uphill"
!
!  The fluxes (y_flux_X) and tendancy (y_d_X) are calculated 
!  for X=H, Q, U and V, for all vertical levels.
!
!****************************************************************************************

        IF (ok_forc_tsurf) THEN
            DO j=1,knon
                ytsurf_new(j)=tg
                y_d_ts(j) = ytsurf_new(j) - yts(j) 
            ENDDO
        ENDIF ! ok_forc_tsurf

        IF (ok_flux_surf) THEN
          IF (prt_level >=10) THEN
           PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
          ENDIF
          y_flux_t1(:) =  fsens
          y_flux_q1(:) =  flat/RLVTT
          yfluxlat(:) =  flat
!
!!  Test sur iflag_split retire le 2/02/2018, sans vraiment comprendre la raison de ce test. (jyg)
!!          IF (iflag_split .eq.0) THEN
             DO j=1,knon
             Kech_h(j) = ycdragh(j) * (1.0+SQRT(yu(j,1)**2+yv(j,1)**2)) * &
                  ypplay(j,1)/(RD*yt(j,1))
             ENDDO
!!          ENDIF ! (iflag_split .eq.0)

          DO j = 1, knon
            yt1_new=(1./RCPD)*(AcoefH(j)+BcoefH(j)*y_flux_t1(j)*dtime)
            ytsurf_new(j)=yt1_new-y_flux_t1(j)/(Kech_h(j)*RCPD)
            ! for cases forced in flux and for which forcing in Ts is needed
            ! to prevent the latter to reach unrealistic value (even if not used,
            ! Ts is calculated and hgardfou can appear during the calculation
            ! of surface saturation humidity for example
            if (ok_forc_tsurf) ytsurf_new(j)=tg
          ENDDO

          DO j=1,knon
          y_d_ts(j) = ytsurf_new(j) - yts(j) 
          ENDDO

        ELSE ! (ok_flux_surf)
          DO j=1,knon
          y_flux_t1(j) =  yfluxsens(j)
          y_flux_q1(j) = -yevap(j)
#ifdef ISO
          y_flux_xt1(:,:) = -yxtevap(:,:)
#endif
          ENDDO
        ENDIF ! (ok_flux_surf)

        ! flux of blowing snow at the first level 
        IF (ok_bs) THEN 
        DO j=1,knon
        y_flux_bs(j)=yfluxbs(j)
        ENDDO
        ENDIF
!
! ------------------------------------------------------------------------------
! 12a)  Splitting
! ------------------------------------------------------------------------------

       IF (iflag_split .GE. 1) THEN
#ifdef ISO
        call abort_physic('pbl_surface_mod 2607','isos pas encore dans iflag_split=1',1)
#endif

         IF (nsrf .ne. is_oce) THEN

!         Compute potential evaporation and aridity factor  (jyg, 20200328)
          ybeta_prev(:) = ybeta(:)
             DO j = 1, knon
               yqa(j) = AcoefQ(j) - BcoefQ(j)*yevap(j)*dtime
             ENDDO

          CALL wx_evappot(knon, yqa, yTsurf_new, yevap_pot)

          ybeta(1:knon) = min(yevap(1:knon)/yevap_pot(1:knon), 1.)
          
          IF (prt_level >=10) THEN
           DO j=1,knon
            print*,'y_flux_t1,yfluxlat,wakes' &
 &                ,  y_flux_t1(j), yfluxlat(j), ywake_s(j)
            print*,'beta_prev, beta, ytsurf_new', ybeta_prev(j), ybeta(j), ytsurf_new(j)
            print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j)
           ENDDO
          ENDIF  ! (prt_level >=10)
!
! Second call to wx_pbl0_merge and wx_pbl_dts_merge in order to take into account 
! the update of the aridity coeficient beta.
!
        CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta,  &
                        BcoefQ_x, BcoefQ_w  &
                        )

        CALL wx_pbl0_merge(knon, ypplay, ypaprs,  &
                          ywake_s, ydTs0, ydqs0, &
                          yt_x, yt_w, yq_x, yq_w, &
                          yu_x, yu_w, yv_x, yv_w, &
                          ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
                          ycdragm_x, ycdragm_w, &
                          AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
                          AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
                          BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
                          BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
                          AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &
                          BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
                          ycdragh, ycdragq, ycdragm, &
                          yt1, yq1, yu1, yv1 &
                          )
          IF (iflag_split .eq. 2) THEN

            CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
                            ywake_s, ybeta, ywake_cstar, ywake_dens, &
                            AcoefH_x, AcoefH_w, &
                            BcoefH_x, BcoefH_w, &
                            AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
                            AcoefH, AcoefQ, BcoefH, BcoefQ,  &
                            HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
                            phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
                            yg_T, yg_Q, &
                            yGamma_dTs_phiT, yGamma_dQs_phiQ, &
                            ydTs_ins, ydqs_ins &
                            )
          ELSE !
            AcoefH(:) = AcoefH_0(:)
            AcoefQ(:) = AcoefQ_0(:)
            BcoefH(:) = BcoefH_0(:)
            BcoefQ(:) = BcoefQ_0(:)
            yg_T(:) = 0.
            yg_Q(:) = 0.
            yGamma_dTs_phiT(:) = 0.
            yGamma_dQs_phiQ(:) = 0.
            ydTs_ins(:) = 0.
            ydqs_ins(:) = 0.
          ENDIF   ! (iflag_split .eq. 2)
!
        ELSE    ! (nsrf .ne. is_oce)
          ybeta(1:knon) = 1.
          yevap_pot(1:knon) = yevap(1:knon)
          AcoefH(:) = AcoefH_0(:)
          AcoefQ(:) = AcoefQ_0(:)
          BcoefH(:) = BcoefH_0(:)
          BcoefQ(:) = BcoefQ_0(:)
          yg_T(:) = 0.
          yg_Q(:) = 0.
          yGamma_dTs_phiT(:) = 0.
          yGamma_dQs_phiQ(:) = 0.
          ydTs_ins(:) = 0.
          ydqs_ins(:) = 0.
        ENDIF   ! (nsrf .ne. is_oce)
 
        CALL wx_pbl_split(knon, nsrf, dtime, ywake_s, ybeta, iflag_split, &
                       yg_T, yg_Q, &
                       yGamma_dTs_phiT, yGamma_dQs_phiQ, &
                       ydTs_ins, ydqs_ins, &
                       y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &
                       phiQ0_b, phiT0_b, &
                       y_flux_t1_x, y_flux_t1_w, &
                       y_flux_q1_x, y_flux_q1_w, &
                       y_flux_u1_x, y_flux_u1_w, &
                       y_flux_v1_x, y_flux_v1_w, &
                       yfluxlat_x, yfluxlat_w, &
                       y_delta_qsats, &
                       y_delta_tsurf_new, y_delta_qsurf &
                       )

         CALL wx_pbl_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
                       yTs, y_delta_tsurf,  & 
                       yqsurf, yTsurf_new,  &
                       y_delta_tsurf_new, y_delta_qsats,  &
                       AcoefH_x, AcoefH_w, &
                       BcoefH_x, BcoefH_w, &
                       AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
                       AcoefH, AcoefQ, BcoefH, BcoefQ,  &
                       y_flux_t1, y_flux_q1,  &
                       y_flux_t1_x, y_flux_t1_w, &
                       y_flux_q1_x, y_flux_q1_w)

         IF (nsrf .ne. is_oce) THEN

           CALL wx_pbl_dts_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
                         yTs, y_delta_tsurf,  & 
                         yqsurf, yTsurf_new,  &
                         y_delta_qsats, y_delta_tsurf_new, y_delta_qsurf,  &
                         AcoefH_x, AcoefH_w, &
                         BcoefH_x, BcoefH_w, &
                         AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
                         AcoefH, AcoefQ, BcoefH, BcoefQ,  &
                         HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
                         phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
                         yg_T, yg_Q, &
                         yGamma_dTs_phiT, yGamma_dQs_phiQ, &
                         ydTs_ins, ydqs_ins, &
                         y_flux_t1, y_flux_q1,  &
                         y_flux_t1_x, y_flux_t1_w, &
                         y_flux_q1_x, y_flux_q1_w )
         ENDIF   ! (nsrf .ne. is_oce)

       ELSE  ! (iflag_split .ge. 1)
         ybeta(1:knon) = 1.
         yevap_pot(1:knon) = yevap(1:knon)
       ENDIF  ! (iflag_split .ge. 1)

       IF (prt_level >= 10) THEN
         print *,'pbl_surface, ybeta , yevap, yevap_pot ', &
                               ybeta(1:knon) , yevap(1:knon), yevap_pot(1:knon)
       ENDIF  ! (prt_level >= 10)

       IF (iflag_split .ge. 1) THEN
       IF (prt_level >=10) THEN
        DO j = 1, knon
         print*,'Chx,Chw,Ch', ycdragh_x(j), ycdragh_w(j), ycdragh(j)
         print*,'Khx,Khw,Kh', Kech_h_x(j), Kech_h_w(j), Kech_h(j)
         print*,'t1x, t1w, t1, t1_ancien', &
 &               yt_x(j,1), yt_w(j,1),  yt(j,1), t(j,1)
         print*,'delta_coef,delta_flux,delta_tsurf,tau', delta_coef(j), y_delta_flux_t1(j), y_delta_tsurf(j), tau_eq(j) 
        ENDDO

        DO j=1,knon
         print*,'fluxT_x, fluxT_w, y_flux_t1, fluxQ_x, fluxQ_w, yfluxlat, wakes' &
 &             , y_flux_t1_x(j), y_flux_t1_w(j), y_flux_t1(j), y_flux_q1_x(j)*RLVTT, y_flux_q1_w(j)*RLVTT, yfluxlat(j), ywake_s(j)
         print*,'beta, ytsurf_new ', ybeta(j), ytsurf_new(j)
         print*,'inertia, facteur, cstar', inertia, facteur,wake_cstar(j)
        ENDDO
       ENDIF  ! (prt_level >=10)

       ENDIF  ! (iflag_split .ge.1)

       IF (iflag_split .eq.0) THEN

        CALL climb_hq_up(knon, ni, dtime, yt, yq, &
            y_flux_q1, y_flux_t1, ypaprs, ypplay, &
            AcoefH, AcoefQ, BcoefH, BcoefQ, &
            CcoefH, CcoefQ, DcoefH, DcoefQ, &
            Kcoef_hq, gama_q, gama_h, &
            y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:) &
#ifdef ISO
        &    ,yxt,y_flux_xt1 &
        &    ,AcoefXT,BcoefXT,CcoefXT,DcoefXT,gama_xt &
        &    ,y_flux_xt(:,:,:),y_d_xt(:,:,:) &
#endif
        &    )    
       ELSE  !(iflag_split .eq.0)

        CALL climb_hq_up(knon, ni, dtime, yt_x, yq_x, &
            y_flux_q1_x, y_flux_t1_x, ypaprs, ypplay, &
            AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x, &
            CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
            Kcoef_hq_x, gama_q_x, gama_h_x, &
            y_flux_q_x(:,:), y_flux_t_x(:,:), y_d_q_x(:,:), y_d_t_x(:,:) &
#ifdef ISO
        &    ,yxt_x,y_flux_xt1_x &
        &    ,AcoefXT_x,BcoefXT_x,CcoefXT_x,DcoefXT_x,gama_xt_x &
        &    ,y_flux_xt_x(:,:,:),y_d_xt_x(:,:,:) &
#endif
        &    )    
!

       CALL climb_hq_up(knon, ni, dtime, yt_w, yq_w, &
            y_flux_q1_w, y_flux_t1_w, ypaprs, ypplay, &
            AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w, &
            CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
            Kcoef_hq_w, gama_q_w, gama_h_w, &
            y_flux_q_w(:,:), y_flux_t_w(:,:), y_d_q_w(:,:), y_d_t_w(:,:) &
#ifdef ISO
        &    ,yxt_w,y_flux_xt1_w &
        &    ,AcoefXT_w,BcoefXT_w,CcoefXT_w,DcoefXT_w,gama_xt_w &
        &    ,y_flux_xt_w(:,:,:),y_d_xt_w(:,:,:) &
#endif
        &    )    
       ENDIF  ! (iflag_split .eq.0)

       IF (iflag_split .eq.0) THEN

        CALL climb_wind_up(knon, ni, dtime, yu, yv, y_flux_u1, y_flux_v1, &
            AcoefU, AcoefV, BcoefU, BcoefV, &
            CcoefU, CcoefV, DcoefU, DcoefV, &
            Kcoef_m, &
            y_flux_u, y_flux_v, y_d_u, y_d_v)
        
     
     y_d_t_diss(:,:)=0.
     IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN

        CALL yamada_c( knon, knon, dtime,ypaprs,ypplay &
    &   ,yu,yv,yt,y_d_u,y_d_v,y_d_t,ycdragm,ytke,ycoefm,ycoefh,ycoefq,y_d_t_diss,yustar &
    &   ,iflag_pbl)
     ENDIF

       ELSE  !(iflag_split .eq.0)

        CALL climb_wind_up( knon, ni, dtime, yu_x, yv_x, y_flux_u1_x, y_flux_v1_x, &
            AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x, &
            CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
            Kcoef_m_x, &
            y_flux_u_x, y_flux_v_x, y_d_u_x, y_d_v_x)

     y_d_t_diss_x(:,:)=0.
     IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN

        CALL yamada_c(knon, knon,dtime,ypaprs,ypplay &
    &   ,yu_x,yv_x,yt_x,y_d_u_x,y_d_v_x,y_d_t_x,ycdragm_x,ytke_x,ycoefm_x,ycoefh_x &
        ,ycoefq_x,y_d_t_diss_x,yustar_x &
    &   ,iflag_pbl)
     ENDIF

        CALL climb_wind_up(knon, ni, dtime, yu_w, yv_w, y_flux_u1_w, y_flux_v1_w, &
            AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w, &
            CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
            Kcoef_m_w, &
            y_flux_u_w, y_flux_v_w, y_d_u_w, y_d_v_w)

     y_d_t_diss_w(:,:)=0.
     IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN

        CALL yamada_c( knon, knon,dtime,ypaprs,ypplay &
    &   ,yu_w,yv_w,yt_w,y_d_u_w,y_d_v_w,y_d_t_w,ycdragm_w,ytke_w,ycoefm_w,ycoefh_w &
        ,ycoefq_w,y_d_t_diss_w,yustar_w &
    &   ,iflag_pbl)
     ENDIF

        IF (prt_level >=10) THEN
         print *, 'After climbing up, lfuxlat_x, fluxlat_w ', &
               yfluxlat_x(1:knon), yfluxlat_w(1:knon)
        ENDIF
!
       ENDIF  ! (iflag_split .eq.0)

       IF (ok_bs) THEN

            CALL climb_qbs_up(knon, ni, dtime, yqbs, &
            y_flux_bs, ypaprs, ypplay, &
            AcoefQBS, BcoefQBS, &
            CcoefQBS, DcoefQBS, &
            Kcoef_qbs, gama_qbs, &
            y_flux_qbs(:,:), y_d_qbs(:,:))
       ENDIF


!****************************************************************************************
! 13) Transform variables for output format : 
!     - Decompress
!     - Multiply with pourcentage of current surface
!     - Cumulate in global variable
!
!****************************************************************************************


       IF (iflag_split.EQ.0) THEN

        DO k = 1, klev
           DO j = 1, knon
             i = ni(j)
             y_d_t_diss(j,k)  = y_d_t_diss(j,k) * ypct(j)
             y_d_t(j,k)  = y_d_t(j,k) * ypct(j)
             y_d_q(j,k)  = y_d_q(j,k) * ypct(j)
             y_d_u(j,k)  = y_d_u(j,k) * ypct(j)
             y_d_v(j,k)  = y_d_v(j,k) * ypct(j)

             IF  (nsrf .EQ. is_ter .and. ifl_pbltree .GE. 1) THEN

               y_d_u(j,k) =y_d_u(j,k) + y_d_u_frein(j,k)*ypct(j)
               y_d_v(j,k) =y_d_v(j,k) + y_d_v_frein(j,k)*ypct(j)
               treedrg(i,k,nsrf)=y_treedrg(j,k)
             ELSE 
               treedrg(i,k,nsrf)=0.
             ENDIF

             flux_t(i,k,nsrf) = y_flux_t(j,k)
             flux_q(i,k,nsrf) = y_flux_q(j,k)
             flux_u(i,k,nsrf) = y_flux_u(j,k)
             flux_v(i,k,nsrf) = y_flux_v(j,k)

#ifdef ISO
             DO ixt=1,ntraciso
                y_d_xt(ixt,j,k)  = y_d_xt(ixt,j,k) * ypct(j)
                flux_xt(ixt,i,k,nsrf) = y_flux_xt(ixt,j,k)
             ENDDO ! DO ixt=1,ntraciso
             h1_diag(i)=h1(j)
#endif

           ENDDO
        ENDDO

#ifdef ISO
#ifdef ISOVERIF
        if (iso_eau.gt.0) then
         call iso_verif_egalite_vect2D( &
                y_d_xt,y_d_q, &
                'pbl_surface_mod 2600',ntraciso,klon,klev)
        endif       
#endif
#endif

       ELSE  !(iflag_split .eq.0)

! Tendances hors poches
        DO k = 1, klev
          DO j = 1, knon
            i = ni(j)
            y_d_t_diss_x(j,k)  = y_d_t_diss_x(j,k) * ypct(j)
            y_d_t_x(j,k)  = y_d_t_x(j,k) * ypct(j)
            y_d_q_x(j,k)  = y_d_q_x(j,k) * ypct(j)
            y_d_u_x(j,k)  = y_d_u_x(j,k) * ypct(j)
            y_d_v_x(j,k)  = y_d_v_x(j,k) * ypct(j)

            flux_t_x(i,k,nsrf) = y_flux_t_x(j,k)
            flux_q_x(i,k,nsrf) = y_flux_q_x(j,k)
            flux_u_x(i,k,nsrf) = y_flux_u_x(j,k)
            flux_v_x(i,k,nsrf) = y_flux_v_x(j,k)

#ifdef ISO
            DO ixt=1,ntraciso
              y_d_xt_x(ixt,j,k)  = y_d_xt_x(ixt,j,k) * ypct(j)
              flux_xt_x(ixt,i,k,nsrf) = y_flux_xt_x(ixt,j,k)
            ENDDO ! DO ixt=1,ntraciso
#endif
          ENDDO
        ENDDO

! Tendances dans les poches
        DO k = 1, klev
          DO j = 1, knon
            i = ni(j)
            y_d_t_diss_w(j,k)  = y_d_t_diss_w(j,k) * ypct(j)
            y_d_t_w(j,k)  = y_d_t_w(j,k) * ypct(j)
            y_d_q_w(j,k)  = y_d_q_w(j,k) * ypct(j)
            y_d_u_w(j,k)  = y_d_u_w(j,k) * ypct(j)
            y_d_v_w(j,k)  = y_d_v_w(j,k) * ypct(j)

            flux_t_w(i,k,nsrf) = y_flux_t_w(j,k)
            flux_q_w(i,k,nsrf) = y_flux_q_w(j,k)
            flux_u_w(i,k,nsrf) = y_flux_u_w(j,k)
            flux_v_w(i,k,nsrf) = y_flux_v_w(j,k)

#ifdef ISO
            DO ixt=1,ntraciso
              y_d_xt_w(ixt,j,k)  = y_d_xt_w(ixt,j,k) * ypct(j)
              flux_xt_w(ixt,i,k,nsrf) = y_flux_xt_w(ixt,j,k)
            ENDDO ! do ixt=1,ntraciso
#endif

          ENDDO
        ENDDO

! Flux, tendances et Tke moyenne dans la maille
        DO k = 1, klev
          DO j = 1, knon
            i = ni(j)
            flux_t(i,k,nsrf) = flux_t_x(i,k,nsrf)+ywake_s(j)*(flux_t_w(i,k,nsrf)-flux_t_x(i,k,nsrf))
            flux_q(i,k,nsrf) = flux_q_x(i,k,nsrf)+ywake_s(j)*(flux_q_w(i,k,nsrf)-flux_q_x(i,k,nsrf))
            flux_u(i,k,nsrf) = flux_u_x(i,k,nsrf)+ywake_s(j)*(flux_u_w(i,k,nsrf)-flux_u_x(i,k,nsrf))
            flux_v(i,k,nsrf) = flux_v_x(i,k,nsrf)+ywake_s(j)*(flux_v_w(i,k,nsrf)-flux_v_x(i,k,nsrf))
#ifdef ISO
            DO ixt=1,ntraciso
              flux_xt(ixt,i,k,nsrf) = flux_xt_x(ixt,i,k,nsrf)+ywake_s(j)*(flux_xt_w(ixt,i,k,nsrf)-flux_xt_x(ixt,i,k,nsrf))
            ENDDO ! do ixt=1,ntraciso
#endif
          ENDDO
        ENDDO
        DO j=1,knon
          yfluxlat(j)=yfluxlat_x(j)+ywake_s(j)*(yfluxlat_w(j)-yfluxlat_x(j))
        ENDDO
        IF (prt_level >=10) THEN
          print *,' nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf) ', &
                    nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf)
        ENDIF

        DO k = 1, klev
          DO j = 1, knon
            y_d_t_diss(j,k) = y_d_t_diss_x(j,k)+ywake_s(j)*(y_d_t_diss_w(j,k) -y_d_t_diss_x(j,k))
            y_d_t(j,k) = y_d_t_x(j,k)+ywake_s(j)*(y_d_t_w(j,k) -y_d_t_x(j,k))
            y_d_q(j,k) = y_d_q_x(j,k)+ywake_s(j)*(y_d_q_w(j,k) -y_d_q_x(j,k))
            y_d_u(j,k) = y_d_u_x(j,k)+ywake_s(j)*(y_d_u_w(j,k) -y_d_u_x(j,k))
            y_d_v(j,k) = y_d_v_x(j,k)+ywake_s(j)*(y_d_v_w(j,k) -y_d_v_x(j,k))
          ENDDO
        ENDDO

       ENDIF  ! (iflag_split .eq.0)


       ! tendencies of blowing snow 
       IF (ok_bs) THEN 
           DO k = 1, klev   
            DO j = 1, knon 
                i = ni(j)
                y_d_qbs(j,k)=y_d_qbs(j,k) * ypct(j)
                flux_qbs(i,k,nsrf) = y_flux_qbs(j,k) 
            ENDDO
          ENDDO
       ENDIF


       DO j = 1, knon
          i = ni(j)
          evap(i,nsrf) = - flux_q(i,1,nsrf)                  !jyg
          if (ok_bs) then ; snowerosion(i,nsrf)=flux_qbs(i,1,nsrf); endif
          beta(i,nsrf) = ybeta(j)                             !jyg
          d_ts(i,nsrf) = y_d_ts(j)
!albedo SB >>>
          DO k=1,nsw
            alb_dir(i,k,nsrf) = yalb_dir_new(j,k)
            alb_dif(i,k,nsrf) = yalb_dif_new(j,k)
          ENDDO
!albedo SB <<<
          snow(i,nsrf) = ysnow(j)  
          qsurf(i,nsrf) = yqsurf(j)
          z0m(i,nsrf) = yz0m(j)
          z0h(i,nsrf) = yz0h(j)
          fluxlat(i,nsrf) = yfluxlat(j)
          agesno(i,nsrf) = yagesno(j)  
          cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
          cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
          dflux_t(i) = dflux_t(i) + y_dflux_t(j)*ypct(j)
          dflux_q(i) = dflux_q(i) + y_dflux_q(j)*ypct(j)
#ifdef ISO
        DO ixt=1,niso
          xtsnow(ixt,i,nsrf) = yxtsnow(ixt,j)  
        ENDDO
        DO ixt=1,ntraciso
          xtevap(ixt,i,nsrf) = - flux_xt(ixt,i,1,nsrf)
          dflux_xt(ixt,i) = dflux_xt(ixt,i) + y_dflux_xt(ixt,j)*ypct(j)
        ENDDO  
        IF (nsrf == is_lic) THEN
          DO ixt=1,niso
            Rland_ice(ixt,i) = yRland_ice(ixt,j)  
          ENDDO
        ENDIF !IF (nsrf == is_lic) THEN     
#ifdef ISOVERIF
        IF (iso_eau.gt.0) THEN  
          call iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, &
     &         'pbl_surf_mod 1230',errmax,errmaxrel)
        ENDIF !if (iso_eau.gt.0) then
#endif        
#endif
       ENDDO

       IF (iflag_split .ge.1) THEN

        DO j = 1, knon
          i = ni(j)
          fluxlat_x(i,nsrf) = yfluxlat_x(j)
          fluxlat_w(i,nsrf) = yfluxlat_w(j)
          delta_tsurf(i,nsrf)=y_delta_tsurf_new(j)
          delta_qsurf(i) = delta_qsurf(i) + y_delta_qsurf(j)*ypct(j)
          cdragh_x(i) = cdragh_x(i) + ycdragh_x(j)*ypct(j)
          cdragh_w(i) = cdragh_w(i) + ycdragh_w(j)*ypct(j)
          cdragm_x(i) = cdragm_x(i) + ycdragm_x(j)*ypct(j)
          cdragm_w(i) = cdragm_w(i) + ycdragm_w(j)*ypct(j)
          kh(i) = kh(i) + Kech_h(j)*ypct(j)
          kh_x(i) = kh_x(i) + Kech_h_x(j)*ypct(j)
          kh_w(i) = kh_w(i) + Kech_h_w(j)*ypct(j)
        ENDDO
       ENDIF  ! (iflag_split .ge.1)

       IF (iflag_split .eq.0) THEN
        wake_dltke(:,:,nsrf) = 0.
        DO k = 1, klev+1
           DO j = 1, knon
              i = ni(j)
              tke_x(i,k,nsrf)    = ytke(j,k)
              tke_x(i,k,is_ave)  = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)
              eps_x(i,k,nsrf)    = yeps(j,k)
              eps_x(i,k,is_ave)  = eps_x(i,k,is_ave) + yeps(j,k)*ypct(j)
           ENDDO
        ENDDO

       ELSE  ! (iflag_split .eq.0)
        DO k = 1, klev+1
          DO j = 1, knon
            i = ni(j)
            wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
            tke_x(i,k,nsrf)   = ytke_x(j,k)
            tke_x(i,k,is_ave)   = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)       
            eps_x(i,k,nsrf)   = yeps_x(j,k)
            eps_x(i,k,is_ave)   = eps_x(i,k,is_ave) + eps_x(i,k,nsrf)*ypct(j) 
            wake_dltke(i,k,is_ave)   = wake_dltke(i,k,is_ave) + wake_dltke(i,k,nsrf)*ypct(j)
          ENDDO
        ENDDO
       ENDIF  ! (iflag_split .eq.0)

       DO k = 2, klev
          DO j = 1, knon
             i = ni(j)
             zcoefh(i,k,nsrf) = ycoefh(j,k)
             zcoefm(i,k,nsrf) = ycoefm(j,k)
             zcoefh(i,k,is_ave) = zcoefh(i,k,is_ave) + ycoefh(j,k)*ypct(j)
             zcoefm(i,k,is_ave) = zcoefm(i,k,is_ave) + ycoefm(j,k)*ypct(j)
          ENDDO
       ENDDO

       IF ( nsrf .EQ. is_ter ) THEN 
          DO j = 1, knon
             i = ni(j)
             qsol(i) = yqsol(j)
#ifdef ISO
             runoff_diag(i)=yrunoff_diag(j)   
             DO ixt=1,niso
               xtsol(ixt,i) = yxtsol(ixt,j)
               xtrunoff_diag(ixt,i)=yxtrunoff_diag(ixt,j)
             ENDDO
#endif
          ENDDO
       ENDIF
       
       DO k = 1, nsoilmx
          DO j = 1, knon
             i = ni(j)
             ftsoil(i, k, nsrf) = ytsoil(j,k)
          ENDDO
       ENDDO

#ifdef ISO
#ifdef ISOVERIF
       DO i = 1, klon
         DO ixt=1,niso
           call iso_verif_noNaN(xtsol(ixt,i),'pbl_surface 1405')
         ENDDO
       ENDDO
#endif
#ifdef ISOVERIF
     IF (iso_eau.gt.0) THEN
        call iso_verif_egalite_vect2D( &
                y_d_xt,y_d_q, &
                'pbl_surface_mod 1261',ntraciso,klon,klev)
     ENDIF !if (iso_eau.gt.0) then
#endif
#endif

       IF (iflag_split .ge.1) THEN

        DO k = 1, klev
          DO j = 1, knon
           i = ni(j)
           d_t_diss_x(i,k) = d_t_diss_x(i,k) + y_d_t_diss_x(j,k)
           d_t_x(i,k) = d_t_x(i,k) + y_d_t_x(j,k)
           d_q_x(i,k) = d_q_x(i,k) + y_d_q_x(j,k)
           d_u_x(i,k) = d_u_x(i,k) + y_d_u_x(j,k)
           d_v_x(i,k) = d_v_x(i,k) + y_d_v_x(j,k)
!
           d_t_diss_w(i,k) = d_t_diss_w(i,k) + y_d_t_diss_w(j,k)
           d_t_w(i,k) = d_t_w(i,k) + y_d_t_w(j,k)
           d_q_w(i,k) = d_q_w(i,k) + y_d_q_w(j,k)
           d_u_w(i,k) = d_u_w(i,k) + y_d_u_w(j,k)
           d_v_w(i,k) = d_v_w(i,k) + y_d_v_w(j,k)
#ifdef ISO
           DO ixt=1,ntraciso
             d_xt_x(ixt,i,k) = d_xt_x(ixt,i,k) + y_d_xt_x(ixt,j,k)
             d_xt_w(ixt,i,k) = d_xt_w(ixt,i,k) + y_d_xt_w(ixt,j,k)
           ENDDO ! DO ixt=1,ntraciso
#endif

          ENDDO
        ENDDO
      ENDIF  ! (iflag_split .ge.1)
       
       DO k = 1, klev
          DO j = 1, knon
             i = ni(j)
             d_t_diss(i,k) = d_t_diss(i,k) + y_d_t_diss(j,k)
             d_t(i,k) = d_t(i,k) + y_d_t(j,k)
             d_q(i,k) = d_q(i,k) + y_d_q(j,k)
#ifdef ISO
             DO ixt=1,ntraciso
               d_xt(ixt,i,k) = d_xt(ixt,i,k) + y_d_xt(ixt,j,k)
             ENDDO !DO ixt=1,ntraciso
#endif
             d_u(i,k) = d_u(i,k) + y_d_u(j,k)
             d_v(i,k) = d_v(i,k) + y_d_v(j,k)
          ENDDO
       ENDDO


       IF (ok_bs) THEN
         DO k = 1, klev
         DO j = 1, knon
         i = ni(j)
         d_qbs(i,k) = d_qbs(i,k) + y_d_qbs(j,k)
         ENDDO
         ENDDO
        ENDIF

#ifdef ISO
#ifdef ISOVERIF
        call iso_verif_noNaN_vect2D( &
     &           d_xt, &
     &           'pbl_surface 1385',ntraciso,klon,klev)  
     IF (iso_eau >= 0) THEN
        call iso_verif_egalite_vect2D( &
                y_d_xt,y_d_q, &
                'pbl_surface_mod 2945',ntraciso,klon,klev)
        call iso_verif_egalite_vect2D( &
                d_xt,d_q, &
                'pbl_surface_mod 1276',ntraciso,klon,klev)
     ENDIF !IF (iso_eau >= 0) THEN
#endif
#endif

       IF (prt_level >=10) THEN
         print *, 'pbl_surface tendencies for w: d_t_w, d_t_x, d_t ', &
          d_t_w(1:knon,1), d_t_x(1:knon,1), d_t(1:knon,1)
       ENDIF

       if (nsrf == is_oce .and. activate_ocean_skin >= 1) then
          delta_sal = missing_val
          ds_ns = missing_val
          dt_ns = missing_val
          delta_sst = missing_val
          dter = missing_val
          dser = missing_val
          tkt = missing_val
          tks = missing_val
          taur = missing_val
          sss = missing_val
          
          delta_sal(ni(:knon)) = ydelta_sal(:knon)
          ds_ns(ni(:knon)) = yds_ns(:knon)
          dt_ns(ni(:knon)) = ydt_ns(:knon)
          delta_sst(ni(:knon)) = ydelta_sst(:knon)
          dter(ni(:knon)) = ydter(:knon)
          dser(ni(:knon)) = ydser(:knon)
          tkt(ni(:knon)) = ytkt(:knon)
          tks(ni(:knon)) = ytks(:knon)
          taur(ni(:knon)) = ytaur(:knon)
          sss(ni(:knon)) = ysss(:knon)

          if (activate_ocean_skin == 2 .and. type_ocean == "couple") then
             dt_ds = missing_val
             dt_ds(ni(:knon)) = ydt_ds(:knon)
          end if
       end if


!****************************************************************************************
! 14) Calculate the temperature and relative humidity at 2m and the wind at 10m 
!     Call HBTM
!
!****************************************************************************************
!!!
!
#undef T2m     
#define T2m     
#ifdef T2m
! Calculations of diagnostic t,q at 2m and u, v at 10m

      IF (iflag_split .eq.0) THEN
        DO j=1, knon
          uzon(j) = yu(j,1) + y_d_u(j,1)
          vmer(j) = yv(j,1) + y_d_v(j,1)
          tair1(j) = yt(j,1) + y_d_t(j,1) + y_d_t_diss(j,1)
          qair1(j) = yq(j,1) + y_d_q(j,1)
          zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
               * (ypaprs(j,1)-ypplay(j,1))
          tairsol(j) = yts(j) + y_d_ts(j)
          qairsol(j) = yqsurf(j)
        ENDDO
       ELSE  ! (iflag_split .eq.0)
        DO j=1, knon
          uzon_x(j) = yu_x(j,1) + y_d_u_x(j,1)
          vmer_x(j) = yv_x(j,1) + y_d_v_x(j,1)
          tair1_x(j) = yt_x(j,1) + y_d_t_x(j,1) + y_d_t_diss_x(j,1)
          qair1_x(j) = yq_x(j,1) + y_d_q_x(j,1)
          zgeo1_x(j) = RD * tair1_x(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
               * (ypaprs(j,1)-ypplay(j,1))
          tairsol(j) = yts(j) + y_d_ts(j)
          tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf_new(j)
          qairsol(j) = yqsurf(j)
        ENDDO
        DO j=1, knon
          uzon_w(j) = yu_w(j,1) + y_d_u_w(j,1)
          vmer_w(j) = yv_w(j,1) + y_d_v_w(j,1)
          tair1_w(j) = yt_w(j,1) + y_d_t_w(j,1) + y_d_t_diss_w(j,1)
          qair1_w(j) = yq_w(j,1) + y_d_q_w(j,1)
          zgeo1_w(j) = RD * tair1_w(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
               * (ypaprs(j,1)-ypplay(j,1))
          tairsol_w(j) = tairsol(j) + (1.- ywake_s(j))*y_delta_tsurf(j)
          qairsol(j) = yqsurf(j)
        ENDDO
      ENDIF  ! (iflag_split .eq.0)

       DO j=1, knon
          psfce(j)=ypaprs(j,1)
          patm(j)=ypplay(j,1)
       ENDDO

       IF (iflag_pbl_surface_t2m_bug==1) THEN
          yz0h_oupas(1:knon)=yz0m(1:knon)
       ELSE
          yz0h_oupas(1:knon)=yz0h(1:knon)
       ENDIF
       

! Calculate the temperature and relative humidity at 2m and the wind at 10m 
       IF (iflag_split .eq.0) THEN
        IF (iflag_new_t2mq2m==1) THEN
           
           CALL stdlevvarn(knon, knon, nsrf, zxli, &
            uzon, vmer, tair1, qair1, zgeo1, &
            tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
            yt2m, yq2m, yt10m, yq10m, yu10m, yustar, &
            yn2mout(:, :, :))
        ELSE

        CALL stdlevvar(knon, knon, nsrf, zxli, &
            uzon, vmer, tair1, qair1, zgeo1, &
            tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
            yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, yzxtsol)
        ENDIF
       ELSE  !(iflag_split .eq.0)
        IF (iflag_new_t2mq2m==1) THEN

         CALL stdlevvarn(knon, knon, nsrf, zxli, &
            uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
            tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
            yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, &
            yn2mout_x(:, :, :))

         CALL stdlevvarn(knon, knon, nsrf, zxli, &
            uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
            tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
            yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, &
            yn2mout_w(:, :, :))
        ELSE

        CALL stdlevvar(knon, knon, nsrf, zxli, &
            uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
            tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
            yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, ypblh_x, rain_f, yzxtsol)

        CALL stdlevvar(knon, knon, nsrf, zxli, &
            uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
            tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
            yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, ypblh_w, rain_f, yzxtsol)
        ENDIF

       ENDIF  ! (iflag_split .eq.0)

       IF (iflag_split .eq.0) THEN
        DO j=1, knon
          i = ni(j)
          t2m(i,nsrf)=yt2m(j)
          q2m(i,nsrf)=yq2m(j)
     ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
          ustar(i,nsrf)=yustar(j)
          u10m(i,nsrf)=(yu10m(j) * uzon(j))/max(SQRT(uzon(j)**2+vmer(j)**2), smallestreal)
          v10m(i,nsrf)=(yu10m(j) * vmer(j))/max(SQRT(uzon(j)**2+vmer(j)**2), smallestreal)

          DO k = 1, 6
           n2mout(i,nsrf,k) = yn2mout(j,nsrf,k)
          END DO  

        ENDDO
       ELSE  !(iflag_split .eq.0)
        DO j=1, knon
          i = ni(j)
          t2m_x(i,nsrf)=yt2m_x(j)
          q2m_x(i,nsrf)=yq2m_x(j)
     ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
          ustar_x(i,nsrf)=yustar_x(j)
          u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/max(SQRT(uzon_x(j)**2+vmer_x(j)**2), smallestreal)
          v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/max(SQRT(uzon_x(j)**2+vmer_x(j)**2), smallestreal)

          DO k = 1, 6
           n2mout_x(i,nsrf,k) = yn2mout_x(j,nsrf,k)
          END DO  

        ENDDO
        DO j=1, knon
          i = ni(j)
          t2m_w(i,nsrf)=yt2m_w(j)
          q2m_w(i,nsrf)=yq2m_w(j)
     ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
          ustar_w(i,nsrf)=yustar_w(j)
          u10m_w(i,nsrf)=(yu10m_w(j) * uzon_w(j))/max(SQRT(uzon_w(j)**2+vmer_w(j)**2), smallestreal)
          v10m_w(i,nsrf)=(yu10m_w(j) * vmer_w(j))/max(SQRT(uzon_w(j)**2+vmer_w(j)**2), smallestreal)

          ustar(i,nsrf) = ustar_x(i,nsrf) + wake_s(i)*(ustar_w(i,nsrf)-ustar_x(i,nsrf))
          u10m(i,nsrf) = u10m_x(i,nsrf) + wake_s(i)*(u10m_w(i,nsrf)-u10m_x(i,nsrf))
          v10m(i,nsrf) = v10m_x(i,nsrf) + wake_s(i)*(v10m_w(i,nsrf)-v10m_x(i,nsrf))

          DO k = 1, 6
           n2mout_w(i,nsrf,k) = yn2mout_w(j,nsrf,k)
          END DO  

        ENDDO

       ENDIF  ! (iflag_split .eq.0)


!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
!IM Ajoute dependance type surface
       IF (thermcep) THEN

       IF (iflag_split .eq.0) THEN
          DO j = 1, knon
             i=ni(j)
             zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))
             zx_qs1  = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)
             zx_qs1  = MIN(0.5,zx_qs1)
             zcor1   = 1./(1.-RETV*zx_qs1)
             zx_qs1  = zx_qs1*zcor1
             
             rh2m(i)   = rh2m(i)   + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
             qsat2m(i) = qsat2m(i) + zx_qs1  * pctsrf(i,nsrf)
          ENDDO
       ELSE  ! (iflag_split .eq.0)
          DO j = 1, knon
             i=ni(j)
             zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_x(j) ))
             zx_qs1  = r2es * FOEEW(yt2m_x(j),zdelta1)/paprs(i,1)
             zx_qs1  = MIN(0.5,zx_qs1)
             zcor1   = 1./(1.-RETV*zx_qs1)
             zx_qs1  = zx_qs1*zcor1
             
             rh2m_x(i)   = rh2m_x(i)   + yq2m_x(j)/zx_qs1 * pctsrf(i,nsrf)
             qsat2m_x(i) = qsat2m_x(i) + zx_qs1  * pctsrf(i,nsrf)
          ENDDO
          DO j = 1, knon
             i=ni(j)
             zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_w(j) ))
             zx_qs1  = r2es * FOEEW(yt2m_w(j),zdelta1)/paprs(i,1)
             zx_qs1  = MIN(0.5,zx_qs1)
             zcor1   = 1./(1.-RETV*zx_qs1)
             zx_qs1  = zx_qs1*zcor1
             
             rh2m_w(i)   = rh2m_w(i)   + yq2m_w(j)/zx_qs1 * pctsrf(i,nsrf)
             qsat2m_w(i) = qsat2m_w(i) + zx_qs1  * pctsrf(i,nsrf)
          ENDDO

       ENDIF  ! (iflag_split .eq.0)

       ENDIF
!
       IF (prt_level >=10) THEN
         print *, 'T2m, q2m, RH2m ', &
          t2m(1:knon,:), q2m(1:knon,:), rh2m(1:knon)
       ENDIF


       IF (iflag_split .eq.0) THEN

        CALL hbtm(knon, ypaprs, ypplay, &
            yt2m,yt10m,yq2m,yq10m,yustar,ywstar, &
            y_flux_t,y_flux_q,yu,yv,yt,yq, &
            ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
            ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
          IF (prt_level >=10) THEN
       print *,' Arg. de HBTM: yt2m ',yt2m(1:knon)
       print *,' Arg. de HBTM: yt10m ',yt10m(1:knon)
       print *,' Arg. de HBTM: yq2m ',yq2m(1:knon)
       print *,' Arg. de HBTM: yq10m ',yq10m(1:knon)
       print *,' Arg. de HBTM: yustar ',yustar(1:knon)
       print *,' Arg. de HBTM: y_flux_t ',y_flux_t(1:knon,:)
       print *,' Arg. de HBTM: y_flux_q ',y_flux_q(1:knon,:)
       print *,' Arg. de HBTM: yu ',yu(1:knon,:)
       print *,' Arg. de HBTM: yv ',yv(1:knon,:)
       print *,' Arg. de HBTM: yt ',yt(1:knon,:)
       print *,' Arg. de HBTM: yq ',yq(1:knon,:)
          ENDIF
       ELSE  ! (iflag_split .eq.0)

        CALL HBTM(knon, ypaprs, ypplay, &
            yt2m_x,yt10m_x,yq2m_x,yq10m_x,yustar_x,ywstar_x, &
            y_flux_t_x,y_flux_q_x,yu_x,yv_x,yt_x,yq_x, &
            ypblh_x,ycapCL_x,yoliqCL_x,ycteiCL_x,ypblT_x, &
            ytherm_x,ytrmb1_x,ytrmb2_x,ytrmb3_x,ylcl_x)
          IF (prt_level >=10) THEN
       print *,' Arg. de HBTM: yt2m_x ',yt2m_x(1:knon)
       print *,' Arg. de HBTM: yt10m_x ',yt10m_x(1:knon)
       print *,' Arg. de HBTM: yq2m_x ',yq2m_x(1:knon)
       print *,' Arg. de HBTM: yq10m_x ',yq10m_x(1:knon)
       print *,' Arg. de HBTM: yustar_x ',yustar_x(1:knon)
       print *,' Arg. de HBTM: y_flux_t_x ',y_flux_t_x(1:knon,:)
       print *,' Arg. de HBTM: y_flux_q_x ',y_flux_q_x(1:knon,:)
       print *,' Arg. de HBTM: yu_x ',yu_x(1:knon,:)
       print *,' Arg. de HBTM: yv_x ',yv_x(1:knon,:)
       print *,' Arg. de HBTM: yt_x ',yt_x(1:knon,:)
       print *,' Arg. de HBTM: yq_x ',yq_x(1:knon,:)
          ENDIF

        CALL HBTM(knon, ypaprs, ypplay, &
            yt2m_w,yt10m_w,yq2m_w,yq10m_w,yustar_w,ywstar_w, &
            y_flux_t_w,y_flux_q_w,yu_w,yv_w,yt_w,yq_w, &
            ypblh_w,ycapCL_w,yoliqCL_w,ycteiCL_w,ypblT_w, &
            ytherm_w,ytrmb1_w,ytrmb2_w,ytrmb3_w,ylcl_w)
      
       ENDIF  ! (iflag_split .eq.0)

       IF (iflag_split .eq.0) THEN

        DO j=1, knon
          i = ni(j)
          pblh(i,nsrf)   = ypblh(j)
          wstar(i,nsrf)  = ywstar(j)
          plcl(i,nsrf)   = ylcl(j)
          capCL(i,nsrf)  = ycapCL(j)
          oliqCL(i,nsrf) = yoliqCL(j)
          cteiCL(i,nsrf) = ycteiCL(j)
          pblT(i,nsrf)   = ypblT(j)
          therm(i,nsrf)  = ytherm(j)
          trmb1(i,nsrf)  = ytrmb1(j)
          trmb2(i,nsrf)  = ytrmb2(j)
          trmb3(i,nsrf)  = ytrmb3(j)
        ENDDO
        IF (prt_level >=10) THEN
          print *, 'After HBTM: pblh ', pblh(1:knon,:)
          print *, 'After HBTM: plcl ', plcl(1:knon,:)
          print *, 'After HBTM: cteiCL ', cteiCL(1:knon,:)
        ENDIF
       ELSE  !(iflag_split .eq.0)
        DO j=1, knon
          i = ni(j)
          pblh_x(i,nsrf)   = ypblh_x(j)
          wstar_x(i,nsrf)  = ywstar_x(j)
          plcl_x(i,nsrf)   = ylcl_x(j)
          capCL_x(i,nsrf)  = ycapCL_x(j)
          oliqCL_x(i,nsrf) = yoliqCL_x(j)
          cteiCL_x(i,nsrf) = ycteiCL_x(j)
          pblT_x(i,nsrf)   = ypblT_x(j)
          therm_x(i,nsrf)  = ytherm_x(j)
          trmb1_x(i,nsrf)  = ytrmb1_x(j)
          trmb2_x(i,nsrf)  = ytrmb2_x(j)
          trmb3_x(i,nsrf)  = ytrmb3_x(j)
        ENDDO
        IF (prt_level >=10) THEN
          print *, 'After HBTM: pblh_x ', pblh_x(1:knon,:)
          print *, 'After HBTM: plcl_x ', plcl_x(1:knon,:)
          print *, 'After HBTM: cteiCL_x ', cteiCL_x(1:knon,:)
        ENDIF
        DO j=1, knon
          i = ni(j)
          pblh_w(i,nsrf)   = ypblh_w(j)
          wstar_w(i,nsrf)  = ywstar_w(j)
          plcl_w(i,nsrf)   = ylcl_w(j)
          capCL_w(i,nsrf)  = ycapCL_w(j)
          oliqCL_w(i,nsrf) = yoliqCL_w(j)
          cteiCL_w(i,nsrf) = ycteiCL_w(j)
          pblT_w(i,nsrf)   = ypblT_w(j)
          therm_w(i,nsrf)  = ytherm_w(j)
          trmb1_w(i,nsrf)  = ytrmb1_w(j)
          trmb2_w(i,nsrf)  = ytrmb2_w(j)
          trmb3_w(i,nsrf)  = ytrmb3_w(j)
        ENDDO
        IF (prt_level >=10) THEN
          print *, 'After HBTM: pblh_w ', pblh_w(1:knon,:)
          print *, 'After HBTM: plcl_w ', plcl_w(1:knon,:)
          print *, 'After HBTM: cteiCL_w ', cteiCL_w(1:knon,:)
        ENDIF

       ENDIF  ! (iflag_split .eq.0)

#else 
! T2m not defined
! No calculation
       PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
#endif

!****************************************************************************************
! 15) End of loop over different surfaces
!
!****************************************************************************************
!    ENDDO loop_nbsrf
  END SUBROUTINE pbl_surface_subsrf

END MODULE pbl_surface_subsrf_mod
