mode_read_extern.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.
!     #####################
MODULE mode_read_extern
!     #####################
!-------------------------------------------------------------------
!
USE modd_surf_par,       ONLY : nundef, xundef
USE modd_data_cover_par, ONLY : ncover, ntype, nvegtype, jpcover, nvegtype_old, nvegtype_ecosg
!
USE mode_read_surf_cov, ONLY : read_surf_cov
!
USE modi_read_lecoclimap
!
USE modi_old_name
USE modi_open_aux_io_surf
USE modi_close_aux_io_surf
USE modi_read_surf
USE modi_make_choice_array
USE modi_abor1_sfx
USE modi_read_pgd_cover_garden
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
CONTAINS
!
!---------------------------------------------------------------------------------------
!
!     #######################
      SUBROUTINE read_extern_depth (U, DTCO, GCP, IO, &
                                    HFILE,HFILETYPE,HFILEPGD,HFILEPGDTYPE,&
                                    KLUOUT,HISBA,HNAT,HFIELD,KNI,KLAYER, &
                                   KPATCH,PSOILGRID,PDEPTH,KVERSION,KWG_LAYER          )
!     #######################
!
USE modd_surf_atm_n, ONLY : surf_atm_t
USE modd_grid_conf_proj_n, ONLY : grid_conf_proj_t
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_isba_options_n, ONLY : isba_options_t
!
USE modd_surfex_mpi, ONLY : nrank,npio
!
USE modd_prep_teb_greenroof,   ONLY : ngrid_level, xgrid_soil
!
USE modi_open_aux_io_surf
USE modi_close_aux_io_surf
USE modi_read_surf_isba_par_n
USE modi_convert_cover_isba
USE modi_garden_soil_depth
USE modi_read_arrange_cover
!
! Modifications :
! P.Marguinaud : 11-09-2012 : shorten field name
! G.Delautier : 24-06-2015 : bug for arome compressed files
!
IMPLICIT NONE
!
!* dummy arguments
!  ---------------
!
TYPE(surf_atm_t), INTENT(INOUT) :: U
TYPE(grid_conf_proj_t),INTENT(INOUT) :: GCP
!
TYPE(data_cover_t), INTENT(INOUT) :: DTCO
TYPE(isba_options_t), INTENT(INOUT) :: IO
!
 CHARACTER(LEN=28),    INTENT(IN)    :: HFILE  ! type of input file
 CHARACTER(LEN=6),     INTENT(IN)    :: HFILETYPE  ! type of input file
 CHARACTER(LEN=28),    INTENT(IN)    :: HFILEPGD  ! type of input file
 CHARACTER(LEN=6),     INTENT(IN)    :: HFILEPGDTYPE  ! type of input file
INTEGER,              INTENT(IN)    :: KLUOUT    ! logical unit of output listing
 CHARACTER(LEN=3),     INTENT(IN)    :: HISBA     ! type of ISBA soil scheme
 CHARACTER(LEN=3),     INTENT(IN)    :: HNAT      ! type of surface (nature, gardens)
 CHARACTER(LEN=7),     INTENT(IN)    :: HFIELD    ! field name
INTEGER,              INTENT(IN)    :: KNI       ! number of points
INTEGER,              INTENT(IN)    :: KLAYER    ! number of layers
INTEGER,              INTENT(IN)    :: KPATCH    ! number of patch
INTEGER,              INTENT(IN)    :: KVERSION  ! surface version
REAL, DIMENSION(:),   INTENT(IN)    :: PSOILGRID !
REAL, DIMENSION(:,:,:), POINTER     :: PDEPTH    ! depth of each layer over each patches
INTEGER, DIMENSION(:,:), INTENT(OUT):: KWG_LAYER
!
!* local variables
!  ---------------
!
 CHARACTER(LEN=4 ) :: YLVL
 CHARACTER(LEN=12) :: YRECFM         ! Name of the article to be read
 CHARACTER(LEN=16) :: YRECFM2
 CHARACTER(LEN=100):: YCOMMENT       ! Comment string
 CHARACTER(LEN=6) :: YSURF
INTEGER           :: IRESP          ! reading return code
INTEGER           :: JL         ! loop counter
INTEGER           :: JP         ! loop counter
INTEGER           :: JJ, JI, IEND
INTEGER           :: IVERSION
INTEGER           :: IBUGFIX
!
LOGICAL, DIMENSION(:), ALLOCATABLE   :: GCOVER ! flag to read the covers
REAL,    DIMENSION(:,:), ALLOCATABLE :: ZCOVER ! cover fractions
REAL,    DIMENSION(:,:), ALLOCATABLE :: ZGROUND_DEPTH ! cover fractions
REAL, DIMENSION(:), ALLOCATABLE      :: ZSOILGRID
REAL,  DIMENSION(KNI)                :: ZHVEG  ! high vegetation fraction
REAL,  DIMENSION(KNI)                :: ZLVEG  ! low  vegetation fraction
REAL,  DIMENSION(KNI)                :: ZNVEG  ! no   vegetation fraction
REAL,  DIMENSION(KNI)                :: ZPERM  ! permafrost distribution
 CHARACTER(LEN=4)                    :: YHVEG  ! type of high vegetation
 CHARACTER(LEN=4)                    :: YLVEG  ! type of low  vegetation
 CHARACTER(LEN=4)                    :: YNVEG  ! type of no   vegetation
INTEGER                              :: INVEGTYPE_SAVE, IJPCOVER_SAVE
LOGICAL                              :: GECOCLIMAP ! T if ecoclimap is used
LOGICAL                              :: GECOSG
LOGICAL                              :: GPAR_GARDEN! T if garden data are used
LOGICAL, DIMENSION(NVEGTYPE_ECOSG)   :: GDATA_DG
LOGICAL, DIMENSION(NVEGTYPE_ECOSG)   :: GDATA_GROUND_DEPTH, GDATA_ROOT_DEPTH
LOGICAL                              :: GPERM
LOGICAL                              :: GREAD_EXT
LOGICAL      :: GREAD_OK, GDIM, GDIM2, GWATER_TO_NATURE, GTOWN_TO_ROCK, GGARDEN
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!
!------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('MODE_READ_EXTERN:READ_EXTERN_DEPTH',0,zhook_handle)
!
 CALL open_aux_io_surf(hfilepgd,hfilepgdtype,'FULL  ')

 yrecfm='VERSION'
 CALL read_surf(hfilepgdtype,yrecfm,iversion,iresp,hdir='-')
