!
! $Id: pbl_surface_mod.F90 5990 2025-12-19 18:14:40Z ymeurdesoif $
!
MODULE pbl_surface_uncompressed_post_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_uncompressed_post

CONTAINS

  SUBROUTINE pbl_surface_uncompressed_post( &
       itap, dtime,         &
       u,        v,        &
       wake_s,                  &
       pctsrf,                  &
       ts,ustar, u10m, v10m,wstar, &
       zu1,    zv1,              &
       zxsens,   zxevap,  zxsnowerosion,      &
       zxtsol,    zxfluxlat, zt2m,     qsat2m, zn2mout,                 &
       zxsens_x,  zxfluxlat_x,zxsens_w,zxfluxlat_w,  &
       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,                          &
       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, bilg_cumul, iflag_split_ref,  &
       & n2mout, n2mout_x, 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, t2m_x, q2m_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, qsat2m_w,  & 
       pblh_w, plcl_w, pblh, plcl, capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3 &        
#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 hbtm_mod, only: hbtm
    USE indice_sol_mod
    USE mod_grid_phy_lmdz,  ONLY : grid1dto2d_glo
#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 dimsoil_mod_h, ONLY: nsoilmx
    USE surf_param_mod, ONLY: eff_surf_param  !AM
    USE yomcst_mod_h
    USE ocean_forced_mod,ONLY : ocean_forced_ice_reset_bilg_cumul
IMPLICIT NONE

    INCLUDE "FCTTRE.h"

!****************************************************************************************
    INTEGER,                      INTENT(IN)        :: itap    ! time step
    REAL,                         INTENT(IN)        :: dtime   ! interval du temps (secondes)
    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: u       ! u speed
    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: v       ! v speed
    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
    REAL, DIMENSION(klon),        INTENT(IN)        :: wake_s    ! Fraction de poches froides

#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(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

! Output variables
!****************************************************************************************
    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),        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)       :: zxtsol     ! temperature at surface, mean for each grid point
    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
#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)       :: 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, 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(IN)       :: solsw      ! net shortwave radiation at surface 
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: sollw      ! net longwave radiation at surface
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: d_ts       ! change in temperature at surface
    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: evap       ! evaporation at surface
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: fluxlat    ! latent flux
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: 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(IN) :: flux_t     ! sensible heat flux (CpT) J/m**2/s (W/m**2)
                                                                  ! positve orientation downwards
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(IN) :: flux_u     ! u wind tension (kg m/s)/(m**2 s) or Pascal
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(IN) :: 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(IN)        :: dflux_t    ! change of sensible heat flux 
    REAL, DIMENSION(klon),       INTENT(IN)        :: 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(IN)        :: q2m        ! water vapour at 2 meter height
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(IN) :: flux_q     ! water vapour flux(latent flux) (kg/m**2/s)
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(IN) :: flux_qbs   ! blowind snow vertical flux (kg/m**2
    REAL, DIMENSION(klon),       INTENT(INOUT)     :: bilg_cumul      ! flux cumulated
    INTEGER,                     INTENT(INOUT)      :: iflag_split_ref 

#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

! Other local variables
!****************************************************************************************
    INTEGER                            :: iflag_split
    INTEGER                            :: i, k, nsrf 
    REAL                               :: amn, amx
    INTEGER, DIMENSION(klon, nbsrf, 6), INTENT(IN) :: n2mout, n2mout_x
    REAL, DIMENSION(klon)              :: meansqT ! mean square deviation of subsurface temperatures
    LOGICAL, PARAMETER                 :: zxli=.FALSE. ! utiliser un jeu de fonctions simples
    LOGICAL, PARAMETER                 :: check=.FALSE.

    REAL, DIMENSION(klon, klev, nbsrf), INTENT(IN) :: flux_t_x, flux_q_x, flux_t_w, flux_q_w
    REAL, DIMENSION(klon, klev, nbsrf), INTENT(IN) :: flux_u_x, flux_v_x, flux_u_w, flux_v_w
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: fluxlat_x, fluxlat_w
    REAL, DIMENSION(klon, klev)        :: zxfluxt_x, zxfluxq_x, zxfluxt_w, zxfluxq_w
    REAL, DIMENSION(klon, klev)        :: zxfluxu_x, zxfluxv_x, zxfluxu_w, zxfluxv_w
#ifdef ISO
    REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN)         :: zxfluxxt_x
    REAL, DIMENSION(ntraciso,klon,klev), INTENT(IN)         :: zxfluxxt_w
    REAL, DIMENSION(ntraciso,klon,klev,nbsrf), INTENT(IN)   :: flux_xt_x, flux_xt_w
