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, DGU, UG, U, IM, DST, 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
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE modd_surfex_n, ONLY : isba_model_t
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_diag_surf_atm_n, ONLY : diag_surf_atm_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_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_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_read_surf
USE modi_read_isba_n
USE modi_init_isba_landuse
USE modi_read_isba_canopy_n
USE modi_init_veg_n
USE modi_averaged_albedo_emis_isba
USE modi_diag_isba_init_n
USE modi_init_surf_topd
USE modi_isba_soc_parameters
!
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 yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments
!              -------------------------
!
!
TYPE(isba_model_t), INTENT(INOUT) :: im
TYPE(data_cover_t), INTENT(INOUT) :: dtco
TYPE(diag_surf_atm_t), INTENT(INOUT) :: dgu
TYPE(surf_atm_grid_t), INTENT(INOUT) :: ug
TYPE(surf_atm_t), INTENT(INOUT) :: u
TYPE(dst_t), INTENT(INOUT) :: dst
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
!              -------------------------------
!
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,IM%I%NPATCH)  :: zwg1 ! work array for surface water content
REAL, DIMENSION(KI,IM%I%NPATCH)  :: ztg1 ! work array for surface temperature
!
REAL, DIMENSION(KI)         :: zm, zwork
REAL, DIMENSION(KI,IM%I%NPATCH)  :: zf, zpertbuf
!
INTEGER :: idim_full, jl
INTEGER           :: jilu     ! loop increment
INTEGER           :: iluout   ! unit of output listing file
INTEGER           :: iresp   ! return code
INTEGER           :: idecade, idecade2  ! decade of simulation
INTEGER           :: jpatch  ! loop counter on tiles
INTEGER           :: infompi
INTEGER           :: isize_lmeb_patch  ! Number of patches with MEB=true
!
LOGICAL                           :: lwork
!
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
!
zf(:,:) = xundef
zm(:)   = xundef
zwork(:)   = xundef
!
!*       2.3    Physiographic data fields from land cover:
!               -----------------------------------------
!
 CALL allocate_physio(im%I, &
                      im%I%CPHOTO, im%I%CISBA, ki, nvegtype, im%I%NGROUND_LAYER, im%I%NPATCH, &
                     im%I%XVEGTYPE, im%I%XLAI, im%I%XVEG, im%I%XZ0, im%I%XEMIS, im%I%XDG, im%I%XD_ICE,      &
                     im%I%XRSMIN, im%I%XGAMMA, im%I%XWRMAX_CF, im%I%XRGL, im%I%XCV,               &
                     im%I%XZ0_O_Z0H, im%I%XALBNIR_VEG, im%I%XALBVIS_VEG, im%I%XALBUV_VEG,    &
                     im%I%XH_TREE, im%I%XRE25, im%I%XLAIMIN, im%I%XBSLAI, im%I%XSEFOLD,           &
                     im%I%XGMES, im%I%XGC, im%I%XF2I, im%I%XDMAX, im%I%LSTRESS,                   &
                     im%I%XCE_NITRO, im%I%XCF_NITRO, im%I%XCNA_NITRO,                   &
                     im%I%TSEED, im%I%TREAP, im%I%XWATSUP, im%I%XIRRIG,                      &
                     im%I%XROOTFRAC, im%I%NWG_LAYER, im%I%XDROOT, im%I%XDG2,                 &
                     im%I%XGNDLITTER,im%I%XRGLGV,im%I%XGAMMAGV,im%I%XRSMINGV,    &
                     im%I%XROOTFRACGV,im%I%XWRMAX_CFGV,im%I%XLAIGV,im%I%XZ0LITTER,im%I%XH_VEG     )
!
IF (im%I%TTIME%TDATE%MONTH /= nundef) THEN
  idecade = 3 * ( im%I%TTIME%TDATE%MONTH - 1 ) + min(im%I%TTIME%TDATE%DAY-1,29) / 10 + 1
ELSE
  idecade = 1
END IF
!
idecade2 = idecade
!
 CALL init_isba_mixpar(dtco, im%DTI, im%IG, im%I, &
                       im%I%CISBA,idecade,idecade2,im%I%XCOVER,im%I%LCOVER,im%I%CPHOTO,'NAT')
