compute_isba_parameters.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!#############################################################
SUBROUTINE compute_isba_parameters (DTCO, OREAD_BUDGETC, UG, U, &
                                    IO, DTI, SB, S, IG, K, NK, NIG, NP, NPE,   &
                                    NAG, NISS, ISS, NCHI, CHI, ID, GB, NGB,    &
                                    NDST, SLT, SV, HPROGRAM,HINIT,OLAND_USE,   &
                                    KI,KSV,KSW,HSV,PCO2,PRHOA,                 &
                                    PZENITH,PSW_BANDS,PDIR_ALB,PSCA_ALB,       &
                                    PEMIS,PTSRAD,PTSURF, HTEST             )  
!#############################################################
!
!!****  *COMPUTE_ISBA_PARAMETERS_n* - routine to initialize ISBA
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!      V. Masson   *Meteo France*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!      Modified by P. Le Moigne (11/2004): miscellaneous diagnostics
!!      Modified by P. Le Moigne (06/2006): seeding and irrigation    
!!      Modified by B. Decharme    (2008) : SGH and Flooding scheme
!!      Modified by B. Decharme  (01/2009): optional deep soil temperature as in Arpege
!!      Modified by R. Hamdi     (01/2009): Cp and L
!!      Modified by B. Decharme  (06/2009): read topographic index statistics
!!      Modified by P. Le Moigne (01/2009): Beljaars sso
!!      Modified by B. Decharme  (08/2009): Active Trip coupling variable if Earth System Model
!!      A.L. Gibelin   04/09 : change BSLAI_NITRO initialisation
!!      A.L. Gibelin   04/09 : modifications for CENTURY model 
!!      A.L. Gibelin   06/09 : soil carbon initialisation
!!      Modified by B. Decharme  (09/2012): Bug in exponential profile calculation with DIF
!!      F. Bouttier    08/13 : apply random perturbation patterns for ensembles
!!      B. Vincendon   03/14 : bug correction for CISBA=3L and CKSAT=EXP (TOPD coupling)
!!      Modified by B. Decharme  (04/2013): Subsurface runoff if SGH (DIF option only)
!!                                          Delete CTOPREG (never used)
!!                                          Delete NWG_LAYER_TOT, NWG_SIZE
!!                                          water table / Surface coupling
!!      P. Samuelsson  02/14 : MEB
!!      B. Decharme    01/16 : Bug when vegetation veg, z0 and emis are imposed whith interactive vegetation
!!      B. Decharme   10/2016  bug surface/groundwater coupling 
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_s_t, isba_p_t, isba_pe_t, isba_k_t, isba_nk_t, &
                        isba_np_t, isba_npe_t
USE modd_data_isba_n, ONLY : data_isba_t
USE modd_sfx_grid_n, ONLY : grid_t, grid_np_t
USE modd_agri_n, ONLY : agri_t, agri_np_t
USE modd_sso_n, ONLY : sso_t, sso_np_t
USE modd_ch_isba_n, ONLY : ch_isba_t, ch_isba_np_t
USE modd_canopy_n, ONLY : canopy_t
USE modd_gr_biog_n, ONLY : gr_biog_t, gr_biog_np_t
USE modd_surfex_n, ONLY : isba_diag_t
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_surf_atm_n, ONLY : surf_atm_t
USE modd_dst_n, ONLY : dst_np_t, dst_t
USE modd_slt_n, ONLY : slt_t
USE modd_sv_n, ONLY : sv_t
!
USE modd_sfx_oasis,  ONLY : lcpl_land, lcpl_flood, lcpl_gw, lcpl_calving
!
!
#ifdef TOPD
USE modd_dummy_exp_profile,ONLY : xc_depth_ratio
#endif
!
USE modd_assim, ONLY : cassim_isba, lassim
!
USE modd_deepsoil,       ONLY : lphysdomc, ldeepsoil, xtdeep_cli, xgammat_cli
USE modd_agri,           ONLY : lagrip, xthreshold
!
!
USE modd_sgh_par,        ONLY : ndimtab, xice_deph_max, xf_decay
!
USE modd_data_cover_par, ONLY : nvegtype
USE modd_surf_par,       ONLY : xundef, nundef
USE modd_snow_par,       ONLY : xemissn
!
USE modd_topd_par, ONLY : nunit
USE modd_topodyn, ONLY : nncat, nmesht
!
USE mode_random
!
USE modi_get_1d_mask
USE modi_get_z0rel
USE modi_get_luout
USE modi_abor1_sfx
USE modi_init_io_surf_n
USE modi_allocate_physio
USE modi_init_isba_mixpar
USE modi_convert_patch_isba
USE modi_init_veg_pgd_n
USE modi_init_top
USE modi_exp_decay_soil_fr
USE modi_carbon_init
USE modi_soiltemp_arp_par
USE modi_end_io_surf_n
!
USE modi_make_choice_array
USE modi_read_surf
USE modi_read_isba_n
USE modi_init_isba_landuse
USE modi_read_sbl_n
USE modi_init_veg_n
USE modi_init_chemical_n
USE modi_open_namelist
USE modi_ch_init_dep_isba_n
USE modi_close_namelist
USE modi_init_dst
USE modi_init_slt
USE modi_averaged_albedo_emis_isba
USE modi_diag_isba_init_n
USE modi_init_surf_topd
USE modi_isba_soc_parameters
USE modi_pack_same_rank
!
USE modi_read_and_send_mpi
USE modi_isba_to_topd
USE modi_open_file
USE modi_close_file
USE modi_fix_meb_veg
USE modi_av_pgd
USE modi_surf_patch
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments
!              -------------------------
!
TYPE(data_cover_t), INTENT(INOUT) :: DTCO
LOGICAL, INTENT(IN) :: OREAD_BUDGETC
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
TYPE(surf_atm_t), INTENT(INOUT) :: U
!
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(data_isba_t), INTENT(INOUT) :: DTI
TYPE(canopy_t), INTENT(INOUT) :: SB
TYPE(isba_s_t), INTENT(INOUT) :: S
TYPE(grid_t), INTENT(INOUT) :: IG
TYPE(isba_k_t), INTENT(INOUT) :: K
TYPE(isba_nk_t), INTENT(INOUT) :: NK
TYPE(grid_np_t), INTENT(INOUT) :: NIG
TYPE(isba_np_t), INTENT(INOUT) :: NP
TYPE(isba_npe_t), INTENT(INOUT) :: NPE
TYPE(agri_np_t), INTENT(INOUT) :: NAG
TYPE(sso_np_t), INTENT(INOUT) :: NISS
TYPE(sso_t), INTENT(INOUT) :: ISS
TYPE(ch_isba_np_t), INTENT(INOUT) :: NCHI
TYPE(ch_isba_t), INTENT(INOUT) :: CHI
TYPE(isba_diag_t), INTENT(INOUT) :: ID
TYPE(gr_biog_t), INTENT(INOUT) :: GB
TYPE(gr_biog_np_t), INTENT(INOUT) :: NGB
!
TYPE(dst_np_t), INTENT(INOUT) :: NDST
TYPE(slt_t), INTENT(INOUT) :: SLT
TYPE(sv_t), INTENT(INOUT) :: SV
!
 CHARACTER(LEN=6),                INTENT(IN)  :: HPROGRAM  ! program calling surf. schemes
 CHARACTER(LEN=3),                INTENT(IN)  :: HINIT     ! choice of fields to initialize
