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, DTV, IG, IO, S, K, ISS, UG, U, USS, HPROGRAM)
!     ##############################################################
!
!!**** *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
!!    B. Decharme    10/2016  bug surface/groundwater coupling
!!
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
#ifdef SFX_OL
USE modn_io_offline,     ONLY : lwr_vegtype
#endif
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_data_isba_n, ONLY : data_isba_t
USE modd_sfx_grid_n, ONLY : grid_t
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_s_t, isba_k_t
USE modd_sso_n, ONLY : sso_t
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_surf_atm_n, ONLY : surf_atm_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, nvt_tebd, nvt_bone, nvt_trbe, &
                                nvt_trbd, nvt_tebe, nvt_tene, nvt_bobd, nvt_bond
USE modd_data_cover,     ONLY : xdata_vegtype
!
USE modd_isba_par,       ONLY : noptimlayer, xoptimgrid
!
USE modi_av_pgd
!
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
!
USE modi_read_namelists_isba
!
#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) :: DTV
TYPE(grid_t), INTENT(INOUT) :: IG
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_s_t), INTENT(INOUT) :: S
TYPE(isba_k_t), INTENT(INOUT) :: K
TYPE(sso_t), INTENT(INOUT) :: ISS
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
TYPE(surf_atm_t), INTENT(INOUT) :: U
TYPE(sso_t), INTENT(INOUT) :: USS
!
 CHARACTER(LEN=6), INTENT(IN)  :: HPROGRAM   ! program calling surf. schemes
!
!
!*    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
INTEGER                  :: JVEGTYPE
 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
 CHARACTER(LEN=4)         :: YALBEDO
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=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
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
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
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(NVEGTYPE)   :: GMEB_PATCH
LOGICAL, DIMENSION(NVEGTYPE)   :: 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, yalbedo, 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,        &                       
                       yrunoffb, yrunoffbfiletype, xunif_runoffb,                &
                       ywdrain,  ywdrainfiletype , xunif_wdrain, zsoilgrid,      &
                       yph, yphfiletype, xunif_ph, yfert, yfertfiletype,         &
                       xunif_fert                          )  
!
io%NPATCH        = ipatch
io%NGROUND_LAYER = iground_layer
io%CISBA         = yisba
io%CPEDOTF       = ypedotf
io%CPHOTO        = yphoto
io%LTR_ML        = gtr_ml
io%CALBEDO       = yalbedo
io%XRM_PATCH     = max(min(zrm_patch,1.),0.)
!
!
 CALL read_namelists_isba(hprogram)
!
!-------------------------------------------------------------------------------
!
!*    2.2      Reading of ISBA MEB namelist
!             -----------------------------
!
IF (io%NPATCH<1 .OR. io%NPATCH>nvegtype) THEN
  WRITE(iluout,*) '*****************************************'
  WRITE(iluout,*) '* Number of patch must be between 1 and ', nvegtype
  WRITE(iluout,*) '* You have chosen NPATCH = ', io%NPATCH
  WRITE(iluout,*) '*****************************************'
  CALL abor1_sfx('PGD_ISBA: NPATCH MUST BE BETWEEN 1 AND NVEGTYPE')
END IF
!
ALLOCATE(io%LMEB_PATCH(io%NPATCH))
!
io%LMEB_PATCH(:) = .false.
io%LFORC_MEASURE = .false.
io%LMEB_LITTER   = .false.
io%LMEB_GNDRES   = .false.

