Skip to content
Snippets Groups Projects
compute_isba_parameters.F90 41.9 KiB
Newer Older
  • Learn to ignore specific revisions
  • !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, MGN, MSF,  ID,  &
    
                                        GB, NGB, NDST, DST, SLT,BLOWSNW, SV,       &
                                        HPROGRAM,HINIT,                            &
                                        KI,KMONTH,KDAY,KSV,KSW,HSV,PCO2,PRHOA,     &
    
                                        PZENITH,PSW_BANDS,PDIR_ALB,PSCA_ALB,       &
    
                                        PEMIS,PTSRAD,PTSURF, HTEST, PMEGAN_FIELDS  )  
    
    !#############################################################
    !
    !!****  *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 
    
    !!      P. Tulet      06/2016 : call init_megan for coupling megan with surfex
    !!      R. Schoetter  01/2018  add UTCI
    !!      A. Druel      02/2019  Deep change to accept patch duplication (for irrigation or agricultural pratices), mainly for ECOCLIMAP-SG
    
    !!    Séférian/Decharme    08/16 : include fire , river-carbon cycle module + new land-use change implementation
    !!      B. Decharme   02/2022 : split init_diag_isba routine
    
    !!
    !-------------------------------------------------------------------------------
    !
    !*       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_BLOWSNW_n, ONLY : BLOWSNW_t
    
    USE MODD_MEGAN_n, ONLY : MEGAN_t
    USE MODD_MEGAN_SURF_FIELDS_n, ONLY : MEGAN_SURF_FIELDS_t
    !
    
    USE MODD_SFX_OASIS,  ONLY : LCPL_LAND, LCPL_FLOOD, LCPL_GW, LCPL_CALVING, &
                                LCPL_RIVCARB
    
    !
    !
    #ifdef TOPD
    USE MODD_DUMMY_EXP_PROFILE,ONLY : XC_DEPTH_RATIO
    
    USE MODD_ASSIM, ONLY : CASSIM_ISBA, LASSIM
    
    !
    USE MODD_DEEPSOIL,       ONLY : LPHYSDOMC, LDEEPSOIL, XTDEEP_CLI, XGAMMAT_CLI
    
    USE MODD_AGRI,           ONLY : LAGRIP, LIRRIGMODE, NVEG_IRR, 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, LEN_HREC
    
    USE MODD_TOPD_PAR, ONLY : NUNIT
    
    USE MODE_RANDOM
    
    USE MODE_BLOWSNW_SEDIM_LKT1D
    
    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_CONVERT_PATCH_ALB_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_READ_ISBA_NUDGING_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_INIT_SURF_TOPD
    USE MODI_ISBA_SOC_PARAMETERS
    
    USE MODI_PACK_SAME_RANK
    
    USE MODI_DIAG_ISBA_INIT_n
    USE MODI_DIAG_MISC_ISBA_INIT_n
    !
    
    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 MODI_INIT_MEGAN_n
    !
    
    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(MEGAN_t),         INTENT(INOUT) :: MGN
    
    TYPE(MEGAN_SURF_FIELDS_t), INTENT(INOUT) :: MSF
    
    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(DST_t),           INTENT(INOUT) :: DST
    TYPE(SLT_t),           INTENT(INOUT) :: SLT
    TYPE(SV_t),            INTENT(INOUT) :: SV
    
    TYPE(BLOWSNW_t), INTENT(INOUT) :: BLOWSNW
    
    CHARACTER(LEN=6),                 INTENT(IN)  :: HPROGRAM  ! program calling surf. schemes
    CHARACTER(LEN=3),                 INTENT(IN)  :: HINIT     ! choice of fields to initialize
    
    INTEGER,                          INTENT(IN)  :: KI        ! number of points
    
    INTEGER,                          INTENT(IN)  :: KMONTH    ! Month
    INTEGER,                          INTENT(IN)  :: KDAY      ! Day
    
    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)
    
    REAL,             DIMENSION(KI,MSF%NMEGAN_NBR),INTENT(IN), OPTIONAL :: PMEGAN_FIELDS
    
    !
     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
    
    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 :: GFIX      ! update fixed fields
    LOGICAL :: GTIME     ! update time varing fields
    LOGICAL :: GMEB      ! update meb fields
    LOGICAL :: GIRR      ! update irrigation fields
    LOGICAL :: GALBSOIL  ! update bare soil albedo fields
    LOGICAL :: GALBVEG   ! update only vegetation albedo
    LOGICAL :: GURBTREE  ! update urban tree fields
    !
    
    LOGICAL :: GDIM, GCAS1, GCAS2, GCAS3 
    INTEGER :: JVEG, IVERSION, IBUGFIX, IMASK, JMAXLOC
    !
    
    CHARACTER(LEN=4)        :: YLVL
    CHARACTER(LEN=LEN_HREC) :: 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 : 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', U%LECOSG)
    
    !
    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'.OR.IO%CRESPSL=='DIF').AND.IO%CPHOTO=='NCB')THEN
      CALL CARBON_INIT
    ENDIF
    
    !
    !----------------------------------------------------------------------------------
    !----------------------------------------------------------------------------------
    !
    !    PART 2 : 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
    !
    
    ! Take care about irrigation fraction if ECOSG and IRRIGATION
    ALLOCATE(S%XVEGTYPE2(KI,NVEGTYPE+NVEG_IRR))
    S%XVEGTYPE2(:,:) = 0.
    IF (DTI%LDATA_MIXPAR) THEN                                                                 ! Probably sufficient case.  
      ALLOCATE(DTI%XPAR_VEGTYPE2(SIZE(DTI%XPAR_VEGTYPE,1),SIZE(DTI%XPAR_VEGTYPE,2)+NVEG_IRR))  ! because if .NOT.(LIRRIGMODE .OR. LAGRIP),
      DTI%XPAR_VEGTYPE2(:,:) = 0.                                                              ! DTI%LDATA_MIXPAR automaticaly = .TRUE.
    ELSE
      ALLOCATE(DTI%XPAR_VEGTYPE2(1,1))
    ENDIF
    !
    IF ( U%LECOSG .AND. ( LIRRIGMODE .OR. LAGRIP ) ) THEN
      !
      IF ( NVEG_IRR == 0 ) CALL ABOR1_SFX('COMPUTE_ISBA_PARAMETERS - WITH LECOSG AND (LIRRIGMODE OR LAGRIP), NVEG_IRR HAVE TO BE >0')
      !
      IF ( NVEG_IRR /= SIZE(DTI%NPAR_VEG_IRR_USE)) &
        CALL ABOR1_SFX('PCOMPUTE_ISBA_PARAMETER: WHEN LIRRIGMODE or LAGRIP+LECOSG, NVEG_IRR AND NPAR_VEG_IRR_USE HAVE TO BE EQUAL')
      !
      ! IF ECOSG and IRRIGATION, compute fraction of non irrigated vegtype + irrigated vegtype
      DO JVEG=1,NVEGTYPE+NVEG_IRR
        IF ( JVEG <= NVEGTYPE ) THEN
          IF ( ANY(DTI%NPAR_VEG_IRR_USE(:)==JVEG ) ) THEN
            WHERE ( DTI%XPAR_IRRIGFRAC(:,JVEG) /= XUNDEF .AND. DTI%XPAR_IRRIGFRAC(:,JVEG) /= 0. )
              S%XVEGTYPE2(:,JVEG) = S%XVEGTYPE(:,JVEG) * (1 - DTI%XPAR_IRRIGFRAC(:,JVEG) )
            ELSEWHERE
              S%XVEGTYPE2(:,JVEG) = S%XVEGTYPE(:,JVEG)
            ENDWHERE
            IF (DTI%LDATA_MIXPAR) THEN
              WHERE ( DTI%XPAR_IRRIGFRAC(:,JVEG) /= XUNDEF .AND. DTI%XPAR_IRRIGFRAC(:,JVEG) /= 0. )
                DTI%XPAR_VEGTYPE2(:,JVEG) = DTI%XPAR_VEGTYPE(:,JVEG) * (1 - DTI%XPAR_IRRIGFRAC(:,JVEG) )
              ELSEWHERE
                DTI%XPAR_VEGTYPE2(:,JVEG) = DTI%XPAR_VEGTYPE(:,JVEG)
              ENDWHERE
            ENDIF
          ELSE
            S%XVEGTYPE2(:,JVEG) = S%XVEGTYPE(:,JVEG)
            IF ( DTI%LDATA_MIXPAR ) DTI%XPAR_VEGTYPE2(:,JVEG) = DTI%XPAR_VEGTYPE(:,JVEG)
          ENDIF
        ELSE
          WHERE ( DTI%XPAR_IRRIGFRAC(:,DTI%NPAR_VEG_IRR_USE(JVEG-NVEGTYPE)) /= XUNDEF .AND.        &
                  DTI%XPAR_IRRIGFRAC(:,DTI%NPAR_VEG_IRR_USE(JVEG-NVEGTYPE)) /= 0. )
            S%XVEGTYPE2(:,JVEG) = S%XVEGTYPE(:,DTI%NPAR_VEG_IRR_USE(JVEG-NVEGTYPE))                & 
                                  * DTI%XPAR_IRRIGFRAC(:,DTI%NPAR_VEG_IRR_USE(JVEG-NVEGTYPE))
          ELSEWHERE
            S%XVEGTYPE2(:,JVEG) = 0.
          ENDWHERE
          IF (DTI%LDATA_MIXPAR) THEN
            WHERE ( DTI%XPAR_IRRIGFRAC(:,DTI%NPAR_VEG_IRR_USE(JVEG-NVEGTYPE)) /= XUNDEF .AND.      &
                    DTI%XPAR_IRRIGFRAC(:,DTI%NPAR_VEG_IRR_USE(JVEG-NVEGTYPE)) /= 0. )
              DTI%XPAR_VEGTYPE2(:,JVEG) = DTI%XPAR_VEGTYPE(:,DTI%NPAR_VEG_IRR_USE(JVEG-NVEGTYPE))  &
                                  * DTI%XPAR_IRRIGFRAC(:,DTI%NPAR_VEG_IRR_USE(JVEG-NVEGTYPE))
            ELSEWHERE
              DTI%XPAR_VEGTYPE2(:,JVEG) = 0.
            ENDWHERE
          ENDIF
        ENDIF
      ENDDO
      !
    ELSE
      !
      IF ( NVEG_IRR /= 0 ) CALL ABOR1_SFX('COMPUTE_ISBA_PARAMETERS - WITHOUT LECOSG AND (LIRRIGMODE OR LAGRIP), NVEG_IRR HAVE TO BE =0')
      !
      S%XVEGTYPE2(:,:) = S%XVEGTYPE(:,:)
      IF (DTI%LDATA_MIXPAR) DTI%XPAR_VEGTYPE2(:,:) = DTI%XPAR_VEGTYPE(:,:)
      !
    ENDIF
    !
    !
    
    ! patches come from vegtypes
    
    !
    ALLOCATE(S%XPATCH(KI, IO%NPATCH))
    ALLOCATE(S%XVEGTYPE_PATCH(KI, NVEGTYPE+NVEG_IRR, IO%NPATCH))
    CALL SURF_PATCH(IO%NPATCH,S%XVEGTYPE2,DTI%NPAR_VEG_IRR_USE,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 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)
    
      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+NVEG_IRR))
    
      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))
    
      !
      !
      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, U%LECOSG)
    
      !-------------------------------------------------------------------------------
      !
      ! Initialize general parameters
      !
      GFIX     = .TRUE.  ! update fixed fields
      GTIME    = .TRUE.  ! update time varing fields
      GMEB     = .TRUE.  ! update meb fields
      GIRR     = .TRUE.  ! update irrigation fields
      GURBTREE = .FALSE. ! update urban tree fields
      !
    
      CALL CONVERT_PATCH_ISBA(DTCO, DTI, IO, KMONTH, KDAY, IDECADE, IDECADE2, S%XCOVER, S%LCOVER, &
    
                              LAGRIP, U%LECOSG, LIRRIGMODE, 'NAT', JP, KK, PK, PEK, GFIX,         &
                              GTIME, GMEB, GIRR, GURBTREE, PSOILGRID=IO%XSOILGRID, PPERM=KK%XPERM )
      !
      !-------------------------------------------------------------------------------
      !
      ! Initialize vegetation albedo only
      !
      GALBSOIL = .FALSE. ! update bare soil albedo fields
      GALBVEG  = .TRUE.  ! update only vegetation albedo
      !
      CALL CONVERT_PATCH_ALB_ISBA(DTCO, DTI, IO, IDECADE, IDECADE2, S%XCOVER, S%LCOVER, &
                                  'NAT',JP, KK, PK, PEK, GALBSOIL, GALBVEG              )
    
      !
      !-------------------------------------------------------------------------------
      !
      ! 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, &
    
                          .TRUE., LDEEPSOIL, LPHYSDOMC, .TRUE., .TRUE., XTDEEP_CLI, XGAMMAT_CLI,   & 
                          LAGRIP,LIRRIGMODE,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
    
        ALLOCATE(KK%XMUF(0))
    
      !
      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
    
    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%SLTI, CHI%DSTI, &
                          CHI%CCH_NAMES, CHI%CAER_NAMES,  &
    
                          HDSTNAMES=CHI%CDSTNAMES, HSLTNAMES=CHI%CSLTNAMES,          &
                          HSNWNAMES=CHI%CSNWNAMES      )         
    
    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, &
               DTI%XPAR_IRRIGFRAC, DTI%NPAR_VEG_IRR_USE)
    
        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_VEG(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%NSNWEQ >=1) THEN
        ALLOCATE (BLOWSNW%XSNW_FSED(KI,CHI%SVI%NSNWEQ+1))  !Output array
        ALLOCATE (BLOWSNW%XSNW_FTURB(KI,CHI%SVI%NSNWEQ+1))  !Output array
        ALLOCATE (BLOWSNW%XSNW_FNET(KI,CHI%SVI%NSNWEQ+1))  !Output array
        ALLOCATE (BLOWSNW%XSNW_FSALT(KI,CHI%SVI%NSNWEQ+1))  !Output array
        ALLOCATE (BLOWSNW%XSFSNW(KI,CHI%SVI%NSNWEQ+1))  !Output array
        ALLOCATE (BLOWSNW%XSNW_SUBL(KI,CHI%SVI%NSNWEQ+1))  !Output array
        BLOWSNW%XSNW_FSED (:,:) = 0.
        BLOWSNW%XSNW_FTURB(:,:) = 0.
        BLOWSNW%XSNW_FNET (:,:) = 0.
        BLOWSNW%XSNW_FSALT(:,:) = 0.
        BLOWSNW%XSNW_SUBL (:,:) = 0.
        BLOWSNW%XSFSNW    (:,:) = 0.
        !Read in look up tables of snow particles properties
        !No arguments, all look up tables are defined in module
        !mode_snowdrift_sedim_lkt
        CALL BLOWSNW_SEDIM_LKT1D_SET
      ELSE
        ALLOCATE(BLOWSNW%XSNW_FSED(0,0))
        ALLOCATE(BLOWSNW%XSNW_FTURB(0,0))
        ALLOCATE(BLOWSNW%XSNW_FSALT(0,0))
        ALLOCATE(BLOWSNW%XSNW_FNET(0,0))
        ALLOCATE(BLOWSNW%XSNW_SUBL(0,0))
        ALLOCATE(BLOWSNW%XSFSNW(0,0))
      END IF 
    
    
    !-------------------------------------------------------------------------------
    
    !        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')
    
      !
      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 )
      !
    
      !  
      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))
      !                  
    
    !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)
    
          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)
    
          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))
    
        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
        !
    
    ! 6.C. Initialize required coupling fields :
    !-------------------------------------------
    
    IO%LCPL_RRM = .FALSE.
    IO%LFLOOD   = .FALSE.
    IO%LWTD     = .FALSE.
    
    IF(LCPL_LAND)THEN
    !
    
    ! * River coupling
    !
      IO%LCPL_RRM = .TRUE.
    
      ALLOCATE(S%XCPL_DRAIN (KI))
      ALLOCATE(S%XCPL_RUNOFF(KI))
    
      ALLOCATE(S%XCPL_TWS   (KI))
    
      S%XCPL_DRAIN (:) = 0.0
      S%XCPL_RUNOFF(:) = 0.0
    
      S%XCPL_TWS   (:) = 0.0
    !
    ! * Calving coupling
    
    !
      IF(IO%LGLACIER)THEN
         ALLOCATE(S%XCPL_ICEFLUX(KI))
         S%XCPL_ICEFLUX(:) = 0.0
    
         ALLOCATE(S%XCPL_ICEFLUX(0))
    
    !
    ! * Groundwater coupling
    !
      IF(LCPL_GW)THEN
        IO%LWTD = .TRUE.
      ENDIF
    !
    ! * Floodplains coupling
    
    !
      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    
    
        ALLOCATE(S%XCPL_EFLOOD(0))
        ALLOCATE(S%XCPL_PFLOOD(0))
        ALLOCATE(S%XCPL_IFLOOD(0))     
    
    ! * Riverine-carbon cycle coupling
    !
      IF(IO%LCLEACH)THEN
         ALLOCATE(S%XCPL_DOCFLUX(KI))
         S%XCPL_DOCFLUX(:) = 0.0
      ELSE
         ALLOCATE(S%XCPL_DOCFLUX(0))
      ENDIF
    !
    
      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))
    
      ALLOCATE(S%XCPL_DOCFLUX (0))
      ALLOCATE(S%XCPL_TWS     (0))
    
    IF (LCPL_LAND) THEN
    
    !
    ! * Groundwater coupling
    !
    
      ALLOCATE(K%XFWTD(KI))
      ALLOCATE(K%XWTD (KI))
      K%XFWTD(:) = 0.0
      K%XWTD (:) = XUNDEF
    
    !
    ! * Floodplains coupling
    !
    
      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
      !
    
      !
      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
    
    IF(LCPL_RIVCARB)THEN
       IF(.NOT.IO%LCLEACH)THEN
         CALL ABOR1_SFX('COMPUTE_ISBA_PARAMETERS: LCLEACH MUST BE ACTIVATED IF LCPL_RIVCARB')
       ENDIF
    ENDIF
    !
    IF(IO%LCLEACH.AND.IO%CPHOTO/='NCB')THEN
      WRITE(ILUOUT,*)'COMPUTE_ISBA_PARAMETERS: Riverine - Carbon cycle coupling '
      WRITE(ILUOUT,*)'COMPUTE_ISBA_PARAMETERS: Please check your prep namelist where CPHOTO must be NCB '
      CALL ABOR1_SFX('COMPUTE_ISBA_PARAMETERS: Withe LCLEACH, NCB vegetation module is required')
    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)
      !