!
! $Id: pbl_surface_mod.F90 5990 2025-12-19 18:14:40Z ymeurdesoif $
!
MODULE pbl_surface_main_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_init, pbl_surface_final, pbl_surface_newfrac, pbl_surface_main
#ifdef ISO
  PUBLIC :: pbl_surface_init_iso
#endif

CONTAINS

!
!****************************************************************************************
!
!GG
!  SUBROUTINE pbl_surface_init(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
  SUBROUTINE pbl_surface_init(fder_rst, snow_rst, qsurf_rst, ftsoil_rst, hice_rst,tice_rst,bilg_cumul_rst)
!GG

! This routine should be called after the restart file has been read.
! This routine initialize the restart variables and does some validation tests
! for the index of the different surfaces and tests the choice of type of ocean.
    USE pbl_surface_data
    USE indice_sol_mod
    USE print_control_mod, ONLY: lunout
    USE ioipsl_getin_p_mod, ONLY : getin_p
    USE dimsoil_mod_h, ONLY: nsoilmx
    USE flux_arp_mod_h
    USE cdrag_mod, ONLY : cdrag_init
    USE climb_hq_mod, ONLY : climb_hq_init
    USE climb_wind_mod, ONLY : climb_wind_init
    USE climb_qbs_mod, ONLY : climb_qbs_init
    USE yamada_c_mod, ONLY : yamada_c_init
    USE soil_mod, ONLY : soil_init
    USE surf_landice_mod, ONLY : surf_landice_init

    IMPLICIT NONE
 
! Input variables
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(IN)                 :: fder_rst
!GG
    REAL, DIMENSION(klon), INTENT(IN)                 :: hice_rst
    REAL, DIMENSION(klon), INTENT(IN)                 :: tice_rst
    REAL, DIMENSION(klon), INTENT(IN)                 :: bilg_cumul_rst
!GG
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: snow_rst
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: qsurf_rst
    REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(IN) :: ftsoil_rst
  
! Local variables
!****************************************************************************************
    INTEGER                       :: ierr
    CHARACTER(len=80)             :: abort_message
    CHARACTER(len = 20)           :: modname = 'pbl_surface_init'

!****************************************************************************************
! Initialize some module variables
!****************************************************************************************    
    smallestreal = tiny(smallestreal)
    
!****************************************************************************************
! Allocate and initialize module variables with fields read from restart file.
!
!****************************************************************************************    

    ALLOCATE(fder(klon), stat=ierr)
    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)

!GG
    ALLOCATE(hice(klon), stat=ierr)
    IF (ierr /= 0) CALL abort_physic('pbl_surface_init hice', 'pb in allocation',1)

    ALLOCATE(tice(klon), stat=ierr)
    IF (ierr /= 0) CALL abort_physic('pbl_surface_init tice', 'pb in allocation',1)

    ALLOCATE(bilg_cumul(klon), stat=ierr)
    IF (ierr /= 0) CALL abort_physic('pbl_surface_init bilg', 'pb in allocation',1)
!GG

    ALLOCATE(snow(klon,nbsrf), stat=ierr)
    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)

    ALLOCATE(qsurf(klon,nbsrf), stat=ierr)
    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)

    ALLOCATE(ftsoil(klon,nsoilmx,nbsrf), stat=ierr)
    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)

    ALLOCATE(ydTs0(klon), stat=ierr)
    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)

    ALLOCATE(ydqs0(klon), stat=ierr)
    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)

    fder(:)       = fder_rst(:)
!GG
    hice(:)       = hice_rst(:)
    tice(:)       = tice_rst(:)
    bilg_cumul(:)       = bilg_cumul_rst(:)
!GG
    snow(:,:)     = snow_rst(:,:)
    qsurf(:,:)    = qsurf_rst(:,:)
    ftsoil(:,:,:) = ftsoil_rst(:,:,:)
    ydTs0(:) = 0.
    ydqs0(:) = 0.