!
yrecfm='BUG'
 CALL read_surf(hfilepgdtype,yrecfm,ibugfix,iresp,hdir='-')
!
gdim = (iversion>8 .OR. iversion==8 .AND. ibugfix>0)
gdim2 = gdim
IF (gdim) CALL read_surf(hfilepgdtype,'SPLIT_PATCH',gdim2,iresp) 
!
 CALL read_lecoclimap(hfilepgdtype,gecoclimap,gecosg,hdir='-')
IF (hnat=='NAT') THEN  
  ysurf = "NATURE"
  CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
ELSE
  CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
  CALL open_aux_io_surf(hfilepgd,hfilepgdtype,'TOWN  ')
  CALL read_surf(hfilepgdtype,'PAR_GARDEN',gpar_garden,iresp,hdir='-')
  gecoclimap = .NOT. gpar_garden
  IF (.NOT.gecoclimap) gecosg = .false.
  ysurf = "TOWN  "
  CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
END IF
!
invegtype_save = nvegtype
ijpcover_save = jpcover
IF (gecosg) THEN
  nvegtype = nvegtype_ecosg
  jpcover = sum(ntype)
ELSE
  nvegtype = nvegtype_old
  jpcover = ncover
ENDIF
!
!------------------------------------------------------------------------------
!
ALLOCATE(pdepth(kni,klayer,kpatch))
pdepth(:,:,:) = xundef
!
kwg_layer(:,:) = nundef
!
 CALL open_aux_io_surf(hfilepgd,hfilepgdtype,ysurf)
!
!-------------------------------------------------------------------
!
gread_ok = .false.
!
!
iend = 1
IF (gdim) iend = nvegtype_ecosg
!
IF (hnat=='NAT' .AND. (iversion>=7 .OR. .NOT.gecoclimap)) THEN
  !
  !* directly read soil layers in the file for nature ISBA soil layers
  !
  gdata_dg(:) = .false.
  IF (.NOT.gecoclimap) gdata_dg(:) = .true.
  IF (iversion>=7) THEN
    yrecfm='L_DG'
    ycomment=yrecfm
    CALL read_surf(hfilepgdtype,yrecfm,gdata_dg(1:iend),iresp,hcomment=ycomment,hdir='-')
  ENDIF
  !
  IF (any(gdata_dg(1:iend))) THEN
    !
    DO jl=1,klayer
      IF (gdim) THEN
        WRITE(yrecfm,fmt='(A6,I2.2)') 'D_DG_L',jl
      ELSE
        IF (jl<10)  WRITE(yrecfm,fmt='(A4,I1.1)') 'D_DG',jl
        IF (jl>=10) WRITE(yrecfm,fmt='(A4,I2.2)') 'D_DG',jl
      ENDIF
      CALL read_surf_isba_par_n(dtco, u, gcp, kpatch, hfilepgdtype, yrecfm, kluout, kni, &
                                iversion, ibugfix, gdata_dg, pdepth(:,jl,:),iresp,hdir='E')
    END DO
    gread_ok = .true.
    !
  ENDIF
  !
  IF (iversion>7 .OR. iversion==7 .AND. ibugfix>=2) THEN
    !
    !cas when root_depth and ground_depth were extrapolated in extrapol_field
    !during pgd step
    IF (.NOT.any(gdata_dg(1:iend)) .AND. hisba=="3-L") THEN
      !
      yrecfm2='L_ROOT_DEPTH'
      ycomment=yrecfm2
      CALL read_surf(hfilepgdtype,yrecfm2,gdata_root_depth(1:iend),iresp,hcomment=ycomment,hdir='-')
      !
      IF (any(gdata_root_depth(1:iend))) THEN
        IF (gdim) THEN
          yrecfm2='D_RTDPT_'
        ELSE
          yrecfm2='D_ROOT_DEPTH'
        ENDIF                   
        CALL read_surf_isba_par_n(dtco, u, gcp, kpatch, hfilepgdtype, yrecfm2, kluout, kni, &
                                  iversion, ibugfix, gdata_root_depth, pdepth(:,2,:),iresp,hdir='E')
      ENDIF
      !
    ENDIF
    !    
    yrecfm2='L_GROUND_DEPTH'
    IF (iversion>7 .OR. iversion==7 .AND. ibugfix>=3) yrecfm2='L_GROUND_DPT'
    ycomment=yrecfm2
    CALL read_surf(hfilepgdtype,yrecfm2,gdata_ground_depth(1:iend),iresp,hcomment=ycomment,hdir='-')
    !
    IF (any(gdata_ground_depth(1:iend))) THEN
      !
      ALLOCATE(zground_depth(kni,kpatch))
      IF (gdim) THEN
        yrecfm2='D_GRDPT_'
      ELSE
        yrecfm2='D_GROUND_DEPTH'  
        IF (iversion>7 .OR. iversion==7 .AND. ibugfix>=3) yrecfm2='D_GROUND_DPT'
      ENDIF            
      CALL read_surf_isba_par_n(dtco, u, gcp, kpatch, hfilepgdtype, yrecfm2, kluout, kni, &
                                iversion, ibugfix, gdata_ground_depth, zground_depth(:,:),iresp,hdir='E')
      !
      IF (.NOT.any(gdata_dg(1:iend))) THEN
        !
        IF (hisba=="2-L") THEN
          !
          pdepth(:,2,:) = zground_depth(:,:)
          pdepth(:,1,:) = xundef
          WHERE (zground_depth(:,:)/=xundef) pdepth(:,1,:) = 0.01
          gread_ok = .true.
          !
        ELSEIF (hisba=="3-L") THEN
          !
          pdepth(:,3,:) = zground_depth(:,:)
          pdepth(:,1,:) = xundef
          WHERE (zground_depth(:,:)/=xundef) pdepth(:,1,:) = 0.01
          IF (any(gdata_root_depth(1:iend))) gread_ok = .true.
          !
        ELSEIF (hisba=="DIF") THEN
          !
          ALLOCATE(zsoilgrid(klayer))
          DO jl=1,klayer
            WRITE(ylvl,'(I4)') jl
            yrecfm2='SOILGRID'//adjustl(ylvl(:len_trim(ylvl)))
            CALL read_surf(hfilepgdtype,yrecfm,zsoilgrid(jl),iresp)
            pdepth(:,jl,:) = zsoilgrid(jl)
          ENDDO
          DEALLOCATE(zsoilgrid)
          gread_ok = .true.
          !
        ENDIF
      ENDIF
      !
      DO jp=1,kpatch
        DO jj=1,kni
          DO jl=1,klayer
            IF ( pdepth(jj,jl,jp) <= zground_depth(jj,jp) .AND. zground_depth(jj,jp) < xundef ) &
                kwg_layer(jj,jp) = jl
          ENDDO
        ENDDO
      ENDDO
      DEALLOCATE(zground_depth)
      !
    ENDIF
    !
  ENDIF
  !