#endif
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: t2m_x
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: q2m_x
    REAL, DIMENSION(klon), INTENT(IN)              :: qsat2m_x
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: u10m_x
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: v10m_x
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: ustar_x
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: wstar_x
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: pblh_x
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: plcl_x
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: capCL_x
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: oliqCL_x
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: cteiCL_x
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: pblt_x
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: therm_x
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: trmb1_x
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: trmb2_x
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: trmb3_x
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: t2m_w
    REAL, DIMENSION(klon), INTENT(IN)              :: qsat2m_w
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: pblh_w
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)       :: plcl_w
    REAL, PARAMETER                    :: facteur = 2. / 1.772  ! ( == 2. / SQRT(3.14))
    REAL, PARAMETER                    :: inertia=2000.

    REAL, DIMENSION(klon,nbsrf), INTENT(IN)        :: pblh         ! height of the planetary boundary layer
    REAL, DIMENSION(klon,nbsrf), INTENT(IN)        :: plcl         ! condensation level
    REAL, DIMENSION(klon,nbsrf), INTENT(IN)        :: capCL
    REAL, DIMENSION(klon,nbsrf), INTENT(IN)        :: oliqCL
    REAL, DIMENSION(klon,nbsrf), INTENT(IN)        :: cteiCL
    REAL, DIMENSION(klon,nbsrf), INTENT(IN)        :: pblT
    REAL, DIMENSION(klon,nbsrf), INTENT(IN)        :: therm
    REAL, DIMENSION(klon,nbsrf), INTENT(IN)        :: trmb1        ! deep cape
    REAL, DIMENSION(klon,nbsrf), INTENT(IN)        :: trmb2        ! inhibition
    REAL, DIMENSION(klon,nbsrf), INTENT(IN)        :: trmb3        ! point Omega
#ifdef ISO
    INTEGER                     :: ixt
#endif

!
!----------------------------------------------------------------------------------------
!   Reset iflag_split 
!
   iflag_split=iflag_split_ref

#ifdef ISO
#ifdef ISOVERIF

    IF (iso_eau >= 0) THEN
        call iso_verif_egalite_vect2D( &
                d_xt,d_q, &
                'pbl_surface_mod 1276',ntraciso,klon,klev)
    ENDIF !IF (iso_eau >= 0) THEN
#endif
#endif

!YM something bad to check 
    CALL ocean_forced_ice_reset_bilg_cumul(itap, dtime, bilg_cumul)
!****************************************************************************************
! 16) Calculate the mean value over all sub-surfaces for some variables
!
!****************************************************************************************
    
    z0m(:,nbsrf+1) = 0.0
    z0h(:,nbsrf+1) = 0.0
    DO nsrf = 1, nbsrf
       DO i = 1, klon
          z0m(i,nbsrf+1) = z0m(i,nbsrf+1) + z0m(i,nsrf)*pctsrf(i,nsrf)
          z0h(i,nbsrf+1) = z0h(i,nbsrf+1) + z0h(i,nsrf)*pctsrf(i,nsrf)
       ENDDO
    ENDDO

    zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
    zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
    zxfluxt_x(:,:) = 0.0 ; zxfluxq_x(:,:) = 0.0
    zxfluxu_x(:,:) = 0.0 ; zxfluxv_x(:,:) = 0.0
    zxfluxt_w(:,:) = 0.0 ; zxfluxq_w(:,:) = 0.0
    zxfluxu_w(:,:) = 0.0 ; zxfluxv_w(:,:) = 0.0
#ifdef ISO
      zxfluxxt(:,:,:) = 0.0 
      zxfluxxt_x(:,:,:) = 0.0
      zxfluxxt_w(:,:,:) = 0.0