LOGICAL,                          INTENT(IN)  :: OLAND_USE !
INTEGER,                          INTENT(IN)  :: KI        ! number of points
INTEGER,                          INTENT(IN)  :: KSV       ! number of scalars
INTEGER,                          INTENT(IN)  :: KSW       ! number of short-wave spectral bands
 CHARACTER(LEN=6), DIMENSION(KSV),INTENT(IN)  :: HSV       ! name of all scalar variables
REAL,             DIMENSION(KI),  INTENT(IN)  :: PCO2      ! CO2 concentration (kg/m3)
REAL,             DIMENSION(KI),  INTENT(IN)  :: PRHOA     ! air density
REAL,             DIMENSION(KI),  INTENT(IN)  :: PZENITH   ! solar zenithal angle
REAL,             DIMENSION(KSW), INTENT(IN)  :: PSW_BANDS ! middle wavelength of each band
REAL,             DIMENSION(KI,KSW),INTENT(OUT) :: PDIR_ALB  ! direct albedo for each band
REAL,             DIMENSION(KI,KSW),INTENT(OUT) :: PSCA_ALB  ! diffuse albedo for each band
REAL,             DIMENSION(KI),  INTENT(OUT) :: PEMIS     ! emissivity
REAL,             DIMENSION(KI),  INTENT(OUT) :: PTSRAD    ! radiative temperature
REAL,             DIMENSION(KI),  INTENT(OUT) :: PTSURF    ! surface effective temperature         (K)
!
 CHARACTER(LEN=2),                 INTENT(IN)  :: HTEST       ! must be equal to 'OK'
!
!
!*       0.2   Declarations of local variables
!              -------------------------------
!
TYPE(grid_t), POINTER :: GK
TYPE(isba_p_t), POINTER :: PK
TYPE(isba_k_t), POINTER :: KK
TYPE(isba_pe_t), POINTER :: PEK
TYPE(agri_t), POINTER :: AGK
TYPE(sso_t), POINTER :: ISSK
TYPE(dst_t), POINTER :: DSTK
!
REAL, DIMENSION(U%NDIM_FULL)   :: ZF_PARAM, ZC_DEPTH_RATIO
!
REAL, DIMENSION(KI)     :: ZTSRAD_NAT !radiative temperature
REAL, DIMENSION(KI)     :: ZTSURF_NAT !effective temperature
REAL, DIMENSION(KI)     :: ZM
!
REAL, DIMENSION(KI)  :: ZWG1 ! work array for surface water content
REAL, DIMENSION(KI,IO%NPATCH)  :: ZTG1 ! work array for surface temperature
REAL, DIMENSION(KI,IO%NPATCH)  :: ZF
!
REAL, DIMENSION(:,:), ALLOCATABLE :: ZWORK
REAL, DIMENSION(:,:), ALLOCATABLE :: ZDG_SOIL, ZDG_SOIL_P
REAL, DIMENSION(:), ALLOCATABLE :: ZSUM_PATCH
!
INTEGER :: ICH     ! unit of input chemistry file
INTEGER           :: JI, JL     ! loop increment
INTEGER           :: ILUOUT   ! unit of output listing file
INTEGER           :: IRESP   ! return code
INTEGER           :: IDECADE, IDECADE2  ! decade of simulation
INTEGER           :: JP  ! loop counter on tiles
INTEGER           :: ISIZE_LMEB_PATCH  ! Number of patches with MEB=true
!
LOGICAL :: GDIM, GCAS1, GCAS2, GCAS3 
INTEGER :: JVEG, IVERSION, IBUGFIX, IMASK, JMAXLOC
!
 CHARACTER(LEN=4)  :: YLVL
 CHARACTER(LEN=12) :: YRECFM
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
!               Initialisation for IO
!
IF (lhook) CALL dr_hook('COMPUTE_ISBA_PARAMETERS',0,zhook_handle)
 CALL get_luout(hprogram,iluout)
!
IF (htest/='OK') THEN
  CALL abor1_sfx('COMPUTE_ISBA_PARAMETERS: FATAL ERROR DURING ARGUMENT TRANSFER')
END IF
!
!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------
!
!    PART 1 : Arrays of vegtypes & patches 
!    -------------------------------------
!
! We need XVEGTYPE, XPATCH and XVEGTYPE_PATCH with dimension "PATCH" for some
! cases: initialized here
!
! Vegtypes first
ALLOCATE(s%XVEGTYPE(ki,nvegtype))
IF (dti%LDATA_VEGTYPE) THEN
  s%XVEGTYPE = dti%XPAR_VEGTYPE
ELSE
  !classical ecoclimap case
  DO jveg=1,nvegtype
    CALL av_pgd(dtco, s%XVEGTYPE(:,jveg),s%XCOVER ,dtco%XDATA_VEGTYPE(:,jveg),'NAT','ARI',s%LCOVER)
  END DO
ENDIF
!
! patches come from vegtypes
ALLOCATE(s%XPATCH(ki,io%NPATCH))
ALLOCATE(s%XVEGTYPE_PATCH(ki,nvegtype,io%NPATCH))
 CALL surf_patch(io%NPATCH,s%XVEGTYPE,s%XPATCH,s%XVEGTYPE_PATCH)