!****************************************************************************************
! Test for sub-surface indices
!
!****************************************************************************************
    IF (is_ter /= 1) THEN 
      WRITE(lunout,*)" *** Warning ***"
      WRITE(lunout,*)" is_ter n'est pas le premier surface, is_ter = ",is_ter
      WRITE(lunout,*)"or on doit commencer par les surfaces continentales"
      abort_message="voir ci-dessus"
      CALL abort_physic(modname,abort_message,1)
    ENDIF

    IF ( is_oce > is_sic ) THEN
      WRITE(lunout,*)' *** Warning ***'
      WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
      WRITE(lunout,*)' l''ocean doit etre traite avant la banquise'
      WRITE(lunout,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic
      abort_message='voir ci-dessus'
      CALL abort_physic(modname,abort_message,1)
    ENDIF

    IF ( is_lic > is_sic ) THEN
      WRITE(lunout,*)' *** Warning ***'
      WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
      WRITE(lunout,*)' la glace contineltalle doit etre traite avant la glace de mer'
      WRITE(lunout,*)' or is_lic = ',is_lic, '> is_sic = ',is_sic
      abort_message='voir ci-dessus'
      CALL abort_physic(modname,abort_message,1)
    ENDIF

!****************************************************************************************
! Validation of ocean mode
!
!****************************************************************************************

    IF (type_ocean /= 'slab  ' .AND. type_ocean /= 'force ' .AND. type_ocean /= 'couple') THEN
       WRITE(lunout,*)' *** Warning ***'
       WRITE(lunout,*)'Option couplage pour l''ocean = ', type_ocean
       abort_message='option pour l''ocean non valable'
       CALL abort_physic(modname,abort_message,1)
    ENDIF

    iflag_pbl_surface_t2m_bug=0
    CALL getin_p('iflag_pbl_surface_t2m_bug',iflag_pbl_surface_t2m_bug)
    WRITE(lunout,*) 'iflag_pbl_surface_t2m_bug=',iflag_pbl_surface_t2m_bug
!FC
!    iflag_frein = 0
!    CALL getin_p('iflag_frein',iflag_frein)
!
!jyg<
!****************************************************************************************
! Allocate variables for pbl splitting
!
!****************************************************************************************

!****************************************************************************************
!   Initialisation and validation tests 
!   moved from only done first time entering this subroutine
!
!****************************************************************************************
    iflag_new_t2mq2m=1
    CALL getin_p('iflag_new_t2mq2m',iflag_new_t2mq2m)
    WRITE(lunout,*) 'pbl_iflag_new_t2mq2m=',iflag_new_t2mq2m

    ok_bug_zg_wk_pbl=.TRUE.
    CALL getin_p('ok_bug_zg_wk_pbl',ok_bug_zg_wk_pbl)
    WRITE(lunout,*) 'ok_bug_zg_wk_pbl=',ok_bug_zg_wk_pbl

    print*,'PBL SURFACE AVEC GUSTINESS'
      
    ! Initialize ok_flux_surf (for 1D model)
    IF (klon_glo>1) ok_flux_surf=.FALSE.
    IF (klon_glo>1) ok_forc_tsurf=.FALSE.

    ! intialize beta_land
    beta_land = 0.5
    call getin_p('beta_land', beta_land)


    CALL wx_pbl_init
!>jyg

    CALL cdrag_init
    CALL climb_hq_init
    CALL climb_wind_init
    CALL climb_qbs_init
    CALL yamada_c_init
    CALL soil_init
    CALL surf_landice_init

  END SUBROUTINE pbl_surface_init

#ifdef ISO
  SUBROUTINE pbl_surface_init_iso(xtsnow_rst,Rland_ice_rst)

! This routine should be called after the restart file has been read.
! This routine initialize the restart variables and does some validation tests
! for the index of the different surfaces and tests the choice of type of ocean.
    USE pbl_surface_data
    USE indice_sol_mod
    USE print_control_mod, ONLY: lunout
#ifdef ISOVERIF
    USE isotopes_mod, ONLY: iso_eau,ridicule
    USE isotopes_verif_mod
#endif
    USE dimsoil_mod_h, ONLY: nsoilmx
    IMPLICIT NONE
 
! Input variables
!****************************************************************************************
    REAL, DIMENSION(niso,klon, nbsrf), INTENT(IN)          :: xtsnow_rst
    REAL, DIMENSION(niso,klon), INTENT(IN)          :: Rland_ice_rst
  
! Local variables
!****************************************************************************************
    INTEGER                       :: ierr
    CHARACTER(len=80)             :: abort_message
    CHARACTER(len = 20)           :: modname = 'pbl_surface_init'
    integer i,ixt
    
!****************************************************************************************
! Allocate and initialize module variables with fields read from restart file.
!
!****************************************************************************************    

    ALLOCATE(xtsnow(niso,klon,nbsrf), stat=ierr)
    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)

    ALLOCATE(Rland_ice(niso,klon), stat=ierr)
    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)

    ALLOCATE(Roce(niso,klon), stat=ierr)
    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)

    xtsnow(:,:,:)  = xtsnow_rst(:,:,:)
    Rland_ice(:,:) = Rland_ice_rst(:,:)
    Roce(:,:)      = 0.0

