MODULE pbl_surface_uncompress_pre_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_uncompress_pre

CONTAINS

  SUBROUTINE pbl_surface_uncompress_pre( &
       itap,          &
       solsw_m,  solswfdiff_m, sollw_m,       &
           paprs,     pctsrf,                  &
       ts,SFRWL,   alb_dir, alb_dif,ustar, u10m, v10m,wstar, &
       cdragh,    cdragm,   zu1,    zv1,              &
       alb_dir_m,    alb_dif_m,  zxsens,   zxevap,  zxsnowerosion,      &
       icesub_ice, icemelt_ice, alb3_lic,  runoff,    snowhgt,   qsnow,     to_ice,    sissnow,  &
       zxtsol,    zxfluxlat, zt2m,     qsat2m, zn2mout,                 &
       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,                             & 
       zxsens_x,  zxfluxlat_x,zxsens_w,zxfluxlat_w,  &
       cdragh_x,cdragh_w,      &
       cdragm_x,cdragm_w,kh,kh_x,kh_w,               &
       zcoefh,    zcoefm,    slab_wfbils,            &
       qsol,    zq2m,      s_pblh,   s_plcl,         &
       s_pblh_x, s_plcl_x,   s_pblh_w, s_plcl_w,     &
       s_capCL,   s_oliqCL,  s_cteiCL, s_pblT,       &
       s_therm,   s_trmb1,   s_trmb2,  s_trmb3,      &
       zustar,zu10m,  zv10m,    fder_print,          &
       zxqsurf, delta_qsurf,                         &
       rh2m,      zxfluxu,  zxfluxv,                 &
       z0m, z0h,     sollw,    solsw,         &
       d_ts,      evap,    fluxlat,   t2m,           &
       wfbils,    wfevap,                            & 
       flux_t,   flux_u, flux_v,                     &
       dflux_t,   dflux_q,   zxsnow,                 &
       zxfluxt,   zxfluxq, zxfluxqbs,   q2m, flux_q, flux_qbs, tke_x, eps_x, &
       wake_dltke, 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, 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, snowerosion, alb &          
#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      
     &   )
!$gpum horizontal 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, carbon_cycle_tr 
    USE carbon_cycle_mod,   ONLY : co2_send, nbcf_out, fields_out, cfname_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
#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:  frac_tersrf, albedo_tersrf !AM 
    USE dimsoil_mod_h, ONLY: nsoilmx
    USE surf_param_mod, ONLY: eff_surf_param  !AM
    USE yomcst_mod_h
    USE phys_local_var_mod, only: l_mixmin, l_mix, wprime
IMPLICIT NONE

    INCLUDE "FCTTRE.h"
!FC

!****************************************************************************************
    INTEGER,                      INTENT(IN)        :: itap    ! time step
    REAL, DIMENSION(klon),        INTENT(IN)        :: solsw_m ! net shortwave radiation at mean surface
    REAL, DIMENSION(klon),        INTENT(IN)        :: solswfdiff_m ! diffuse fraction fordownward shortwave radiation at mean surface
    REAL, DIMENSION(klon),        INTENT(IN)        :: sollw_m ! net longwave radiation at mean surface
    REAL, DIMENSION(klon,klev+1), INTENT(IN)        :: paprs   ! pression between layers (Pa) 
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)        :: pctsrf  ! sub-surface fraction
#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
#ifdef ISO
    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(IN)        :: wake_dlxt   