!
! removing little fractions of patches must be done of the XPATCH with dimension
! "PATCH"
IF (io%XRM_PATCH/=0.) THEN
  !
  WRITE(iluout,*) " REMOVE PATCH below 5 % add to dominant patch " 
  ! remove small fraction of PATCHES and add to MAIN PATCH
  DO ji = 1,ki
    !1) find most present patch maximum value 
    jmaxloc = maxval(maxloc(s%XPATCH(ji,:)))
    !2) FIND small value of cover 
    DO jp = 1,io%NPATCH
      IF ( s%XPATCH(ji,jp)<io%XRM_PATCH ) THEN
        s%XPATCH(ji,jmaxloc) = s%XPATCH(ji,jmaxloc) + s%XPATCH(ji,jp)
        s%XPATCH(ji,jp) = 0.0
       ENDIF
    ENDDO
  ENDDO
  !
ENDIF
!
!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------
!
!    PART 2 : Things depending only on options and / or needed first
!    --------------------------------------------------------------

!*       Physiographic data fields from land cover:
!         -----------------------------------------
!
IF (s%TTIME%TDATE%MONTH /= nundef) THEN
  idecade = 3 * ( s%TTIME%TDATE%MONTH - 1 ) + min(s%TTIME%TDATE%DAY-1,29) / 10 + 1
ELSE
  idecade = 1
END IF
!
idecade2 = idecade
!
! concern DATA_ISBA, so no dependence on patches
 CALL init_isba_mixpar(dtco, dti, ig%NDIM, io, idecade, idecade2, s%XCOVER, s%LCOVER, 'NAT')
!
isize_lmeb_patch=count(io%LMEB_PATCH(:))
IF (isize_lmeb_patch>0)  THEN
  CALL fix_meb_veg(dti, ig%NDIM, io%LMEB_PATCH, io%NPATCH)
ENDIF
!
!
!*       Soil carbon
!        -----------                        
!
IF (hinit == 'ALL' .AND. io%CRESPSL=='CNT' .AND. io%CPHOTO == 'NCB') CALL carbon_init
!
!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------
!
!    PART 3 : Loop on patches for general initialization
!    --------------------------------------------------
!
! loop on patches
DO jp = 1, io%NPATCH
  !
  kk => nk%AL(jp)  
  pk => np%AL(jp)
  pek => npe%AL(jp)
  agk => nag%AL(jp)
  issk => niss%AL(jp)
  !
  ! dimension of the patch 
  pk%NSIZE_P = count(s%XPATCH(:,jp) > 0.0)
  !
  ! mask of the patch in tile nature
  ALLOCATE(pk%NR_P    (pk%NSIZE_P))
  CALL get_1d_mask(pk%NSIZE_P, ki, s%XPATCH(:,jp), pk%NR_P) 
  !
  ! the array of vegtypes, patches and vegtypes by patches reduced on this patches
  ALLOCATE(kk%XVEGTYPE(pk%NSIZE_P,nvegtype))
  CALL pack_same_rank(pk%NR_P,s%XVEGTYPE,kk%XVEGTYPE)
  !
  ALLOCATE(pk%XPATCH(pk%NSIZE_P))
  ALLOCATE(pk%XVEGTYPE_PATCH (pk%NSIZE_P,nvegtype))
  CALL pack_same_rank(pk%NR_P,s%XPATCH(:,jp),pk%XPATCH)
  CALL pack_same_rank(pk%NR_P,s%XVEGTYPE_PATCH(:,:,jp),pk%XVEGTYPE_PATCH)
  !
  !
  ! soon needed packed fields
  !
  IF (io%LPERM) THEN
    ALLOCATE(kk%XPERM(pk%NSIZE_P))
    CALL pack_same_rank(pk%NR_P, k%XPERM, kk%XPERM)
  ELSE
    ALLOCATE(kk%XPERM(0))
  ENDIF
  !
  !
  ALLOCATE(kk%XSAND(pk%NSIZE_P,io%NGROUND_LAYER))
  ALLOCATE(kk%XCLAY(pk%NSIZE_P,io%NGROUND_LAYER))
  !
  ALLOCATE(issk%XAOSIP(pk%NSIZE_P))
  ALLOCATE(issk%XAOSIM(pk%NSIZE_P))
  ALLOCATE(issk%XAOSJP(pk%NSIZE_P))
  ALLOCATE(issk%XAOSJM(pk%NSIZE_P))
  ALLOCATE(issk%XHO2IP(pk%NSIZE_P))
  ALLOCATE(issk%XHO2IM(pk%NSIZE_P))
  ALLOCATE(issk%XHO2JP(pk%NSIZE_P))
  ALLOCATE(issk%XHO2JM(pk%NSIZE_P))
  !
  !
  CALL pack_same_rank(pk%NR_P, k%XSAND, kk%XSAND)
  CALL pack_same_rank(pk%NR_P, k%XCLAY, kk%XCLAY)
  !
  CALL pack_same_rank(pk%NR_P,iss%XAOSIP,issk%XAOSIP)
  CALL pack_same_rank(pk%NR_P,iss%XAOSIM,issk%XAOSIM)
  CALL pack_same_rank(pk%NR_P,iss%XAOSJP,issk%XAOSJP)
  CALL pack_same_rank(pk%NR_P,iss%XAOSJM,issk%XAOSJM)
  CALL pack_same_rank(pk%NR_P,iss%XHO2IP,issk%XHO2IP)
  CALL pack_same_rank(pk%NR_P,iss%XHO2IM,issk%XHO2IM)
  CALL pack_same_rank(pk%NR_P,iss%XHO2JP,issk%XHO2JP)
  CALL pack_same_rank(pk%NR_P,iss%XHO2JM,issk%XHO2JM)
  !
  !
  !*       2.5    Physiographic fields
  !               --------------------
  !
  CALL allocate_physio(io, kk, pk, pek, nvegtype  )
  !
  CALL convert_patch_isba(dtco, dti, io, idecade, idecade2, s%XCOVER, s%LCOVER, &
                          lagrip, 'NAT', jp, kk, pk, pek, &
                          .true., .true., .true., .true., .false., .false., &
                          psoilgrid=io%XSOILGRID, pperm=kk%XPERM  )
  !
  !-------------------------------------------------------------------------------
  !
  ! in init_veg_pgd_n, things needed also by garden and greenroof
  CALL init_veg_pgd_n(issk, dti, io, s, k, kk, pk, pek, agk, ki,  &
                      hprogram, 'NATURE', iluout, pk%NSIZE_P, s%TTIME%TDATE%MONTH, &
                      ldeepsoil, lphysdomc, xtdeep_cli, xgammat_cli,   & 
                      lagrip, xthreshold, hinit, pco2, prhoa        )
  !
  !-------------------------------------------------------------------------------
  !
  ! Other fields needed to be initialized for isba only
  !
  !Rainfall spatial distribution
  !CRAIN used in HYDRO_VEG and HYDRO_SGH and VEG_SGH_UPDATE
  IF(io%CRAIN=='SGH')THEN
    ALLOCATE(kk%XMUF(pk%NSIZE_P))
    kk%XMUF(:)=0.0
  ELSE
    ALLOCATE(kk%XMUF(0))
  ENDIF
  !
  ALLOCATE(kk%XFSAT(pk%NSIZE_P))  
  kk%XFSAT(:) = 0.0
  !
  ! * Initialize flood scheme :
  !
  ALLOCATE(kk%XFFLOOD (pk%NSIZE_P))
  ALLOCATE(kk%XPIFLOOD(pk%NSIZE_P))
  ALLOCATE(kk%XFF     (pk%NSIZE_P))
  ALLOCATE(kk%XFFG    (pk%NSIZE_P))
  ALLOCATE(kk%XFFV    (pk%NSIZE_P))
  ALLOCATE(kk%XFFROZEN(pk%NSIZE_P))
  ALLOCATE(kk%XALBF   (pk%NSIZE_P))
  ALLOCATE(kk%XEMISF  (pk%NSIZE_P)) 
  kk%XFFLOOD       = 0.0
  kk%XPIFLOOD      = 0.0
  kk%XFF           = 0.0
  kk%XFFG          = 0.0
  kk%XFFV          = 0.0
  kk%XFFROZEN      = 0.0
  kk%XALBF         = 0.0
  kk%XEMISF        = 0.0  
  !
