subroutine soil_settings(nid,ngrid,nsoil,tsurf,tsoil)
      implicit none

!======================================================================
!  Author: Ehouarn Millour (07/2006)
!
!  Purpose: Read and/or initialise soil depths and properties
!
!  This subroutine reads from a NetCDF file (opened by the caller)
!  of "startfi.nc" format.
!  The various actions and variable read/initialized are:
!  1. Check out the number of soil layers (in datafile); if it isn't equal
!     to nsoil, then some interpolation will be required
!     Also check if data in file "startfi.nc" is in older format (ie:
!     thermal inertia was depth-independent; and there was no "depth"
!     coordinate.
!     Read/build layer (and midlayer) depths
!  2. Read volumetric specific heat (or initialise it to default value)
!  3. Read Thermal inertia
!  4. Read soil temperatures
!  5. Interpolate thermal inertia and temperature on the new grid, if
!     necessary
!======================================================================

#include "dimensions.h"
#include "dimphys.h"

#include "comsoil.h"
#include "netcdf.inc"
!======================================================================
!  arguments
!  ---------
!  inputs:
      integer nid	! Input Netcdf file ID 
      integer ngrid	! # of horizontal grid points
      integer nsoil	! # of soil layers
      real tsurf(ngrid) ! surface temperature
!  output:
      real tsoil(ngridmx,nsoilmx)	! soil temperature

!======================================================================
! local variables:
      integer ierr	! status (returned by NetCDF functions)
      integer nvarid	! ID of NetCDF variable
      integer dimid	! ID of NetCDF dimension
      integer dimlen	! length along the "depth" dimension
      integer ndims	! # of dimensions of read <inertiedat> data
      integer ig,iloop	! loop counters

      logical :: olddepthdef=.false. ! flag
      logical :: interpol=.false. ! flag: true if interpolation will be requiered

      ! to store "old" values
      real,dimension(:),allocatable :: surfinertia !surface thermal inertia
      real,dimension(:),allocatable :: oldmlayer
      real,dimension(:,:),allocatable :: oldinertiedat
      real,dimension(:,:),allocatable :: oldtsoil
      
      ! for interpolation
      real,dimension(:),allocatable :: oldgrid
      real,dimension(:),allocatable :: oldval
      real,dimension(:),allocatable :: newval

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

      real,parameter :: default_volcapa=1.e6

!======================================================================

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

! 1.1 Start by reading how many layers of soil there are

	ierr=NF_INQ_DIMID(nid,"subsurface_layers",dimid)
	if (ierr.ne.NF_NOERR) then
	 write(*,*)'soil_settings: Problem reading <subsurface_layers>'
	 call abort
	endif

	ierr=NF_INQ_DIMLEN(nid,dimid,dimlen)
	if (ierr.ne.NF_NOERR) then
	 write(*,*)'soil_settings: Problem getting <subsurface_layers>',
     &             'length'
	 call abort
	endif

        if (dimlen.ne.nsoil) then
          write(*,*)'soil_settings: Interpolation of soil temperature ',
     &              'and thermal inertia will be required!'
	! if dimlen doesn't match nsoil, then interpolation of
	! soil temperatures and thermal inertia will be requiered
	  interpol=.true.
	endif

! 1.2 Find out the # of dimensions <inertiedat> was defined as using
!     (in ye old days, thermal inertia was only given at the "surface")
      ! Look for thermal inertia data
      ierr=NF_INQ_VARID(nid, "inertiedat", nvarid)
      if (ierr.NE.NF_NOERR) then
         write(*,*)'soil_settings: Field <inertiedat> not found!'
         call abort
      endif

      ! Read the # of dimensions <inertidat> was defined as using
      ierr=NF_INQ_VARNDIMS(nid,nvarid,ndims)
      ! if (ndims.eq.1) then we have the "old 2D-surface" format

! 1.3 Read depths values or set olddepthdef flag and values
      if (ndims.eq.1) then ! we know that there is none
        write(*,*)'soil_settings: no <soildepth> field expected'
	write(*,*)'building one mimicking old definitions'
        olddepthdef=.true.
	interpol=.true.
        ! allocate oldmlayer
        if (.not.allocated(oldmlayer)) then
          allocate(oldmlayer(dimlen),stat=ierr)
          if (ierr.ne.0) then
            write(*,*) 'soil_settings: failed allocation of oldmlayer!'
            stop
          endif
        endif
	do iloop=1,dimlen
	  oldmlayer(iloop)=sqrt(887.75/3.14)*((2.**(iloop-0.5))-1.)
	enddo
      else ! Look for depth 
        ierr=NF_INQ_VARID(nid,"soildepth",nvarid)
        if (ierr.ne.NF_NOERR) then
          write(*,*)'soil_settings: Field <soildepth> not found!'
	  call abort
        endif
        ! read <depth> coordinate
        if (interpol) then !put values in oldmlayer
#ifdef NC_DOUBLE
           ierr = NF_GET_VAR_DOUBLE(nid,nvarid,oldmlayer)
#else
           ierr = NF_GET_VAR_REAL(nid,nvarid,oldmlayer)
#endif
          if (ierr.ne.NF_NOERR) then
           write(*,*)'soil_settings: Problem while reading <soildepth>'
           call abort
          endif
        else ! put values in mlayer
#ifdef NC_DOUBLE
           ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayer)