ELSE IF (hnat=='GD' .AND. .NOT.gecoclimap ) THEN
  !
  !* computes soil layers from vegetation fractions read in the file
  !
  CALL read_surf(hfilepgdtype,'D_TYPE_HVEG',yhveg,iresp,hdir='E')
  CALL read_surf(hfilepgdtype,'D_TYPE_LVEG',ylveg,iresp,hdir='E')
  CALL read_surf(hfilepgdtype,'D_TYPE_NVEG',ynveg,iresp,hdir='E')
  CALL read_surf(hfilepgdtype,'D_FRAC_HVEG',zhveg,iresp,hdir='E')
  CALL read_surf(hfilepgdtype,'D_FRAC_LVEG',zlveg,iresp,hdir='E')
  CALL read_surf(hfilepgdtype,'D_FRAC_NVEG',znveg,iresp,hdir='E')
  ! Ground layers
  CALL garden_soil_depth(ynveg,ylveg,yhveg,znveg,zlveg,zhveg,pdepth)
  !
ELSEIF (hnat=='GR' .AND. .NOT.gecoclimap ) THEN
! Depth of greenroof ground layers
  pdepth(:, 1,:) = xgrid_soil(ngrid_level - 5)
  pdepth(:, 2,:) = xgrid_soil(ngrid_level - 4)
  pdepth(:, 3,:) = xgrid_soil(ngrid_level - 3)
  pdepth(:, 4,:) = xgrid_soil(ngrid_level - 2)
  pdepth(:, 5,:) = xgrid_soil(ngrid_level - 1)
  pdepth(:, 6,:) = xgrid_soil(ngrid_level - 0)
END IF
!
 CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
!
IF (gecoclimap .AND. .NOT.gread_ok ) THEN
  !
  IF (iversion>8 .OR. (iversion==8 .AND. ibugfix>=1)) THEN
    CALL open_aux_io_surf(hfile,hfiletype,'FULL  ')   
    CALL read_surf(hfiletype,'WRITE_EXT  ',gread_ext,iresp,hdir='-')
    CALL close_aux_io_surf(hfile,hfiletype)
  ELSE
    gread_ext = .false.
  ENDIF
  !  
  IF (gread_ext) THEN
    !
    CALL open_aux_io_surf(hfile,hfiletype,ysurf)
    yrecfm='VERSION'
    CALL read_surf(hfiletype,yrecfm,iversion,iresp)
    yrecfm='BUG'
    CALL read_surf(hfiletype,yrecfm,ibugfix,iresp)
    gdim = (iversion>8 .OR. iversion==8 .AND. ibugfix>0)
    gdim2 = gdim
    IF (gdim) CALL read_surf(hfiletype,'SPLIT_PATCH',gdim2,iresp) 
    DO jl=1,klayer
      WRITE(ylvl,'(I4)') jl
      yrecfm='DG'//adjustl(ylvl(:len_trim(ylvl)))
      CALL make_choice_array(hfiletype, kpatch, gdim2, yrecfm, pdepth(:,jl,:),hdir='E')   
    END DO     
    CALL close_aux_io_surf(hfile,hfiletype)
    !
  ELSE
    !
    CALL open_aux_io_surf(hfilepgd,hfilepgdtype,'FULL  ')
    !
    !* reading of the cover to obtain the depth of inter-layers
    !
    IF (gdim.AND.gecosg) THEN
      ALLOCATE(gcover(sum(ntype)))
    ELSE
      ALLOCATE(gcover(ncover))
    ENDIF
    !
    CALL old_name(hfilepgdtype,'COVER_LIST      ',yrecfm,hdir='-')
    CALL read_surf(hfilepgdtype,yrecfm,gcover(:),iresp,hdir='-')
    !
    IF (nrank==npio) THEN 
      ALLOCATE(zcover(kni,count(gcover)))
    ELSE
      ALLOCATE(zcover(0,0))
    ENDIF
    yrecfm='COVER'
    CALL read_surf_cov(hfilepgdtype,yrecfm,zcover(:,:),gcover(:),iresp,hdir='E')
    !
    CALL close_aux_io_surf(hfilepgd,hfilepgdtype) 
    !
    !* computes soil layers
    !
    !* permafrost distribution for soil depth
    !
    gperm   =.false.
    zperm(:)=0.0
    !
    IF (hnat=='NAT'.AND.(iversion>7 .OR. (iversion==7 .AND. ibugfix>3)))THEN
      CALL open_aux_io_surf(hfilepgd,hfilepgdtype,'NATURE')        
      yrecfm='PERMAFROST'
      CALL read_surf(hfilepgdtype,yrecfm,gperm,iresp,hdir='-')
      IF(gperm)THEN
        yrecfm='PERM'
        CALL read_surf(hfilepgdtype,yrecfm,zperm(:),iresp,hdir='E')           
      ENDIF
      CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
    ENDIF
    !
    IF (SIZE(gcover)/=jpcover) THEN
      CALL open_aux_io_surf(hfilepgd,hfilepgdtype,'NATURE') 
      CALL read_pgd_cover_garden(hfilepgdtype,ggarden)
      CALL read_arrange_cover(hfilepgdtype,gwater_to_nature,gtown_to_rock,'A')
      CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
    ELSE
      ggarden = .false.
      gwater_to_nature = .false.
      gtown_to_rock    = .false.
    ENDIF
    !
    IF (nrank==npio) THEN
      CALL convert_cover_isba(dtco, io%CALBEDO, &
                              hisba,io%LTR_ML,1,zcover,gcover,'   ',hnat,psoilgrid=psoilgrid, &
                              pperm=zperm,pdg=pdepth,kwg_layer=kwg_layer, &
                              owater_to_nature=gwater_to_nature, otown_to_rock=gtown_to_rock, &
                              ogarden=ggarden )
    ENDIF
    !
    DEALLOCATE(gcover,zcover)
    !
  ENDIF
  !
