subroutine soil_settings_PEM(ngrid,nslope,nsoil_PEM,nsoil_GCM,
     &   TI_GCM,TI_PEM)


!      use netcdf
      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, 
     &   depth_breccia,depth_bedrock,index_breccia,index_bedrock  
      use comsoil_h, only: inertiedat,layer,mlayer, volcapa                       


      use iostart, only: inquire_field_ndims, get_var, get_field,
     &                   inquire_field, inquire_dimension_length

      implicit none

!======================================================================
!  Author: LL, based on work by  Ehouarn Millour (07/2006)
!
!  Purpose: Read and/or initialise soil depths and properties
!
! Modifications: Aug.2010 EM : use NetCDF90 to load variables (enables using
!                      r4 or r8 restarts independently of having compiled
!                      the GCM in r4 or r8)
!                June 2013 TN : Possibility to read files with a time axis
!
!
!  The various actions and variable read/initialized are:
!  1. Read/build layer (and midlayer) depth
!  2. Interpolate thermal inertia and temperature on the new grid, if
!     necessary
!======================================================================

!======================================================================
!  arguments
!  ---------
!  inputs:
      integer,intent(in) :: ngrid	! # of horizontal grid points
      integer,intent(in) :: nslope	! # of subslope wihtin the mesh
      integer,intent(in) :: nsoil_PEM	! # of soil layers in the PEM
      integer,intent(in) :: nsoil_GCM	! # of soil layers in the GCM
      real,intent(in) :: TI_GCM(ngrid,nsoil_GCM,nslope)	!  Thermal inertia  in the GCM [SI]
      real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope)	! Thermal inertia   in the PEM [SI]

!======================================================================
! local variables:
      integer ig,iloop,islope	! loop counters
      logical found

      real alpha,lay1 ! coefficients for building layers
      real xmin,xmax ! to display min and max of a field
      
!======================================================================

! 1. Depth coordinate
! -------------------

! 1.4 Build mlayer(), if necessary
      ! default mlayer distribution, following a power law:
      !  mlayer(k)=lay1*alpha**(k-1/2)
        lay1=2.e-4
	alpha=2
        do iloop=0,nsoil_PEM-1
	  mlayer_PEM(iloop)=lay1*(alpha**(iloop-0.5))
	enddo

! 1.5 Build layer(); following the same law as mlayer()
      ! Assuming layer distribution follows mid-layer law:
      ! layer(k)=lay1*alpha**(k-1)
      lay1=sqrt(mlayer_PEM(0)*mlayer_PEM(1))
      alpha=mlayer_PEM(1)/mlayer_PEM(0)
      do iloop=1,nsoil_PEM
        layer_PEM(iloop)=lay1*(alpha**(iloop-1))
      enddo

! 2. Thermal inertia (note: it is declared in comsoil_h)
! ------------------

      do ig = 1,ngrid
        do islope = 1,nslope
          do iloop = 1,nsoil_GCM
            TI_PEM(ig,iloop,islope) = TI_GCM(ig,iloop,islope)
          enddo
          if(nsoil_PEM.gt.nsoil_GCM) then
           do iloop = nsoil_GCM+1,nsoil_PEM
             TI_PEM(ig,iloop,islope) = TI_GCM(ig,nsoil_GCM,islope)
           enddo
          endif
        enddo
      enddo


! 3. Index for subsurface layering
! ------------------
      index_breccia= 1
      do iloop = 1,nsoil_PEM-1
        if (depth_breccia.ge.layer_PEM(iloop)) then
          index_breccia=iloop 
        else
         exit
        endif
      enddo

      index_bedrock= 1
      do iloop = 1,nsoil_PEM-1
        if (depth_bedrock.ge.layer_PEM(iloop)) then
          index_bedrock=iloop 
        else
         exit
        endif
      enddo


      end subroutine soil_settings_PEM