#else
           ierr = NF_GET_VAR_REAL(nid,nvarid,mlayer)
#endif
          if (ierr.ne.NF_NOERR) then
           write(*,*)'soil_settings: Problem while reading <soildepth>'
           call abort
          endif
        endif !of if (interpol)
      endif !of if (ndims.eq.1)

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

! 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(0)*mlayer(1))
      alpha=mlayer(1)/mlayer(0)
      do iloop=1,nsoil
        layer(iloop)=lay1*(alpha**(iloop-1))
      enddo

! 2. Volumetric heat capacity (note: it is declared in comsoil.h)
! ---------------------------
! "volcapa" is (so far) 0D and written in "controle" table of startfi file 
! volcapa is read or set when "controle" is read (see tabfi.F)
! Just in case, we check here that it is not zero. If it is, we
! set it to "default_volcapa"

      if (volcapa.le.0.0) then
        write(*,*)'soil_settings: Warning, volcapa = ',volcapa
	write(*,*)'               That doesn t seem right'
        write(*,*)'        Initializing Volumetric heat capacity to ',
     &             default_volcapa
	volcapa=default_volcapa
      endif
! Look for it
!      ierr=NF_INQ_VARID(nid,"volcapa",nvarid)
!      if (ierr.NE.NF_NOERR) then
!         write(*,*)'soil_settings: Field <volcapa> not found!'
!         write(*,*)'Initializing Volumetric heat capacity to ',
!     &             default_volcapa
!         volcapa=default_volcapa
!      else
!#ifdef NC_DOUBLE
!       ierr = NF_GET_VAR_DOUBLE(nid,nvarid,volcapa)
!#else
!       ierr = NF_GET_VAR_REAL(nid,nvarid,volcapa)
!#endif
!        if (ierr.ne.NF_NOERR) then
!	  write(*,*)'soil_settings: Problem while reading <volcapa>'
!         call abort
!	endif
!      endif

! 3. Thermal inertia (note: it is declared in comsoil.h)
! ------------------

! 3.1 Look (again) for thermal inertia data (to reset nvarid)
      ierr=NF_INQ_VARID(nid, "inertiedat", nvarid)
      if (ierr.NE.NF_NOERR) then
         write(*,*)'soil_settings: Field <inertiedat> not found!'
         call abort
      endif

! 3.2 Knowing the # of dimensions <inertidat> was defined as using,
!     read/build thermal inertia

      if (ndims.eq.1) then ! "old 2D-surface" format
       write(*,*)'soil_settings: Thermal inertia is only given as surfac
     &e data!'
       ! Read Surface thermal inertia
       allocate(surfinertia(ngrid))
#ifdef NC_DOUBLE
       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, surfinertia)
#else
       ierr = NF_GET_VAR_REAL(nid, nvarid, surfinertia)
#endif
        if (ierr.NE.NF_NOERR) then
         write(*,*)'soil_settings: Problem while reading <inertiedat>'
         call abort
        endif
       
       write(*,*)' => Building soil thermal inertia (using reference sur
     &face thermal inertia)'
       do iloop=1,nsoil
         inertiedat(:,iloop)=surfinertia(:)
       enddo
       deallocate(surfinertia)

      else ! "3D surface+depth" format
       if (interpol) then ! put values in oldinertiedat
         if (.not.allocated(oldinertiedat)) then
           allocate(oldinertiedat(ngrid,dimlen),stat=ierr)
           if (ierr.ne.0) then
            write(*,*) 'soil_settings: failed allocation of ',
     &                 'oldinertiedat!'
            stop
           endif
         endif
#ifdef NC_DOUBLE
        ierr = NF_GET_VAR_DOUBLE(nid,nvarid,oldinertiedat)
#else
        ierr = NF_GET_VAR_REAL(nid,nvarid,oldinertiedat)
#endif
        if (ierr.NE.NF_NOERR) then
         write(*,*)'soil_settings: Problem while reading <inertiedat>'
         call abort
        endif
       else ! put values in therm_i
#ifdef NC_DOUBLE
        ierr = NF_GET_VAR_DOUBLE(nid,nvarid,inertiedat)
#else
        ierr = NF_GET_VAR_REAL(nid,nvarid,inertiedat)
#endif
        if (ierr.NE.NF_NOERR) then
         write(*,*)'soil_settings: Problem while reading <inertiedat>'
         call abort
        endif
       endif ! of if (interpol)
      endif ! of if (ndims.eq.1)
      