ENDIF
!
nvegtype = invegtype_save
jpcover = ijpcover_save
!-------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('MODE_READ_EXTERN:READ_EXTERN_DEPTH',1,zhook_handle)
!-------------------------------------------------------------------
!
END SUBROUTINE read_extern_depth
!
!
!-------------------------------------------------------------------
!---------------------------------------------------------------------------------------
!
!     #######################
      SUBROUTINE read_extern_isba (U, DTCO, GCP, IO, &
                                   HFILE,HFILETYPE,HFILEPGD,HFILEPGDTYPE,&
                                  KLUOUT,KNI,HFIELD,HNAME,PFIELD,PDEPTH,OKEY)
!     #######################
!
!
USE modd_isba_n, ONLY : isba_np_t, isba_k_t, isba_p_t, isba_np_init, isba_k_init
USE modd_surf_atm_n, ONLY : surf_atm_t
USE modd_grid_conf_proj_n, ONLY : grid_conf_proj_t
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_isba_options_n, ONLY : isba_options_t
!
USE modd_surfex_mpi, ONLY : nrank
USE modd_isba_par,    ONLY : xoptimgrid
!
USE mode_soil
USE modi_read_surf
USE modi_isba_soc_parameters
USE modi_pack_same_rank
!
IMPLICIT NONE
!
!* dummy arguments
!  ---------------
!
TYPE(surf_atm_t), INTENT(INOUT) :: U
TYPE(data_cover_t), INTENT(INOUT) :: DTCO
TYPE(grid_conf_proj_t),INTENT(INOUT) :: GCP
!
TYPE(isba_options_t), INTENT(INOUT) :: IO
!
 CHARACTER(LEN=28),    INTENT(IN)  :: HFILE     ! name of file
 CHARACTER(LEN=6),     INTENT(IN)  :: HFILETYPE ! type of input file
 CHARACTER(LEN=28),    INTENT(IN)  :: HFILEPGD     ! name of file
 CHARACTER(LEN=6),     INTENT(IN)  :: HFILEPGDTYPE ! type of input file
INTEGER,              INTENT(IN)  :: KLUOUT    ! logical unit of output listing
INTEGER,              INTENT(IN)  :: KNI       ! number of points
 CHARACTER(LEN=7),     INTENT(IN)  :: HFIELD    ! field name
 CHARACTER(LEN=*),     INTENT(IN)  :: HNAME     ! field name in the file
REAL, DIMENSION(:,:,:), POINTER   :: PFIELD    ! field to initialize
REAL, DIMENSION(:,:,:), POINTER   :: PDEPTH    ! depth of each inter-layer
LOGICAL, OPTIONAL,  INTENT(INOUT) :: OKEY
!
!
!* local variables
!  ---------------
!
TYPE(isba_k_t) :: YK
TYPE(isba_np_t) :: YNP
TYPE(isba_p_t), POINTER :: PK
!
 CHARACTER(LEN=12) :: YRECFM         ! Name of the article to be read
#ifdef MNH_PARALLEL
 CHARACTER(LEN=8)  :: YPATCH
#endif 
 CHARACTER(LEN=4)  :: YLVL
 CHARACTER(LEN=4)  :: YPEDOTF        ! type of pedo-transfert function
 CHARACTER(LEN=3)  :: YISBA          ! type of ISBA soil scheme
 CHARACTER(LEN=3)  :: YNAT           ! type of surface (nature, garden)
!
INTEGER           :: IRESP          ! reading return code
INTEGER           :: ILAYER, ILAYER_SAVE         ! number of layers
INTEGER           :: JL         ! loop counter
INTEGER           :: IPATCH         ! number of patch
INTEGER           :: JP         ! loop counter
INTEGER           :: JVEGTYPE       ! loop counter
LOGICAL           :: GTEB           ! TEB field
LOGICAL           :: GGD 
INTEGER           :: IWORK          ! work integer
INTEGER           :: JI
!
REAL,  DIMENSION(:,:),   ALLOCATABLE   :: ZCLAY     ! clay fraction
REAL,  DIMENSION(:,:),   ALLOCATABLE   :: ZSAND     ! sand fraction
REAL,  DIMENSION(:,:),   ALLOCATABLE   :: ZCONDSAT  ! sand fraction
REAL,  DIMENSION(:,:),   ALLOCATABLE   :: ZMPOTSAT  ! sand fraction
REAL,  DIMENSION(:,:),   ALLOCATABLE   :: ZBCOEF    ! sand fraction
REAL,  DIMENSION(:,:),   ALLOCATABLE   :: ZSOC_GR   ! sand fraction
REAL,  DIMENSION(:),   ALLOCATABLE   :: ZSOILGRID
REAL,  DIMENSION(:),   ALLOCATABLE   :: ZNAT      ! natural surface fraction
!
REAL,  DIMENSION(:,:),   ALLOCATABLE :: ZWWILT    ! wilting point
REAL,  DIMENSION(:,:),   ALLOCATABLE :: ZWFC      ! field capacity
REAL,  DIMENSION(:,:),   ALLOCATABLE :: ZWSAT     ! saturation
REAL,  DIMENSION(:,:),   ALLOCATABLE :: ZPATCH
!
REAL,  DIMENSION(KNI,2)              :: ZSOC
!
REAL, DIMENSION(:,:), ALLOCATABLE :: ZWORK, ZFRACSOC
INTEGER, DIMENSION(:,:), ALLOCATABLE :: IWG_LAYER
!
LOGICAL                              :: GTEMP_ARP ! Arpege soil temperature profile
LOGICAL                              :: GSOC_DATA ! Soil organic carbon (data in pgd)
LOGICAL                              :: GSOC      ! Soil organic carbon (physical option)
!
REAL, PARAMETER   :: ZWSAT_OM      = 0.9       ! Porosity of OM (m3/m3)
REAL, PARAMETER   :: ZMPOT_WWILT   = -150.     ! Matric potential at wilting point (m)
REAL, PARAMETER   :: ZHYDCOND_WFC  = 1.157e-9  ! Hydraulic conductivity at field capacity (m/s)
!
INTEGER :: IVERSION   ! surface version
INTEGER :: IBUGFIX, ISIZE
LOGICAL :: GDIM, GDIM2
LOGICAL :: GDATA_WSAT, GDATA_WWILT, GDATA_WFC, GDATA_CONDSAT, GDATA_MPOTSAT, GDATA_BCOEF, GCALC
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MODE_READ_EXTERN:READ_EXTERN_ISBA',0,zhook_handle)
WRITE  (kluout,*) ' | Reading ',hfield,' in externalized file'
!
!------------------------------------------------------------------------------
! Init
!
gteb = (hname(1:3)=='TWN' .OR. hname(1:3)=='GD_' .OR. hname(1:3)=='GR_' &
        .OR. hname(4:6)=='GD_' .OR. hname(4:6)=='GR_')