#endif


       IF (iflag_split .ge.1) THEN

        DO nsrf = 1, nbsrf
          DO k = 1, klev
            DO i = 1, klon
              zxfluxt_x(i,k) = zxfluxt_x(i,k) + flux_t_x(i,k,nsrf) * pctsrf(i,nsrf)
              zxfluxq_x(i,k) = zxfluxq_x(i,k) + flux_q_x(i,k,nsrf) * pctsrf(i,nsrf)
              zxfluxu_x(i,k) = zxfluxu_x(i,k) + flux_u_x(i,k,nsrf) * pctsrf(i,nsrf)
              zxfluxv_x(i,k) = zxfluxv_x(i,k) + flux_v_x(i,k,nsrf) * pctsrf(i,nsrf)
              zxfluxt_w(i,k) = zxfluxt_w(i,k) + flux_t_w(i,k,nsrf) * pctsrf(i,nsrf)
              zxfluxq_w(i,k) = zxfluxq_w(i,k) + flux_q_w(i,k,nsrf) * pctsrf(i,nsrf)
              zxfluxu_w(i,k) = zxfluxu_w(i,k) + flux_u_w(i,k,nsrf) * pctsrf(i,nsrf)
              zxfluxv_w(i,k) = zxfluxv_w(i,k) + flux_v_w(i,k,nsrf) * pctsrf(i,nsrf)
#ifdef ISO
              DO ixt=1,ntraciso
                zxfluxxt_x(ixt,i,k) = zxfluxxt_x(ixt,i,k) + flux_xt_x(ixt,i,k,nsrf) * pctsrf(i,nsrf)
                zxfluxxt_w(ixt,i,k) = zxfluxxt_w(ixt,i,k) + flux_xt_w(ixt,i,k,nsrf) * pctsrf(i,nsrf)
              ENDDO ! DO ixt=1,ntraciso
#endif
            ENDDO
          ENDDO
        ENDDO

    DO i = 1, klon
      zxsens_x(i) = - zxfluxt_x(i,1)
      zxsens_w(i) = - zxfluxt_w(i,1)
    ENDDO
!!!
       ENDIF  ! (iflag_split .ge.1)
!!!

    DO nsrf = 1, nbsrf
       DO k = 1, klev
          DO i = 1, klon
             zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
             zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
             zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
             zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
#ifdef ISO
             DO ixt=1,niso
               zxfluxxt(ixt,i,k) = zxfluxxt(ixt,i,k) + flux_xt(ixt,i,k,nsrf) * pctsrf(i,nsrf)
             ENDDO ! DO ixt=1,niso
#endif
          ENDDO
       ENDDO
    ENDDO

    DO i = 1, klon
       zxsens(i)     = - zxfluxt(i,1) ! flux de chaleur sensible au sol
       zxevap(i)     = - zxfluxq(i,1) ! flux d'evaporation au sol
       fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
    ENDDO

    ! if blowing snow
    if (ok_bs) then  
       DO nsrf = 1, nbsrf 
       DO k = 1, klev
       DO i = 1, klon 
         zxfluxqbs(i,k) = zxfluxqbs(i,k) + flux_qbs(i,k,nsrf) * pctsrf(i,nsrf)
       ENDDO
       ENDDO
       ENDDO

       DO i = 1, klon
        zxsnowerosion(i)     = zxfluxqbs(i,1) ! blowings snow flux at the surface 
       END DO
    endif

#ifdef ISO
    DO i = 1, klon
      DO ixt=1,ntraciso
        zxxtevap(ixt,i)     = - zxfluxxt(ixt,i,1)
      ENDDO
    ENDDO
#endif

