pgd_isba.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 pgd_isba (DTCO, DTI, DGU, IG, I, UG, U, USS, &
                           hprogram,oecoclimap)
!     ##############################################################
!
!!**** *PGD_ISBA* monitor for averaging and interpolations of ISBA physiographic fields
!!
!!    PURPOSE
!!    -------
!!
!!    METHOD
!!    ------
!!   
!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!
!!    V. Masson        Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    10/12/97
!!    P. Le Moigne  12/2004 : add type of photosynthesis and correct computation
!!                            of ground layers number in diffusion case
!!    P. Le Moigne  09/2005 : AGS modifs of L. Jarlan
!!    B. Decharme      2008 :  XWDRAIN
!!    E. Martin     12/2008 : files of data for runoffb and wdrain
!!    B. Decharme   06/2009 : files of data for topographic index
!!    A.L. Gibelin  04/2009 : dimension NBIOMASS for ISBA-A-gs
!!    R. Alkama     05/2012 : npatch must be 12 or 19 if CPHOTO/='NON'
!!    B. Decharme   11/2013 : groundwater distribution for water table/surface coupling
!!    P. Samuelsson 02/2012 : MEB
!!
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_data_isba_n, ONLY : data_isba_t
USE modd_diag_surf_atm_n, ONLY : diag_surf_atm_t
USE modd_isba_grid_n, ONLY : isba_grid_t
USE modd_isba_n, ONLY : isba_t
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_surf_atm_n, ONLY : surf_atm_t
USE modd_surf_atm_sso_n, ONLY : surf_atm_sso_t
!
USE modd_surf_par,       ONLY : xundef, nundef
USE modd_pgd_grid,       ONLY : nl
USE modd_pgdwork,        ONLY : catype
USE modd_data_cover_par, ONLY : nvegtype, jpcover
!
USE modd_isba_par,       ONLY : noptimlayer, xoptimgrid
!
USE modi_get_luout
USE modi_read_nam_pgd_isba
USE modi_read_nam_pgd_isba_meb
USE modi_pgd_field
USE modi_test_nam_var_surf
!
USE modi_get_aos_n
USE modi_get_sso_n
USE modi_get_surf_size_n
USE modi_pack_pgd_isba
USE modi_pack_pgd
USE modi_write_cover_tex_isba
USE modi_write_cover_tex_isba_par
USE modi_pgd_topo_index
USE modi_open_namelist
USE modi_close_namelist
USE modi_pgd_isba_par
USE modi_pgd_topd
USE mode_pos_surf
!
USE modi_read_surf
USE modi_init_io_surf_n
USE modi_end_io_surf_n
!
#ifdef SFX_ASC
USE modd_io_surf_asc, ONLY : cfilein
#endif
#ifdef SFX_FA
USE modd_io_surf_fa,  ONLY : cfilein_fa
#endif
#ifdef SFX_LFI
USE modd_io_surf_lfi, ONLY : cfilein_lfi
#endif
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
!
TYPE(data_cover_t), INTENT(INOUT) :: dtco
TYPE(data_isba_t), INTENT(INOUT) :: dti
TYPE(diag_surf_atm_t), INTENT(INOUT) :: dgu
TYPE(isba_grid_t), INTENT(INOUT) :: ig
TYPE(isba_t), INTENT(INOUT) :: i
TYPE(surf_atm_grid_t), INTENT(INOUT) :: ug
TYPE(surf_atm_t), INTENT(INOUT) :: u
TYPE(surf_atm_sso_t), INTENT(INOUT) :: uss
!
 CHARACTER(LEN=6), INTENT(IN)  :: hprogram   ! program calling surf. schemes