!
isize_lmeb_patch=count(im%I%LMEB_PATCH(:))
IF (isize_lmeb_patch>0)  THEN
  CALL fix_meb_veg(im%DTI, im%IG, im%I, &
                   im%I%NPATCH)
ENDIF
!
 CALL convert_patch_isba(dtco, im%DTI, im%I, &
                         im%I%CISBA,idecade,idecade2,im%I%XCOVER,im%I%LCOVER,im%I%CPHOTO,lagrip,   &
                        im%I%LPERM,im%I%LTR_ML,'NAT',pveg=im%I%XVEG,plai=im%I%XLAI,                &
                        prsmin=im%I%XRSMIN,pgamma=im%I%XGAMMA,pwrmax_cf=im%I%XWRMAX_CF,       &
                        prgl=im%I%XRGL,pcv=im%I%XCV,psoilgrid=im%I%XSOILGRID,                 &
                        pdg=im%I%XDG,kwg_layer=im%I%NWG_LAYER,pdroot=im%I%XDROOT,pdg2=im%I%XDG2,   &
                        pz0=im%I%XZ0,pz0_o_z0h=im%I%XZ0_O_Z0H,                           &
                        palbnir_veg=im%I%XALBNIR_VEG,palbvis_veg=im%I%XALBVIS_VEG,       &
                        palbuv_veg=im%I%XALBUV_VEG,pemis_eco=im%I%XEMIS,                 &
                        pvegtype=im%I%XVEGTYPE,prootfrac=im%I%XROOTFRAC,                 &
                        pgmes=im%I%XGMES,pbslai=im%I%XBSLAI,plaimin=im%I%XLAIMIN,             &
                        psefold=im%I%XSEFOLD,pgc=im%I%XGC,                               &
                        pdmax=im%I%XDMAX,pf2i=im%I%XF2I,ostress=im%I%LSTRESS,ph_tree=im%I%XH_TREE, &
                        pre25=im%I%XRE25,pce_nitro=im%I%XCE_NITRO,pcf_nitro=im%I%XCF_NITRO,   &
                        pcna_nitro=im%I%XCNA_NITRO,pd_ice=im%I%XD_ICE,tpseed=im%I%TSEED,      &
                        tpreap=im%I%TREAP,pwatsup=im%I%XWATSUP,pirrig=im%I%XIRRIG,            &
                        pgndlitter=im%I%XGNDLITTER,                                           &
                        prglgv=im%I%XRGLGV,   &
                        pgammagv=im%I%XGAMMAGV,prsmingv=im%I%XRSMINGV,                   &
                        prootfracgv=im%I%XROOTFRACGV,pwrmax_cfgv=im%I%XWRMAX_CFGV,       &
                        plaigv=im%I%XLAIGV,pz0litter=im%I%XZ0LITTER,ph_veg=im%I%XH_VEG        )