ENDDO
!
IF (dti%LDATA_CONDSAT) DEALLOCATE(dti%XPAR_CONDSAT)
!
!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------
!
!    PART 4 : Initialization not depending on patches 
!    ------------------------------------------------
!
! Fields needed also unpacked 
!
IF(io%CRAIN=='SGH')THEN
  ALLOCATE(k%XMUF(ki))
  k%XMUF(:)=0.0
ENDIF
!
!
ALLOCATE(iss%XZ0REL(ki))
 CALL get_z0rel(iss)
!
!-------------------------------------------------------------------------------
!
!        PART 5:  Initialize Chemical Deposition
!            -----------------------------------
!
!        3.1 Chemical gazes
!            --------------
!
    !* for the time being, chemistry on vegetation works only for
    ! ISBA on nature tile (not for gardens), because subroutine INIT_CHEMICAL_n
    ! contains explicitely modules from ISBAn. It should be cleaned in a future
    ! version.
 CALL init_chemical_n(iluout, ksv, hsv, chi%SVI, chi%CCH_NAMES, chi%CAER_NAMES,  &
                      hdstnames=chi%CDSTNAMES, hsltnames=chi%CSLTNAMES        )
!
IF (ksv /= 0) THEN
  !
  IF (chi%SVI%NBEQ > 0) THEN
    !* for the time being, chemistry deposition on vegetation works only for
    ! ISBA on nature tile (not for gardens), because subroutine CH_INIT_DEP_ISBA_n
    ! contains explicitely modules from ISBAn. It should be cleaned in a future
    ! version.
    CALL open_namelist(hprogram, ich, hfile=chi%CCHEM_SURF_FILE)
    CALL ch_init_dep_isba_n(chi, nchi, np, dtco, io%NPATCH, s%LCOVER, s%XCOVER, ich, iluout, ki)
    CALL close_namelist(hprogram, ich)
  END IF
  !
  DO jp = 1,io%NPATCH
    !
    dstk => ndst%AL(jp)
    !
    IF (chi%SVI%NDSTEQ >=1) THEN
      !
      ALLOCATE (dstk%XSFDST (pk%NSIZE_P, chi%SVI%NDSTEQ))  !Output array
      ALLOCATE (dstk%XSFDSTM(pk%NSIZE_P, chi%SVI%NDSTEQ))  !Output array
      dstk%XSFDST (:,:)  = 0.
      dstk%XSFDSTM(:,:) = 0.     
      CALL init_dst(dstk, u, hprogram, pk%NSIZE_P, pk%NR_P, pk%XVEGTYPE_PATCH)    
    ELSE
      ALLOCATE(dstk%XSFDST (0,0))
      ALLOCATE(dstk%XSFDSTM(0,0))
    END IF
    !
  ENDDO
  !
  IF (chi%SVI%NSLTEQ >=1) THEN
    CALL init_slt(slt, hprogram)
  END IF
  !
ENDIF
!
!-------------------------------------------------------------------------------
!
!        PART 6:  Specific options
!        --------------------------

!6.A. DIF option :
!---------------
!    Anisotropy coeficient for hydraulic conductivity for topmodel drainage (Fan et al. 2006)
!    Soil organic matter effect and/or Exponential decay for DIF option
!    Must be call before INIT_TOP
!
!
IF(io%CISBA=='DIF' .AND. io%CKSAT=='SGH') THEN
  !
  WRITE(iluout,*)'THE KSAT EXP PROFILE WITH ISBA-DF IS NOT PHYSIC AND HAS BEEN REMOVED FOR NOW' 
  WRITE(iluout,*)'A NEW PHYSICAL APPROACH WILL BE DEVELLOPED ACCOUNTING FOR COMPACTION IN ALL ' 
  WRITE(iluout,*)'HYDRODYNAMIC PARAMETERS (WSAT, PSISAT, KSAT, B) AND NOT ONLY IN KSAT        ' 
  CALL abor1_sfx('CKSAT=SGH is not physic with ISBA-DF and has been removed for now')
  !    