LOGICAL,          INTENT(IN)  :: oecoclimap ! T if parameters are computed with ecoclimap
!                                           ! F if all parameters must be specified
!
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
INTEGER                           :: iluout    ! output listing logical unit
INTEGER                           :: jlayer    ! loop counter
INTEGER                           :: ilu       ! number of points
INTEGER                           :: ilunam    ! namelist file logical unit
REAL, DIMENSION(NL)               :: zaosip    ! A/S i+ on all surface points
REAL, DIMENSION(NL)               :: zaosim    ! A/S i- on all surface points
REAL, DIMENSION(NL)               :: zaosjp    ! A/S j+ on all surface points
REAL, DIMENSION(NL)               :: zaosjm    ! A/S j- on all surface points
REAL, DIMENSION(NL)               :: zho2ip    ! h/2 i+ on all surface points
REAL, DIMENSION(NL)               :: zho2im    ! h/2 i- on all surface points
REAL, DIMENSION(NL)               :: zho2jp    ! h/2 j+ on all surface points
REAL, DIMENSION(NL)               :: zho2jm    ! h/2 j- on all surface points
REAL, DIMENSION(NL)               :: zsso_slope! subgrid slope on all surface points
INTEGER                           :: iresp     ! error code
LOGICAL                           :: gmeb      ! Multi-energy balance (MEB)
LOGICAL                           :: gfound    ! flag when namelist is present
!
!*    0.3    Declaration of namelists
!            ------------------------
!
!
INTEGER                  :: ipatch           ! number of patches
INTEGER                  :: iground_layer    ! number of soil layers
 CHARACTER(LEN=3)         :: yisba            ! ISBA option
 CHARACTER(LEN=4)         :: ypedotf          ! Pedo transfert function for DIF
 CHARACTER(LEN=3)         :: yphoto           ! photosynthesis option
LOGICAL                  :: gtr_ml           ! new radiative transfert
REAL                     :: zrm_patch        ! threshold to remove little fractions of patches
 CHARACTER(LEN=28)        :: ysand            ! file name for sand fraction
 CHARACTER(LEN=28)        :: yclay            ! file name for clay fraction
 CHARACTER(LEN=28)        :: ysoc_top         ! file name for organic carbon top soil
 CHARACTER(LEN=28)        :: ysoc_sub         ! file name for organic carbon sub soil
 CHARACTER(LEN=28)        :: ycti             ! file name for topographic index
 CHARACTER(LEN=28)        :: yrunoffb         ! file name for runoffb parameter
 CHARACTER(LEN=28)        :: ywdrain          ! file name for wdrain parameter
 CHARACTER(LEN=28)        :: yperm            ! file name for permafrost distribution
 CHARACTER(LEN=28)        :: ygw               ! file name for groundwater map
 CHARACTER(LEN=6)         :: ysandfiletype    ! sand data file type
 CHARACTER(LEN=6)         :: yclayfiletype    ! clay data file type
 CHARACTER(LEN=6)         :: ysocfiletype     ! organic carbon data file type
 CHARACTER(LEN=6)         :: yctifiletype     ! topographic index data file type
 CHARACTER(LEN=6)         :: yrunoffbfiletype ! subgrid runoff data file type
 CHARACTER(LEN=6)         :: ywdrainfiletype  ! subgrid drainage data file type
 CHARACTER(LEN=6)         :: ypermfiletype    ! permafrost distribution data file type
 CHARACTER(LEN=6)         :: ygwfiletype      ! groundwater distribution data file type
REAL                     :: xunif_sand       ! uniform value of sand fraction  (-)
REAL                     :: xunif_clay       ! uniform value of clay fraction  (-)
REAL                     :: xunif_soc_top    ! uniform value of organic carbon top soil (kg/m2)
REAL                     :: xunif_soc_sub    ! uniform value of organic carbon sub soil (kg/m2)
REAL                     :: xunif_runoffb    ! uniform value of subgrid runoff coefficient
REAL                     :: xunif_wdrain     ! uniform subgrid drainage parameter
REAL                     :: xunif_perm       ! uniform permafrost distribution
REAL                     :: xunif_gw         ! uniform groundwater distribution
LOGICAL                  :: limp_sand        ! Imposed maps of Sand
LOGICAL                  :: limp_clay        ! Imposed maps of Clay
LOGICAL                  :: limp_soc         ! Imposed maps of organic carbon
LOGICAL                  :: limp_cti         ! Imposed maps of topographic index statistics
LOGICAL                  :: limp_perm        ! Imposed maps of permafrost distribution
LOGICAL                  :: limp_gw          ! Imposed maps of groundwater distribution
REAL, DIMENSION(150)     :: zsoilgrid        ! Soil grid reference for DIF
 CHARACTER(LEN=28)        :: yph           ! file name for pH
 CHARACTER(LEN=28)        :: yfert         ! file name for fertilisation rate
 CHARACTER(LEN=6)         :: yphfiletype   ! pH data file type
 CHARACTER(LEN=6)         :: yfertfiletype ! fertilisation data file type