#endif
! Input/Output variables
!****************************************************************************************
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ts      ! temperature at surface (K)
    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(OUT)   :: eps_x      ! TKE dissipation rate

    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh     ! drag coefficient for T and Q
    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragm     ! drag coefficient for wind
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zu1        ! u wind speed in first layer
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zv1        ! v wind speed in first layer
    REAL, DIMENSION(klon, nsw),   INTENT(OUT)       :: alb_dir_m,alb_dif_m
    REAL, DIMENSION(klon),        INTENT(OUT)       :: alb3_lic
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsens     ! sensible heat flux at surface with inversed sign 
                                                                  ! (=> positive sign upwards)
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxevap     ! water vapour flux at surface, positiv upwards
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsnowerosion     ! blowing snow flux at surface
    REAL, DIMENSION(klon),        INTENT(OUT)       :: icesub_ice ! ice (no snow!) sublimation flux over iced surfaces [kg/m2/s]
    REAL, DIMENSION(klon),        INTENT(OUT)       :: icemelt_ice ! ice (no snow!) melting flux over iced surfaces [kg/m2/s]
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxtsol     ! temperature at surface, mean for each grid point
    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_t_w      !   !
    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_q_w      !      !  Tendances dans les poches
    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_t_x      !   !
    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_q_x      !      !  Tendances hors des poches
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat  ! latent flux, mean for each grid point
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zt2m       ! temperature at 2m, mean for each grid point
    INTEGER, DIMENSION(klon, 6),  INTENT(OUT)       :: zn2mout    ! number of times the 2m temperature is out of the [tsol,temp]
    REAL, DIMENSION(klon),        INTENT(OUT)       :: qsat2m
    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_t        ! change in temperature 
    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_t_diss       ! change in temperature 
    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_q        ! change in water vapour
    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_u        ! change in u speed
    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_v        ! change in v speed
    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_qbs        ! change in blowing snow specific content
    REAL, INTENT(OUT):: 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(OUT):: 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(OUT)       :: zxxtevap     ! water vapour flux at surface, positiv upwards
    REAL, DIMENSION(ntraciso,klon, klev),  INTENT(OUT)       :: d_xt        ! change in water vapour
    REAL, DIMENSION(klon),                 INTENT(OUT)       :: runoff_diag
    REAL, DIMENSION(niso,klon),            INTENT(OUT)       :: xtrunoff_diag
    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)       :: d_xt_w
    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)       :: d_xt_x
#endif
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsens_x   ! Flux sensible hors poche
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsens_w   ! Flux sensible dans la poche
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat_x! Flux latent hors poche
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat_w! Flux latent dans la poche
! Output only for diagnostics
    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh_x
    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh_w
    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragm_x
    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragm_w
    REAL, DIMENSION(klon),        INTENT(OUT)       :: kh
    REAL, DIMENSION(klon),        INTENT(OUT)       :: kh_x
    REAL, DIMENSION(klon),        INTENT(OUT)       :: kh_w
    REAL, DIMENSION(klon),        INTENT(OUT)       :: slab_wfbils! heat balance at surface only for slab at ocean points
    REAL, DIMENSION(klon),        INTENT(OUT)       :: qsol     ! water height in the soil (mm)
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zq2m       ! water vapour at 2m, mean for each grid point
    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblh     ! height of the planetary boundary layer(HPBL)
    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblh_x   ! height of the PBL in the off-wake region
    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblh_w   ! height of the PBL in the wake region
    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_plcl     ! condensation level
    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_plcl_x   ! condensation level in the off-wake region
    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_plcl_w   ! condensation level in the wake region
    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_capCL    ! CAPE of PBL
    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_oliqCL   ! liquid water intergral of PBL
    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_cteiCL   ! cloud top instab. crit. of PBL
    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblT     ! temperature at PBLH
    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_therm    ! thermal virtual temperature excess
    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb1    ! deep cape, mean for each grid point
    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb2    ! inhibition, mean for each grid point
    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb3    ! point Omega, mean for each grid point
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zustar     ! u*
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zu10m      ! u speed at 10m, mean for each grid point
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zv10m      ! v speed at 10m, mean for each grid point
    REAL, DIMENSION(klon),        INTENT(OUT)       :: fder_print ! fder for printing (=fder(i) + dflux_t(i) + dflux_q(i))
    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxqsurf    ! humidity at surface, mean for each grid point
    REAL, DIMENSION(klon),        INTENT(OUT)       :: delta_qsurf! humidity difference at surface, mean for each grid point
    REAL, DIMENSION(klon),        INTENT(OUT)       :: rh2m       ! relative humidity at 2m
    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zxfluxu    ! u wind tension, mean for each grid point
    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zxfluxv    ! v wind tension, mean for each grid point
    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(OUT)       :: solsw      ! net shortwave radiation at surface 
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: sollw      ! net longwave radiation at surface
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: d_ts       ! change in temperature at surface
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: evap       ! evaporation at surface
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: fluxlat    ! latent flux
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: t2m        ! temperature at 2 meter height
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: wfbils     ! heat balance at surface
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: wfevap     ! water balance (evap) at surface weighted by srf
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t     ! sensible heat flux (CpT) J/m**2/s (W/m**2)
                                                                  ! positve orientation downwards
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u     ! u wind tension (kg m/s)/(m**2 s) or Pascal
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_v     ! v wind tension (kg m/s)/(m**2 s) or Pascal
#ifdef ISO        
    REAL, DIMENSION(niso,klon),   INTENT(OUT)       :: xtsol      ! water height in the soil (mm)
    REAL, DIMENSION(ntraciso,klon, nbsrf)           :: xtevap     ! evaporation at surface
    REAL, DIMENSION(klon),        INTENT(OUT)       :: h1_diag    ! just diagnostic, not useful