ENDIF
!  
IF(io%CISBA=='DIF' .AND. io%LSOC)THEN   
  !
  IF(.NOT.io%LSOCP)THEN
    WRITE(iluout,*)'LSOC = T can be activated only if SOC data given in PGD fields'
    CALL abor1_sfx('LSOC = T can be activated only if SOC data given in PGD fields')
  ENDIF
  !
  ALLOCATE(s%XFRACSOC(ki,io%NGROUND_LAYER))
  CALL isba_soc_parameters(io%CRUNOFF, s%XSOC, k, np, s%XFRACSOC, &
                           k%XWSAT, k%XWFC, k%XWWILT, io%NPATCH )
  !
ELSE
  ALLOCATE(s%XFRACSOC(0,0))
ENDIF
!
!
!6.B. Topmodel
!--------------
!
zf(:,:) = xundef
zm(:)   = xundef
!
!CRUNOFF used in hydro_sgh and isba_sgh_update
IF( io%CRUNOFF=='SGH '.AND. hinit/='PRE' .AND. .NOT.lassim ) THEN 
  !
  ! Subsurface flow by layer (m/s)
  DO jp = 1,io%NPATCH
    pk => np%AL(jp)
    IF(io%CISBA=='DIF') THEN
      ALLOCATE(pk%XTOPQS(pk%NSIZE_P,io%NGROUND_LAYER))
      pk%XTOPQS(:,:) = 0.0
    ELSE
      ALLOCATE(pk%XTOPQS(0,0))
    ENDIF
  ENDDO
  !  
  ALLOCATE(s%XTAB_FSAT(ki,ndimtab))
  ALLOCATE(s%XTAB_WTOP(ki,ndimtab))
  ALLOCATE(s%XTAB_QTOP(ki,ndimtab))
  s%XTAB_FSAT(:,:) = 0.0
  s%XTAB_WTOP(:,:) = 0.0
  s%XTAB_QTOP(:,:) = 0.0
  !
  WHERE(k%XCLAY(:,1)==xundef.AND.s%XTI_MEAN(:)/=xundef) s%XTI_MEAN(:)=xundef
  CALL init_top(io, s, k, nk, np, iluout, zm )
  !
ELSE                  
  !  
  DO jp = 1,io%NPATCH
    pk => np%AL(jp)
    ALLOCATE(pk%XTOPQS(0,0))
  ENDDO  
  !  
  ALLOCATE(s%XTAB_FSAT(0,0))
  ALLOCATE(s%XTAB_WTOP(0,0))
  ALLOCATE(s%XTAB_QTOP(0,0))
  !                  
ENDIF  
!
!
!Exponential decay for ISBA-FR option
!CKSAT used in hydro_soil.F90 and soil.F90
IF ( io%CISBA/='DIF' .AND. hinit/='PRE' .AND. .NOT.lassim ) THEN 
  !
  gcas1 = (io%CKSAT=='EXP' .AND. io%CISBA=='3-L')
  gcas2 = (io%CKSAT=='SGH')
  gcas3 = (hprogram/='AROME ' .AND. hprogram/='MESONH ')
  !
  IF ( gcas1 .OR. gcas2 ) THEN
    !
    ALLOCATE(s%XF_PARAM (ki))
    s%XF_PARAM(:) = xundef
    !
    IF ( gcas1 .AND. gcas3 ) THEN
      !
      !reading of XF_PARAM in external file
      CALL open_file('ASCII ',nunit,hfile='carte_f_dc.txt',hform='FORMATTED',haction='READ ')
      DO ji = 1,u%NDIM_FULL
        READ(nunit,*) zf_param(ji), zc_depth_ratio(ji)
      ENDDO
      CALL close_file('ASCII ',nunit)
      CALL read_and_send_mpi(zf_param,s%XF_PARAM,u%NR_NATURE)
#ifdef TOPD
      IF (.NOT.ALLOCATED(xc_depth_ratio))  ALLOCATE(xc_depth_ratio(ki))
      xc_depth_ratio(:) = xundef
      CALL read_and_send_mpi(zc_depth_ratio,xc_depth_ratio,u%NR_NATURE)
#endif
      !
    ELSEIF ( gcas1 ) THEN
      WRITE(iluout,*) "COMPUTE_ISBA_PARAMETERS: WITH CKSAT=EXP, IN NOT OFFLINE "//&
                      "MODE, TOPMODEL FILE FOR F_PARAM IS NOT READ "
    ENDIF
    !
    ! definition of ZF functions of options
    !
    ! Exponential decay factor calculate using soil properties 
    ! (eq. 11, Decharme et al., J. Hydrometeor, 2006)
    DO jp = 1,io%NPATCH
      pk => np%AL(jp)
      !
      DO ji = 1,pk%NSIZE_P
        imask = pk%NR_P(ji)

        IF ( gcas2 .AND. io%CRUNOFF=='SGH' .AND. zm(imask)/=xundef ) THEN
          zf(ji,jp) = (k%XWSAT(imask,1)-k%XWD0(imask,1)) / zm(imask)
        ELSEIF ( gcas1 ) THEN
          zf(ji,jp) = s%XF_PARAM(imask)
         ENDIF
      ENDDO 
    ENDDO
    !
    DO jp = 1,io%NPATCH
      pk => np%AL(jp)
      !
      WHERE ( zf(1:pk%NSIZE_P,jp)==xundef.AND.pk%XDG(:,2)/=xundef ) 
        zf(1:pk%NSIZE_P,jp) = 4.0/pk%XDG(:,2)
      ENDWHERE
      zf(1:pk%NSIZE_P,jp) = min(zf(1:pk%NSIZE_P,jp),xf_decay)
      !
      zc_depth_ratio(1:pk%NSIZE_P) = 1.
#ifdef TOPD      
      IF (ALLOCATED(xc_depth_ratio)) THEN
        CALL pack_same_rank(pk%NR_P,xc_depth_ratio,zc_depth_ratio(1:pk%NSIZE_P))
      ENDIF
#endif      
      CALL exp_decay_soil_fr(io%CISBA, zf(1:pk%NSIZE_P,jp), pk, zc_depth_ratio(1:pk%NSIZE_P))
    ENDDO
    !
    IF ( gcas2 ) THEN 
      !
      DO ji = 1,np%AL(1)%NSIZE_P
        imask = np%AL(1)%NR_P(ji)
        s%XF_PARAM(imask) = zf(ji,1)
      ENDDO
      !
    ENDIF
    !
  ENDIF
  ! 