ggd = .false.
IF (gteb) ggd = (hname(1:3)=='TWN' .OR. hname(1:3)=='GD_' .OR. hname(4:6)=='GD_')
!
gtemp_arp = .false.
gsoc      = .false.
gsoc_data = .false.
!
!------------------------------------------------------------------------------
!
 CALL open_aux_io_surf(hfilepgd,hfilepgdtype,'FULL  ')
yrecfm='VERSION'
 CALL read_surf(hfilepgdtype,yrecfm,iversion,iresp,hdir='-')
yrecfm='BUG'
 CALL read_surf(hfilepgdtype,yrecfm,ibugfix,iresp,hdir='-')
 CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
gdim = (iversion>8 .OR. iversion==8 .AND. ibugfix>0)

!
IF (gteb) THEN
  CALL open_aux_io_surf(hfilepgd,hfilepgdtype,'TOWN  ')
ELSE
  CALL open_aux_io_surf(hfilepgd,hfilepgdtype,'NATURE')
ENDIF
!
!* Read number of soil layers
!
yrecfm='GROUND_LAYER'
IF (gteb) THEN 
  yrecfm='TWN_LAYER'
  IF (iversion>7 .OR. iversion==7 .AND. ibugfix>=3) THEN
    IF (ggd) THEN 
      yrecfm='GD_LAYER'
    ELSE
      yrecfm='GR_LAYER'
    ENDIF
  ENDIF
ENDIF
 CALL read_surf(hfilepgdtype,yrecfm,ilayer,iresp,hdir='-')
!
!* number of tiles
!
ipatch=1
IF (.NOT. gteb) THEN
  yrecfm='PATCH_NUMBER'
  CALL read_surf(hfilepgdtype,yrecfm,ipatch,iresp,hdir='-')
END IF
!
!* soil scheme
!
yrecfm='ISBA'
IF (gteb) THEN 
  yrecfm='TWN_ISBA'
  IF (iversion>7 .OR. iversion==7 .AND. ibugfix>=3) THEN
    IF (ggd) THEN 
      yrecfm='GD_ISBA'
    ELSE
      yrecfm='GR_ISBA'
    ENDIF
  ENDIF
ENDIF
 CALL read_surf(hfilepgdtype,yrecfm,yisba,iresp,hdir='-')
IF(yisba=='DIF'.AND.PRESENT(okey))THEN
  okey=.false.
ENDIF
!
IF (iversion>=7) THEN
  !
  !* Pedo-transfert function
  !
  yrecfm='PEDOTF'
  IF (gteb) THEN 
    yrecfm='TWN_PEDOTF'
    IF (iversion>7 .OR. iversion==7 .AND. ibugfix>=3) yrecfm='GD_PEDOTF'
  ENDIF
  CALL read_surf(hfilepgdtype,yrecfm,ypedotf,iresp,hdir='-')
  !
ELSE
  ypedotf = 'CH78'
ENDIF
!
!Only Brook and Corey with Force-Restore scheme
IF(yisba/='DIF')THEN
  ypedotf='CH78'