#endif

! Output not needed
    REAL, DIMENSION(klon),       INTENT(OUT)        :: dflux_t    ! change of sensible heat flux 
    REAL, DIMENSION(klon),       INTENT(OUT)        :: dflux_q    ! change of water vapour flux
    REAL, DIMENSION(klon),       INTENT(OUT)        :: zxsnow     ! snow at surface, mean for each grid point
    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxt    ! sensible heat flux, mean for each grid point
    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxq    ! water vapour flux, mean for each grid point
    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxqbs    ! blowing snow flux, mean for each grid point
    REAL, DIMENSION(klon, nbsrf),INTENT(OUT)        :: q2m        ! water vapour at 2 meter height
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q     ! water vapour flux(latent flux) (kg/m**2/s)
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_qbs   ! blowind snow vertical flux (kg/m**2

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

    REAL, DIMENSION(klon),       INTENT(OUT)        :: qsnow      ! snow water content
    REAL, DIMENSION(klon),       INTENT(OUT)        :: snowhgt    ! snow height
    REAL, DIMENSION(klon),       INTENT(OUT)        :: to_ice     ! snow passed to ice
    REAL, DIMENSION(klon),       INTENT(OUT)        :: sissnow    ! snow in snow model
    REAL, DIMENSION(klon),       INTENT(OUT)        :: runoff     ! runoff on land ice
    INTEGER,                     INTENT(INOUT)      :: iflag_split_ref 
! Other local variables
!****************************************************************************************
    INTEGER                            :: n 
    INTEGER                            :: iflag_split
    INTEGER                            :: i, k, nsrf 
    REAL                               :: f1 ! fraction de longeurs visibles parmi tout SW intervalle
    REAL, DIMENSION(klon)              :: r_co2_ppm     ! taux CO2 atmosphere
    REAL, DIMENSION(klon)              :: ztsol
    REAL, DIMENSION(klon)              :: meansqT ! mean square deviation of subsurface temperatures
    REAL, DIMENSION(klon)              :: alb_m  ! mean albedo for whole SW interval
    REAL, DIMENSION(klon,klev), INTENT(OUT)         :: delp
    LOGICAL, PARAMETER                 :: zxli=.FALSE. ! utiliser un jeu de fonctions simples
    LOGICAL, PARAMETER                 :: check=.FALSE.
    REAL, DIMENSION(klon,klev), INTENT(OUT)         :: d_t_diss_x, d_t_diss_w
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t_x, flux_q_x, flux_t_w, flux_q_w
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u_x, flux_v_x, flux_u_w, flux_v_w
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: fluxlat_x, fluxlat_w
    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: d_u_x
    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: d_u_w
    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: d_v_x
    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: d_v_w 
#ifdef ISO
    REAL, DIMENSION(ntraciso,klon,klev,nbsrf), INTENT(OUT)   :: flux_xt_x, flux_xt_w
#endif
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: windsp
!
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: t2m_x
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: q2m_x
    REAL, DIMENSION(klon), INTENT(OUT)              :: rh2m_x
    REAL, DIMENSION(klon), INTENT(OUT)              :: qsat2m_x
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: u10m_x
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: v10m_x
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: ustar_x
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: wstar_x
!              
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: pblh_x
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: plcl_x
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: capCL_x
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: oliqCL_x
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: cteiCL_x
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: pblt_x
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: therm_x
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: trmb1_x
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: trmb2_x
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: trmb3_x
!
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: t2m_w
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: q2m_w
    REAL, DIMENSION(klon) , INTENT(OUT)             :: rh2m_w
    REAL, DIMENSION(klon), INTENT(OUT)              :: qsat2m_w
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: u10m_w
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: v10m_w
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: ustar_w
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: wstar_w
!                           
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: pblh_w
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: plcl_w
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: capCL_w
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: oliqCL_w
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: cteiCL_w
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: pblt_w
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: therm_w
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: trmb1_w
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: trmb2_w
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: trmb3_w

    REAL, PARAMETER                    :: facteur = 2. / 1.772  ! ( == 2. / SQRT(3.14))
    REAL, PARAMETER                    :: inertia=2000.
    REAL, DIMENSION(klon,nbsrf), INTENT(OUT)        :: pblh         ! height of the planetary boundary layer
    REAL, DIMENSION(klon,nbsrf), INTENT(OUT)        :: plcl         ! condensation level
    REAL, DIMENSION(klon,nbsrf), INTENT(OUT)        :: capCL
    REAL, DIMENSION(klon,nbsrf), INTENT(OUT)        :: oliqCL
    REAL, DIMENSION(klon,nbsrf), INTENT(OUT)        :: cteiCL
    REAL, DIMENSION(klon,nbsrf), INTENT(OUT)        :: pblT
    REAL, DIMENSION(klon,nbsrf), INTENT(OUT)        :: therm
    REAL, DIMENSION(klon,nbsrf), INTENT(OUT)        :: trmb1        ! deep cape
    REAL, DIMENSION(klon,nbsrf), INTENT(OUT)        :: trmb2        ! inhibition
    REAL, DIMENSION(klon,nbsrf), INTENT(OUT)        :: trmb3        ! point Omega
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: alb          ! mean albedo for whole SW interval
    REAL, DIMENSION(klon,nbsrf), INTENT(OUT)        :: snowerosion    
    REAL, DIMENSION(klon) ::  albedo_eff
#ifdef ISO
    INTEGER                     :: ixt
#endif

!****************************************************************************************
! End of declarations
!****************************************************************************************

      IF (prt_level >=10) print *,' -> pbl_surface, itap ',itap
!
!!jyg      iflag_split = mod(iflag_pbl_split,2)
!!jyg      iflag_split = mod(iflag_pbl_split,10)
!
! Flags controlling the splitting of the turbulent boundary layer:
!   iflag_split_ref = 0  ==> no splitting
!                   = 1  ==> splitting without coupling with surface temperature
!                   = 2  ==> splitting with coupling with surface temperature over land
!                   = 3  ==> splitting over ocean; no splitting over land
!   iflag_split: actual flag controlling the splitting.
!   iflag_split = iflag_split_ref outside the sub-surface loop
!               = iflag_split_ref if iflag_split_ref = 0, 1, or 2
!               = 0 over land  if iflga_split_ref = 3
!               = 1 over ocean if iflga_split_ref = 3

      iflag_split_ref = mod(iflag_pbl_split,10)
      iflag_split = iflag_split_ref

#ifdef ISO      
#ifdef ISOVERIF
      DO i=1,klon
        DO ixt=1,niso
          CALL iso_verif_noNaN(xtsol(ixt,i),'pbl_surface 608')
        ENDDO
      ENDDO
#endif
#ifdef ISOVERIF
      DO i=1,klon  
        IF (iso_eau >= 0) THEN  
          CALL iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, &
     &         'pbl_surf_mod 585',errmax,errmaxrel)
          CALL iso_verif_egalite_choix(xtsnow_f(iso_eau,i),snow_f(i), &
     &         'pbl_surf_mod 594',errmax,errmaxrel)
          IF (iso_verif_egalite_choix_nostop(xtsol(iso_eau,i),qsol(i), &
     &         'pbl_surf_mod 596',errmax,errmaxrel) == 1) THEN
                WRITE(*,*) 'i=',i
                STOP
          ENDIF
          DO nsrf=1,nbsrf
            CALL iso_verif_egalite_choix(xtsnow(iso_eau,i,nsrf),snow(i,nsrf), &
     &         'pbl_surf_mod 598',errmax,errmaxrel)
          ENDDO
        ENDIF !IF (iso_eau >= 0) THEN   
      ENDDO !DO i=1,knon  
      DO k=1,klev
        DO i=1,klon  
          IF (iso_eau >= 0) THEN  
            CALL iso_verif_egalite_choix(xt(iso_eau,i,k),q(i,k), &
     &           'pbl_surf_mod 595',errmax,errmaxrel)
          ENDIF !IF (iso_eau >= 0) THEN  
        ENDDO !DO i=1,knon  
      ENDDO !DO k=1,klev
