prep_hor_isba_field.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 prep_hor_isba_field (DTCO, UG, U, USS, GCP, IG, IO, S, NK, NP, NPE, TPTIME,  &
                                HPROGRAM,HSURF,HATMFILE,HATMFILETYPE,HPGDFILE,HPGDFILETYPE,YDCTL,OKEY)
!     #################################################################################
!
!!****  *PREP_HOR_ISBA_FIELD* - reads, interpolates and prepares an ISBA field
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     V. Masson 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!      P. Le Moigne 10/2005, Phasage Arome
!!      P. Le Moigne 03/2007, Ajout initialisation par ascllv
!!      B. Decharme  01/2009, Optional Arpege deep soil temperature initialization
!!      M. Lafaysse  07/2012, allow netcdf input files
!!      B. Decharme  07/2012, Bug init uniform snow
!!      M. Lafaysse 11/2012,  snow liquid water content
!!      B. Decharme  03/2014, external init with FA files
!!                            new vertical interpolation
!!      P Samuelsson 10/2014, MEB
!!      P. Marguinaud10/2014, Support for a 2-part PREP
!!------------------------------------------------------------------
!
USE modd_type_snow
USE modd_type_date_surf, ONLY : date_time
USE modd_surfex_mpi, ONLY : nrank, npio, ncomm, nproc
!
USE modd_data_cover_n, ONLY : data_cover_t
!
USE modd_sfx_grid_n, ONLY : grid_t
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_nk_t, isba_np_t, isba_npe_t, isba_k_t, isba_s_t, &
                        isba_pe_t, isba_p_t
!
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_surf_atm_n, ONLY : surf_atm_t
USE modd_sso_n, ONLY : sso_t
USE modd_grid_conf_proj_n, ONLY : grid_conf_proj_t
!
USE modd_prep,     ONLY : cingrid_type, cinterp_type, xzs_ls, &
                          xlat_out, xlon_out, xx_out, xy_out, &
                          linterp, cmask
USE modd_grid_grib, ONLY : cinmodel  
!
USE modd_prep_isba, ONLY : xgrid_soil, ngrid_level, lsnow_ideal,    &
                           xwsnow, xrsnow, xtsnow, xlwcsnow, xasnow,  &
                           xsg1snow, xsg2snow, xhistsnow, xagesnow


USE modd_isba_par,       ONLY : xwgmin
USE modd_data_cover_par, ONLY : nvegtype
USE modd_surf_par,       ONLY : xundef,nundef
!
USE mode_prep_ctl, ONLY : prep_ctl, prep_ctl_int_part2, prep_ctl_int_part4
!
USE modi_prep_grib_grid
USE modi_read_prep_isba_conf
USE modi_read_prep_isba_snow
USE modi_prep_isba_ascllv
USE modi_prep_isba_grib
USE modi_prep_isba_unif
USE modi_prep_isba_buffer
USE modi_abor1_sfx
USE modi_hor_interpol
USE modi_put_on_all_vegtypes
USE modi_vegtype_grid_to_patch_grid
USE modi_prep_hor_snow_fields
USE modi_get_luout
USE modi_prep_isba_extern
USE modi_prep_isba_netcdf
USE modi_pack_same_rank
USE modi_allocate_gr_snow
USE modi_get_prep_interp
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
#ifdef SFX_MPI
include "mpif.h"
#endif
!
!*      0.1    declarations of arguments
!
!
TYPE(data_cover_t), INTENT(INOUT) :: DTCO
!
TYPE(grid_t), INTENT(INOUT) :: IG
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_s_t), INTENT(INOUT) :: S
TYPE(isba_nk_t), INTENT(INOUT) :: NK
TYPE(isba_np_t), INTENT(INOUT) :: NP
TYPE(isba_npe_t), INTENT(INOUT) :: NPE
TYPE(date_time), INTENT(IN) :: TPTIME
!
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
TYPE(surf_atm_t), INTENT(INOUT) :: U
TYPE(sso_t), INTENT(INOUT) :: USS
TYPE(grid_conf_proj_t),INTENT(INOUT) :: GCP
!
type(prep_ctl),    INTENT(INOUT) :: ydctl
!
 CHARACTER(LEN=6),   INTENT(IN)  :: HPROGRAM  ! program calling surf. schemes
 CHARACTER(LEN=7),   INTENT(IN)  :: HSURF     ! type of field
 CHARACTER(LEN=28),  INTENT(IN)  :: HATMFILE    ! name of the Atmospheric file
 CHARACTER(LEN=6),   INTENT(IN)  :: HATMFILETYPE! type of the Atmospheric file
 CHARACTER(LEN=28),  INTENT(IN)  :: HPGDFILE    ! name of the Atmospheric file
 CHARACTER(LEN=6),   INTENT(IN)  :: HPGDFILETYPE! type of the Atmospheric file