ENDIF
!
IF (hfield=='WG    ' .OR. hfield=='WGI   ') THEN
  !
  ALLOCATE (zwfc(kni,ilayer))
  ALLOCATE (zwwilt(kni,ilayer))
  ALLOCATE (zwsat(kni,ilayer))
  IF (gteb.AND..NOT.ggd) THEN
    ALLOCATE(zbcoef(kni,ilayer))
    ALLOCATE(zmpotsat(kni,ilayer))
    ALLOCATE(zcondsat(kni,ilayer))    
  ENDIF 
  !
  IF (gdim) THEN
    yrecfm='L_WFC'
    CALL read_surf(hfilepgdtype,yrecfm,gdata_wfc,iresp,hdir='E')
    IF (gdata_wfc) THEN
      DO jl=1,ilayer
        WRITE(yrecfm,fmt='(A9,I2.2)') 'D_WFC_L',jl
        CALL read_surf(hfilepgdtype,yrecfm,zwfc(:,jl),iresp,hdir='E')
      ENDDO      
    ENDIF     
    yrecfm='L_WWILT'
    CALL read_surf(hfilepgdtype,yrecfm,gdata_wwilt,iresp,hdir='E')
    IF (gdata_wwilt) THEN
      DO jl=1,ilayer
        WRITE(yrecfm,fmt='(A9,I2.2)') 'D_WWILT_L',jl
        CALL read_surf(hfilepgdtype,yrecfm,zwwilt(:,jl),iresp,hdir='E')
      ENDDO  
    ENDIF   
    yrecfm='L_WSAT'
    CALL read_surf(hfilepgdtype,yrecfm,gdata_wsat,iresp,hdir='E')
    IF (gdata_wsat) THEN
      DO jl=1,ilayer
        WRITE(yrecfm,fmt='(A9,I2.2)') 'D_WSAT_L',jl
        CALL read_surf(hfilepgdtype,yrecfm,zwsat(:,jl),iresp,hdir='E')
      ENDDO 
    ENDIF
    IF (gteb.AND..NOT.ggd) THEN         
      yrecfm='L_CONDSAT'
      CALL read_surf(hfilepgdtype,yrecfm,gdata_condsat,iresp,hdir='E')
      IF (gdata_condsat) THEN
        DO jl=1,ilayer
          WRITE(yrecfm,fmt='(A9,I2.2)') 'D_CNDSAT_L',jl
          CALL read_surf(hfilepgdtype,yrecfm,zcondsat(:,jl),iresp,hdir='E')
        ENDDO 
      ENDIF
      yrecfm='L_MPOTSAT'
      CALL read_surf(hfilepgdtype,yrecfm,gdata_mpotsat,iresp,hdir='E')
      IF (gdata_mpotsat) THEN
        DO jl=1,ilayer
          WRITE(yrecfm,fmt='(A9,I2.2)') 'D_MPTSAT_L',jl
          CALL read_surf(hfilepgdtype,yrecfm,zmpotsat(:,jl),iresp,hdir='E')
        ENDDO 
      ENDIF
      yrecfm='L_BCOEF'
      CALL read_surf(hfilepgdtype,yrecfm,gdata_bcoef,iresp,hdir='E')
      IF (gdata_bcoef) THEN
        DO jl=1,ilayer
          WRITE(yrecfm,fmt='(A9,I2.2)') 'D_BCOEF_L',jl
          CALL read_surf(hfilepgdtype,yrecfm,zbcoef(:,jl),iresp,hdir='E')
        ENDDO 
      ENDIF
    ENDIF      
  ELSE
    gdata_wwilt   = .false.
    gdata_wfc     = .false.
    gdata_wsat    = .false.  
    gdata_condsat = .false.
    gdata_mpotsat = .false.
    gdata_bcoef   = .false.
  ENDIF
  !
  gcalc = .false.
  IF (.NOT.(gdata_wwilt.AND.gdata_wfc.AND.gdata_wsat)) THEN
    IF ((.NOT.gteb.OR.ggd).OR..NOT.(gdata_condsat.AND.gdata_mpotsat.AND.gdata_bcoef)) gcalc = .true.
  ENDIF
  !
  IF (gcalc) THEN
    !
    !-------------------------------------------------------------------------------
    !
    ! *.  Read clay fraction
    !     ------------------
    !
    ALLOCATE(zclay(kni,ilayer))
    yrecfm='CLAY'
    IF (gteb) THEN 
      yrecfm='TWN_CLAY'
      IF (iversion>7 .OR. iversion==7 .AND. ibugfix>=3) yrecfm='GD_CLAY'
    ENDIF
    IF (.NOT.gteb.OR.ggd) THEN
      CALL read_surf(hfilepgdtype,yrecfm,zclay(:,1),iresp,hdir='E')
      DO jl=2,ilayer
        zclay(:,jl) = zclay(:,1)
      ENDDO
    ELSE
      DO jl=1,ilayer
        WRITE(yrecfm,fmt='(A9,I2.2)') 'D_CLAY_GR',jl
        CALL read_surf(hfilepgdtype,yrecfm,zclay(:,jl),iresp,hdir='E')
      ENDDO
    ENDIF  
    !
    !-------------------------------------------------------------------------------
    !
    ! *.  Read sand fraction
    !     ------------------
    !
    ALLOCATE(zsand(kni,ilayer))
    yrecfm='SAND'
    IF (gteb) THEN
      yrecfm='TWN_SAND'
      IF (iversion>7 .OR. iversion==7 .AND. ibugfix>=3) yrecfm='GD_SAND'
    ENDIF
    IF (.NOT.gteb.OR.ggd) THEN
      CALL read_surf(hfilepgdtype,yrecfm,zsand(:,1),iresp,hdir='E')
      DO jl=2,ilayer
        zsand(:,jl) = zsand(:,1)
      ENDDO
    ELSE
      DO jl=1,ilayer
        WRITE(yrecfm,fmt='(A9,I2.2)') 'D_SAND_GR',jl
        CALL read_surf(hfilepgdtype,yrecfm,zsand(:,jl),iresp,hdir='E')
      ENDDO
    ENDIF
    !
    IF (gteb.AND..NOT.ggd) THEN
      ALLOCATE(zsoc_gr(kni,ilayer))
      DO jl=1,ilayer
        WRITE(yrecfm,fmt='(A7,I2.2)') 'D_OM_GR',jl
        CALL read_surf(hfilepgdtype,yrecfm,zsoc_gr(:,jl),iresp,hdir='E')
      ENDDO      
    ENDIF
    !
  ENDIF
  !
ENDIF
!-------------------------------------------------------------------------------
!
!
! *.  Soil organic carbon profile
!     ---------------------------
!
IF ( (.NOT.gteb).AND.(iversion>7.OR.(iversion==7.AND.ibugfix>3)) &
     .AND.(yisba=='DIF').AND.(hfield=='WG    '.OR.hfield=='WGI   ') ) THEN
   yrecfm='SOCP'
   CALL read_surf(hfilepgdtype,yrecfm,gsoc_data,iresp)
   IF(gsoc_data)THEN
     yrecfm='SOC_TOP'
     CALL read_surf(hfilepgdtype,yrecfm,zsoc(:,1),iresp,hdir='E')
     yrecfm='SOC_SUB'
     CALL read_surf(hfilepgdtype,yrecfm,zsoc(:,2),iresp,hdir='E')
     WHERE(zsoc(:,:)==xundef)zsoc(:,:)=0.0
   ENDIF
ENDIF
!
!-------------------------------------------------------------------------------
!
! *.  Read soil grid
!     --------------
!
!* Reference grid for DIF
!
IF(yisba=='DIF') THEN
  ALLOCATE(zsoilgrid(ilayer))
  zsoilgrid=xundef
  IF (.NOT.gteb.OR.ggd) THEN
    IF (iversion>=8) THEN
      DO jl=1,ilayer
        WRITE(ylvl,'(I4)') jl
        yrecfm='SOILGRID'//adjustl(ylvl(:len_trim(ylvl)))
        IF (gteb) THEN 
          yrecfm='GD_SGRID'//adjustl(ylvl(:len_trim(ylvl)))
        ENDIF
        CALL read_surf(hfilepgdtype,yrecfm,zsoilgrid(jl),iresp,hdir='A')
      ENDDO
    ELSEIF (iversion==7 .AND. ibugfix>=2.AND.ggd) THEN
      yrecfm='SOILGRID'
      IF (gteb) THEN 
        yrecfm='TWN_SOILGRID'
        IF (iversion>7 .OR. iversion==7 .AND. ibugfix>=3) yrecfm='GD_SOILGRID'
      ENDIF
      CALL read_surf(hfilepgdtype,yrecfm,zsoilgrid,iresp,hdir='A')
    ELSE
      zsoilgrid(1:ilayer) = xoptimgrid(1:ilayer)
    ENDIF
  ELSE
    zsoilgrid(1:ilayer) = xoptimgrid(1:ilayer)
  ENDIF