#endif
#endif


         
!****************************************************************************************
! Force soil water content to qsol0 if qsol0>0 and VEGET=F (use bucket
! instead of ORCHIDEE)
    IF (qsol0>=0.) THEN
      PRINT*,'WARNING : On impose qsol=',qsol0
      qsol(:)=qsol0
#ifdef ISO
      DO ixt=1,niso
        xtsol(ixt,:)=qsol0*Rdefault(ixt)
      ENDDO
#ifdef ISOTRAC      
      DO ixt=1+niso,ntraciso
        xtsol(ixt,:)=qsol0*Rdefault(index_iso(ixt))
      ENDDO
#endif       
#endif
    ENDIF
!****************************************************************************************

!****************************************************************************************
! 2) Initialization to zero 
!****************************************************************************************
!
! 2a) Initialization of all argument variables with INTENT(OUT)
!****************************************************************************************
 cdragh(:)=0. ; cdragm(:)=0.
 zu1(:)=0. ; zv1(:)=0.
!albedo SB >>>
  alb_dir_m=0. ; alb_dif_m=0. ; alb3_lic(:)=0.
!albedo SB <<<
 zxsens(:)=0. ; zxevap(:)=0. ; zxtsol(:)=0. ; zxsnowerosion(:)=0.
 d_t_w(:,:)=0. ; d_q_w(:,:)=0. ; d_t_x(:,:)=0. ; d_q_x(:,:)=0.
 zxfluxlat(:)=0.
 zt2m(:)=0. ; zq2m(:)=0. ; qsat2m(:)=0. ; rh2m(:)=0.
 zn2mout(:,:)=0 ;
 d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_qbs(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0.
 zcoefh(:,:,:)=0. ; zcoefm(:,:,:)=0.
 zxsens_x(:)=0. ; zxsens_w(:)=0. ; zxfluxlat_x(:)=0. ; zxfluxlat_w(:)=0.
 cdragh_x(:)=0. ; cdragh_w(:)=0. ; cdragm_x(:)=0. ; cdragm_w(:)=0.
 kh(:)=0. ; kh_x(:)=0. ; kh_w(:)=0.
 slab_wfbils(:)=0.
 s_pblh(:)=0. ; s_pblh_x(:)=0. ; s_pblh_w(:)=0.
 s_plcl(:)=0. ; s_plcl_x(:)=0. ; s_plcl_w(:)=0.
 s_capCL(:)=0. ; s_oliqCL(:)=0. ; s_cteiCL(:)=0. ; s_pblT(:)=0.
 s_therm(:)=0.
 s_trmb1(:)=0. ; s_trmb2(:)=0. ; s_trmb3(:)=0.
 zustar(:)=0.
 zu10m(:)=0. ; zv10m(:)=0.
 fder_print(:)=0.
 zxqsurf(:)=0.
 delta_qsurf(:) = 0.
 zxfluxu(:,:)=0. ; zxfluxv(:,:)=0.
 solsw(:,:)=0. ; sollw(:,:)=0.
 d_ts(:,:)=0.
 evap(:,:)=0.
 snowerosion(:,:)=0. 
 fluxlat(:,:)=0.
 wfbils(:,:)=0. ; wfevap(:,:)=0. ;
 flux_t(:,:,:)=0. ; flux_q(:,:,:)=0. ; flux_u(:,:,:)=0. ; flux_v(:,:,:)=0.
 flux_qbs(:,:,:)=0.
 dflux_t(:)=0. ; dflux_q(:)=0.
 zxsnow(:)=0.
 zxfluxt(:,:)=0. ; zxfluxq(:,:)=0.; zxfluxqbs(:,:)=0.
 qsnow(:)=0. ; snowhgt(:)=0. ; to_ice(:)=0. ; sissnow(:)=0.
 runoff(:)=0. ; icesub_ice(:)=0. ; icemelt_ice(:)=0.
 l_mixmin(:,:,:)=0.
 l_mix(:,:,:) = 0.
 wprime(:,:,:) = 0 
#ifdef ISO
zxxtevap(:,:)=0.
 d_xt(:,:,:)=0. 
 d_xt_x(:,:,:)=0.
 d_xt_w(:,:,:)=0.
 flux_xt(:,:,:,:)=0. 
! xtsnow(:,:,:)=0.! attention, xtsnow est l'équivalent de snow et non de qsnow
 xtevap(:,:,:)=0.
#endif
    IF (iflag_pbl<20.or.iflag_pbl>=30) THEN
       zcoefh(:,:,:) = 0.0
       zcoefh(:,1,:) = 999999. ! zcoefh(:,k=1) should never be used
       zcoefm(:,:,:) = 0.0
       zcoefm(:,1,:) = 999999. !
    ELSE
      zcoefm(:,:,is_ave)=0.
      zcoefh(:,:,is_ave)=0.
    ENDIF
!!
!  The components "is_ave" of tke_x and wake_deltke are "OUT" variables
!jyg<
!!    tke(:,:,is_ave)=0.
    tke_x(:,:,is_ave)=0.
    eps_x(:,:,is_ave)=0.

    wake_dltke(:,:,is_ave)=0.
!>jyg
!!! jyg le 23/02/2013
    t2m(:,:)       = 999999.     ! t2m and q2m are meaningfull only over sub-surfaces
    q2m(:,:)       = 999999.     ! actually present in the grid cell.
!!!
    rh2m(:) = 0. ; qsat2m(:) = 0.
!!!
!!! jyg le 10/02/2012
    rh2m_x(:) = 0. ; qsat2m_x(:) = 0. ; rh2m_w(:) = 0. ; qsat2m_w(:) = 0.


#ifdef ISO
   dflux_xt=0.0 
#endif 

! 2c) Initialization of all local variables computed within the subsurface loop and used later on
!****************************************************************************************
    d_t_diss_x(:,:) = 0. ;        d_t_diss_w(:,:) = 0.
    d_u_x(:,:)=0. ;               d_u_w(:,:)=0. 
    d_v_x(:,:)=0. ;               d_v_w(:,:)=0.
    flux_t_x(:,:,:)=0. ;          flux_t_w(:,:,:)=0. 
    flux_q_x(:,:,:)=0. ;          flux_q_w(:,:,:)=0.