!
LOGICAL, OPTIONAL,  INTENT(INOUT):: OKEY
!
!*      0.2    declarations of local variables
!
 CHARACTER(LEN=6)              :: YFILETYPE ! type of input file
 CHARACTER(LEN=28)             :: YFILE     ! name of file
 CHARACTER(LEN=6)              :: YFILETYPE_SNOW ! type of input file
 CHARACTER(LEN=28)             :: YFILE_SNOW     ! name of file
 CHARACTER(LEN=6)              :: YFILEPGDTYPE_SNOW ! type of input file
 CHARACTER(LEN=28)             :: YFILEPGD_SNOW     ! name of file 
 CHARACTER(LEN=6)              :: YFILEPGDTYPE ! type of input file
 CHARACTER(LEN=28)             :: YFILEPGD     ! name of file
REAL, POINTER, DIMENSION(:,:,:) :: ZFIELDIN   ! field to interpolate horizontally
REAL, POINTER, DIMENSION(:,:,:) :: ZFIELDOUTP ! field interpolated   horizontally
REAL, POINTER, DIMENSION(:,:,:) :: ZFIELDOUTV !
!
TYPE(nsurf_snow) :: TNPSNOW
!
TYPE(isba_k_t), POINTER :: KK
TYPE(isba_p_t), POINTER :: PK
TYPE(isba_pe_t), POINTER :: PEK
!
TYPE fout
  REAL, DIMENSION(:,:), ALLOCATABLE :: ZOUT
END TYPE fout
TYPE nfout
  TYPE(fout), DIMENSION(:), ALLOCATABLE :: AL
END TYPE nfout
type(nfout) :: zw
TYPE(nfout) :: ZF
!
REAL, ALLOCATABLE, DIMENSION(:,:) :: ZDG       ! out T grid (x, output soil grid, npatch)
type(date_time)                :: tztime_grib    ! current date and time
!
 CHARACTER(LEN=3) :: YSNOW_SCHEME
INTEGER :: ISNOW_NLAYER
!
INTEGER, DIMENSION(IO%NPATCH) :: ISIZE_P
INTEGER, DIMENSION(SIZE(IG%XLAT),IO%NPATCH) :: IR_P
!
INTEGER                       :: ILUOUT    ! output listing logical unit
!
LOGICAL                       :: GUNIF     ! flag for prescribed uniform field
LOGICAL                       :: GUNIF_SNOW! flag for prescribed uniform field
INTEGER                       :: JP    ! loop on patches
INTEGER                       :: JVEG  ! loop on vegtypes
INTEGER                       :: INI, INL, INP, JI, JL! Work integer
INTEGER                       :: INFOMPI
INTEGER, DIMENSION(SIZE(IG%XLAT)) :: IWORK
!
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZVEGTYPE_PATCH
REAL, DIMENSION(:,:), ALLOCATABLE :: ZPATCH
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------------
!
!
!*      1.     Reading of input file name and type
!
IF (lhook) CALL dr_hook('PREP_HOR_ISBA_FIELD',0,zhook_handle)
 CALL get_luout(hprogram,iluout)