REAL                     :: xunif_ph      ! uniform value of pH
REAL                     :: xunif_fert    ! uniform value of fertilisation rate
LOGICAL, DIMENSION(19)   :: gmeb_patch
LOGICAL, DIMENSION(19)   :: gmeb_patch_rec ! Recommended MEB patch settings
!
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('PGD_ISBA',0,zhook_handle)
!
 CALL get_luout(hprogram,iluout)
!
!-------------------------------------------------------------------------------
!
!*    2.1      Reading of ISBA namelist
!             -------------------------
!
 CALL read_nam_pgd_isba(hprogram, ipatch, iground_layer,                          &
                       yisba,  ypedotf, yphoto, gtr_ml, zrm_patch,               &
                       yclay, yclayfiletype, xunif_clay, limp_clay,              &
                       ysand, ysandfiletype, xunif_sand, limp_sand,              &
                       ysoc_top, ysoc_sub, ysocfiletype, xunif_soc_top,          &
                       xunif_soc_sub, limp_soc, ycti, yctifiletype, limp_cti,    &
                       yperm, ypermfiletype, xunif_perm, limp_perm, gmeb,        &                       
                       ygw, ygwfiletype, xunif_gw, limp_gw,                      &                       
                       yrunoffb, yrunoffbfiletype, xunif_runoffb,                &
                       ywdrain,  ywdrainfiletype , xunif_wdrain, zsoilgrid,      &
                       yph, yphfiletype, xunif_ph, yfert, yfertfiletype,         &
                       xunif_fert                          )  
!
i%NPATCH        = ipatch
i%NGROUND_LAYER = iground_layer
i%CISBA         = yisba
i%CPEDOTF       = ypedotf
i%CPHOTO        = yphoto
i%LTR_ML        = gtr_ml
i%XRM_PATCH     = max(min(zrm_patch,1.),0.)
!
!
!-------------------------------------------------------------------------------
!
!*    2.2      Reading of ISBA MEB namelist
!             -----------------------------
!
IF (i%NPATCH<1 .OR. i%NPATCH>nvegtype) THEN
  WRITE(iluout,*) '*****************************************'
  WRITE(iluout,*) '* Number of patch must be between 1 and ', nvegtype
  WRITE(iluout,*) '* You have chosen NPATCH = ', i%NPATCH
  WRITE(iluout,*) '*****************************************'
  CALL abor1_sfx('PGD_ISBA: NPATCH MUST BE BETWEEN 1 AND NVEGTYPE')
END IF
!
ALLOCATE(i%LMEB_PATCH(i%NPATCH))
!
i%LMEB_PATCH(:) = .false.
i%LFORC_MEASURE = .false.
i%LMEB_LITTER   = .false.
i%LMEB_GNDRES   = .false.

IF(gmeb)THEN

  i%LTR_ML      = .true. ! Always use this SW radiative transfer option with MEB

  CALL read_nam_pgd_isba_meb(hprogram,iluout,gmeb_patch,i%LFORC_MEASURE,i%LMEB_LITTER,i%LMEB_GNDRES)