!
!jyg<
    flux_u_x(:,:,:)=0. ;          flux_u_w(:,:,:)=0.
    flux_v_x(:,:,:)=0. ;          flux_v_w(:,:,:)=0.
    fluxlat_x(:,:)=0. ;           fluxlat_w(:,:)=0. 
!>jyg
#ifdef ISO
    flux_xt_x(:,:,:,:)=0. ;          flux_xt_w(:,:,:,:)=0.
#endif
!
!jyg<
! pblh,plcl,capCL,cteiCL ... are meaningfull only over sub-surfaces
! actually present in the grid cell  ==> value set to 999999.
!                           
!jyg<
       ustar(:,:)   = 999999.
       wstar(:,:)   = 999999.
       windsp(:,:)  = SQRT(u10m(:,:)**2 + v10m(:,:)**2 )
       u10m(:,:)    = 999999. 
       v10m(:,:)    = 999999. 
!>jyg
!
       pblh(:,:)   = 999999.        ! Hauteur de couche limite
       plcl(:,:)   = 999999.        ! Niveau de condensation de la CLA
       capCL(:,:)  = 999999.        ! CAPE de couche limite
       oliqCL(:,:) = 999999.        ! eau_liqu integree de couche limite
       cteiCL(:,:) = 999999.        ! cloud top instab. crit. couche limite
       pblt(:,:)   = 999999.        ! T a la Hauteur de couche limite
       therm(:,:)  = 999999.
       trmb1(:,:)  = 999999.        ! deep_cape
       trmb2(:,:)  = 999999.        ! inhibition 
       trmb3(:,:)  = 999999.        ! Point Omega