!
! Incrementer la temperature du sol
!
    zxtsol(:) = 0.0  ; zxfluxlat(:) = 0.0
    zt2m(:) = 0.0    ; zq2m(:) = 0.0 ; zn2mout(:,:) = 0
    zustar(:)=0.0 ; zu10m(:) = 0.0   ; zv10m(:) = 0.0
    s_pblh(:) = 0.0  ; s_plcl(:) = 0.0 

     s_pblh_x(:) = 0.0  ; s_plcl_x(:) = 0.0 
     s_pblh_w(:) = 0.0  ; s_plcl_w(:) = 0.0 

    s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
    s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
    s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
    s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
    wstar(:,is_ave)=0.
    
    zxfluxlat_x(:) = 0.0  ;  zxfluxlat_w(:) = 0.0
    
    DO nsrf = 1, nbsrf
       DO i = 1, klon          
          ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
          
          wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
               + flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)

          wfevap(i,nsrf) = evap(i,nsrf)*pctsrf(i,nsrf)

          zxtsol(i)    = zxtsol(i)    + ts(i,nsrf)      * pctsrf(i,nsrf)
          zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
       ENDDO
    ENDDO
!
!<al1 order 2 correction to zxtsol, for radiation computations (main atm effect of Ts)
   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)-zxtsol(i))**2 *pctsrf(i,nsrf)
     ENDDO
    ENDDO
    zxtsol(:) = zxtsol(:)+1.5*meansqT(:)/zxtsol(:)
   ENDIF   ! iflag_order2_sollw == 1

       IF (iflag_split .eq.0) THEN
        DO nsrf = 1, nbsrf
         DO i = 1, klon          
          zt2m(i)  = zt2m(i)  + t2m(i,nsrf)  * pctsrf(i,nsrf)
          zq2m(i)  = zq2m(i)  + q2m(i,nsrf)  * pctsrf(i,nsrf)
!
          DO k = 1, 6
           zn2mout(i,k)  = zn2mout(i,k)  + n2mout(i,nsrf,k)  * pctsrf(i,nsrf)
          ENDDO  
!
          zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)
          wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf)
          zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
          zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)

          s_pblh(i)   = s_pblh(i)   + pblh(i,nsrf)  * pctsrf(i,nsrf)
          s_plcl(i)   = s_plcl(i)   + plcl(i,nsrf)  * pctsrf(i,nsrf)
          s_capCL(i)  = s_capCL(i)  + capCL(i,nsrf) * pctsrf(i,nsrf)
          s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
          s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
          s_pblT(i)   = s_pblT(i)   + pblT(i,nsrf)  * pctsrf(i,nsrf)
          s_therm(i)  = s_therm(i)  + therm(i,nsrf) * pctsrf(i,nsrf)
          s_trmb1(i)  = s_trmb1(i)  + trmb1(i,nsrf) * pctsrf(i,nsrf)
          s_trmb2(i)  = s_trmb2(i)  + trmb2(i,nsrf) * pctsrf(i,nsrf)
          s_trmb3(i)  = s_trmb3(i)  + trmb3(i,nsrf) * pctsrf(i,nsrf)
         ENDDO
        ENDDO
       ELSE  !(iflag_split .eq.0)
        DO nsrf = 1, nbsrf
         DO i = 1, klon          
          zxfluxlat_x(i) = zxfluxlat_x(i) + fluxlat_x(i,nsrf) * pctsrf(i,nsrf)
          zxfluxlat_w(i) = zxfluxlat_w(i) + fluxlat_w(i,nsrf) * pctsrf(i,nsrf)
!!!
!!! jyg le 08/02/2012
!!  Pour le moment, on sort les valeurs dans (x) et (w) de pblh et de plcl ; 
!!  pour zt2m, on fait la moyenne surfacique sur les sous-surfaces ;
!!  pour qsat2m, on fait la moyenne surfacique sur (x) et (w) ;
!!  pour les autres variables, on sort les valeurs de la region (x).
          zt2m(i)  = zt2m(i)  + (t2m_x(i,nsrf)+wake_s(i)*(t2m_w(i,nsrf)-t2m_x(i,nsrf))) * pctsrf(i,nsrf)
          zq2m(i)  = zq2m(i)  + q2m_x(i,nsrf)  * pctsrf(i,nsrf)
!
          DO k = 1, 6
           zn2mout(i,k)  = zn2mout(i,k)  + n2mout_x(i,nsrf,k)  * pctsrf(i,nsrf)
          ENDDO