ENDIF
!
!
! 6.C. Initialize required coupling fields :
!-------------------------------------------
!
io%LCPL_RRM = .false.
io%LFLOOD   = .false.
io%LWTD     = .false.
!
IF(lcpl_land)THEN
!    
  io%LCPL_RRM = .true.
!
  IF(lcpl_gw)THEN
    io%LWTD = .true.
  ENDIF
!
  ALLOCATE(s%XCPL_DRAIN (ki))
  ALLOCATE(s%XCPL_RUNOFF(ki))
  s%XCPL_DRAIN (:) = 0.0
  s%XCPL_RUNOFF(:) = 0.0
!
  IF(io%LGLACIER)THEN
     ALLOCATE(s%XCPL_ICEFLUX(ki))
     s%XCPL_ICEFLUX(:) = 0.0
  ELSE
     ALLOCATE(s%XCPL_ICEFLUX(0))
  ENDIF
!
  IF(lcpl_flood)THEN
     io%LFLOOD = .true.
     ALLOCATE(s%XCPL_EFLOOD(ki))
     ALLOCATE(s%XCPL_PFLOOD(ki))
     ALLOCATE(s%XCPL_IFLOOD(ki))
     s%XCPL_EFLOOD(:)= 0.0
     s%XCPL_PFLOOD(:)= 0.0
     s%XCPL_IFLOOD(:)= 0.0    
  ELSE
    ALLOCATE(s%XCPL_EFLOOD(0))
    ALLOCATE(s%XCPL_PFLOOD(0))
    ALLOCATE(s%XCPL_IFLOOD(0))     
  ENDIF     
!
ELSE
!
  ALLOCATE(s%XCPL_RUNOFF  (0))
  ALLOCATE(s%XCPL_DRAIN   (0))
  ALLOCATE(s%XCPL_ICEFLUX (0))
  ALLOCATE(s%XCPL_EFLOOD  (0))
  ALLOCATE(s%XCPL_PFLOOD  (0))
  ALLOCATE(s%XCPL_IFLOOD  (0))
!
ENDIF
!
!
IF (lcpl_land) THEN
  !
  ALLOCATE(k%XFWTD(ki))
  ALLOCATE(k%XWTD (ki))
  k%XFWTD(:) = 0.0
  k%XWTD (:) = xundef
  !
  IF(lcpl_flood)THEN
    ALLOCATE(k%XFFLOOD (ki))
    ALLOCATE(k%XPIFLOOD(ki))
    k%XFFLOOD (:) = 0.0
    k%XPIFLOOD(:) = 0.0
    !
  ELSE
    !
    ALLOCATE(k%XFFLOOD (0))
    ALLOCATE(k%XPIFLOOD(0))
    !   
  ENDIF
  !
ELSE
  !
  ALLOCATE(k%XFWTD(0))
  ALLOCATE(k%XWTD (0))
  ALLOCATE(k%XFFLOOD (0))
  ALLOCATE(k%XPIFLOOD(0))
  !   
ENDIF
!
! * Check some key :
!
IF(lcpl_calving)THEN
   IF(.NOT.io%LGLACIER)THEN
     CALL abor1_sfx('COMPUTE_ISBA_PARAMETERS: LGLACIER MUST BE ACTIVATED IF LCPL_CALVING')
   ENDIF
ENDIF
!
!-------------------------------------------------------------------------------
!
!*   6.D. ISBA time-varying deep force-restore temperature initialization
!    --------------------------------------------------------------------
!
 CALL soiltemp_arp_par(io, hprogram)
!
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!
!        PART 7:  We packed needed fields and free unless ones
!        -----------------------------------------------------
!
! 
DO jp = 1,io%NPATCH
  !
  kk => nk%AL(jp)
  pk => np%AL(jp)
  issk => niss%AL(jp)
  gk => nig%AL(jp)
  !
  ALLOCATE(kk%XMPOTSAT(pk%NSIZE_P,io%NGROUND_LAYER))
  ALLOCATE(kk%XBCOEF  (pk%NSIZE_P,io%NGROUND_LAYER))
  ! needed to be written as diagnostics, so not free
  ALLOCATE(kk%XWWILT  (pk%NSIZE_P,io%NGROUND_LAYER))
  ALLOCATE(kk%XWFC    (pk%NSIZE_P,io%NGROUND_LAYER))
  ALLOCATE(kk%XWSAT   (pk%NSIZE_P,io%NGROUND_LAYER))
  !
  CALL pack_same_rank(pk%NR_P,k%XMPOTSAT,kk%XMPOTSAT)
  CALL pack_same_rank(pk%NR_P,k%XBCOEF,kk%XBCOEF)
  !
  CALL pack_same_rank(pk%NR_P,k%XWWILT,kk%XWWILT)
  CALL pack_same_rank(pk%NR_P,k%XWFC,kk%XWFC)
  CALL pack_same_rank(pk%NR_P,k%XWSAT,kk%XWSAT)
  !  
  IF (io%CISBA=='2-L' .OR. io%CISBA=='3-L') THEN
    ALLOCATE(kk%XCGSAT(pk%NSIZE_P))          
    ALLOCATE(kk%XC4B  (pk%NSIZE_P))
    ALLOCATE(kk%XACOEF(pk%NSIZE_P))
    ALLOCATE(kk%XPCOEF(pk%NSIZE_P))
    CALL pack_same_rank(pk%NR_P,k%XCGSAT,kk%XCGSAT)    
    CALL pack_same_rank(pk%NR_P,k%XC4B,  kk%XC4B)
    CALL pack_same_rank(pk%NR_P,k%XACOEF,kk%XACOEF)
    CALL pack_same_rank(pk%NR_P,k%XPCOEF,kk%XPCOEF)
  ENDIF
  !
  IF (io%CSCOND=='PL98'.OR.io%CISBA=='DIF') THEN
    ALLOCATE(kk%XHCAPSOIL(pk%NSIZE_P,io%NGROUND_LAYER))
    ALLOCATE(kk%XCONDDRY (pk%NSIZE_P,io%NGROUND_LAYER))
    ALLOCATE(kk%XCONDSLD (pk%NSIZE_P,io%NGROUND_LAYER))
    CALL pack_same_rank(pk%NR_P,k%XHCAPSOIL,kk%XHCAPSOIL)
    CALL pack_same_rank(pk%NR_P,k%XCONDDRY ,kk%XCONDDRY)
    CALL pack_same_rank(pk%NR_P,k%XCONDSLD ,kk%XCONDSLD)
  ENDIF
  !
  ALLOCATE(kk%XWDRAIN (pk%NSIZE_P))
  ALLOCATE(kk%XRUNOFFB(pk%NSIZE_P))
  CALL pack_same_rank(pk%NR_P,k%XWDRAIN,kk%XWDRAIN)
  CALL pack_same_rank(pk%NR_P,k%XRUNOFFB,kk%XRUNOFFB)
  !
  ! needed to be written as diagnostics, so not free
  ALLOCATE(issk%XZ0REL    (pk%NSIZE_P))
  ALLOCATE(issk%XSSO_SLOPE(pk%NSIZE_P))
  !
  CALL pack_same_rank(pk%NR_P,iss%XZ0REL,issk%XZ0REL)
  CALL pack_same_rank(pk%NR_P,iss%XSSO_SLOPE,issk%XSSO_SLOPE)
  !
  ALLOCATE(gk%XLAT(pk%NSIZE_P))
  ALLOCATE(gk%XLON(pk%NSIZE_P))
  !
  CALL pack_same_rank(pk%NR_P,ig%XLAT,gk%XLAT)
  CALL pack_same_rank(pk%NR_P,ig%XLON,gk%XLON)
  !