!
       t2m_x(:,:)    = 999999. 
       q2m_x(:,:)    = 999999. 
       ustar_x(:,:)   = 999999.
       wstar_x(:,:)   = 999999.
       u10m_x(:,:)   = 999999. 
       v10m_x(:,:)   = 999999. 
!                           
       pblh_x(:,:)   = 999999.      ! Hauteur de couche limite
       plcl_x(:,:)   = 999999.      ! Niveau de condensation de la CLA
       capCL_x(:,:)  = 999999.      ! CAPE de couche limite
       oliqCL_x(:,:) = 999999.      ! eau_liqu integree de couche limite
       cteiCL_x(:,:) = 999999.      ! cloud top instab. crit. couche limite
       pblt_x(:,:)   = 999999.      ! T a la Hauteur de couche limite
       therm_x(:,:)  = 999999.      
       trmb1_x(:,:)  = 999999.      ! deep_cape
       trmb2_x(:,:)  = 999999.      ! inhibition 
       trmb3_x(:,:)  = 999999.      ! Point Omega
!
       t2m_w(:,:)    = 999999. 
       q2m_w(:,:)    = 999999. 
       ustar_w(:,:)   = 999999.
       wstar_w(:,:)   = 999999.
       u10m_w(:,:)   = 999999. 
       v10m_w(:,:)   = 999999. 
                           
       pblh_w(:,:)   = 999999.      ! Hauteur de couche limite
       plcl_w(:,:)   = 999999.      ! Niveau de condensation de la CLA
       capCL_w(:,:)  = 999999.      ! CAPE de couche limite
       oliqCL_w(:,:) = 999999.      ! eau_liqu integree de couche limite
       cteiCL_w(:,:) = 999999.      ! cloud top instab. crit. couche limite
       pblt_w(:,:)   = 999999.      ! T a la Hauteur de couche limite
       therm_w(:,:)  = 999999.      
       trmb1_w(:,:)  = 999999.      ! deep_cape
       trmb2_w(:,:)  = 999999.      ! inhibition 
       trmb3_w(:,:)  = 999999.      ! Point Omega