ELSE
  ALLOCATE(zsoilgrid(0))
ENDIF
!
ALLOCATE(iwg_layer(kni,ipatch))
iwg_layer(:,:) = nundef        
!
! *.  Read fraction of nature
!     --------------
!
ALLOCATE(znat(kni))
IF (iversion>=7) THEN
  IF (gteb) THEN
    CALL read_surf(hfilepgdtype,'FRAC_TOWN',znat,iresp,hdir='E')
  ELSE
    CALL read_surf(hfilepgdtype,'FRAC_NATURE',znat,iresp,hdir='E')
  ENDIF
ELSE
  znat=1.0
ENDIF
!
!
 CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
!
!
IF (.NOT.gteb .AND. hfield=='TG    ' .AND. (yisba=='2-L' .OR. yisba=='3-L') ) THEN
  IF (iversion>7) THEN
     yrecfm='TEMPARP'
     CALL open_aux_io_surf(hfile,hfiletype,'NATURE')
     CALL read_surf(hfiletype,yrecfm,gtemp_arp,iresp,hdir='-')
     IF(gtemp_arp)THEN
       yrecfm = 'NTEMPLARP'
       CALL read_surf(hfiletype,yrecfm,ilayer,iresp,hdir='-')     
     ENDIF
     CALL close_aux_io_surf(hfile,hfiletype)
   ENDIF
ENDIF  
!
!
IF ((hfield=='TG    ') .AND. (yisba=='2-L' .OR. yisba=='3-L')) THEN
  ALLOCATE(pdepth(kni,ilayer,nvegtype))
  DO jvegtype=1,nvegtype
    pdepth(:,1,jvegtype) = 0.01
    pdepth(:,2,jvegtype) = 0.40
    IF (ilayer==3) pdepth(:,3,jvegtype) = 5.00
!   GTEMP_ARP case
    IF (gtemp_arp) THEN
       pdepth(:,3,jvegtype) = 1.0
       DO jl=4,ilayer
          pdepth(:,jl,jvegtype) = pdepth(:,jl-1,jvegtype)+1.
       ENDDO
    ENDIF    
  END DO
ELSE
  ynat='NAT'
  IF (gteb) THEN
    IF (ggd) THEN
      ynat='GD'
    ELSE
      ynat='GR'
    ENDIF
  ENDIF
  !
  CALL read_extern_depth(u, dtco,  gcp, io,                      &
                         hfile,hfiletype,hfilepgd,hfilepgdtype,  &
                         kluout,yisba,ynat,hfield,kni,           &
                         ilayer,ipatch,zsoilgrid,pdepth,iversion,iwg_layer)
  !
END IF
!
DEALLOCATE(zsoilgrid)
!
!
! *.  Read soil variable profile
!     --------------------------
!
IF (gteb) THEN
  CALL open_aux_io_surf(hfile,hfiletype,'TOWN  ')
ELSE
  CALL open_aux_io_surf(hfile,hfiletype,'NATURE')
ENDIF
!
yrecfm='VERSION'
 CALL read_surf(hfiletype,yrecfm,iversion,iresp)
yrecfm='BUG'
 CALL read_surf(hfiletype,yrecfm,ibugfix,iresp)
gdim = (iversion>8 .OR. iversion==8 .AND. ibugfix>0)
gdim2 = gdim
IF (gdim) CALL read_surf(hfiletype,'SPLIT_PATCH',gdim2,iresp)
!
iwork=ilayer
IF(yisba=='2-L'.OR.yisba=='3-L') THEN
  SELECT CASE(hfield)
    CASE('TG    ')
      IF(gtemp_arp)THEN
        iwork=ilayer
      ELSE
        iwork=2
      ENDIF
    CASE('WGI   ')
      iwork=2
  END SELECT
ENDIF
!
ALLOCATE(pfield(kni,ilayer,ipatch))
!
DO jl=1,iwork
  WRITE(ylvl,'(I4)') jl
  yrecfm=trim(hname)//adjustl(ylvl(:len_trim(ylvl)))
  IF (gteb) THEN
    CALL make_choice_array(hfiletype, ipatch, gdim2, yrecfm, pfield(:,jl,:),hdir='E',kpatch=0)
  ELSE
    CALL make_choice_array(hfiletype, ipatch, gdim2, yrecfm, pfield(:,jl,:),hdir='E')
  ENDIF
END DO
!
DO jp=1,ipatch
  DO jl=1,iwork
    WHERE (znat(:)==0.) pfield(:,jl,jp) = xundef
  ENDDO
END DO
!
DEALLOCATE (znat)
!
IF(yisba=='3-L') THEN
  SELECT CASE(hfield)
    CASE('TG    ')
      IF(.NOT.gtemp_arp)pfield(:,3,:)=pfield(:,2,:)
    CASE('WGI   ')       
     pfield(:,3,:)=pfield(:,2,:)
  END SELECT
ENDIF
!
 CALL close_aux_io_surf(hfile,hfiletype)