!
          zustar(i) = zustar(i) + ustar_x(i,nsrf) * pctsrf(i,nsrf)
          wstar(i,is_ave)=wstar(i,is_ave)+wstar_x(i,nsrf)*pctsrf(i,nsrf)
          zu10m(i) = zu10m(i) + u10m_x(i,nsrf) * pctsrf(i,nsrf)
          zv10m(i) = zv10m(i) + v10m_x(i,nsrf) * pctsrf(i,nsrf)
!
          s_pblh(i)     = s_pblh(i)     + pblh_x(i,nsrf)  * pctsrf(i,nsrf)
          s_pblh_x(i)   = s_pblh_x(i)   + pblh_x(i,nsrf)  * pctsrf(i,nsrf)
          s_pblh_w(i)   = s_pblh_w(i)   + pblh_w(i,nsrf)  * pctsrf(i,nsrf)
!
          s_plcl(i)     = s_plcl(i)     + plcl_x(i,nsrf)  * pctsrf(i,nsrf)
          s_plcl_x(i)   = s_plcl_x(i)   + plcl_x(i,nsrf)  * pctsrf(i,nsrf)
          s_plcl_w(i)   = s_plcl_w(i)   + plcl_w(i,nsrf)  * pctsrf(i,nsrf)
!
          s_capCL(i)  = s_capCL(i)  + capCL_x(i,nsrf) * pctsrf(i,nsrf)
          s_oliqCL(i) = s_oliqCL(i) + oliqCL_x(i,nsrf)* pctsrf(i,nsrf)
          s_cteiCL(i) = s_cteiCL(i) + cteiCL_x(i,nsrf)* pctsrf(i,nsrf)
          s_pblT(i)   = s_pblT(i)   + pblT_x(i,nsrf)  * pctsrf(i,nsrf)
          s_therm(i)  = s_therm(i)  + therm_x(i,nsrf) * pctsrf(i,nsrf)
          s_trmb1(i)  = s_trmb1(i)  + trmb1_x(i,nsrf) * pctsrf(i,nsrf)
          s_trmb2(i)  = s_trmb2(i)  + trmb2_x(i,nsrf) * pctsrf(i,nsrf)
          s_trmb3(i)  = s_trmb3(i)  + trmb3_x(i,nsrf) * pctsrf(i,nsrf)
         ENDDO
        ENDDO
        DO i = 1, klon          
          qsat2m(i)= qsat2m_x(i)+ wake_s(i)*(qsat2m_x(i)-qsat2m_w(i))
        ENDDO
!!!
       ENDIF  ! (iflag_split .eq.0)
!!!

    IF (check) THEN
       amn=MIN(ts(1,is_ter),1000.)
       amx=MAX(ts(1,is_ter),-1000.)
       DO i=2, klon
          amn=MIN(ts(i,is_ter),amn)
          amx=MAX(ts(i,is_ter),amx)
       ENDDO
       PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
    ENDIF

    DO i = 1, klon
       fder(i) = - 4.0*RSIGMA*zxtsol(i)**3 
    ENDDO
    
    zxqsurf(:) = 0.0
    zxsnow(:)  = 0.0
#ifdef ISO
    zxxtsnow(:,:)  = 0.0
#endif

    DO nsrf = 1, nbsrf
       DO i = 1, klon
          zxqsurf(i) = zxqsurf(i) + MAX(qsurf(i,nsrf),0.0) * pctsrf(i,nsrf)
          zxsnow(i)  = zxsnow(i)  + snow(i,nsrf)  * pctsrf(i,nsrf)
#ifdef ISO
          DO ixt=1,niso
            zxxtsnow(ixt,i)  = zxxtsnow(ixt,i)  + xtsnow(ixt,i,nsrf)  * pctsrf(i,nsrf)
          ENDDO ! DO ixt=1,niso
#endif
       ENDDO
    ENDDO

! Premier niveau de vent sortie dans physiq.F
    zu1(:) = u(:,1)
    zv1(:) = v(:,1)
    
   END SUBROUTINE pbl_surface_uncompressed_post
!  
!****************************************************************************************
!  
END MODULE pbl_surface_uncompressed_post_mod