!!!      
!
!!!
!****************************************************************************************
! 3) - Calculate pressure thickness of each layer
!    - Calculate the wind at first layer
!    - Mean calculations of albedo
!    - Calculate net radiance at sub-surface
!****************************************************************************************
    DO k = 1, klev
       DO i = 1, klon
          delp(i,k) = paprs(i,k)-paprs(i,k+1)
       ENDDO
    ENDDO

!****************************************************************************************
! Test for rugos........ from physiq.. A la fin plutot???
!
!****************************************************************************************

    DO nsrf = 1, nbsrf
       DO i = 1, klon
          z0m(i,nsrf) = MAX(z0m(i,nsrf),z0min)
          z0h(i,nsrf) = MAX(z0h(i,nsrf),z0min)
       ENDDO
    ENDDO

    ! AM heterogeneous continental subsurfaces
    ! compute time-independent effective surface parameters
    IF (iflag_hetero_surf .GT. 0) THEN
      CALL eff_surf_param(klon, nbtersrf, albedo_tersrf, frac_tersrf, 'ARI', albedo_eff)
    ENDIF

! Mean calculations of albedo
!
! * alb  : mean albedo for whole SW interval
!
! Mean albedo for grid point
! * alb_m  : mean albedo at whole SW interval

    alb_dir_m(:,:) = 0.0
    alb_dif_m(:,:) = 0.0
    DO k = 1, nsw
     DO nsrf = 1, nbsrf
       DO i = 1, klon
          ! AM heterogeneous continental sub-surfaces
          IF (nsrf .EQ. is_ter .AND. iflag_hetero_surf .GT. 0) THEN
            alb_dir(i,k,nsrf) = albedo_eff(i)
            alb_dif(i,k,nsrf) = albedo_eff(i)
          ENDIF
          !
          alb_dir_m(i,k) = alb_dir_m(i,k) + alb_dir(i,k,nsrf) * pctsrf(i,nsrf)
          alb_dif_m(i,k) = alb_dif_m(i,k) + alb_dif(i,k,nsrf) * pctsrf(i,nsrf)
       ENDDO
     ENDDO
    ENDDO