ENDDO
!
! Useledd fields from now on
iss%XAOSIP => null()
iss%XAOSIM => null()
iss%XAOSJP => null()
iss%XAOSJM => null()
iss%XHO2IP => null()
iss%XHO2IM => null()
iss%XHO2JP => null()
iss%XHO2JM => null()
!
k%XMPOTSAT => null()
k%XBCOEF   => null()
!
k%XCGSAT => null()
k%XC4B   => null()
k%XACOEF => null()
k%XPCOEF => null()  
!
k%XHCAPSOIL => null()
k%XCONDDRY  => null()
k%XCONDSLD  => null()
!
k%XWDRAIN  => null()
k%XRUNOFFB => null()
!
!-------------------------------------------------------------------------------
!
!* if only physiographic fields are to be initialized, stop here.
!
IF (hinit/='ALL' .AND. hinit/='SOD') THEN
  IF (lhook) CALL dr_hook('COMPUTE_ISBA_PARAMETERS',1,zhook_handle)
  RETURN
END IF
!
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!
!        PART 8: Reading of prognostic variables
!        ----------------------------------------
!
IF (cassim_isba=="ENKF ") CALL init_random_seed()
!
!
 CALL init_io_surf_n(dtco, u, hprogram,'NATURE','ISBA  ','READ ')
!
!*      10.     Prognostic and semi-prognostic fields
!               -------------------------------------
!
 CALL read_isba_n(dtco, io, s, np, npe, k%XCLAY, u, hprogram)
!
IF (hinit/='ALL') THEN
  CALL end_io_surf_n(hprogram)
  IF (lhook) CALL dr_hook('COMPUTE_ISBA_PARAMETERS',1,zhook_handle)
  RETURN
END IF
!
IF (hinit=='PRE' .AND. npe%AL(1)%TSNOW%SCHEME.NE.'3-L' .AND. &
        npe%AL(1)%TSNOW%SCHEME.NE.'CRO' .AND. io%CISBA=='DIF') &
   CALL abor1_sfx("INIT_ISBAN: WITH CISBA = DIF, CSNOW MUST BE 3-L OR CRO")
!
!
!*       Extrapolation of the prognostic and semi-prognostic fields
!                           LAND USE case 
!               -------------------------------------
!
IF (oland_use) THEN
  !
  CALL read_surf(hprogram,'VERSION',iversion,iresp)
  CALL read_surf(hprogram,'BUG',ibugfix,iresp)
  gdim = (iversion>8 .OR. iversion==8 .AND. ibugfix>0)
  IF (gdim) CALL read_surf(hprogram,'SPLIT_PATCH',gdim,iresp)
  !  
  ALLOCATE(zwork(ki,io%NPATCH))
  !
  !* read old patch fraction
  !
  DO jp = 1,io%NPATCH
    ALLOCATE(np%AL(jp)%XPATCH_OLD(np%AL(jp)%NSIZE_P))
  ENDDO
  !
  CALL make_choice_array(hprogram, io%NPATCH, gdim, 'PATCH', zwork)
  DO jp = 1,io%NPATCH
    CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),np%AL(jp)%XPATCH_OLD(:))
  ENDDO
  !
  !* read old soil layer thicknesses (m)
  !
  DO jp = 1,io%NPATCH
    ALLOCATE(np%AL(jp)%XDG_OLD(np%AL(jp)%NSIZE_P,io%NGROUND_LAYER))
  ENDDO
  !
  DO jl=1,io%NGROUND_LAYER
    WRITE(ylvl,'(I4)') jl
    yrecfm='OLD_DG'//adjustl(ylvl(:len_trim(ylvl)))
    CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
    DO jp = 1,io%NPATCH
      CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),np%AL(jp)%XDG_OLD(:,jl))
    ENDDO
  END DO
  DEALLOCATE(zwork)
  !
   CALL init_isba_landuse(dtco, ug, u, io, nk, np, npe, ig%XMESH_SIZE, &
                          hprogram)  
END IF
!
!
!*      12.     Canopy air fields:
!               -----------------
!
 CALL read_sbl_n(dtco, u, sb, io%LCANOPY, hprogram, "NATURE")