!
!-------------------------------------------------------------------------------
!
 CALL init_veg_pgd_n(im%CHI, dtco, dst, im%I, slt, u, &
                    hprogram, 'NATURE', iluout, ki, im%I%NPATCH, im%I%NGROUND_LAYER,      &
                    im%I%TTIME%TDATE%MONTH,                                          &
                    im%I%XVEGTYPE, im%I%XPATCH, im%I%XVEGTYPE_PATCH, im%I%NSIZE_NATURE_P,           &
                    im%I%NR_NATURE_P, im%I%XRM_PATCH,                                     &
                    ldeepsoil, lphysdomc, xtdeep_cli, xgammat_cli, im%I%XTDEEP,      &
                    im%I%XGAMMAT, lagrip, xthreshold, im%AG%NIRRINUM, im%AG%LIRRIDAY, im%AG%LIRRIGATE, &
                    im%AG%XTHRESHOLDSPT,                                              &
                    im%I%CPHOTO, hinit, im%I%LTR_ML, im%I%NNBIOMASS, pco2, prhoa, im%I%XABC, im%I%XPOI,  &
                    im%I%XGMES, im%I%XGC, im%I%XDMAX, im%I%XANMAX, im%I%XFZERO, im%I%XEPSO, im%I%XGAMM, im%I%XQDGAMM,   & 
                    im%I%XQDGMES, im%I%XT1GMES, im%I%XT2GMES, im%I%XAMAX, im%I%XQDAMAX, im%I%XT1AMAX, im%I%XT2AMAX,&
                    im%I%XAH, im%I%XBH, im%I%XTAU_WOOD, im%I%XINCREASE, im%I%XTURNOVER,                  &
                    ksv, hsv, im%CHI%SVI, im%CHI%CCH_NAMES, im%CHI%CAER_NAMES, &
                    im%CHI%CDSTNAMES, im%CHI%CSLTNAMES, im%CHI%CCHEM_SURF_FILE,                      &
                    dst%XSFDST, dst%XSFDSTM, slt%XSFSLT,                                    &
                    im%I%XAOSIP, im%I%XAOSIM, im%I%XAOSJP, im%I%XAOSJM, im%I%XHO2IP, im%I%XHO2IM, im%I%XHO2JP,     &
                    im%I%XHO2JM, im%I%XZ0, im%I%XZ0EFFIP, im%I%XZ0EFFIM, im%I%XZ0EFFJP, im%I%XZ0EFFJM, im%I%XZ0REL,&
                    im%I%XCLAY, im%I%XSAND, im%I%CPEDOTF,                                      &
                    im%I%XCONDSAT, im%I%XMPOTSAT, im%I%XBCOEF, im%I%XWWILT, im%I%XWFC, im%I%XWSAT, im%I%XWD0,      &
                    im%I%XKANISO, im%I%CRUNOFF,                                           &
                    im%I%XTAUICE, im%I%XCGSAT, im%I%XC1SAT, im%I%XC2REF, im%I%XC3, im%I%XC4B, im%I%XACOEF, im%I%XPCOEF, &
                    im%I%XC4REF, im%I%XPCPS, im%I%XPLVTT, im%I%XPLSTT,                              &
                    im%I%CSCOND, im%I%CISBA, im%I%XHCAPSOIL, im%I%XCONDDRY, im%I%XCONDSLD, im%I%CCPSURF,      &
                    im%I%XDG, im%I%XDROOT, im%I%XDG2, im%I%XROOTFRAC, im%I%XRUNOFFD, im%I%XDZG, im%I%XDZDIF,       &
                    im%I%XSOILWGHT, im%I%NWG_LAYER, im%I%NLAYER_HORT, im%I%NLAYER_DUN, im%I%XD_ICE,      &
                    im%I%XKSAT_ICE, im%I%XALBNIR_DRY, im%I%XALBVIS_DRY, im%I%XALBUV_DRY,            &
                    im%I%XALBNIR_WET, im%I%XALBVIS_WET, im%I%XALBUV_WET, im%I%XBSLAI_NITRO,         &
                    im%I%XCE_NITRO, im%I%XCNA_NITRO, im%I%XCF_NITRO, im%I%XFWTD, im%I%XWTD               )  