#ifdef ISOVERIF 
      IF (iso_eau >= 0) THEN
         CALL iso_verif_egalite_vect2D( &
     &           xtsnow,snow, &
     &           'pbl_surface_mod 170',niso,klon,nbsrf)
         DO i=1,klon  
            IF (iso_eau >= 0) THEN  
              CALL iso_verif_egalite(Rland_ice(iso_eau,i),1.0, &
     &         'pbl_surf_mod 177')
            ENDIF
         ENDDO
      ENDIF
#endif

  END SUBROUTINE pbl_surface_init_iso
#endif
!
!****************************************************************************************
!
  SUBROUTINE pbl_surface_final(fder_rst, snow_rst, qsurf_rst, ftsoil_rst &
#ifdef ISO
       ,xtsnow_rst,Rland_ice_rst &
#endif       
       )
    USE pbl_surface_data
    USE indice_sol_mod
#ifdef ISO
#ifdef ISOVERIF
    USE isotopes_mod, ONLY: iso_eau,ridicule
    USE isotopes_verif_mod, ONLY: errmax,errmaxrel
#endif    
#endif
    USE dimsoil_mod_h, ONLY: nsoilmx

! Ouput variables
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)                 :: fder_rst
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: snow_rst
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: qsurf_rst
    REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
#ifdef ISO
    REAL, DIMENSION(niso,klon, nbsrf), INTENT(OUT)     :: xtsnow_rst
    REAL, DIMENSION(niso,klon), INTENT(OUT)            :: Rland_ice_rst
#endif

 
!****************************************************************************************
! Return module variables for writing to restart file
!
!****************************************************************************************    
    fder_rst(:)       = fder(:)
    snow_rst(:,:)     = snow(:,:)
    qsurf_rst(:,:)    = qsurf(:,:)
    ftsoil_rst(:,:,:) = ftsoil(:,:,:)
#ifdef ISO
    xtsnow_rst(:,:,:)  = xtsnow(:,:,:) 
    Rland_ice_rst(:,:) = Rland_ice(:,:)
#endif

!****************************************************************************************
! Deallocate module variables
!
!****************************************************************************************
!   DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
    IF (ALLOCATED(fder)) DEALLOCATE(fder)
    IF (ALLOCATED(hice)) DEALLOCATE(hice)
    IF (ALLOCATED(tice)) DEALLOCATE(tice)
    IF (ALLOCATED(bilg_cumul)) DEALLOCATE(bilg_cumul)
    IF (ALLOCATED(snow)) DEALLOCATE(snow)
    IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
    IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
    IF (ALLOCATED(ydTs0)) DEALLOCATE(ydTs0)
    IF (ALLOCATED(ydqs0)) DEALLOCATE(ydqs0)
#ifdef ISO
    IF (ALLOCATED(xtsnow)) DEALLOCATE(xtsnow)
    IF (ALLOCATED(Rland_ice)) DEALLOCATE(Rland_ice)
    IF (ALLOCATED(Roce)) DEALLOCATE(Roce)
#endif

!jyg<
!****************************************************************************************
! Deallocate variables for pbl splitting
!
!****************************************************************************************

    CALL wx_pbl_final
!>jyg

  END SUBROUTINE pbl_surface_final
!  
!****************************************************************************************
! 

!albedo SB >>>
  SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, &
       evap, z0m, z0h, agesno,                                  &
       tsurf,alb_dir,alb_dif, ustar, u10m, v10m, tke &
#ifdef ISO
      ,xtevap  &
#endif
&      )  
    !albedo SB <<<
    ! Give default values where new fraction has appread
    USE pbl_surface_data
    USE compbl_mod_h
    USE clesphys_mod_h
    USE indice_sol_mod
    use phys_state_var_mod, only: delta_sal, ds_ns, dt_ns, delta_sst, dter, &
         dser, dt_ds
    use config_ocean_skin_m, only: activate_ocean_skin

! Input variables
!****************************************************************************************
    INTEGER, INTENT(IN)                     :: itime
    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old

! InOutput variables
!****************************************************************************************
    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: tsurf
!albedo SB >>>
    REAL, DIMENSION(klon,nsw,nbsrf), INTENT(INOUT)       :: alb_dir, alb_dif 
    INTEGER :: k