! Current recommendation is to use MEB for tree patches only.
! Here follows a test in which non-tree patches in LMEB_PATCH are set to FALSE.
! Thus, if you wish to test MEB for non-tree patches you can set 
! GMEB_PATCH_REC(:)=.TRUE.
! in the following line:

  gmeb_patch_rec(:)=.false.

  IF(i%NPATCH==1 .AND. gmeb_patch(1))THEN
    WRITE(iluout,*) '*****************************************'
    WRITE(iluout,*) '* WARNING!'
    WRITE(iluout,*) '* Using MEB for one patch only is not recommended.'
    WRITE(iluout,*) '* LMEB_PATCH(1) has been set to .FALSE.'
    WRITE(iluout,*) '*****************************************'
  ELSEIF(i%NPATCH>=2 .AND. i%NPATCH<=6)THEN
    gmeb_patch_rec(2)=.true.  ! Only the tree patch (number 2) is allowed to be TRUE
  ELSEIF(i%NPATCH>=7 .AND. i%NPATCH<=8)THEN
    gmeb_patch_rec(3)=.true.  ! Only the tree patch (number 3) is allowed to be TRUE
  ELSEIF(i%NPATCH==9)THEN
    gmeb_patch_rec(3:4)=(/.true.,.true./)  ! Only the tree patches (numbers 3-4) are allowed to be TRUE
  ELSEIF(i%NPATCH==10)THEN
    gmeb_patch_rec(3:5)=(/.true.,.true.,.true./)  ! Only the tree patches (numbers 3-5) are allowed to be TRUE
  ELSEIF(i%NPATCH>=11 .AND. i%NPATCH<=12)THEN
    gmeb_patch_rec(4:6)=(/.true.,.true.,.true./)  ! Only the tree patches (numbers 4-6) are allowed to be TRUE
  ELSEIF(i%NPATCH==19)THEN
    gmeb_patch_rec(4:6)=(/.true.,.true.,.true./)  ! The "old" tree patches (numbers 4-6) are allowed to be TRUE
    gmeb_patch_rec(13:17)=(/.true.,.true.,.true.,.true.,.true./)  ! The "new" tree patches (numbers 13-17) are allowed to be TRUE
  ENDIF

  IF(count(.NOT.gmeb_patch_rec(:) .AND. gmeb_patch(:))>0)THEN
    WRITE(iluout,*) '*****************************************'
    WRITE(iluout,*) '* WARNING!'
    WRITE(iluout,*) '* Using MEB for non-tree patches is not yet recommended.'
    WRITE(iluout,*) '* Therefor, LMEB_PATCH for non-tree patches has been set to .FALSE.'
    WRITE(iluout,*) '* The final LMEB_PATCH vector becomes:'
    WRITE(iluout,*) gmeb_patch(1:i%NPATCH).AND.gmeb_patch_rec(1:i%NPATCH)
    WRITE(iluout,*) '*****************************************'
  ENDIF
  gmeb_patch(:)=gmeb_patch(:).AND.gmeb_patch_rec(:)

  i%LMEB_PATCH(1:i%NPATCH) = gmeb_patch(1:i%NPATCH)
!
  IF (i%LMEB_LITTER) THEN
   i%LMEB_GNDRES = .false.
  ENDIF
!
ENDIF
!
!-------------------------------------------------------------------------------
!
!*    3.      Coherence of options
!             --------------------
!
 CALL test_nam_var_surf(iluout,'CISBA',i%CISBA,'2-L','3-L','DIF')
 CALL test_nam_var_surf(iluout,'CPEDOTF',i%CPEDOTF,'CH78','CO84')
 CALL test_nam_var_surf(iluout,'CPHOTO',i%CPHOTO,'NON','AGS','LAI','AST','LST','NIT','NCB')
!
SELECT CASE (i%CISBA)
!
  CASE ('2-L')
!          
    i%NGROUND_LAYER = 2
    i%CPEDOTF       ='CH78'   
    WRITE(iluout,*) '*****************************************'
    WRITE(iluout,*) '* With option CISBA = ',i%CISBA,'         *'
    WRITE(iluout,*) '* the number of soil layers is set to 2 *'
    WRITE(iluout,*) '* Pedo transfert function = CH78        *'    
    WRITE(iluout,*) '*****************************************'
!    
  CASE ('3-L')
!          
    i%NGROUND_LAYER = 3
    i%CPEDOTF       ='CH78'    
    WRITE(iluout,*) '*****************************************'
    WRITE(iluout,*) '* With option CISBA = ',i%CISBA,'         *'
    WRITE(iluout,*) '* the number of soil layers is set to 3 *'
    WRITE(iluout,*) '* Pedo transfert function = CH78        *'    
    WRITE(iluout,*) '*****************************************'
!    
  CASE ('DIF')
!          
    IF(i%NGROUND_LAYER==nundef)THEN
      IF(oecoclimap)THEN
        i%NGROUND_LAYER=noptimlayer
      ELSE
        WRITE(iluout,*) '****************************************'
        WRITE(iluout,*) '* Number of ground layer not specified *'
        WRITE(iluout,*) '****************************************'
        CALL abor1_sfx('PGD_ISBA: NGROUND_LAYER MUST BE DONE IN NAM_ISBA')
      ENDIF
    ENDIF