!
 CALL read_prep_isba_conf(hprogram,hsurf,yfile,yfiletype,yfilepgd,yfilepgdtype,&
                         hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,iluout,gunif)
!
cmask = 'NATURE'
!
ini=SIZE(ig%XLAT)
!
!-------------------------------------------------------------------------------------
!
!*      2.     Snow variables case?
!
IF (hsurf=='SN_VEG ') THEN
  CALL read_prep_isba_snow(hprogram, ysnow_scheme, isnow_nlayer, yfile_snow, yfiletype_snow,&
                           yfilepgd_snow, yfilepgdtype_snow, gunif_snow)
  !
  DO jp = 1,io%NPATCH
    npe%AL(jp)%TSNOW%SCHEME = ysnow_scheme
    npe%AL(jp)%TSNOW%NLAYER = isnow_nlayer
    isize_p(jp) = np%AL(jp)%NSIZE_P
    ir_p(:,jp) = 0
    ir_p(1:isize_p(jp),jp) = np%AL(jp)%NR_P
  ENDDO
  !
  IF(.NOT.gunif_snow.AND.len_trim(yfile_snow)==0.AND.len_trim(yfiletype_snow)==0)THEN
    IF(len_trim(yfile)/=0.AND.len_trim(yfiletype)/=0)THEN
      yfile_snow        = yfile
      yfiletype_snow    = yfiletype
      yfilepgd_snow     = yfilepgd
      yfilepgdtype_snow = yfilepgdtype       
    ELSE
      gunif_snow=.true.
      IF(all(xwsnow==xundef))xwsnow=0.0
    ENDIF
  ENDIF
  !
  ALLOCATE(tnpsnow%AL(io%NPATCH))
  ALLOCATE(zvegtype_patch(SIZE(s%XVEGTYPE_PATCH,1),SIZE(s%XVEGTYPE_PATCH,2),SIZE(s%XVEGTYPE_PATCH,3)))
  ALLOCATE(zpatch(SIZE(s%XPATCH,1),SIZE(s%XPATCH,2)))
  zvegtype_patch(:,:,:) = 0.
  zpatch(:,:) = 0.
  DO jp = 1,io%NPATCH
    CALL pack_same_rank(np%AL(jp)%NR_P,s%XVEGTYPE_PATCH(:,:,jp),zvegtype_patch(1:np%AL(jp)%NSIZE_P,:,jp))
    CALL pack_same_rank(np%AL(jp)%NR_P,s%XPATCH(:,jp),zpatch(1:np%AL(jp)%NSIZE_P,jp))
    tnpsnow%AL(jp)%SCHEME = npe%AL(jp)%TSNOW%SCHEME
    tnpsnow%AL(jp)%NLAYER = npe%AL(jp)%TSNOW%NLAYER
  ENDDO
  !
  CALL prep_hor_snow_fields(dtco, ig, u, gcp, hprogram, hsurf,   &
                            yfile_snow, yfiletype_snow,          &
                            yfilepgd_snow, yfilepgdtype_snow,    &
                            iluout, gunif_snow, io%NPATCH, 1,    &
                            ini,tnpsnow, tptime,                 &
                            xwsnow, xrsnow, xtsnow, xlwcsnow,    &
                            xasnow, lsnow_ideal, xsg1snow,       &
                            xsg2snow, xhistsnow, xagesnow, ydctl,&
                            pvegtype_patch=zvegtype_patch,       &
                            ppatch=zpatch, ksize_p=isize_p,      &
                            kr_p=ir_p, okey=okey   )
  !
  DEALLOCATE(zpatch)
  DEALLOCATE(zvegtype_patch)
  !
  DO jp = 1,io%NPATCH
    pek => npe%AL(jp)
    CALL allocate_gr_snow(pek%TSNOW,np%AL(jp)%NSIZE_P)
    pek%TSNOW%WSNOW = tnpsnow%AL(jp)%WSNOW
    pek%TSNOW%RHO   = tnpsnow%AL(jp)%RHO
    pek%TSNOW%ALB   = tnpsnow%AL(jp)%ALB
    IF (pek%TSNOW%SCHEME/='D95') pek%TSNOW%HEAT = tnpsnow%AL(jp)%HEAT
    IF (pek%TSNOW%SCHEME=='CRO'.OR.pek%TSNOW%SCHEME=='3-L') &
      pek%TSNOW%AGE = tnpsnow%AL(jp)%AGE
    IF (pek%TSNOW%SCHEME=='CRO') THEN
      pek%TSNOW%GRAN1 = tnpsnow%AL(jp)%GRAN1
      pek%TSNOW%GRAN2 = tnpsnow%AL(jp)%GRAN2
      pek%TSNOW%HIST = tnpsnow%AL(jp)%HIST
    ENDIF
    !
    CALL type_snow_init(tnpsnow%AL(jp))
  ENDDO
  DEALLOCATE(tnpsnow%AL)
  !
  DEALLOCATE(xwsnow)
  DEALLOCATE(xrsnow)
  DEALLOCATE(xtsnow)
  DEALLOCATE(xlwcsnow)
  DEALLOCATE(xsg1snow)
  DEALLOCATE(xsg2snow)
  DEALLOCATE(xhistsnow)
  DEALLOCATE(xagesnow)
  IF (lhook) CALL dr_hook('PREP_HOR_ISBA_FIELD',1,zhook_handle)
  RETURN