! 4. Read soil temperatures
! -------------------------

      ierr=NF_INQ_VARID(nid,"tsoil",nvarid)
      if (ierr.ne.NF_NOERR) then
        write(*,*)'soil_settings: Field <tsoil> not found!'
	write(*,*)' => Building <tsoil> from surface values <tsurf>'
	do iloop=1,nsoil
	  tsoil(:,iloop)=tsurf(:)
	enddo
      else ! <tsoil> found
       if (interpol) then ! put values in oldtsoil
         if (.not.allocated(oldtsoil)) then
           allocate(oldtsoil(ngrid,dimlen),stat=ierr)
           if (ierr.ne.0) then
             write(*,*) 'soil_settings: failed allocation of ',
     &                  'oldtsoil!'
             stop
           endif
         endif
#ifdef NC_DOUBLE
        ierr=NF_GET_VAR_DOUBLE(nid,nvarid,oldtsoil)
#else
        ierr=NF_GET_VAR_REAL(nid,nvarid,oldtsoil)
#endif
        if (ierr.ne.NF_NOERR) then
	 write(*,*)'soil_settings: Problem while reading <tsoil>'
	 call abort
	endif
       else ! put values in tsoil
#ifdef NC_DOUBLE
        ierr=NF_GET_VAR_DOUBLE(nid,nvarid,tsoil)
#else
        ierr=NF_GET_VAR_REAL(nid,nvarid,tsoil)
#endif
        if (ierr.ne.NF_NOERR) then
	 write(*,*)'soil_settings: Problem while reading <tsoil>'
	 call abort
	endif
       endif ! of if (interpol) 
      endif

! 5. If necessary, interpolate soil temperatures and thermal inertias
! -------------------------------------------------------------------

      if (olddepthdef) then
      ! tsoil needs to be interpolated, but not therm_i
        allocate(oldgrid(dimlen+1))
        allocate(oldval(dimlen+1))
	allocate(newval(nsoil))

        do ig=1,ngrid
	  ! copy values
	  oldval(1)=tsurf(ig)
	  oldval(2:dimlen+1)=oldtsoil(ig,1:dimlen)
	  ! build vertical coordinate
	  oldgrid(1)=0. ! ground
	  oldgrid(2:dimlen+1)=oldmlayer(1:dimlen)*
     &                (inertiedat(ig,1)/volcapa)
	  ! interpolate
	  call interp_line(oldgrid,oldval,dimlen+1,mlayer,newval,nsoil)
	  ! copy result in tsoil
	  tsoil(ig,:)=newval(:)
	enddo

        ! cleanup
	deallocate(oldgrid)
	deallocate(oldval)
	deallocate(newval)
	interpol=.false. ! no need for interpolation any more      
      endif !of if (olddepthdef)

      if (interpol) then
      write(*,*)'soil_settings: Vertical interpolation along new grid'
      ! interpolate soil temperatures and thermal inertias
        if (.not.allocated(oldgrid)) then
          allocate(oldgrid(dimlen+1))
          allocate(oldval(dimlen+1))
	  allocate(newval(nsoil))
        endif

      ! thermal inertia
        do ig=1,ngrid
	  ! copy data
	  oldval(1:dimlen)=oldinertiedat(ig,dimlen)
	  ! interpolate
	  call interp_line(oldmlayer,oldval,dimlen,mlayer,newval,nsoil)
	  !copy result in inertiedat
	  inertiedat(ig,:)=newval(:)
	enddo
        
      ! soil temperature
        ! vertical coordinate
	oldgrid(1)=0.0
	oldgrid(2:dimlen+1)=oldmlayer(1:dimlen)
        do ig=1,ngrid
	  ! data
	  oldval(1)=tsurf(ig)
	  oldval(2:dimlen+1)=oldtsoil(ig,1:dimlen)
	  ! interpolate
	  call interp_line(oldgrid,oldval,dimlen+1,mlayer,newval,nsoil)
	  ! copy result in inertiedat
	  tsoil(ig,:)=newval(:)
	enddo
	
	!cleanup
        deallocate(oldgrid)
	deallocate(oldval)
	deallocate(newval)
	deallocate(oldinertiedat)
	deallocate(oldtsoil)
      endif ! of if (interpol)
      
! 6. Report min and max values of soil temperatures and thermal inertias
! ----------------------------------------------------------------------

      write(*,*)
      write(*,*)'Soil volumetric heat capacity:',volcapa

      xmin = MINVAL(inertiedat)
      xmax = MAXVAL(inertiedat)
      write(*,*)'Soil thermal inertia <inertiedat>:',xmin,xmax

      xmin = 1.0E+20
      xmax = -1.0E+20
      xmin = MINVAL(tsoil)
      xmax = MAXVAL(tsoil)
      write(*,*)'Soil temperature <tsoil>:',xmin,xmax

      end