!albedo SB <<<
    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: ustar,u10m, v10m
    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: evap, agesno
    REAL, DIMENSION(klon,nbsrf+1), INTENT(INOUT)        :: z0m,z0h
    REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT) :: tke
#ifdef ISO
    REAL, DIMENSION(ntraciso,klon,nbsrf), INTENT(INOUT)        :: xtevap
#endif

! Local variables
!****************************************************************************************
    INTEGER           :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
    CHARACTER(len=80) :: abort_message
    CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
    INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
#ifdef ISO
    INTEGER           :: ixt
#endif
!
! All at once !! 
!****************************************************************************************
    
    DO nsrf = 1, nbsrf
       ! First decide complement sub-surfaces
       SELECT CASE (nsrf)
       CASE(is_oce)
          nsrf_comp1=is_sic
          nsrf_comp2=is_ter
          nsrf_comp3=is_lic
       CASE(is_sic)
          nsrf_comp1=is_oce
          nsrf_comp2=is_ter
          nsrf_comp3=is_lic
       CASE(is_ter)
          nsrf_comp1=is_lic
          nsrf_comp2=is_oce
          nsrf_comp3=is_sic
       CASE(is_lic)
          nsrf_comp1=is_ter
          nsrf_comp2=is_oce
          nsrf_comp3=is_sic
       END SELECT

       ! Initialize all new fractions
       DO i=1, klon
          IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
             
             IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
                ! Use the complement sub-surface, keeping the continents unchanged
                qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
                evap(i,nsrf)  = evap(i,nsrf_comp1)
                z0m(i,nsrf) = z0m(i,nsrf_comp1)
                z0h(i,nsrf) = z0h(i,nsrf_comp1)
                tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
!albedo SB >>>
                DO k=1,nsw
                 alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp1)
                 alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp1)
                ENDDO
!albedo SB <<<
                ustar(i,nsrf)  = ustar(i,nsrf_comp1)
                u10m(i,nsrf)  = u10m(i,nsrf_comp1)
                v10m(i,nsrf)  = v10m(i,nsrf_comp1)
#ifdef ISO
                DO ixt=1,ntraciso
                  xtevap(ixt,i,nsrf) = xtevap(ixt,i,nsrf_comp1)       
                ENDDO       