END IF
!
!-------------------------------------------------------------------------------------
!
!*      3.     Reading of input  configuration (Grid and interpolation type)
!
NULLIFY (zfieldin, zfieldoutp, zfieldoutv)
!
IF (ydctl%LPART1) THEN
!
  IF (gunif) THEN
    CALL prep_isba_unif(iluout,hsurf,zfieldin)
  ELSE IF (yfiletype=='ASCLLV') THEN
    CALL prep_isba_ascllv(dtco, ug, u, uss, hprogram,hsurf,iluout,zfieldin)
  ELSE IF (yfiletype=='GRIB  ') THEN
    CALL prep_grib_grid(yfile,iluout,cinmodel,cingrid_type,cinterp_type,tztime_grib)        
    IF (nrank==npio) CALL prep_isba_grib(hprogram,hsurf,yfile,iluout,zfieldin)        
  ELSE IF (yfiletype=='MESONH' .OR. yfiletype=='ASCII ' .OR. yfiletype=='LFI   '&
          .OR.yfiletype=='FA    '.OR. yfiletype=='AROME '.OR.yfiletype=='NC    ') THEN
    CALL prep_isba_extern(dtco, io, u, gcp, hprogram,hsurf,yfile,yfiletype,yfilepgd,yfilepgdtype,iluout,zfieldin,okey)
  ELSE IF (yfiletype=='BUFFER') THEN
    CALL prep_isba_buffer(ig, u, hprogram,hsurf,iluout,zfieldin)
  ELSE IF (yfiletype=='NETCDF') THEN
    CALL prep_isba_netcdf(dtco, u, hprogram,hsurf,yfile,iluout,zfieldin)
  ELSE
    CALL abor1_sfx('PREP_HOR_ISBA_FIELD: data file type not supported : '//yfiletype)
  END IF
!
  inl = SIZE(zfieldin,2)
  inp = SIZE(zfieldin,3)
!
ENDIF
!
!-------------------------------------------------------------------------------------
!
!*      5.     Horizontal interpolation
!
 CALL prep_ctl_int_part2 (ydctl, hsurf, cmask, 'NATURE', zfieldin)
!
IF (ydctl%LPART3) THEN
!
  IF (nrank==npio) THEN
    inl = SIZE(zfieldin,2)
    inp = SIZE(zfieldin,3)
  ELSEIF (.NOT.ASSOCIATED(zfieldin)) THEN
    ALLOCATE(zfieldin(0,0,0))
  ENDIF
!
  IF (nproc>1) THEN
#ifdef SFX_MPI
    CALL mpi_bcast(inl,kind(inl)/4,mpi_integer,npio,ncomm,infompi)
    CALL mpi_bcast(inp,kind(inp)/4,mpi_integer,npio,ncomm,infompi)
#endif
  ENDIF
  ALLOCATE(zfieldoutp(ini,inl,inp))
!
! ZPATCH is the array of output patches put on the input patches
  ALLOCATE(zpatch(ini,inp))
  zpatch(:,:) = 0.
!
  CALL get_prep_interp(inp,io%NPATCH,s%XVEGTYPE,s%XPATCH,zpatch)
!
  DO jp = 1, inp
  ! we interpolate each point the output patch is present
    linterp(:) = (zpatch(:,jp) > 0.)
    CALL hor_interpol(dtco, u, gcp, iluout,zfieldin(:,:,jp),zfieldoutp(:,:,jp))
    linterp = .true.
  END DO
!
  DEALLOCATE(zfieldin,zpatch)
!
ENDIF
!
 CALL prep_ctl_int_part4 (ydctl, hsurf, 'NATURE', cmask, zfieldin, zfieldoutp)
!
IF (ydctl%LPART5) THEN
!
  inl = SIZE (zfieldoutp,2)
  inp = SIZE (zfieldoutp,3)
!
  IF (trim(hsurf)/="ZS") THEN
  !
    ALLOCATE(zw%AL(io%NPATCH))
  !
    IF (io%NPATCH/=inp) THEN
    !
      ALLOCATE(zfieldoutv(ini,inl,nvegtype))
      CALL put_on_all_vegtypes(ini,inl,inp,nvegtype,zfieldoutp,zfieldoutv)
    !
    !*      6.     Transformation from vegtype grid to patch grid
    !
      DEALLOCATE(zfieldoutp)
    !
      DO jp = 1,io%NPATCH
        pk => np%AL(jp)
        !
        ALLOCATE(zw%AL(jp)%ZOUT(pk%NSIZE_P,inl))
        !
        CALL vegtype_grid_to_patch_grid(jp, io%NPATCH, pk%XVEGTYPE_PATCH, pk%XPATCH,&
                                      pk%NR_P, zfieldoutv, zw%AL(jp)%ZOUT)
      ENDDO
      !
      DEALLOCATE(zfieldoutv)
    !
    ELSE
    !
      DO jp = 1,io%NPATCH
      !
        pk => np%AL(jp)
      !
        ALLOCATE(zw%AL(jp)%ZOUT(pk%NSIZE_P,inl))
      !
        CALL pack_same_rank(pk%NR_P,zfieldoutp(:,:,jp),zw%AL(jp)%ZOUT)
      !
      ENDDO
    !
      DEALLOCATE(zfieldoutp)
    !
    ENDIF
  !
  ENDIF
!
!
!-------------------------------------------------------------------------------------
!
!*      7.     Return to historical variable
!
!
  SELECT CASE (hsurf)
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
    CASE('ZS     ') 
      ALLOCATE(xzs_ls(ini))
      xzs_ls(:) = zfieldoutp(:,1,1)
      DEALLOCATE(zfieldoutp)
    !
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    CASE('WG     ')
    !
      ALLOCATE(zf%AL(io%NPATCH))
    !
      DO jp = 1,io%NPATCH
        kk => nk%AL(jp)
        pk => np%AL(jp)
        pek => npe%AL(jp)
      !
        ALLOCATE(zf%AL(jp)%ZOUT(pk%NSIZE_P,io%NGROUND_LAYER))
      !
      !* interpolates on output levels
        CALL init_from_ref_grid(xgrid_soil,zw%AL(jp)%ZOUT,pk%XDG,zf%AL(jp)%ZOUT)
      !
      !* retrieves soil water content from soil relative humidity
        ALLOCATE(pek%XWG(pk%NSIZE_P,io%NGROUND_LAYER))
        pek%XWG(:,:)=xundef
        IF(io%CISBA=='DIF')THEN
          iwork(1:pk%NSIZE_P)=pk%NWG_LAYER(:)
        ELSE
          iwork(1:pk%NSIZE_P)=SIZE(pek%XWG,2)
        ENDIF
        DO ji=1,pk%NSIZE_P
          IF(iwork(ji)==nundef)cycle
          inl=iwork(ji)
          DO jl=1,inl
            pek%XWG(ji,jl) = kk%XWWILT(ji,jl) + zf%AL(jp)%ZOUT(ji,jl) * (kk%XWFC(ji,jl)-kk%XWWILT(ji,jl))
            pek%XWG(ji,jl) = max(min(pek%XWG(ji,jl),kk%XWSAT(ji,jl)),xwgmin)
          ENDDO
        ENDDO
        !
        WHERE(zf%AL(jp)%ZOUT(:,:)==xundef) pek%XWG(:,:)=xundef
        !
        DEALLOCATE(zf%AL(jp)%ZOUT)
      ENDDO
      !
      DEALLOCATE(zf%AL)
    !
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    CASE('WGI    ')
      ALLOCATE(zf%AL(io%NPATCH))
      !
      DO jp = 1,io%NPATCH
        kk => nk%AL(jp)
        pk => np%AL(jp)
        pek => npe%AL(jp)
        !
        ALLOCATE(zf%AL(jp)%ZOUT(pk%NSIZE_P,io%NGROUND_LAYER))
        !
        !* interpolates on output levels
        CALL init_from_ref_grid(xgrid_soil,zw%AL(jp)%ZOUT,pk%XDG,zf%AL(jp)%ZOUT)       
        !
        !* retrieves soil ice content from soil relative humidity
        ALLOCATE(pek%XWGI(pk%NSIZE_P,io%NGROUND_LAYER))
        pek%XWGI(:,:)=0.0
        IF(io%CISBA=='DIF')THEN
          iwork(1:pk%NSIZE_P)=pk%NWG_LAYER(:)
        ELSE
          iwork(1:pk%NSIZE_P)=2
        ENDIF  
        DO ji=1,pk%NSIZE_P
          IF(iwork(ji)==nundef)cycle
          inl=iwork(ji)
          DO jl=1,inl
            pek%XWGI(ji,jl) = zf%AL(jp)%ZOUT(ji,jl) * kk%XWSAT(ji,jl)
            pek%XWGI(ji,jl) = max(min(pek%XWGI(ji,jl),kk%XWSAT(ji,jl)),0.)
          ENDDO
        END DO
        !
        WHERE(zf%AL(jp)%ZOUT(:,:)==xundef ) pek%XWGI(:,:)=xundef
        WHERE(pek%XWGI(:,:)<=1.0e-10)pek%XWGI(:,:)=0.0
        !
        DEALLOCATE(zf%AL(jp)%ZOUT)
        !
      ENDDO
      !
      DEALLOCATE(zf%AL)
    !
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    CASE('TG     ') 
      IF(io%LTEMP_ARP)THEN
        inl=io%NTEMPLAYER_ARP
      ELSE
        inl=io%NGROUND_LAYER
      ENDIF
      !
      DO jp = 1,io%NPATCH
        !
        pk => np%AL(jp)
        pek => npe%AL(jp)
        !
        ALLOCATE(pek%XTG(pk%NSIZE_P,inl))
        !
        ALLOCATE(zdg(SIZE(pk%XDG,1),inl))
        IF (io%CISBA=='2-L'.OR.io%CISBA=='3-L') THEN
          zdg(:,1) = 0.01
          zdg(:,2) = 0.40                    ! deep temperature for force-restore taken at 20cm
          IF(io%CISBA=='3-L') zdg(:,3) = 5.00   ! climatological temperature, usually not used       
          IF(io%LTEMP_ARP)THEN
            zdg(:,3) = 1.0
            DO jl=4,inl
              zdg(:,jl) = zdg(:,jl-1)+1.0
            ENDDO
          ENDIF
        ELSE
          !* diffusion method, the soil grid is the same as for humidity
          zdg(:,:) = pk%XDG(:,:)
        END IF
        CALL init_from_ref_grid(xgrid_soil,zw%AL(jp)%ZOUT,zdg,pek%XTG)
        DEALLOCATE(zdg)        
        !
      ENDDO
      !
      !
      !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    CASE('WR     ')
      DO jp = 1,io%NPATCH
        ALLOCATE(npe%AL(jp)%XWR(np%AL(jp)%NSIZE_P))
        npe%AL(jp)%XWR(:) = zw%AL(jp)%ZOUT(:,1)
        
      ENDDO
    !
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    CASE('WRL    ') 
      DO jp = 1,io%NPATCH
        ALLOCATE(npe%AL(jp)%XWRL(np%AL(jp)%NSIZE_P))
        npe%AL(jp)%XWRL(:) = zw%AL(jp)%ZOUT(:,1)
      ENDDO
     !
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    CASE('WRLI   ') 
      DO jp = 1,io%NPATCH
        ALLOCATE(npe%AL(jp)%XWRLI(np%AL(jp)%NSIZE_P))
        npe%AL(jp)%XWRLI(:) = zw%AL(jp)%ZOUT(:,1)
      ENDDO
    !
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    CASE('WRVN   ') 
      DO jp = 1,io%NPATCH
        ALLOCATE(npe%AL(jp)%XWRVN(np%AL(jp)%NSIZE_P))
        npe%AL(jp)%XWRVN(:) = zw%AL(jp)%ZOUT(:,1)
      ENDDO
    !
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    CASE('TV     ') 
      DO jp = 1,io%NPATCH
        ALLOCATE(npe%AL(jp)%XTV(np%AL(jp)%NSIZE_P))
        npe%AL(jp)%XTV(:) = zw%AL(jp)%ZOUT(:,1)
      ENDDO
    !
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    CASE('TL     ') 
      DO jp = 1,io%NPATCH
        ALLOCATE(npe%AL(jp)%XTL(np%AL(jp)%NSIZE_P))
        npe%AL(jp)%XTL(:) = zw%AL(jp)%ZOUT(:,1)
      ENDDO
    !
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    CASE('TC     ') 
      DO jp = 1,io%NPATCH
        ALLOCATE(npe%AL(jp)%XTC(np%AL(jp)%NSIZE_P))
        npe%AL(jp)%XTC(:) = zw%AL(jp)%ZOUT(:,1)
      ENDDO
    !
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    CASE('QC     ') 
      DO jp = 1,io%NPATCH
        ALLOCATE(npe%AL(jp)%XQC(np%AL(jp)%NSIZE_P))
        npe%AL(jp)%XQC(:) = zw%AL(jp)%ZOUT(:,1)
      ENDDO
    !
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    CASE('LAI    ') 
    !* LAI is updated only if present and pertinent (evolutive LAI) in input file
      DO jp = 1,io%NPATCH
        IF (any(zw%AL(jp)%ZOUT(:,:)/=xundef)) THEN
          ALLOCATE(npe%AL(jp)%XLAI(np%AL(jp)%NSIZE_P))
          npe%AL(jp)%XLAI(:) = zw%AL(jp)%ZOUT(:,1)
        ENDIF
      ENDDO
    !
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    CASE('ICE_STO') 
      DO jp = 1,io%NPATCH
        ALLOCATE(npe%AL(jp)%XICE_STO(np%AL(jp)%NSIZE_P))
        npe%AL(jp)%XICE_STO(:) = zw%AL(jp)%ZOUT(:,1)
      ENDDO
    !
  END SELECT
  !
  IF (trim(hsurf)/="ZS") THEN
    DO jp = 1,io%NPATCH
      DEALLOCATE(zw%AL(jp)%ZOUT)
    ENDDO
    DEALLOCATE(zw%AL)
  ENDIF
!
ENDIF
!-------------------------------------------------------------------------------------
!
!*      8.     Deallocations
!
IF (lhook) CALL dr_hook('PREP_HOR_ISBA_FIELD',1,zhook_handle)
!
!-------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------
!
CONTAINS
!
!-------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------
!
SUBROUTINE init_from_ref_grid(PGRID1,PT1,PD2,PT2)
!
USE modi_interp_grid_nat
!
REAL, DIMENSION(:,:), INTENT(IN)  :: PT1    ! variable profile
REAL, DIMENSION(:),   INTENT(IN)  :: PGRID1 ! normalized grid
REAL, DIMENSION(:,:), INTENT(IN)  :: PD2    ! output layer thickness
REAL, DIMENSION(:,:), INTENT(OUT) :: PT2    ! variable profile
!
INTEGER                                  :: JI,JL  ! loop counter
REAL, DIMENSION(SIZE(PT1,1),SIZE(PT1,2)) :: ZD1 ! input grid
REAL, DIMENSION(SIZE(PD2,1),SIZE(PD2,2)) :: ZD2 ! output grid
!
INTEGER :: ILAYER1, ILAYER2
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
IF (lhook) CALL dr_hook('INIT_FROM_REF_GRID',0,zhook_handle)
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
IF (SIZE(pt1,2)==3) THEN
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!* 1. case with only 3 input levels (typically coming from 'UNIF')
!     -----------------------------
!
  IF (io%CISBA=='2-L' .OR. io%CISBA=='3-L') THEN

    !* Possible LTEMP_ARP case
    IF(SIZE(pt2,2)>3)THEN
      ilayer1=3
      ilayer2=SIZE(pt2,2)
    ELSE
      ilayer1=SIZE(pt2,2)
      ilayer2=0
    ENDIF
    !* historical 2L or 3L ISBA version
    pt2(:,1:ilayer1) = pt1(:,1:ilayer1) 
    !* Possible LTEMP_ARP case
    IF(ilayer2>0)THEN
      DO jl=ilayer1+1,ilayer2
        pt2(:,jl) = pt2(:,ilayer1)
      ENDDO
    ENDIF

  ELSEIF(io%CISBA=='DIF')THEN

    !surface layer (generally 0.01m imposed)
    pt2(:,1) = pt1(:,1) 
    !second layer
    pt2(:,2) = 0.25*pt1(:,1)+0.75*pt1(:,2)
    !others layers
    DO ji=1,SIZE(pt1,1)
      DO jl=3,io%NGROUND_LAYER
        IF(pd2(ji,jl)<=pk%XDG2(ji))THEN 
          !root layers
          pt2(ji,jl) = pt1(ji,2)
        ELSE
          !deep layers
          pt2(ji,jl) = pt1(ji,3)
        ENDIF
      END DO
    END DO 
  END IF    
!  
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
ELSE
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!* 2. case with fine grid as input (general case)
!     ----------------------------
!
  DO jl=1,SIZE(pt1,2)
     zd1(:,jl) = pgrid1(jl)
  ENDDO
!
  zd2(:,:) = pd2(:,:)
  CALL interp_grid_nat(zd1,pt1(:,:),zd2,pt2(:,:))
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
ENDIF
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
IF (lhook) CALL dr_hook('INIT_FROM_REF_GRID',1,zhook_handle)
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
END SUBROUTINE init_from_ref_grid
!-------------------------------------------------------------------------------------
!
END SUBROUTINE prep_hor_isba_field