!
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!
!        PART 9: initialize radiative and physical properties
!        ----------------------------------------------------
!
DO jp=1,io%NPATCH
  pk => np%AL(jp)
  kk => nk%AL(jp)
  pek => npe%AL(jp)
  !
  ALLOCATE(kk%XDIR_ALB_WITH_SNOW(pk%NSIZE_P,ksw))
  ALLOCATE(kk%XSCA_ALB_WITH_SNOW(pk%NSIZE_P,ksw))
  kk%XDIR_ALB_WITH_SNOW = 0.0
  kk%XSCA_ALB_WITH_SNOW = 0.0
  !
  CALL init_veg_n(io, kk, pk, pek, dti, id%DM%LSURF_DIAG_ALBEDO, pdir_alb, psca_alb, pemis, ptsrad )
  !
  zwg1(1:pk%NSIZE_P)    = pek%XWG(:,1)
  ztg1(1:pk%NSIZE_P,jp) = pek%XTG(:,1)
  !
  CALL convert_patch_isba(dtco, dti, io, idecade, idecade2, s%XCOVER, s%LCOVER,&
                          lagrip, 'NAT', jp, kk, pk, pek, &
                          .false., .false., .false., .false., .true., .false., &
                          pwg1=zwg1(1:pk%NSIZE_P), pwsat=kk%XWSAT)
  !
ENDDO
!
!
! Load randomly perturbed fields. Perturbation ratios are saved in case fields are reset later.
IF(io%LPERTSURF) THEN
  !
  CALL read_surf(hprogram,'VERSION',iversion,iresp)
  CALL read_surf(hprogram,'BUG',ibugfix,iresp)
  gdim = (iversion>8 .OR. iversion==8 .AND. ibugfix>0)
  !
  ALLOCATE(zwork(ki,io%NPATCH))
  !
  CALL make_choice_array(hprogram, io%NPATCH, gdim, 'VEG', zwork)
  ALLOCATE(s%XPERTVEG(ki))
  s%XPERTVEG(:)=zwork(:,1)
!
  CALL make_choice_array(hprogram, io%NPATCH, gdim, 'LAI', zwork)
  ALLOCATE(s%XPERTLAI(ki))
  s%XPERTLAI(:)=zwork(:,1)
!
  CALL make_choice_array(hprogram, io%NPATCH, gdim, 'CV', zwork)
  ALLOCATE(s%XPERTCV(ki))
  s%XPERTCV(:)=zwork(:,1)
!
  CALL make_choice_array(hprogram, io%NPATCH, gdim, 'PERTALB', zwork)
  ALLOCATE(s%XPERTALB(ki))
  s%XPERTALB(:)=zwork(:,1)

  pek => npe%AL(1)
  issk => niss%AL(1)

  WHERE(pek%XALBNIR_VEG (:)/=xundef)  pek%XALBNIR_VEG(:) = pek%XALBNIR_VEG (:) *( 1.+ s%XPERTALB(:) )
  WHERE(pek%XALBVIS_VEG (:)/=xundef)  pek%XALBVIS_VEG(:) = pek%XALBVIS_VEG (:) *( 1.+ s%XPERTALB(:) )
  WHERE(pek%XALBUV_VEG  (:)/=xundef)  pek%XALBUV_VEG (:) = pek%XALBUV_VEG  (:) *( 1.+ s%XPERTALB(:) )
  WHERE(pek%XALBNIR_SOIL(:)/=xundef) pek%XALBNIR_SOIL(:) = pek%XALBNIR_SOIL(:) *( 1.+ s%XPERTALB(:) )
  WHERE(pek%XALBVIS_SOIL(:)/=xundef) pek%XALBVIS_SOIL(:) = pek%XALBVIS_SOIL(:) *( 1.+ s%XPERTALB(:) )
  WHERE(pek%XALBUV_SOIL (:)/=xundef) pek%XALBUV_SOIL (:) = pek%XALBUV_SOIL (:) *( 1.+ s%XPERTALB(:) )
!
  CALL make_choice_array(hprogram, io%NPATCH, gdim, 'PERTZ0LAND', zwork)
  ALLOCATE(s%XPERTZ0(ki))
  s%XPERTZ0(:)=zwork(:,1)
  WHERE(pek%XZ0(:)/=xundef)      pek%XZ0(:)        = pek%XZ0(:)      *( 1.+ s%XPERTZ0(:) )
  WHERE(issk%XZ0EFFIP(:)/=xundef) issk%XZ0EFFIP(:) = issk%XZ0EFFIP(:)*( 1.+ s%XPERTZ0(:) )
  WHERE(issk%XZ0EFFIM(:)/=xundef) issk%XZ0EFFIM(:) = issk%XZ0EFFIM(:)*( 1.+ s%XPERTZ0(:) )
  WHERE(issk%XZ0EFFJP(:)/=xundef) issk%XZ0EFFJP(:) = issk%XZ0EFFJP(:)*( 1.+ s%XPERTZ0(:) )
  WHERE(issk%XZ0EFFJM(:)/=xundef) issk%XZ0EFFJM(:) = issk%XZ0EFFJM(:)*( 1.+ s%XPERTZ0(:) )
!
ENDIF
!
!-------------------------------------------------------------------------------
!
!*       14.    Output radiative fields
!               -----------------------
!
ALLOCATE(s%XEMIS_NAT   (ki))
s%XEMIS_NAT (:) = xundef
!
 CALL averaged_albedo_emis_isba(io, s, nk, np, npe,                           &
                                pzenith, ztg1, psw_bands, pdir_alb, psca_alb, &
                                s%XEMIS_NAT, ztsrad_nat, ztsurf_nat        )  
!
pemis  = s%XEMIS_NAT
ptsrad = ztsrad_nat
ptsurf = ztsurf_nat
!
!-------------------------------------------------------------------------------
!
!*      15.     ISBA diagnostics initialization
!               -------------------------------
!
IF(io%NPATCH<=1) id%O%LPATCH_BUDGET=.false.
!
 CALL diag_isba_init_n(chi, id%DE, id%DEC, id%NDE, id%NDEC, id%O, &
                       id%D, id%DC, id%ND, id%NDC, id%DM, id%NDM, &
                       oread_budgetc, ngb, gb, io, np, npe%AL(1)%TSNOW%SCHEME, &
                       npe%AL(1)%TSNOW%NLAYER, SIZE(s%XABC), hprogram,ki,ksw)
!
!-------------------------------------------------------------------------------
!
 CALL init_surf_topd(id%DEC, io, s, k, np, npe, ug, u, hprogram, u%NDIM_FULL)
!
!-------------------------------------------------------------------------------
!
!         End of IO
!
 CALL end_io_surf_n(hprogram)
!
IF (lhook) CALL dr_hook('COMPUTE_ISBA_PARAMETERS',1,zhook_handle)
!
END SUBROUTINE compute_isba_parameters