#endif
                IF (iflag_pbl > 1) THEN
                 tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
                ENDIF
                mfois(nsrf) = mfois(nsrf) + 1
                ! F. Codron sensible default values for ocean and sea ice
                IF (nsrf.EQ.is_oce) THEN
                   tsurf(i,nsrf) = 271.35
                   ! (temperature of sea water under sea ice, so that
                   ! is also the temperature of appearing sea water)
                   DO k=1,nsw
                      alb_dir(i,k,nsrf) = 0.06 ! typical Ocean albedo
                      alb_dif(i,k,nsrf) = 0.06
                   ENDDO
                   if (activate_ocean_skin >= 1) then
                      if (activate_ocean_skin == 2 &
                           .and. type_ocean == "couple") then
                         delta_sal(i) = 0.
                         delta_sst(i) = 0.
                         dter(i) = 0.
                         dser(i) = 0.
                         dt_ds(i) = 0.
                      end if
                      
                      ds_ns(i) = 0.
                      dt_ns(i) = 0.
                   end if
                ELSE IF (nsrf.EQ.is_sic) THEN
                   tsurf(i,nsrf) = 271.35
                   ! (Temperature at base of sea ice. Surface
                   ! temperature could be higher, up to 0 Celsius
                   ! degrees. We set it to -1.8 Celsius degrees for
                   ! consistency with the ocean slab model.)
                   DO k=1,nsw
                      alb_dir(i,k,nsrf) = 0.3 ! thin ice
                      alb_dif(i,k,nsrf) = 0.3
                   ENDDO
                ENDIF
             ELSE
                ! The continents have changed. The new fraction receives the mean sum of the existent fractions
                qsurf(i,nsrf) = qsurf(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
                evap(i,nsrf)  = evap(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
                z0m(i,nsrf) = z0m(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + z0m(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
                z0h(i,nsrf) = z0h(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + z0h(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
                tsurf(i,nsrf) = tsurf(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
!albedo SB >>>
                DO k=1,nsw
                 alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
                                        alb_dir(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
                 alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
                                        alb_dif(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
                ENDDO
!albedo SB <<<
                ustar(i,nsrf)  = ustar(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + ustar(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
                u10m(i,nsrf)  = u10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
                v10m(i,nsrf)  = v10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
#ifdef ISO
                DO ixt=1,ntraciso
                  xtevap(ixt,i,nsrf) = xtevap(ixt,i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) &
                                     + xtevap(ixt,i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3) 
                ENDDO       
#endif
                IF (iflag_pbl > 1) THEN
                 tke(i,:,nsrf) = tke(i,:,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
                ENDIF
            
                ! Security abort. This option has never been tested. To test, comment the following line.
!                abort_message='The fraction of the continents have changed!'
!                CALL abort_physic(modname,abort_message,1)
                nfois(nsrf) = nfois(nsrf) + 1
             ENDIF
             snow(i,nsrf)     = 0.
             agesno(i,nsrf)   = 0.
             ftsoil(i,:,nsrf) = tsurf(i,nsrf)
#ifdef ISO            
             xtsnow(:,i,nsrf) = 0.
#endif
          ELSE
             pfois(nsrf) = pfois(nsrf)+ 1
          ENDIF
       ENDDO
       
    ENDDO

  END SUBROUTINE pbl_surface_newfrac

  SUBROUTINE pbl_surface_precall
  USE surf_landice_mod, ONLY : surf_landice_precall
  IMPLICIT NONE
    CALL surf_landice_precall
  END SUBROUTINE pbl_surface_precall

  SUBROUTINE pbl_surface_postcall
  USE surf_landice_mod, ONLY : surf_landice_postcall
  IMPLICIT NONE
    CALL surf_landice_postcall
  END SUBROUTINE pbl_surface_postcall


  SUBROUTINE pbl_surface_main( &
       dtime,     date0,     itap,     jour,          &
       debut,     lafin,                              &
       rlon,      rlat,      rugoro,   rmu0,          &
       lwdown_m,  pphi, cldt,          &
       rain_f,    snow_f,    bs_f, solsw_m,  solswfdiff_m, sollw_m,       &
       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,   zu1,    zv1,              &
       beta, &
       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,  &
       delta_tsurf,wake_dens,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,   agesno,  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,                                     &
       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      
     &   )

!****************************************************************************************
! 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 pbl_surface_uncompress_pre_mod, ONLY : pbl_surface_uncompress_pre
  USE pbl_surface_subsrf_mod, ONLY : pbl_surface_subsrf
  USE pbl_surface_uncompressed_post_mod, ONLY : pbl_surface_uncompressed_post

  USE dimphy, ONLY : klon, klev
  USE indice_sol_mod, ONLY : nbsrf, is_ter, is_oce, is_sic, is_lic
  USE clesphys_mod_h, ONLY : nsw
  USE dimsoil_mod_h, ONLY : nsoilmx
#ifdef ISO
  USE infotrac_phy, ONLY: ntraciso=>ntiso    
#endif
  USE print_control_mod,  ONLY : prt_level
  USE mod_phys_lmdz_para, ONLY : is_master
  USE print_control_mod, ONLY: lunout 
IMPLICIT NONE


!****************************************************************************************
    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),        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),   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
                                                                   !wake and off-wake regions
!albedo SB >>>
    REAL, DIMENSION(6), intent(in) :: SFRWL
    REAL, DIMENSION(klon, nsw, nbsrf), INTENT(INOUT)     :: alb_dir,alb_dif
!albedo SB <<<
!jyg Pourquoi ustar et wstar sont-elles INOUT ?
    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!) meltin 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
    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(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

! Martin
! inlandsis
    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
    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) :: n2mout
    INTEGER, DIMENSION(klon, nbsrf, 6) :: n2mout_x
    INTEGER, DIMENSION(klon, nbsrf, 6) :: n2mout_w
    REAL, DIMENSION(klon, klev)        :: d_u_x
    REAL, DIMENSION(klon, klev)        :: d_u_w
    REAL, DIMENSION(klon, klev)        :: d_v_x
    REAL, DIMENSION(klon, klev)        :: d_v_w 
    REAL, DIMENSION(klon, nbsrf)       :: windsp
    REAL, DIMENSION(klon, nbsrf)       :: t2m_x
    REAL, DIMENSION(klon, nbsrf)       :: q2m_x
    REAL, DIMENSION(klon)              :: rh2m_x
    REAL, DIMENSION(klon)              :: qsat2m_x
    REAL, DIMENSION(klon, nbsrf)       :: u10m_x
    REAL, DIMENSION(klon, nbsrf)       :: v10m_x
    REAL, DIMENSION(klon, nbsrf)       :: ustar_x
    REAL, DIMENSION(klon, nbsrf)       :: wstar_x
    REAL, DIMENSION(klon, nbsrf)       :: pblh_x
    REAL, DIMENSION(klon, nbsrf)       :: plcl_x
    REAL, DIMENSION(klon, nbsrf)       :: capCL_x
    REAL, DIMENSION(klon, nbsrf)       :: oliqCL_x
    REAL, DIMENSION(klon, nbsrf)       :: cteiCL_x
    REAL, DIMENSION(klon, nbsrf)       :: pblt_x
    REAL, DIMENSION(klon, nbsrf)       :: therm_x
    REAL, DIMENSION(klon, nbsrf)       :: trmb1_x
    REAL, DIMENSION(klon, nbsrf)       :: trmb2_x
    REAL, DIMENSION(klon, nbsrf)       :: trmb3_x
    REAL, DIMENSION(klon, nbsrf)       :: t2m_w
    REAL, DIMENSION(klon, nbsrf)       :: q2m_w
    REAL, DIMENSION(klon)              :: rh2m_w
    REAL, DIMENSION(klon)              :: qsat2m_w
    REAL, DIMENSION(klon, nbsrf)       :: u10m_w
    REAL, DIMENSION(klon, nbsrf)       :: v10m_w
    REAL, DIMENSION(klon, nbsrf)       :: ustar_w
    REAL, DIMENSION(klon, nbsrf)       :: wstar_w
!                           
    REAL, DIMENSION(klon, nbsrf)       :: pblh_w
    REAL, DIMENSION(klon, nbsrf)       :: plcl_w
    REAL, DIMENSION(klon, nbsrf)       :: capCL_w
    REAL, DIMENSION(klon, nbsrf)       :: oliqCL_w
    REAL, DIMENSION(klon, nbsrf)       :: cteiCL_w
    REAL, DIMENSION(klon, nbsrf)       :: pblt_w
    REAL, DIMENSION(klon, nbsrf)       :: therm_w
    REAL, DIMENSION(klon, nbsrf)       :: trmb1_w
    REAL, DIMENSION(klon, nbsrf)       :: trmb2_w
    REAL, DIMENSION(klon, nbsrf)       :: trmb3_w
!
    REAL, DIMENSION(klon,nbsrf)        :: pblh         ! height of the planetary boundary layer
    REAL, DIMENSION(klon,nbsrf)        :: plcl         ! condensation level
    REAL, DIMENSION(klon,nbsrf)        :: capCL
    REAL, DIMENSION(klon,nbsrf)        :: oliqCL
    REAL, DIMENSION(klon,nbsrf)        :: cteiCL
    REAL, DIMENSION(klon,nbsrf)        :: pblT
    REAL, DIMENSION(klon,nbsrf)        :: therm
    REAL, DIMENSION(klon,nbsrf)        :: trmb1        ! deep cape
    REAL, DIMENSION(klon,nbsrf)        :: trmb2        ! inhibition
    REAL, DIMENSION(klon,nbsrf)        :: trmb3        ! point Omega
    REAL, DIMENSION(klon, nbsrf)       :: alb          ! mean albedo for whole SW interval
    REAL, DIMENSION(klon,nbsrf)        :: snowerosion
    REAL, DIMENSION(klon,klev)         :: delp
    REAL, DIMENSION(klon,klev)         :: d_t_diss_x, d_t_diss_w
    REAL, DIMENSION(klon, klev, nbsrf) :: flux_t_x, flux_q_x, flux_t_w, flux_q_w
    REAL, DIMENSION(klon, klev, nbsrf) :: flux_u_x, flux_v_x, flux_u_w, flux_v_w
    REAL, DIMENSION(klon, nbsrf)       :: fluxlat_x, fluxlat_w
    


    INTEGER :: iflag_split_ref
    INTEGER :: nsrf
    INTEGER :: i
    INTEGER :: knon
    INTEGER :: ni(klon)
    
    
    CALL pbl_surface_precall

    CALL 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      
     &   )
    
  DO nsrf = 1, nbsrf                                        !<<<<<<<<<<<<<
!    IF (nsrf/=is_ter) CYCLE                                                                      !<<<<<<<<<<<<<
    IF (prt_level >=10) print *,' Loop nsrf ',nsrf

    ! 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

    
    CALL pbl_surface_subsrf( nsrf, knon, ni(1:knon),  &
       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)
  
  ENDDO
  

  CALL 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      
     &   )

      CALL pbl_surface_postcall

END SUBROUTINE pbl_surface_main

!  
END MODULE pbl_surface_main_mod