! 
    ALLOCATE(i%XSOILGRID(i%NGROUND_LAYER))
    i%XSOILGRID(:)=xundef
    i%XSOILGRID(:)=zsoilgrid(1:i%NGROUND_LAYER) 
    IF (all(zsoilgrid(:)==xundef)) THEN
      IF(oecoclimap) i%XSOILGRID(1:i%NGROUND_LAYER)=xoptimgrid(1:i%NGROUND_LAYER)
    ELSEIF (count(i%XSOILGRID/=xundef)/=i%NGROUND_LAYER) THEN
      WRITE(iluout,*) '********************************************************'
      WRITE(iluout,*) '* Soil grid reference values /= number of ground layer *'
      WRITE(iluout,*) '********************************************************'
      CALL abor1_sfx('PGD_ISBA: XSOILGRID must be coherent with NGROUND_LAYER in NAM_ISBA') 
    ELSEIF (i%XSOILGRID(1).GT.0.01) THEN
      CALL abor1_sfx('PGD_ISBA: First layer of XSOILGRID must be lower than 1cm')
    ENDIF
!
    WRITE(iluout,*) '*****************************************'
    WRITE(iluout,*) '* Option CISBA            = ',i%CISBA
    WRITE(iluout,*) '* Pedo transfert function = ',i%CPEDOTF    
    WRITE(iluout,*) '* Number of soil layers   = ',i%NGROUND_LAYER
    IF(oecoclimap)THEN
      WRITE(iluout,*) '* Soil layers grid (m)    = ',i%XSOILGRID(1:i%NGROUND_LAYER)
    ENDIF
    WRITE(iluout,*) '*****************************************'
!    
END SELECT
!
SELECT CASE (i%CPHOTO)
  CASE ('AGS','LAI','AST','LST')
    i%NNBIOMASS = 1
  CASE ('NIT')
    i%NNBIOMASS = 3
  CASE ('NCB')
    i%NNBIOMASS = 6