!
!
! *.  Compute relative humidity from units kg/m^2 (SWI)
!     ------------------------------------------------
!
IF (hfield=='WG    ' .OR. hfield=='WGI   ') THEN
  !
  IF (gcalc) THEN
    !
    ! Compute ISBA model constants
    !
    DO jl=1,ilayer
      zwsat(:,jl) = wsat_func(zclay(:,jl),zsand(:,jl),ypedotf)
      zwwilt(:,jl) = wwilt_func(zclay(:,jl),zsand(:,jl),ypedotf)
      IF(yisba=='DIF')THEN
        zwfc(:,jl) = w33_func(zclay(:,jl),zsand(:,jl),ypedotf)
      ELSE
        zwfc(:,jl) = wfc_func(zclay(:,jl),zsand(:,jl),ypedotf)
      ENDIF
    ENDDO
    !
    IF (gteb.AND..NOT.ggd) THEN
      !
      DO jl=1,ilayer 
        zbcoef(:,jl) = bcoef_func(zclay(:,jl),zsand(:,jl),ypedotf)
        zcondsat(:,jl) = hydcondsat_func(zclay(:,jl),zsand(:,jl),ypedotf) 
        zmpotsat(:,jl) = matpotsat_func(zclay(:,jl),zsand(:,jl),ypedotf)
     
        WHERE (zwsat(:,jl)/=xundef) 
          zwsat(:,jl) =  zsoc_gr(:,jl)* zwsat_om +(1-zsoc_gr(:,jl))* zwsat(:,jl)
          zwwilt(:,jl) = exp(((log(-1*zmpot_wwilt)-log(-1*zmpotsat(:,jl)))   &
                          / (-1*zbcoef(:,jl)))+log(zwsat(:,jl)))
          zwfc(:,jl) = exp(((log(zhydcond_wfc)-log(zcondsat(:,jl)))      &
                          / (2*zbcoef(:,jl)+3))+log(zwsat(:,jl)))
        END WHERE
      END DO
      DEALLOCATE(zbcoef,zmpotsat,zcondsat)
      !
      ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      ! Validation case : experimental values for Nancy 2011 case
      ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      ! Substrate layer
      DO jl=1,4
        zwsat(:,jl) = 0.674     ! Value tested
        zwwilt(:,jl) = 0.15      ! from OBS-NANCY      
        zwfc(:,jl) = 0.37      ! from OBS-NANCY
      ENDDO
      ! Drainage layer
      DO jl=5,6
        zwsat(:,jl) = 0.9       ! Value tested
        zwwilt(:,jl) = 0.15      ! sert à initialiser le WG ds la couche
        zwfc(:,jl) = 0.37      ! sert à initialiser le WG ds la couche
      ENDDO

    ENDIF
    !
    DEALLOCATE (zsand)
    DEALLOCATE (zclay)
    !
    IF(gsoc_data)THEN
      !
      ALLOCATE(zwork(kni,ipatch))
      !
      CALL open_aux_io_surf(hfile,hfiletype,'NATURE')
      ! 
      yrecfm='SOC'
      CALL read_surf(hfiletype,yrecfm,gsoc,iresp,hdir='-')
      yrecfm='PATCH'
      CALL make_choice_array(hfiletype, ipatch, gdim2, yrecfm, zwork,hdir='E')
      WHERE(zwork(:,:)==xundef)zwork(:,:)=0.0
      CALL close_aux_io_surf(hfile,hfiletype)    
      !
      IF(gsoc)THEN
        !
        ALLOCATE(ynp%AL(ipatch))
        DO jp = 1,ipatch
          pk => ynp%AL(jp)
          pk%NSIZE_P = kni
          ALLOCATE(pk%NR_P(kni))
          DO ji = 1,kni
            pk%NR_P(ji) = ji
          ENDDO
          ALLOCATE(pk%XPATCH(kni))
          pk%XPATCH = zwork(:,jp)
          ALLOCATE(pk%XDG(kni,ilayer))
          pk%XDG(:,:) = pdepth(:,:,jp)
          ALLOCATE(pk%XCONDSAT(kni,ilayer))
          pk%XCONDSAT(:,:) = 0.0
        ENDDO
        ALLOCATE(yk%XBCOEF   (kni,ilayer))
        ALLOCATE(yk%XMPOTSAT (kni,ilayer))
        ALLOCATE(yk%XHCAPSOIL(kni,ilayer))
        ALLOCATE(yk%XCONDDRY (kni,ilayer))
        ALLOCATE(yk%XCONDSLD (kni,ilayer))
        ALLOCATE(yk%XWD0     (kni,ilayer))
        ALLOCATE(yk%XKANISO  (kni,ilayer))
        ALLOCATE(zfracsoc(kni,ilayer))
        yk%XBCOEF   (:,:)=0.0
        yk%XMPOTSAT (:,:)=0.0
        yk%XHCAPSOIL(:,:)=0.0
        yk%XCONDDRY (:,:)=xundef
        yk%XCONDSLD (:,:)=xundef
        yk%XWD0     (:,:)=0.0
        yk%XKANISO   (:,:)=0.0
        zfracsoc(:,:)=0.0
        CALL isba_soc_parameters ('NONE', zsoc, yk, ynp, zfracsoc, &
                                  zwsat, zwfc, zwwilt, ipatch)
        CALL isba_k_init(yk)
        CALL isba_np_init(ynp,ipatch)
        !
      ENDIF
      !
    ENDIF
    !
  ENDIF
  !
  IF(yisba=='DIF')THEN
    !     
    ! extrapolation of deep layers
    DO jp=1,ipatch
      DO ji=1,kni
        iwork=iwg_layer(ji,jp)
        IF(iwork<ilayer)THEN
          DO jl=iwork+1,ilayer
            pfield(ji,jl,jp)=pfield(ji,iwork,jp)
          ENDDO
        ENDIF
      ENDDO
    ENDDO
  ENDIF
  !
  IF (hfield=='WG    ') THEN
    DO jp=1,ipatch
      DO jl=1,ilayer
        DO ji=1,kni
          IF(pfield(ji,jl,jp)/=xundef)THEN
            pfield(ji,jl,jp) = max(min(pfield(ji,jl,jp),zwsat(ji,jl)),0.)
            !
            pfield(ji,jl,jp) = (pfield(ji,jl,jp) - zwwilt(ji,jl)) / &
                                (zwfc(ji,jl) - zwwilt(ji,jl))
          ENDIF
        END DO
      END DO
    END DO
  ELSE IF (hfield=='WGI   ') THEN
    DO jp=1,ipatch
      DO jl=1,ilayer
        WHERE(pfield(:,jl,jp)/=xundef)
          pfield(:,jl,jp) = pfield(:,jl,jp) / zwsat(:,jl) 
        END WHERE
      END DO
    END DO
  END IF
!
  DEALLOCATE (zwsat)
  DEALLOCATE (zwwilt)
  DEALLOCATE (zwfc)
!
END IF
!
DEALLOCATE(iwg_layer)
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('MODE_READ_EXTERN:READ_EXTERN_ISBA',1,zhook_handle)
!
!------------------------------------------------------------------------------
!
END SUBROUTINE read_extern_isba
!
!------------------------------------------------------------------------------
!
END MODULE mode_read_extern