! We here suppose the fraction f1 of incoming radiance of visible radiance 
! as a fraction of all shortwave radiance 
    f1 = 0.5
!    f1 = 1    ! put f1=1 to recreate old calculations

!f1 is already included with SFRWL values in each surf files
    alb=0.0
    DO k=1,nsw
      DO nsrf = 1, nbsrf
        DO i = 1, klon
            alb(i,nsrf) = alb(i,nsrf) + alb_dir(i,k,nsrf)*SFRWL(k)
        ENDDO
      ENDDO
    ENDDO

    alb_m=0.0
    DO k = 1,nsw
      DO i = 1, klon
        alb_m(i) = alb_m(i) + alb_dir_m(i,k)*SFRWL(k)
      ENDDO
    ENDDO
!albedo SB <<<



! Calculation of mean temperature at surface grid points
    ztsol(:) = 0.0
    DO nsrf = 1, nbsrf
       DO i = 1, klon
          ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)
       ENDDO
    ENDDO

! Linear distrubution on sub-surface of long- and shortwave net radiance
    DO nsrf = 1, nbsrf
       DO i = 1, klon
          sollw(i,nsrf) = sollw_m(i) + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ts(i,nsrf))
!--OB this line is not satisfactory because alb is the direct albedo not total albedo
          solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
       ENDDO
    ENDDO
!
!<al1: second order corrections
!- net = dwn -up; up=sig( T4 + 4sum%T3T' + 6sum%T2T'2 +...)
   IF (iflag_order2_sollw == 1) THEN
    meansqT(:) = 0. ! as working buffer
    DO nsrf = 1, nbsrf
     DO i = 1, klon
      meansqT(i) = meansqT(i)+(ts(i,nsrf)-ztsol(i))**2 *pctsrf(i,nsrf)
     ENDDO
    ENDDO
    DO nsrf = 1, nbsrf
     DO i = 1, klon
      sollw(i,nsrf) = sollw(i,nsrf) &
                + 6.0*RSIGMA*ztsol(i)**2 *(meansqT(i)-(ztsol(i)-ts(i,nsrf))**2)
     ENDDO
    ENDDO
   ENDIF   ! iflag_order2_sollw == 1
!>al1

!--OB add diffuse fraction of SW down
   DO n=1,nbcf_out
       IF (cfname_out(n) == "swdownfdiff" ) fields_out(:,n) = solswfdiff_m(:)
   ENDDO
! >> PC
   IF (carbon_cycle_cpl .AND. carbon_cycle_tr .AND. nbcf_out.GT.0 ) THEN
       r_co2_ppm(1:klon) = co2_send(1:klon)
       DO n=1,nbcf_out
           IF (cfname_out(n) == "atmco2" ) fields_out(1:klon,n) = co2_send(1:klon)
       ENDDO
   ENDIF

   IF ( .NOT. carbon_cycle_tr .AND. nbcf_out.GT.0 ) THEN
       r_co2_ppm(1:klon) = co2_ppm     ! Constant field
       DO n=1,nbcf_out
           IF (cfname_out(n) == "atmco2" ) fields_out(1:klon,n) = co2_ppm
       ENDDO
   ENDIF
   
END SUBROUTINE pbl_surface_uncompress_pre

END MODULE pbl_surface_uncompress_pre_mod