END SELECT
WRITE(iluout,*) '*****************************************'
WRITE(iluout,*) '* With option CPHOTO = ',i%CPHOTO,'               *'
WRITE(iluout,*) '* the number of biomass pools is set to ', i%NNBIOMASS
WRITE(iluout,*) '*****************************************'
!
IF ( i%CPHOTO/='NON' .AND. i%NPATCH/=12 .AND. i%NPATCH/=19 ) THEN
  WRITE(iluout,*) '*****************************************'
  WRITE(iluout,*) '* With option CPHOTO = ', i%CPHOTO
  WRITE(iluout,*) '* Number of patch must be equal to 12 or 19'
  WRITE(iluout,*) '* But you have chosen NPATCH = ', i%NPATCH
  WRITE(iluout,*) '*****************************************'
  CALL abor1_sfx('PGD_ISBA: CPHOTO='//i%CPHOTO//' REQUIRES NPATCH=12 or 19')
END IF
!
IF ( i%CPHOTO=='NON' .AND. i%LTR_ML .AND. .NOT. gmeb) THEN
  WRITE(iluout,*) '*****************************************'
  WRITE(iluout,*) '* With option CPHOTO == NON      '
  WRITE(iluout,*) '* And With MEB = F               '
  WRITE(iluout,*) '* New radiative transfert TR_ML  '
  WRITE(iluout,*) '* cant be used '
  WRITE(iluout,*) '*****************************************'
  CALL abor1_sfx('PGD_ISBA: WITH CPHOTO= NON LTR_ML MUST BE FALSE')
END IF
!
!-------------------------------------------------------------------------------
!
!*    4.      Number of points and packing of general fields
!             ----------------------------------------------
!
 CALL get_surf_size_n(dtco, u, &
                      'NATURE',ilu)
!
ALLOCATE(i%LCOVER     (jpcover))
ALLOCATE(i%XZS        (ilu))
ALLOCATE(ig%XLAT       (ilu))
ALLOCATE(ig%XLON       (ilu))
ALLOCATE(ig%XMESH_SIZE (ilu))
ALLOCATE(i%XZ0EFFJPDIR(ilu))
!
 CALL pack_pgd(dtco, u, &
               hprogram, 'NATURE',                    &
                ig%CGRID,  ig%XGRID_PAR,                     &
                i%LCOVER, i%XCOVER, i%XZS,                   &
                ig%XLAT, ig%XLON, ig%XMESH_SIZE, i%XZ0EFFJPDIR    )  
!
!-------------------------------------------------------------------------------
!
!*    5.      Packing of ISBA specific fields
!             -------------------------------
!
 CALL get_aos_n(uss, &
                hprogram,nl,zaosip,zaosim,zaosjp,zaosjm,zho2ip,zho2im,zho2jp,zho2jm)
 CALL get_sso_n(uss, &
                hprogram,nl,zsso_slope)
!
 CALL pack_pgd_isba(dtco, ig, i, u, &
                    hprogram,                                    &
                     zaosip, zaosim, zaosjp, zaosjm,              &
                     zho2ip, zho2im, zho2jp, zho2jm,              &
                     zsso_slope                                   )  
!
!-------------------------------------------------------------------------------
!
!*    6.      Topographic index for TOPMODEL
!             ------------------------------
!
 CALL pgd_topo_index(dgu, dtco, ug, u, uss, i, &
                     hprogram,ilu,ycti,yctifiletype,limp_cti)
!
!-------------------------------------------------------------------------------
!
!*    7.      Sand fraction
!             -------------
!
 catype='ARI'
!
ALLOCATE(i%XSAND(ilu,i%NGROUND_LAYER))
!
IF(limp_sand)THEN
!
  IF(ysandfiletype=='NETCDF')THEN
     CALL abor1_sfx('Use another format than netcdf for sand input file with LIMP_SAND')
  ELSE
#ifdef SFX_ASC
     cfilein     = adjustl(adjustr(ysand)//'.txt')
#endif
#ifdef SFX_FA
     cfilein_fa  = adjustl(adjustr(ysand)//'.fa')
#endif
#ifdef SFX_LFI
     cfilein_lfi = adjustl(ysand)
#endif
 CALL init_io_surf_n(dtco, dgu, u, &
                         ysandfiletype,'NATURE','ISBA  ','READ ')
  ENDIF     
!   
  CALL read_surf(&
                 ysandfiletype,'SAND',i%XSAND(:,1),iresp) 
!
  CALL end_io_surf_n(ysandfiletype)
!
ELSE
   CALL pgd_field(dtco, ug, u, uss, &
                  hprogram,'sand fraction','NAT',ysand,ysandfiletype,xunif_sand,i%XSAND(:,1))
ENDIF
!
DO jlayer=1,i%NGROUND_LAYER
  i%XSAND(:,jlayer) = i%XSAND(:,1)
END DO
!-------------------------------------------------------------------------------
!
!*    8.      Clay fraction
!             -------------
!
ALLOCATE(i%XCLAY(ilu,i%NGROUND_LAYER))
!
IF(limp_clay)THEN
!
  IF(yclayfiletype=='NETCDF')THEN
     CALL abor1_sfx('Use another format than netcdf for clay input file with LIMP_CLAY')
  ELSE
#ifdef SFX_ASC
     cfilein     = adjustl(adjustr(yclay)//'.txt')
#endif
#ifdef SFX_FA
     cfilein_fa  = adjustl(adjustr(yclay)//'.fa')
#endif
#ifdef SFX_LFI
     cfilein_lfi = adjustl(yclay)
#endif
 CALL init_io_surf_n(dtco, dgu, u, &
                         yclayfiletype,'NATURE','ISBA  ','READ ')
  ENDIF     
!   
  CALL read_surf(&
                 yclayfiletype,'CLAY',i%XCLAY(:,1),iresp) 
!
  CALL end_io_surf_n(yclayfiletype)
!
ELSE
  CALL pgd_field(dtco, ug, u, uss, &
                  hprogram,'clay fraction','NAT',yclay,yclayfiletype,xunif_clay,i%XCLAY(:,1))
ENDIF
!
DO jlayer=1,i%NGROUND_LAYER
  i%XCLAY(:,jlayer) = i%XCLAY(:,1)
END DO
!
!-------------------------------------------------------------------------------
!
!*    9.      organic carbon profile
!             ----------------------
!
IF(len_trim(ysocfiletype)/=0.OR.(xunif_soc_top/=xundef.AND.xunif_soc_sub/=xundef))THEN
!
  ALLOCATE(i%XSOC(ilu,i%NGROUND_LAYER))
!
  i%LSOCP=.true.
!
  IF((len_trim(ysoc_top)==0.AND.len_trim(ysoc_sub)/=0).OR.(len_trim(ysoc_top)/=0.AND.len_trim(ysoc_sub)==0))THEN
    WRITE(iluout,*) ' '
    WRITE(iluout,*) '***********************************************************'
    WRITE(iluout,*) '* Error in soil organic carbon preparation                *'
    WRITE(iluout,*) '* If used, sub and top soil input file must be given      *'
    WRITE(iluout,*) '***********************************************************'
    WRITE(iluout,*) ' '
    CALL abor1_sfx('PGD_ISBA: TOP AND SUB SOC INPUT FILE REQUIRED')        
  ENDIF
!
  IF(limp_soc)THEN
!
!   Topsoil
!
    IF(ysocfiletype=='NETCDF')THEN
       CALL abor1_sfx('Use another format than netcdf for organic carbon input file with LIMP_SOC')
    ELSE
#ifdef SFX_ASC
       cfilein     = adjustl(adjustr(ysoc_top)//'.txt')
#endif
#ifdef SFX_FA
       cfilein_fa  = adjustl(adjustr(ysoc_top)//'.fa')
#endif
#ifdef SFX_LFI
       cfilein_lfi = adjustl(ysoc_top)
#endif
 CALL init_io_surf_n(dtco, dgu, u, &
                         ysocfiletype,'NATURE','ISBA  ','READ ')
    ENDIF     
!   
    CALL read_surf(&
                 ysocfiletype,'SOC_TOP',i%XSOC(:,1),iresp) 
!
    CALL end_io_surf_n(ysocfiletype)
!
!   Subsoil
!
    IF(ysocfiletype=='NETCDF')THEN
       CALL abor1_sfx('Use another format than netcdf for organic carbon input file with LIMP_SOC')
    ELSE
#ifdef SFX_ASC
       cfilein     = adjustl(adjustr(ysoc_sub)//'.txt')
#endif
#ifdef SFX_FA
       cfilein_fa  = adjustl(adjustr(ysoc_sub)//'.fa')
#endif
#ifdef SFX_LFI
       cfilein_lfi = adjustl(ysoc_sub)
#endif
 CALL init_io_surf_n(dtco, dgu, u, &
                         ysocfiletype,'NATURE','ISBA  ','READ ')
    ENDIF     
!   
    CALL read_surf(&
                 ysocfiletype,'SOC_SUB',i%XSOC(:,2),iresp) 
!
    CALL end_io_surf_n(ysocfiletype)
!
  ELSE
    CALL pgd_field(dtco, ug, u, uss, &
                  hprogram,'organic carbon','NAT',ysoc_top,ysocfiletype,xunif_soc_top,i%XSOC(:,1))
    CALL pgd_field(dtco, ug, u, uss, &
                  hprogram,'organic carbon','NAT',ysoc_sub,ysocfiletype,xunif_soc_sub,i%XSOC(:,2))
  ENDIF
!
  DO jlayer=2,i%NGROUND_LAYER
    i%XSOC(:,jlayer) = i%XSOC(:,2)
  END DO
!
ELSE
!
  i%LSOCP=.false.
  ALLOCATE(i%XSOC(0,0))
!
ENDIF
!
!*    10.     Permafrost distribution
!             -----------------------
!
IF(len_trim(yperm)/=0.OR.xunif_perm/=xundef)THEN
!
  ALLOCATE(i%XPERM(ilu))
!
  i%LPERM=.true.
!
  IF(limp_perm)THEN
!
    IF(ypermfiletype=='NETCDF')THEN
       CALL abor1_sfx('Use another format than netcdf for permafrost input file with LIMP_PERM')
    ELSE
#ifdef SFX_ASC
       cfilein     = adjustl(adjustr(yperm)//'.txt')
#endif
#ifdef SFX_FA
       cfilein_fa  = adjustl(adjustr(yperm)//'.fa')
#endif
#ifdef SFX_LFI
       cfilein_lfi = adjustl(yperm)
#endif
 CALL init_io_surf_n(dtco, dgu, u, &
                         ypermfiletype,'NATURE','ISBA  ','READ ')
    ENDIF     
!   
    CALL read_surf(&
                 ypermfiletype,'PERM',i%XPERM(:),iresp) 
!
    CALL end_io_surf_n(ypermfiletype)
  ELSE
    CALL pgd_field(dtco, ug, u, uss, &
                  hprogram,'permafrost mask','NAT',yperm,ypermfiletype,xunif_perm,i%XPERM(:))
  ENDIF
!
ELSE
!
  i%LPERM=.false.  
  ALLOCATE(i%XPERM(0))
!
ENDIF
!
!*    11.     Groundwater bassin distribution
!             -----------------------
!
IF(len_trim(ygw)/=0.OR.xunif_gw/=xundef)THEN
!
  ALLOCATE(i%XGW(ilu))
!
  i%LGW=.true.
!
  IF(limp_gw)THEN
!
    IF(ygwfiletype=='NETCDF')THEN
       CALL abor1_sfx('Use another format than netcdf for groundwater input file with LIMP_GW')
    ELSE
#ifdef SFX_ASC
       cfilein     = adjustl(adjustr(ygw)//'.txt')
#endif
#ifdef SFX_FA
       cfilein_fa  = adjustl(adjustr(ygw)//'.fa')
#endif
#ifdef SFX_LFI
       cfilein_lfi = adjustl(ygw)
#endif
 CALL init_io_surf_n(dtco, dgu, u, &
                         ygwfiletype,'NATURE','ISBA  ','READ ')
    ENDIF     
!   
    CALL read_surf(&
                 ygwfiletype,'GW',i%XGW(:),iresp) 
!
    CALL end_io_surf_n(ygwfiletype)
  ELSE
    CALL pgd_field(dtco, ug, u, uss, &
                  hprogram,'Groundwater bassin','NAT',ygw,ygwfiletype,xunif_gw,i%XGW(:))
  ENDIF
!
ELSE
!
  i%LGW=.false.  
  ALLOCATE(i%XGW(0))
!
ENDIF
!
!-------------------------------------------------------------------------------
!
!*    12.  pH and fertlisation data
!             --------------------------
!
IF((len_trim(yphfiletype)/=0.OR.xunif_ph/=xundef) .AND. (len_trim(yfertfiletype)/=0.OR.xunif_fert/=xundef)) THEN
  !
  ALLOCATE(i%XPH(ilu))
  ALLOCATE(i%XFERT(ilu))
  !
  i%LNOF = .true.
  !
  CALL pgd_field(dtco, ug, u, uss, &
                  hprogram,'pH value','NAT',yph,yphfiletype,xunif_ph,i%XPH(:))
  CALL pgd_field(dtco, ug, u, uss, &
                  hprogram,'fertilisation','NAT',yfert,yfertfiletype,xunif_fert,i%XFERT(:))
  !
ENDIF
!
!-------------------------------------------------------------------------------
!
!*    13.      Subgrid runoff 
!             --------------
!
ALLOCATE(i%XRUNOFFB(ilu))
 CALL pgd_field(dtco, ug, u, uss, &
                hprogram,'subgrid runoff','NAT',yrunoffb,yrunoffbfiletype,xunif_runoffb,i%XRUNOFFB(:))  
!
!-------------------------------------------------------------------------------
!
!*    14.     Drainage coefficient
!             --------------------
!
ALLOCATE(i%XWDRAIN(ilu))
 CALL pgd_field(dtco, ug, u, uss, &
                hprogram,'subgrid drainage','NAT',ywdrain,ywdrainfiletype,xunif_wdrain,i%XWDRAIN(:))
!
!-------------------------------------------------------------------------------
!
!*   15.      ISBA specific fields
!             --------------------
!
i%LECOCLIMAP = oecoclimap
!
 CALL pgd_isba_par(dtco, dgu, ug, u, uss, dti, i, ig, &
                   hprogram)
!
!-------------------------------------------------------------------------------
!
 CALL pgd_topd(i, ug, u, uss, &
               hprogram)
!
!-------------------------------------------------------------------------------
!
!*   16.     Prints of cover parameters in a tex file
!            ----------------------------------------
!
IF (oecoclimap) THEN
  CALL write_cover_tex_isba(i%NPATCH,i%NGROUND_LAYER,i%CISBA)
  CALL write_cover_tex_isba_par(dtco, i, &
                                i%NPATCH,i%NGROUND_LAYER,i%CISBA,i%CPHOTO,i%XSOILGRID)
END IF
IF (lhook) CALL dr_hook('PGD_ISBA',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE pgd_isba