IF(gmeb)THEN

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

  CALL read_nam_pgd_isba_meb(hprogram,iluout,gmeb_patch,io%LFORC_MEASURE,io%LMEB_LITTER,io%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.

  gmeb_patch_rec(:)=.true.
  gmeb_patch_rec(1:3)=.false.
  

  IF(io%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(io%NPATCH>=2 .AND. io%NPATCH<=6)THEN
    gmeb_patch_rec(2)=.true.  ! Only the tree patch (number 2) is allowed to be TRUE
  ELSEIF(io%NPATCH>=7 .AND. io%NPATCH<=8)THEN
    gmeb_patch_rec(3)=.true.  ! Only the tree patch (number 3) is allowed to be TRUE
  ELSEIF(io%NPATCH==9)THEN
    gmeb_patch_rec(3:4)=(/.true.,.true./)  ! Only the tree patches (numbers 3-4) are allowed to be TRUE
  ELSEIF(io%NPATCH==10)THEN
    gmeb_patch_rec(3:5)=(/.true.,.true.,.true./)  ! Only the tree patches (numbers 3-5) are allowed to be TRUE
  ELSEIF(io%NPATCH>=11 .AND. io%NPATCH<=12)THEN
    gmeb_patch_rec(4:6)=(/.true.,.true.,.true./)  ! Only the tree patches (numbers 4-6) are allowed to be TRUE
  ELSEIF(io%NPATCH==nvegtype)THEN
    ! The "old" tree patches (numbers 4-6) are allowed to be TRUE
    gmeb_patch_rec(nvt_tebd) = .true.
    gmeb_patch_rec(nvt_bone) = .true.
    gmeb_patch_rec(nvt_trbe) = .true.
    ! The "new" tree patches (numbers 13-17) are allowed to be TRUE
    gmeb_patch_rec(nvt_trbd) = .true.
    gmeb_patch_rec(nvt_tebe) = .true.
    gmeb_patch_rec(nvt_tene) = .true.
    gmeb_patch_rec(nvt_bobd) = .true.
    gmeb_patch_rec(nvt_bond) = .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:io%NPATCH).AND.gmeb_patch_rec(1:io%NPATCH)
    WRITE(iluout,*) '*****************************************'
  ENDIF
  gmeb_patch(:)=gmeb_patch(:).AND.gmeb_patch_rec(:)

  io%LMEB_PATCH(1:io%NPATCH) = gmeb_patch(1:io%NPATCH)

  IF (io%LMEB_LITTER) io%LMEB_GNDRES = .false.

ENDIF
!
!-------------------------------------------------------------------------------
!
!*    3.      Coherence of options
!             --------------------
!
 CALL test_nam_var_surf(iluout,'CISBA',io%CISBA,'2-L','3-L','DIF')
 CALL test_nam_var_surf(iluout,'CPEDOTF',io%CPEDOTF,'CH78','CO84')
 CALL test_nam_var_surf(iluout,'CPHOTO',io%CPHOTO,'NON','AST','NIT','NCB')
 CALL test_nam_var_surf(iluout,'CALBEDO',io%CALBEDO,'EVOL','DRY ','WET ','MEAN','USER','CM13')
!
SELECT CASE (io%CISBA)
!
  CASE ('2-L')
!          
    io%NGROUND_LAYER = 2
    io%CPEDOTF       ='CH78'
    ALLOCATE(io%XSOILGRID(0)) 
    WRITE(iluout,*) '*****************************************'
    WRITE(iluout,*) '* With option CISBA = ',io%CISBA,'         *'
    WRITE(iluout,*) '* the number of soil layers is set to 2 *'
    WRITE(iluout,*) '* Pedo transfert function = CH78        *'    
    WRITE(iluout,*) '*****************************************'
!    
  CASE ('3-L')
!          
    io%NGROUND_LAYER = 3
    io%CPEDOTF       ='CH78'    
    ALLOCATE(io%XSOILGRID(0))
    WRITE(iluout,*) '*****************************************'
    WRITE(iluout,*) '* With option CISBA = ',io%CISBA,'         *'
    WRITE(iluout,*) '* the number of soil layers is set to 3 *'
    WRITE(iluout,*) '* Pedo transfert function = CH78        *'    
    WRITE(iluout,*) '*****************************************'
!    
  CASE ('DIF')
!          
    IF(io%NGROUND_LAYER==nundef)THEN
      IF(u%LECOCLIMAP)THEN
        io%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(io%XSOILGRID(io%NGROUND_LAYER))
    io%XSOILGRID(:)=0.
    io%XSOILGRID(:)=zsoilgrid(1:io%NGROUND_LAYER) 
    IF (all(zsoilgrid(:)==xundef)) THEN
      IF(u%LECOCLIMAP) io%XSOILGRID(1:io%NGROUND_LAYER)=xoptimgrid(1:io%NGROUND_LAYER)
    ELSEIF (count(io%XSOILGRID/=xundef)/=io%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 (io%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            = ',io%CISBA
    WRITE(iluout,*) '* Pedo transfert function = ',io%CPEDOTF    
    WRITE(iluout,*) '* Number of soil layers   = ',io%NGROUND_LAYER
    IF(u%LECOCLIMAP)THEN
      WRITE(iluout,*) '* Soil layers grid (m)    = ',io%XSOILGRID(1:io%NGROUND_LAYER)
    ENDIF
    WRITE(iluout,*) '*****************************************'
!    
END SELECT
!
SELECT CASE (io%CPHOTO)
  CASE ('AST')
    io%NNBIOMASS = 1
  CASE ('NIT')
    io%NNBIOMASS = 3
  CASE ('NCB')
    io%NNBIOMASS = 6
END SELECT
WRITE(iluout,*) '*****************************************'
WRITE(iluout,*) '* With option CPHOTO = ',io%CPHOTO,'               *'
WRITE(iluout,*) '* the number of biomass pools is set to ', io%NNBIOMASS
WRITE(iluout,*) '*****************************************'
!
IF ( io%CPHOTO/='NON' .AND. io%NPATCH/=12 .AND. io%NPATCH/=nvegtype ) THEN
  WRITE(iluout,*) '*****************************************'
  WRITE(iluout,*) '* With option CPHOTO = ', io%CPHOTO
  WRITE(iluout,*) '* Number of patch must be equal to 12 or NVEGTYPE'
  WRITE(iluout,*) '* But you have chosen NPATCH = ', io%NPATCH
  WRITE(iluout,*) '*****************************************'
  CALL abor1_sfx('PGD_ISBA: CPHOTO='//io%CPHOTO//' REQUIRES NPATCH=12 or NVEGTYPE')
END IF
!
IF ( io%CPHOTO=='NON' .AND. io%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(s%LCOVER      (jpcover))
ALLOCATE(s%XZS         (ilu))
ALLOCATE(ig%XLAT       (ilu))
ALLOCATE(ig%XLON       (ilu))
ALLOCATE(ig%XMESH_SIZE (ilu))
ALLOCATE(iss%XZ0EFFJPDIR(ilu))
!
 CALL pack_pgd(dtco, u, hprogram, 'NATURE',  ig, s%LCOVER, s%XCOVER, s%XZS  )  
!
!-------------------------------------------------------------------------------
!
!*    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%NDIM, iss, u, hprogram,              &
                     zaosip, zaosim, zaosjp, zaosjm,              &
                     zho2ip, zho2im, zho2jp, zho2jm,              &
                     zsso_slope                                   )  
!
!-------------------------------------------------------------------------------
!
!*   15.      ISBA specific fields
!             --------------------
!
io%LECOCLIMAP = u%LECOCLIMAP
!
 CALL pgd_isba_par(dtco, ug, u, uss, dtv, io, s, ig%NDIM, hprogram)
!
!-------------------------------------------------------------------------------
!
#ifdef SFX_OL
IF (lwr_vegtype) THEN
  ALLOCATE(k%XVEGTYPE(ilu,nvegtype))
  IF (dtv%LDATA_VEGTYPE) THEN
    k%XVEGTYPE(:,:) = dtv%XPAR_VEGTYPE(:,:)
  ELSE
    DO jvegtype=1,nvegtype
      CALL av_pgd(dtco,k%XVEGTYPE(:,jvegtype),s%XCOVER,xdata_vegtype(:,jvegtype),'NAT','ARI',s%LCOVER)
    ENDDO
  ENDIF
ENDIF
#endif
!
DEALLOCATE(s%XCOVER)
!
!-------------------------------------------------------------------------------
!
!*    6.      Topographic index for TOPMODEL
!             ------------------------------
!
 CALL pgd_topo_index(dtco, ug, u, uss, s, io%LCTI, &
                     hprogram,ilu,ycti,yctifiletype,limp_cti)
!
!-------------------------------------------------------------------------------
!
!*    7.      Sand fraction
!             -------------
!
CATYPE='ARI'
!
ALLOCATE(k%XSAND(ilu,io%NGROUND_LAYER))
!
 CALL get_field(ysandfiletype,ysand,"SAND",limp_sand,xunif_sand,k%XSAND(:,1))
!
DO jlayer=1,io%NGROUND_LAYER
  k%XSAND(:,jlayer) = k%XSAND(:,1)
END DO
!-------------------------------------------------------------------------------
!
!*    8.      Clay fraction
!             -------------
!
ALLOCATE(k%XCLAY(ilu,io%NGROUND_LAYER))
!
 CALL get_field(yclayfiletype,yclay,"CLAY",limp_clay,xunif_clay,k%XCLAY(:,1))
!
DO jlayer=1,io%NGROUND_LAYER
  k%XCLAY(:,jlayer) = k%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(s%XSOC(ilu,io%NGROUND_LAYER))
!
  io%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
!
 CALL get_field(ysocfiletype,ysoc_top,"SOC_TOP",limp_soc,xunif_soc_top,s%XSOC(:,1))
!
 CALL get_field(ysocfiletype,ysoc_sub,"SOC_SUB",limp_soc,xunif_soc_sub,s%XSOC(:,2))
!
  DO jlayer=2,io%NGROUND_LAYER
    s%XSOC(:,jlayer) = s%XSOC(:,2)
  END DO
!
ELSE
!
  io%LSOCP=.false.
  ALLOCATE(s%XSOC(0,0))
!
ENDIF
!
!*    10.     Permafrost distribution
!             -----------------------
!
IF(len_trim(yperm)/=0.OR.xunif_perm/=xundef)THEN
!
  ALLOCATE(k%XPERM(ilu))
!
  io%LPERM=.true.
! 
 CALL get_field(ypermfiletype,yperm,"PERM",limp_perm,xunif_perm,k%XPERM(:))
!
ELSE
!
  io%LPERM=.false.  
  ALLOCATE(k%XPERM(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(s%XPH(ilu))
  ALLOCATE(s%XFERT(ilu))
  !
  io%LNOF = .true.
  !
  CALL pgd_field(dtco, ug, u, uss, &
                  hprogram,'pH value','NAT',yph,yphfiletype,xunif_ph,s%XPH(:))
  CALL pgd_field(dtco, ug, u, uss, &
                  hprogram,'fertilisation','NAT',yfert,yfertfiletype,xunif_fert,s%XFERT(:))
  !
ENDIF
!
!-------------------------------------------------------------------------------
!
!*    13.      Subgrid runoff 
!             --------------
!
ALLOCATE(k%XRUNOFFB(ilu))
 CALL pgd_field(dtco, ug, u, uss, hprogram,'subgrid runoff','NAT',yrunoffb,yrunoffbfiletype,&
                        xunif_runoffb,k%XRUNOFFB(:))  
!
!-------------------------------------------------------------------------------
!
!*    14.     Drainage coefficient
!             --------------------
!
ALLOCATE(k%XWDRAIN(ilu))
 CALL pgd_field(dtco, ug, u, uss, hprogram,'subgrid drainage','NAT',ywdrain,ywdrainfiletype,&
                        xunif_wdrain,k%XWDRAIN(:))
!
!-------------------------------------------------------------------------------
!
 CALL pgd_topd(io%CISBA, ug%G%CGRID, ug%G%XGRID_PAR, u%NDIM_FULL, uss%XSSO_SLOPE, hprogram)
!
!-------------------------------------------------------------------------------
!
!*   16.     Prints of cover parameters in a tex file
!            ----------------------------------------
!
IF (u%LECOCLIMAP) THEN
  CALL write_cover_tex_isba    (io%NPATCH,io%NGROUND_LAYER,io%CISBA)
  CALL write_cover_tex_isba_par(dtco, io%CALBEDO, io%LTR_ML, &
                                io%NPATCH,io%NGROUND_LAYER,io%CISBA,io%CPHOTO,io%XSOILGRID)
END IF
IF (lhook) CALL dr_hook('PGD_ISBA',1,zhook_handle)
!
CONTAINS
!
SUBROUTINE get_field(HFILETYPE,HFILE,HFIELD,OIMP,PUNIF,PFIELD)
!
IMPLICIT NONE
!
 CHARACTER(LEN=*), INTENT(INOUT) :: HFILETYPE
 CHARACTER(LEN=*), INTENT(IN) :: HFILE
 CHARACTER(LEN=*), INTENT(IN) :: HFIELD
LOGICAL, INTENT(IN) :: OIMP
REAL, INTENT(IN) :: PUNIF
REAL, DIMENSION(:), INTENT(OUT) :: PFIELD
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('PGD_ISBA:GET_FIELD',0,zhook_handle)
!
IF(oimp)THEN
!
  IF(hfiletype=='NETCDF')THEN
     CALL abor1_sfx('Use another format than netcdf for '//trim(hfield)//' input file with LIMP ')
  ELSE
#ifdef SFX_ASC
     cfilein     = adjustl(adjustr(hfile)//'.txt')
#endif
#ifdef SFX_FA
     cfilein_fa  = adjustl(adjustr(hfile)//'.fa')
#endif
#ifdef SFX_LFI
     cfilein_lfi = adjustl(hfile)
#endif
 CALL init_io_surf_n(dtco, u, hfiletype,'NATURE','ISBA  ','READ ')
  ENDIF     
!   
  CALL read_surf(hfiletype,trim(hfield),pfield,iresp) 
!
  CALL end_io_surf_n(hfiletype)
!
ELSE
   CALL pgd_field(dtco, ug, u, uss, &
                  hprogram,hfield,'NAT',hfile,hfiletype,punif,pfield)
ENDIF
!
IF (lhook) CALL dr_hook('PGD_ISBA:GET_FIELD',1,zhook_handle)
!
END SUBROUTINE get_field
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE pgd_isba