!
!-------------------------------------------------------------------------------
!
!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(im%I%CISBA=='DIF') THEN
  !
  IF( im%I%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(im%I%LSOC)THEN   
    IF(.NOT.im%I%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(im%I%XFRACSOC(ki,im%I%NGROUND_LAYER))
    im%I%XFRACSOC(:,:)=0.0
    CALL isba_soc_parameters(im%I%CRUNOFF,im%I%XPATCH,im%I%XDG,im%I%XSOC,im%I%XBCOEF,im%I%XMPOTSAT,   &
                             im%I%XCONDSAT,im%I%XWSAT,im%I%XHCAPSOIL,im%I%XCONDDRY,         &
                             im%I%XCONDSLD,im%I%XWFC,im%I%XWWILT,im%I%XWD0,im%I%XKANISO,im%I%XFRACSOC )
  ELSE
    ALLOCATE(im%I%XFRACSOC(0,0))
  ENDIF
!
ELSE
  ALLOCATE(im%I%XFRACSOC(0,0))
ENDIF
!
!Topmodel
!
!CRUNOFF used in hydro_sgh and isba_sgh_update
IF( im%I%CRUNOFF=='SGH ') THEN 
!
  ALLOCATE(im%I%XTAB_FSAT(ki,ndimtab))
  ALLOCATE(im%I%XTAB_WTOP(ki,ndimtab))
  ALLOCATE(im%I%XTAB_QTOP(ki,ndimtab))
!
  im%I%XTAB_FSAT(:,:) = 0.0
  im%I%XTAB_WTOP(:,:) = 0.0
  im%I%XTAB_QTOP(:,:) = 0.0
!
  IF(hinit/='PRE' .AND. .NOT.lassim)THEN
!
    WHERE(im%I%XCLAY(:,1)==xundef.AND.im%I%XTI_MEAN(:)/=xundef) im%I%XTI_MEAN(:)=xundef
!
    CALL init_top(im%I, &
                   im%I%CISBA, iluout, im%I%XPATCH, im%I%XRUNOFFD,          &
                   im%I%XWD0, im%I%XWSAT, im%I%XTI_MIN,                     &
                   im%I%XTI_MAX, im%I%XTI_MEAN, im%I%XTI_STD, im%I%XTI_SKEW,     &
                   im%I%XSOILWGHT, im%I%XTAB_FSAT, im%I%XTAB_WTOP,          &
                   im%I%XTAB_QTOP, zm                             )  
!
!
    IF (im%I%CKSAT=='SGH' .AND. im%I%CISBA/='DIF') THEN
!     Exponential decay factor calculate using soil properties 
!     (eq. 11, Decharme et al., J. Hydrometeor, 2006)
      DO jilu=1,ki
        IF (zm(jilu)/=xundef) zf(jilu,:) = (im%I%XWSAT(jilu,1)-im%I%XWD0(jilu,1))/zm(jilu)
      ENDDO
!       
    ENDIF
!
  ENDIF
!
! Subsurface flow by layer (m/s)
  IF(im%I%CISBA=='DIF') THEN
    ALLOCATE(im%I%XTOPQS(ki,im%I%NGROUND_LAYER,im%I%NPATCH))
    im%I%XTOPQS(:,:,:)=0.0
  ELSE
    ALLOCATE(im%I%XTOPQS(0,0,0))
  ENDIF
!
ELSE                  
!  
  ALLOCATE(im%I%XTAB_FSAT(0,0))
  ALLOCATE(im%I%XTAB_WTOP(0,0))
  ALLOCATE(im%I%XTAB_QTOP(0,0))
  ALLOCATE(im%I%XTOPQS(0,0,0))  
!                  
ENDIF  
!
!Exponential decay for ISBA-FR option
!CKSAT used in hydro_soil.F90 and soil.F90
IF(hinit/='PRE'.AND.im%I%CISBA/='DIF')THEN 
  !
  IF(im%I%CKSAT=='SGH') THEN
    !
    WHERE(zf(:,:)==xundef.AND.im%I%XDG(:,2,:)/=xundef) 
      zf(:,:) = 4.0/im%I%XDG(:,2,:)
    ENDWHERE
    zf(:,:) = min(zf(:,:),xf_decay)
    !
    ALLOCATE(im%I%XF_PARAM (ki))
    im%I%XF_PARAM(:) = zf(:,1)
    !
    DO jpatch=1,im%I%NPATCH
      IF (im%I%NSIZE_NATURE_P(jpatch) == 0 ) cycle
        CALL exp_decay_soil_fr(im%I%CISBA, zf(:,jpatch),im%I%XC1SAT(:,jpatch),im%I%XC2REF(:,jpatch),   &
                                im%I%XDG(:,:,jpatch),im%I%XD_ICE(:,jpatch),im%I%XC4REF(:,jpatch),      &
                                im%I%XC3(:,:,jpatch),im%I%XCONDSAT(:,:,jpatch),im%I%XKSAT_ICE(:,jpatch))  
    ENDDO                       
    !
  ELSEIF ( im%I%CKSAT=='EXP' .AND. im%I%CISBA=='3-L' ) THEN
    !
    ALLOCATE(im%I%XF_PARAM (ki))
    im%I%XF_PARAM(:) = xundef
    !
    IF (hprogram/='AROME ' .AND. hprogram/='MESONH ') THEN
      !
      CALL open_file('ASCII ',nunit,hfile='carte_f_dc.txt',hform='FORMATTED',haction='READ ')
      DO jilu=1,u%NDIM_FULL
        READ(nunit,*) zf_param(jilu), zc_depth_ratio(jilu)
      ENDDO
      CALL close_file('ASCII ',nunit)
      CALL read_and_send_mpi(zf_param,im%I%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
      !
    ELSE
      WRITE(iluout,*) "COMPUTE_ISBA_PARAMETERS: WITH CKSAT=EXP, IN NOT OFFLINE "//&
                      "MODE, TOPMODEL FILE FOR F_PARAM IS NOT READ "
    ENDIF
    !
    DO jpatch=1,im%I%NPATCH
      WHERE (im%I%XF_PARAM(:)==xundef.AND.im%I%XDG(:,2,jpatch)/=xundef)
        zf(:,jpatch) = 4.0/im%I%XDG(:,2,jpatch)
      ELSEWHERE
        zf(:,jpatch) = im%I%XF_PARAM(:)
      ENDWHERE
    ENDDO
     zf(:,:) = min(zf(:,:),xf_decay)
    !
    DO jpatch=1,im%I%NPATCH
      CALL exp_decay_soil_fr(im%I%CISBA, zf(:,jpatch),im%I%XC1SAT(:,jpatch),im%I%XC2REF(:,jpatch), &
                             im%I%XDG(:,:,jpatch),im%I%XD_ICE(:,jpatch),im%I%XC4REF(:,jpatch),   &
                             im%I%XC3(:,:,jpatch),im%I%XCONDSAT(:,:,jpatch),                &
                             im%I%XKSAT_ICE(:,jpatch))  
    ENDDO    
    !    
  ENDIF
  ! 
ENDIF
!
!
!*       2.10   Soil carbon
!               -----------                        
!
IF (hinit == 'ALL' .AND. im%I%CRESPSL=='CNT' .AND. im%I%CPHOTO == 'NCB') THEN
  CALL carbon_init(im%I%NNBIOMASS, im%I%NNLITTER, im%I%NNLITTLEVS, im%I%NNSOILCARB)
ENDIF
!
!Rainfall spatial distribution
!CRAIN used in HYDRO_VEG and HYDRO_SGH and ISBA_SGH_UPDATE
IF(im%I%CRAIN=='SGH')THEN
  ALLOCATE(im%I%XMUF(ki))
  im%I%XMUF(:)=0.0
ELSE
  ALLOCATE(im%I%XMUF(0))
ENDIF
!
ALLOCATE(im%I%XFSAT(ki))  
im%I%XFSAT(:) = 0.0
!
!-------------------------------------------------------------------------------
!
!*       6.2    Initialize of SFX - RRM coupling:
!               ---------------------------------
!
! * Check some key :
!
IF(lcpl_calving)THEN
   IF(.NOT.im%I%LGLACIER)THEN
     CALL abor1_sfx('COMPUTE_ISBA_PARAMETERS: LGLACIER MUST BE ACTIVATED IF LCPL_CALVING')
   ENDIF
ENDIF
!
! * Initialize required coupling fields :
!
im%I%LCPL_RRM = .false.
im%I%LFLOOD   = .false.
im%I%LWTD     = .false.
!
IF(lcpl_land)THEN
!    
  im%I%LCPL_RRM = .true.
!
  ALLOCATE(im%I%XCPL_DRAIN (ki))
  ALLOCATE(im%I%XCPL_RUNOFF(ki))
  im%I%XCPL_DRAIN (:) = 0.0
  im%I%XCPL_RUNOFF(:) = 0.0
!
  IF(im%I%LGLACIER)THEN
     ALLOCATE(im%I%XCPL_ICEFLUX(ki))
     im%I%XCPL_ICEFLUX(:) = 0.0
  ELSE
     ALLOCATE(im%I%XCPL_ICEFLUX(0))
  ENDIF
!
  IF(lcpl_gw)THEN
    im%I%LWTD = .true.
    ALLOCATE(im%I%XCPL_RECHARGE(ki))
    im%I%XCPL_RECHARGE(:) = 0.0
  ELSE
    ALLOCATE(im%I%XCPL_RECHARGE(0))
  ENDIF
!
  IF(lcpl_flood)THEN
     im%I%LFLOOD = .true.
     ALLOCATE(im%I%XCPL_EFLOOD(ki))
     ALLOCATE(im%I%XCPL_PFLOOD(ki))
     ALLOCATE(im%I%XCPL_IFLOOD(ki))
     im%I%XCPL_EFLOOD(:)= 0.0
     im%I%XCPL_PFLOOD(:)= 0.0
     im%I%XCPL_IFLOOD(:)= 0.0    
  ELSE
    ALLOCATE(im%I%XCPL_EFLOOD(0))
    ALLOCATE(im%I%XCPL_PFLOOD(0))
    ALLOCATE(im%I%XCPL_IFLOOD(0))     
  ENDIF     
!
ELSE
!
  ALLOCATE(im%I%XCPL_RUNOFF  (0))
  ALLOCATE(im%I%XCPL_DRAIN   (0))
  ALLOCATE(im%I%XCPL_ICEFLUX (0))
  ALLOCATE(im%I%XCPL_RECHARGE(0))
  ALLOCATE(im%I%XCPL_EFLOOD  (0))
  ALLOCATE(im%I%XCPL_PFLOOD  (0))
  ALLOCATE(im%I%XCPL_IFLOOD  (0))
!
ENDIF
!
IF(im%I%LWTD.AND..NOT.im%I%LGW)THEN
  WRITE(iluout,*)'COMPUTE_ISBA_PARAMETERS: Groundwater map is required by SFX - Groundwater coupling '
  WRITE(iluout,*)'COMPUTE_ISBA_PARAMETERS: Please check your pgd namelist where this map must be     '
  WRITE(iluout,*)'COMPUTE_ISBA_PARAMETERS: specified (YGW and YGWFILETYPE, or XUNIF_GW, or LIMP_GW)  '
  CALL abor1_sfx('COMPUTE_ISBA_PARAMETERS: Groundwater map is required by SFX - Groundwater coupling')
ENDIF
!
! * Initialize flood scheme :
!
IF(im%I%LFLOOD)THEN
  ALLOCATE(im%I%XFFLOOD (ki))
  ALLOCATE(im%I%XPIFLOOD(ki))
  ALLOCATE(im%I%XFF     (ki,im%I%NPATCH))
  ALLOCATE(im%I%XFFG    (ki,im%I%NPATCH))
  ALLOCATE(im%I%XFFV    (ki,im%I%NPATCH))
  ALLOCATE(im%I%XFFROZEN(ki,im%I%NPATCH))
  ALLOCATE(im%I%XALBF   (ki,im%I%NPATCH))
  ALLOCATE(im%I%XEMISF  (ki,im%I%NPATCH)) 
  im%I%XFFLOOD       = 0.0
  im%I%XPIFLOOD      = 0.0
  im%I%XFF           = 0.0
  im%I%XFFG          = 0.0
  im%I%XFFV          = 0.0
  im%I%XFFROZEN      = 0.0
  im%I%XALBF         = 0.0
  im%I%XEMISF        = 0.0  
ELSE
  ALLOCATE(im%I%XFFLOOD   (0))
  ALLOCATE(im%I%XPIFLOOD  (0))
  ALLOCATE(im%I%XFF     (0,0))
  ALLOCATE(im%I%XFFG    (0,0))
  ALLOCATE(im%I%XFFV    (0,0))
  ALLOCATE(im%I%XFFROZEN(0,0))
  ALLOCATE(im%I%XALBF   (0,0))
  ALLOCATE(im%I%XEMISF  (0,0))
ENDIF
!
!-------------------------------------------------------------------------------
!
!*      7.     ISBA time-varying deep force-restore temperature initialization
!              ---------------------------------------------------------------
!
 CALL soiltemp_arp_par(im%I, &
                       hprogram,im%I%LTEMP_ARP,im%I%NTEMPLAYER_ARP)
!
!-------------------------------------------------------------------------------
!
!*       9.     Prints of cover parameters in a tex file
!               ----------------------------------------
!
!* 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
!
!-------------------------------------------------------------------------------
!
IF (cassim_isba=="ENKF ") THEN
  !
  CALL init_random_seed()
  !
ENDIF
!
 CALL init_io_surf_n(dtco, dgu, u, &
                    hprogram,'NATURE','ISBA  ','READ ')
!
!*      10.     Prognostic and semi-prognostic fields
!               -------------------------------------
!
 CALL read_isba_n(dtco, im%I, 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. im%I%TSNOW%SCHEME.NE.'3-L' .AND. im%I%TSNOW%SCHEME.NE.'CRO' .AND. im%I%CISBA=='DIF') THEN
    CALL abor1_sfx("INIT_ISBAN: WITH CISBA = DIF, CSNOW MUST BE 3-L OR CRO")
ENDIF
!
!-------------------------------------------------------------------------------
!
!*      11.  Extrapolation of the prognostic and semi-prognostic fields
!                           LAND USE case 
!               -------------------------------------
!
IF (oland_use) THEN
   CALL init_isba_landuse(dtco, im%IG, im%I, ug, u, &
                          hprogram)  
END IF
!
!-------------------------------------------------------------------------------
!
!*      12.     Canopy air fields:
!               -----------------
!
 CALL read_isba_canopy_n(dtco, im%ICP, im%I, u, &
                         hprogram)
!
!-------------------------------------------------------------------------------
!
!*      13.     initialize radiative and physical properties
!               --------------------------------------------
!
ALLOCATE(im%I%XDIR_ALB_WITH_SNOW(ki,ksw,im%I%NPATCH))
ALLOCATE(im%I%XSCA_ALB_WITH_SNOW(ki,ksw,im%I%NPATCH))
im%I%XDIR_ALB_WITH_SNOW = 0.0
im%I%XSCA_ALB_WITH_SNOW = 0.0
!
 CALL init_veg_n(im%I%NPATCH, ki, im%I%LCANOPY, im%I%CROUGH, im%I%LAGRI_TO_GRASS, im%I%TSNOW, &
                 im%I%CPHOTO, im%DTI%LIMP_VEG, im%DTI%LIMP_Z0, im%DTI%LIMP_EMIS, &
                 im%I%XLAIMIN, im%I%XH_TREE, im%I%XVEGTYPE_PATCH, im%I%XLAI, &
                 im%I%XZ0, im%I%XVEG, im%I%XEMIS, &
                 im%I%LTR_ML, im%I%XFAPARC, im%I%XFAPIRC, im%I%XLAI_EFFC, im%I%XMUS, &
                 im%I%XALBNIR_SOIL, im%I%XALBVIS_SOIL, im%I%XALBUV_SOIL, im%I%XALBNIR, &
                 im%I%XALBVIS, im%I%XALBUV, &
                 im%DGMI%LSURF_DIAG_ALBEDO, im%I%XPSN, im%I%XPSNG, im%I%XPSNV, im%I%XPSNV_A, &
                 pdir_alb, psca_alb, pemis, ptsrad )
!
DO jpatch=1,im%I%NPATCH
  zwg1(:,jpatch) = im%I%XWG(:,1,jpatch)
  ztg1(:,jpatch) = im%I%XTG(:,1,jpatch)
END DO
!
 CALL convert_patch_isba(dtco, im%DTI, im%I, &
                         im%I%CISBA,idecade,idecade2,im%I%XCOVER,im%I%LCOVER,&
                          im%I%CPHOTO,lagrip,im%I%LPERM,im%I%LTR_ML,'NAT',   &
                          pwg1 = zwg1,               &
                          palbnir_soil=im%I%XALBNIR_SOIL, &
                          palbvis_soil=im%I%XALBVIS_SOIL, &
                          palbuv_soil=im%I%XALBUV_SOIL )
!
! Load randomly perturbed fields. Perturbation ratios are saved in case fields are reset later.
IF(im%I%LPERTSURF) THEN
!
  CALL read_surf(&
                 hprogram,'VEG',im%I%XVEG(:,:),iresp)
  ALLOCATE(im%I%XPERTVEG(ki))
  im%I%XPERTVEG(:)=im%I%XVEG(:,1)
!
  CALL read_surf(&
                 hprogram,'LAI',im%I%XLAI(:,:),iresp)
  ALLOCATE(im%I%XPERTLAI(ki))
  im%I%XPERTLAI(:)=im%I%XLAI(:,1)
!
  CALL read_surf(&
                 hprogram,'CV',im%I%XCV(:,:),iresp)
  ALLOCATE(im%I%XPERTCV(ki))
  im%I%XPERTCV(:)=im%I%XCV(:,1)
!
  CALL read_surf(&
                 hprogram,'PERTALB',zpertbuf(:,:),iresp)
  ALLOCATE(im%I%XPERTALB(ki))
  im%I%XPERTALB(:)=zpertbuf(:,1)
  WHERE(im%I%XALBNIR_VEG(:,1)/=xundef)  im%I%XALBNIR_VEG(:,1) = im%I%XALBNIR_VEG(:,1) *( 1.+ im%I%XPERTALB(:) )
  WHERE(im%I%XALBVIS_VEG(:,1)/=xundef)  im%I%XALBVIS_VEG(:,1) = im%I%XALBVIS_VEG(:,1) *( 1.+ im%I%XPERTALB(:) )
  WHERE(im%I%XALBUV_VEG(:,1)/=xundef)   im%I%XALBUV_VEG(:,1)  = im%I%XALBUV_VEG(:,1)  *( 1.+ im%I%XPERTALB(:) )
  WHERE(im%I%XALBNIR_SOIL(:,1)/=xundef) im%I%XALBNIR_SOIL(:,1)= im%I%XALBNIR_SOIL(:,1)*( 1.+ im%I%XPERTALB(:) )
  WHERE(im%I%XALBVIS_SOIL(:,1)/=xundef) im%I%XALBVIS_SOIL(:,1)= im%I%XALBVIS_SOIL(:,1)*( 1.+ im%I%XPERTALB(:) )
  WHERE(im%I%XALBUV_SOIL(:,1)/=xundef)  im%I%XALBUV_SOIL(:,1) = im%I%XALBUV_SOIL(:,1) *( 1.+ im%I%XPERTALB(:) )
!
  CALL read_surf(&
                 hprogram,'PERTZ0LAND',zpertbuf(:,:),iresp)
  ALLOCATE(im%I%XPERTZ0(ki))
  im%I%XPERTZ0(:)=zpertbuf(:,1)
  WHERE(im%I%XZ0(:,1)/=xundef)      im%I%XZ0(:,1)     =im%I%XZ0(:,1)     *( 1.+ im%I%XPERTZ0(:) )
  WHERE(im%I%XZ0EFFIP(:,1)/=xundef) im%I%XZ0EFFIP(:,1)=im%I%XZ0EFFIP(:,1)*( 1.+ im%I%XPERTZ0(:) )
  WHERE(im%I%XZ0EFFIM(:,1)/=xundef) im%I%XZ0EFFIM(:,1)=im%I%XZ0EFFIM(:,1)*( 1.+ im%I%XPERTZ0(:) )
  WHERE(im%I%XZ0EFFJP(:,1)/=xundef) im%I%XZ0EFFJP(:,1)=im%I%XZ0EFFJP(:,1)*( 1.+ im%I%XPERTZ0(:) )
  WHERE(im%I%XZ0EFFJM(:,1)/=xundef) im%I%XZ0EFFJM(:,1)=im%I%XZ0EFFJM(:,1)*( 1.+ im%I%XPERTZ0(:) )
!
ENDIF
!
!-------------------------------------------------------------------------------
!
!*       14.    Output radiative fields
!               -----------------------
!
ALLOCATE(im%I%XEMIS_NAT   (ki))
im%I%XEMIS_NAT (:) = xundef
!
 CALL averaged_albedo_emis_isba(im%I, &
                                im%I%LFLOOD, im%I%CALBEDO, pzenith,                &
                                 im%I%XVEG,im%I%XZ0,im%I%XLAI,                          &
                                 im%I%LMEB_PATCH,im%I%XGNDLITTER,im%I%XZ0LITTER,im%I%XLAIGV, &
                                 im%I%XH_VEG, im%I%XTV,               &
                                 ztg1,                                   &
                                 im%I%XPATCH,                                 &
                                 psw_bands,                              &
                                 im%I%XALBNIR_VEG,im%I%XALBVIS_VEG,im%I%XALBUV_VEG,     &
                                 im%I%XALBNIR_SOIL,im%I%XALBVIS_SOIL,im%I%XALBUV_SOIL,  &
                                 im%I%XEMIS,                                  &
                                 im%I%TSNOW,                                  &
                                 im%I%XALBNIR,im%I%XALBVIS,im%I%XALBUV,                 &
                                 pdir_alb, psca_alb,                     &
                                 im%I%XEMIS_NAT,ztsrad_nat,ztsurf_nat         )  
!
pemis  = im%I%XEMIS_NAT
ptsrad = ztsrad_nat
ptsurf = ztsurf_nat
!
!-------------------------------------------------------------------------------
!
!*      15.     ISBA diagnostics initialization
!               -------------------------------
!
IF(im%I%NPATCH<=1) im%DGI%LPATCH_BUDGET=.false.
!
 CALL diag_isba_init_n(&
                       im%CHI, im%DGEI, im%DGI, im%DGMI, dgu, im%GB, im%I, &
                       hprogram,ki,ksw)
!
!-------------------------------------------------------------------------------
!
 CALL init_surf_topd(im%DGEI, im%I, 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