!======================================================================
! Program for definition of physical parameters of soil and initial condition
! using observation data in the known observation point

! The measurements of the Falkenberg station have been
! performed at the Falkenberg Boundary Layer Field Site (in German)
! of the Meteorological Observatory Lindenberg - Richard-Aßmann Observatory (MOL-RAO).
! The co-ordinates of the GM Falkenberg are given by:
! 52° 10' 01" N = 52.17° N
! 14° 07' 27" E = 14.12° E
!
!The soil type distribution in the area around Lindenberg is dominated by sandy soils. In the
!forested parts west of Lindenberg, the sand reaches a depth of several
!meters. Dominating soil reference groups are brown soil - Cambic Arenosol, and Ferric
!Podzol. At the GM Falkenberg, sandy soils (pale soil - Eutric Podzoluvisol, brown soil -
!Cambic Arenosol) cover a layer of loam, which can be typically found at a depth of between
!50 cm and 80 cm, locally even below.

program def_soil_ini_soil

implicit none
integer, parameter :: nlon=1, nlat=1, nlevel_soil=22, nst=31, nlevel_snow=10, nlevel_soil_inp_t=14, nlevel_soil_inp_q=7

real, dimension(nlevel_soil) :: level_soil=(/0.02, 0.05, 0.08, 0.10, 0.15, 0.20, 0.30, 0.45, 0.50, &
 0.60, 0.90, 1.00, 1.20, 1.50, 2.21, 3.48, 5.17, 7.34, 10.05, 13.36, &
 17.33, 22.02/), &
 xout, fout

! Input data !!! --->

! 2009/08/23 06:00

real, dimension(nlevel_soil_inp_t) :: &
 ts_inp=(/14.54, 14.37, 15.13, 15.92, 16.73, 17.55, 17.91, 18.03, 17.83, 17.22, &
 16.99, 16.64, 16.07, -9999./), &
 level_soil_inp_t=(/0.00, 0.05, 0.10, 0.15, 0.20, 0.30, 0.45, 0.50, 0.60, 0.90, 1.00, 1.20, 1.50, 13.36/), &
 xinp, finp

real, dimension(nlevel_soil_inp_q) :: &
 qs_inp=(/ 5.80,  5.60,  6.20, 11.47, 18.67, 26.50, -9999./), &
 level_soil_inp_q=(/0.08, 0.15, 0.30, 0.45, 0.60, 0.90, 10.05/)

! depth  8 cm: Max and min values:    33.43       3.87
! depth 15 cm: Max and min values:    30.00       4.70
! depth 30 cm: Max and min values:    25.43       5.40
! depth 45 cm: Max and min values:    26.00       9.90
! depth 60 cm: Max and min values:    30.93      16.43
! depth 90 cm: Max and min values:    35.30      23.80

! <---

real, dimension(nlon,nlat,nst+1) :: soil_fao
real, dimension(nlon,nlat,nlevel_soil) :: soil_qmax, soil_qmin, soil_c, soil_rho, &
 soil_psi, soil_k, soil_par_b, soil_par_c, soil_qrel_wilt, soil_qrel_ref
real, dimension(nlon,nlat) :: soil_water_table, soil_albedo_dry, soil_albedo_wet, &
 soil_emiss_1_dry, soil_emiss_1_wet, soil_emiss_2_dry, soil_emiss_2_wet

real :: t_bottom_0=9.5, water_table_0=10., qrel_bottom_0=1., veg_root_0=0.3, veg_lai_0=1.5, veg_frac_0=0.3, &
 tzer=273.15, valmiss=-9999., val_lai_max=7.5
integer :: index_lev_tbottom=20, index_lev_qbottom=19

real, dimension(nlon,nlat,nlevel_soil) :: ts, qs
real, dimension(nlon,nlat) :: tsurf, qssurf, snow, veg_lai, veg_frac, veg_root
real, dimension(nlon,nlat,0:nlevel_snow) :: snow_lev, snow_t, snow_age, snow_melt_age

integer :: unit_out=31, itype, jlon=1, jlat=1, k
character (len=30) :: file_name_out="observ_initial.dat2"

! Soil type definition
! Cambisols (type 2) + Podzoluvisols (4) 

 soil_fao(:,:,:) = 0.

! itype=2
! soil_fao(jlon,jlat, 1)=float(itype)
! soil_fao(jlon,jlat,itype+1)=0.6
! itype=4
! soil_fao(jlon,jlat,itype+1)=0.4

 itype=2
 soil_fao(jlon,jlat, 1)=float(itype)
 soil_fao(jlon,jlat,itype+1)=1.0 ! see Cambisols in subroutine definition_soil_phys_param

 call definition_soil_phys_param(soil_fao,nlevel_soil,level_soil,nlon,nlat, &
 soil_qmax, soil_qmin, soil_c, soil_rho, soil_psi, soil_k, soil_par_b, soil_par_c, soil_qrel_wilt, soil_qrel_ref, &
 soil_water_table, soil_albedo_dry, soil_albedo_wet, soil_emiss_1_dry, soil_emiss_1_wet, soil_emiss_2_dry, soil_emiss_2_wet)

 print *,' depth  qmax  qmin  qrel_wilt  qrel_ref'
 do k=1,nlevel_soil
  print '(i3,10f8.2)', k, &
 level_soil(k),soil_qmax(jlon,jlat,k), soil_qmin(jlon,jlat,k), soil_qrel_wilt(jlon,jlat,k), soil_qrel_ref(jlon,jlat,k)
 enddo

 print *,' depth  c  rho  psi  k  parb'
 do k=1,nlevel_soil
  print '(i3,f8.2,2f8.0,f9.3,e11.2,f6.2)',k, &
 level_soil(k),soil_c(jlon,jlat,k), soil_rho(jlon,jlat,k), soil_psi(jlon,jlat,k), soil_k(jlon,jlat,k), &
 soil_par_b(jlon,jlat,k)
 enddo

! print *
! print *,soil_water_table

! Vertical interpolation of soil temperature and moisture from input soil levels to model levels

 jlon=1; jlat=1

 xout(1:nlevel_soil) = sqrt(level_soil(1:nlevel_soil))

! Soil temperature

 ts(jlon,jlat,1:nlevel_soil)=t_bottom_0+tzer

 xinp(1:nlevel_soil_inp_t) = sqrt(level_soil_inp_t(1:nlevel_soil_inp_t))
 finp(1:nlevel_soil_inp_t-1) = ts_inp(1:nlevel_soil_inp_t-1)+tzer
 finp(nlevel_soil_inp_t) = ts(jlon,jlat,nlevel_soil)

 call interp(1., 0.0, 0.5, nlevel_soil_inp_t, xinp(1:nlevel_soil_inp_t), finp(1:nlevel_soil_inp_t), &
 xout(1:index_lev_tbottom), fout(1:index_lev_tbottom), index_lev_tbottom)

 ts(jlon,jlat,1:index_lev_tbottom)=fout(1:index_lev_tbottom)

! Soil moisture

!!! qs(jlon,jlat,1:nlevel_soil)=qrel_bottom_0
 qs(jlon,jlat,1:nlevel_soil)=qrel_bottom_0* &
 (soil_qmax(jlon,jlat,index_lev_qbottom)-soil_qmin(jlon,jlat,index_lev_qbottom))+ &
 soil_qmin(jlat,jlon,index_lev_qbottom)

 xinp(1:nlevel_soil_inp_q) = sqrt(level_soil_inp_q(1:nlevel_soil_inp_q))
 finp(1:nlevel_soil_inp_q-1) = qs_inp(1:nlevel_soil_inp_q-1)*1.e-2
 finp(nlevel_soil_inp_q) = qs(jlon,jlat,nlevel_soil)

 call interp(1., 0.0, 0.5, nlevel_soil_inp_q, xinp(1:nlevel_soil_inp_q), finp(1:nlevel_soil_inp_q), &
 xout(1:index_lev_qbottom), fout(1:index_lev_qbottom), index_lev_qbottom)

 qs(jlon,jlat,1:index_lev_qbottom)=fout(1:index_lev_qbottom)

 do k=1,nlevel_soil
   qs(jlon,jlat,k) = max( min( qs(jlon,jlat,k), soil_qmax(jlon,jlat,k)), soil_qmin(jlon,jlat,k))
 enddo

!--------------------------------------------------------------------------

 print *,' depth qs qsrel ts'
 do k=1,nlevel_soil
  print '(i3,10f8.2)', k, &
 level_soil(k),qs(jlon,jlat,k), &
 (qs(jlon,jlat,k)-soil_qmin(jlon,jlat,k))/(soil_qmax(jlon,jlat,k)-soil_qmin(jlon,jlat,k)), &
 ts(jlon,jlat,k)-tzer
 enddo

!--------------------------------------------------------------------------

 jlon=1; jlat=1

 tsurf(jlon,jlat) = ts_inp(1)+tzer
 qssurf(jlon,jlat) = qs(jlon,jlat,1)
 veg_lai(jlon,jlat) = veg_lai_0
 veg_frac(jlon,jlat) = veg_frac_0
 veg_root(jlon,jlat) = veg_root_0
 snow(jlon,jlat) = 0.
 snow_lev(jlon,jlat,:) = valmiss
 snow_t(jlon,jlat,:) = valmiss
 snow_age(jlon,jlat,:) = valmiss
 snow_melt_age(jlon,jlat,:) = valmiss

 open (unit_out, file=file_name_out, form="formatted", status="unknown")

 write (unit_out,101) nlevel_soil, index_lev_tbottom, index_lev_qbottom

 do k=1,nlevel_soil
   write (unit_out,102) &
     level_soil(k), ts(jlon,jlat,k), qs(jlon,jlat,k), &
 soil_qmax(jlon,jlat,k), soil_qmin(jlon,jlat,k), &
 soil_c(jlon,jlat,k), soil_rho(jlon,jlat,k), soil_psi(jlon,jlat,k), &
 soil_k(jlon,jlat,k), soil_par_b(jlon,jlat,k), soil_par_c(jlon,jlat,k), &
 soil_qrel_wilt(jlon,jlat,k), soil_qrel_ref(jlon,jlat,k)
 enddo

 write (unit_out,102) &
 tsurf(jlon,jlat), qssurf(jlon,jlat), snow(jlon,jlat), &
 veg_lai(jlon,jlat), val_lai_max, veg_frac(jlon,jlat), veg_root(jlon,jlat)

 do k=0,nlevel_snow
   write (unit_out,102) &
 snow_lev(jlon,jlat,k), snow_t(jlon,jlat,k), snow_age(jlon,jlat,k), snow_melt_age(jlon,jlat,k)
 enddo

 write (unit_out,102) valmiss 

 close (unit_out)

 101 format (100i3)
 102 format (100e16.8)

stop
end program def_soil_ini_soil

!*********************************************************************

!======================================================================
subroutine definition_soil_phys_param(soil_fao,nlevel_soil,level_soil,nlon,nlat, &
 soil_qmax, soil_qmin, soil_c, soil_rho, soil_psi, soil_k, soil_par_b, soil_par_c, soil_qrel_wilt, soil_qrel_ref, &
 soil_water_table, soil_albedo_dry, soil_albedo_wet, soil_emiss_1_dry, soil_emiss_1_wet, soil_emiss_2_dry, soil_emiss_2_wet)

! Procedure defines physical parameters of soil using soil types (classification FAO)
! defined on the grid points

! Internal basic variables:

!      nst: total number soil soil types;

! Input variables:

!      soil_fao: soil types (classification FAO, 31 types) defined on the grid points, 3D array,
!                1-st and 2-nd array index are grid point index, 3-rd index is soil type number plus 1,
!                if 3-rd index has value 1, then the variable means dominate soil type,
!                values 2-32 of 3-rd index means fraction (of area, proportion) of the each of 31 soil types;

!      nlevel_soil: number of soil levels;
!      level_soil: soil levels values (m);
!      nlon: grid point number along axis of rotated longitude;
!      nlat: grid point number along axis of rotated latidude.

! Output variables:

!      soil_qmax: maximum specific volumetric soil water content (m^3/m^3) at all soil levels, 3D array;
!      soil_qmin: minimum specific volumetric soil water content (m^3/m^3) at all soil levels, 3D array;
!      soil_c: dry soil thermal capacity (J/K/m^3) at all soil levels, 3D array;
!      soil_rho: dry soil density (kg/m^3) at all soil levels, 3D array;
!      soil_psi: idraulic potential of saturated soil (m) at all soil levels, 3D array;
!      soil_k: idraulic conductivity of saturated soil (m/s) at all soil levels, 3D array;
!      soil_par_b: parameter in soil water transport equation (dimensionless) at all soil levels, 3D array;
!      soil_par_c: parameter in soil water transport equation (dimensionless) at all soil levels, 3D array;
!      soil_qrel_wilt: relative soil water content (proportion) at which the vegetation wilt at all soil levels, 3D array;
!      soil_qrel_ref: relative soil water content (proportion) at which the vegetation stop evapotraspiration increase at all soil levels, 3D array;
!      soil_water_table: approximate depth of water table top (m), 2D array;
!      soil_albedo_dry: dry soil albedo (proportion), 2D array;
!      soil_albedo_wet: wet soil albedo (proportion), 2D array;
!      soil_emiss_1_dry: dry soil emissivity in broadband window (proportion), 2D array;
!      soil_emiss_1_wet: wet soil emissivity in broadband window (proportion), 2D array;
!      soil_emiss_2_dry: dry soil emissivity in 8-12 micron window (proportion), 2D array;
!      soil_emiss_2_wet: wet soil emissivity in 8-12 micron window (proportion), 2D array.

implicit none

integer, parameter :: nst=31

! Input:

integer :: nlon, nlat, nlevel_soil
real, dimension(nlon,nlat,nst+1) :: soil_fao
real, dimension(nlevel_soil) :: level_soil

! Ouput:

real, dimension(nlon,nlat,nlevel_soil) :: soil_qmax, soil_qmin, soil_c, soil_rho, &
 soil_psi, soil_k, soil_par_b, soil_par_c, soil_qrel_wilt, soil_qrel_ref
real, dimension(nlon,nlat) :: soil_water_table, soil_albedo_dry, soil_albedo_wet, &
 soil_emiss_1_dry, soil_emiss_1_wet, soil_emiss_2_dry, soil_emiss_2_wet

! Work:

integer, parameter :: soil_horizon_number=5, ival_missing=-9999
real, dimension(soil_horizon_number,nst) :: horizon_bottom, &
 qgmax, qgmin, rogdr, cgdry, psig0, ykghy, bghyd, cghyd, qg_rel_wilt, qg_rel_refw
real, dimension(nst) :: water_table_depth, soildalbed, soilwalbed, &
 soildemis1, soilwemis1, soildemis2, soilwemis2
real, dimension(nlevel_soil,nst) :: qgmax_lev, qgmin_lev, rogdr_lev, cgdry_lev, &
 psig0_lev, ykghy_lev, bghyd_lev, cghyd_lev, qg_rel_wilt_lev, qg_rel_refw_lev
real, dimension(nlon,nlat) :: zzz1, zzz2
integer :: itype, nsmooth, ilev, ihorizon, i, ip, jp, k, n11
real :: zfrac

! Definition of phisical parameters for various soil type at all soil horizons,
! various soil type may have variuos horizon number (maximum 5) with
! various thikness

! Depth (m) of bottom of soil horisons
!                                Horizon
!                   1       2       3       4       5          Soil type
data &
 horizon_bottom/   0.10,   0.40,   0.90,   1.50,   100., &  !   1 Acrisols
!!!                   0.05,   0.30,   0.80,   1.50,   100., &  !   2 Cambisols
                   0.05,   0.30,   0.60,   1.20,   100., &  !   2 Cambisols
                   0.05,   0.80,   1.50,   2.00,   100., &  !   3 Chernozems
                   0.10,   0.30,   0.70,   1.40,   100., &  !   4 Podzoluvisols
                   0.10,   0.50,   100., -9999., -9999., &  !   5 Rendzinans
                   0.70,   2.00,    80.,   100., -9999., &  !   6 Ferrasols
                   0.20,   0.50,   1.50,   100., -9999., &  !   7 Gleysols
                   0.50,   1.00,   2.00,   100., -9999., &  !   8 Phaeozems
                   0.10,   100., -9999., -9999., -9999., &  !   9 Lithosols
                   0.20,   1.50,   100., -9999., -9999., &  !  10 Fluvisols
                   0.20,   0.60,   1.60,   100., -9999., &  !  11 Kastanozems
                   0.05,   0.20,   0.80,   1.50,   100., &  !  12 Luvisols
                   0.20,   0.80,   1.50,   100., -9999., &  !  13 Greyzems
                   0.40,   1.30,   2.00,   100., -9999., &  !  14 Nitosols
                   0.80,   1.20,   2.00,   100., -9999., &  !  15 Histosols
                   0.10,   0.40,   0.80,   1.40,   100., &  !  16 Podzols
                   0.80,   2.00,   100., -9999., -9999., &  !  17 Arenosols
                   0.10,   1.00,   100., -9999., -9999., &  !  18 Regosols
                   0.20,   0.50,   1.00,   100., -9999., &  !  19 Solonetz
                   0.30,   1.50,   100., -9999., -9999., &  !  20 Andosols
                   0.20,   0.50,   100., -9999., -9999., &  !  21 Rankers
                   0.20,   0.50,   1.00,   100., -9999., &  !  22 Vertisols
                   0.10,   2.00,   100., -9999., -9999., &  !  23 Planosols
                   0.50,   0.80,   1.50,   100., -9999., &  !  24 Xerosols
                   0.20,   1.00,   100., -9999., -9999., &  !  25 Yermosols
                   0.30,   1.00,   100., -9999., -9999., &  !  26 Solonchaks
                   2.00,   100., -9999., -9999., -9999., &  !  27 Dunes or shifting sands
                   1.00,   100., -9999., -9999., -9999., &  !  28 Salt flats
                   1.00,   100., -9999., -9999., -9999., &  !  29 Rock debris or desert detritus
                   100., -9999., -9999., -9999., -9999., &  !  30 Glaciers
                   100., -9999., -9999., -9999., -9999./    !  31 Water body
                   

! Maximum volumetric water content of the soil (m**3/m**3) (qgmax*row is soil porosity)
!                                Horizon
!                   1       2       3       4       5          Soil type
data qgmax/       0.600,  0.480,  0.470,  0.490,  0.420, &  !   1 Acrisols
!!!                  0.550,  0.400,  0.450,  0.490,  0.400, &  !   2 Cambisols
                  0.550,  0.370,  0.360,  0.360,  0.400, &  !   2 Cambisols
                  0.550,  0.500,  0.450,  0.400,  0.400, &  !   3 Chernozems
                  0.600,  0.400,  0.450,  0.420,  0.400, &  !   4 Podzoluvisols
                  0.450,  0.400,  0.400, -9999., -9999., &  !   5 Rendzinans
                  0.500,  0.500,  0.450,  0.420, -9999., &  !   6 Ferrasols
                  0.450,  0.400,  0.480,  0.420, -9999., &  !   7 Gleysols
                  0.450,  0.480,  0.490,  0.400, -9999., &  !   8 Phaeozems
                  0.450,  0.420, -9999., -9999., -9999., &  !   9 Lithosols
                  0.500,  0.450,  0.420, -9999., -9999., &  !  10 Fluvisols
                  0.420,  0.450,  0.400,  0.400, -9999., &  !  11 Kastanozems
                  0.550,  0.400,  0.480,  0.450,  0.400, &  !  12 Luvisols
                  0.450,  0.480,  0.490,  0.400, -9999., &  !  13 Greyzems
                  0.480,  0.490,  0.490,  0.420, -9999., &  !  14 Nitosols
                  0.900,  0.490,  0.490,  0.420, -9999., &  !  15 Histosols
                  0.600,  0.400,  0.480,  0.420,  0.400, &  !  16 Podzols
                  0.400,  0.430,  0.400, -9999., -9999., &  !  17 Arenosols
                  0.430,  0.460,  0.400, -9999., -9999., &  !  18 Regosols
                  0.410,  0.470,  0.490,  0.420, -9999., &  !  19 Solonetz
                  0.500,  0.700,  0.400, -9999., -9999., &  !  20 Andosols
                  0.450,  0.400,  0.400, -9999., -9999., &  !  21 Rankers
                  0.350,  0.400,  0.450,  0.420, -9999., &  !  22 Vertisols
                  0.450,  0.470,  0.420, -9999., -9999., &  !  23 Planosols
                  0.450,  0.480,  0.460,  0.400, -9999., &  !  24 Xerosols
                  0.450,  0.420,  0.400, -9999., -9999., &  !  25 Yermosols
                  0.450,  0.480,  0.420, -9999., -9999., &  !  26 Solonchaks
                  0.400,  0.400, -9999., -9999., -9999., &  !  27 Dunes or shifting sands
                  0.420,  0.400, -9999., -9999., -9999., &  !  28 Salt flats
                  0.420,  0.400, -9999., -9999., -9999., &  !  29 Rock debris or desert detritus
                  1.000, -9999., -9999., -9999., -9999., &  !  30 Glaciers
                  1.000, -9999., -9999., -9999., -9999./    !  31 Water body

! Minimum volumetric water content of the soil (m**3/m**3)
!                                Horizon
!                   1       2       3       4       5          Soil type
data qgmin/       0.150,  0.145,  0.140,  0.142,  0.025, &  !   1 Acrisols
!!!                  0.100,  0.030,  0.060,  0.120,  0.020, &  !   2 Cambisols
                  0.005,  0.004,  0.008,  0.016,  0.020, &  !   2 Cambisols
                  0.100,  0.100,  0.080,  0.100,  0.080, &  !   3 Chernozems
                  0.150,  0.050,  0.130,  0.110,  0.050, &  !   4 Podzoluvisols
                  0.050,  0.020,  0.050, -9999., -9999., &  !   5 Rendzinans
                  0.100,  0.150,  0.140,  0.100, -9999., &  !   6 Ferrasols
                  0.040,  0.020,  0.140,  0.100, -9999., &  !   7 Gleysols
                  0.060,  0.100,  0.110,  0.050, -9999., &  !   8 Phaeozems
                  0.030,  0.020, -9999., -9999., -9999., &  !   9 Lithosols
                  0.100,  0.140,  0.025, -9999., -9999., &  !  10 Fluvisols
                  0.080,  0.130,  0.110,  0.050, -9999., &  !  11 Kastanozems
                  0.100,  0.050,  0.130,  0.120,  0.050, &  !  12 Luvisols
                  0.060,  0.100,  0.110,  0.050, -9999., &  !  13 Greyzems
                  0.135,  0.145,  0.140,  0.025, -9999., &  !  14 Nitosols
                  0.100,  0.140,  0.140,  0.025, -9999., &  !  15 Histosols
                  0.150,  0.050,  0.140,  0.110,  0.050, &  !  16 Podzols
                  0.020,  0.120,  0.020, -9999., -9999., &  !  17 Arenosols
                  0.100,  0.120,  0.020, -9999., -9999., &  !  18 Regosols
                  0.050,  0.130,  0.140,  0.050, -9999., &  !  19 Solonetz
                  0.020,  0.010,  0.020, -9999., -9999., &  !  20 Andosols
                  0.050,  0.020,  0.050, -9999., -9999., &  !  21 Rankers
                  0.100,  0.130,  0.140,  0.100, -9999., &  !  22 Vertisols
                  0.080,  0.140,  0.100, -9999., -9999., &  !  23 Planosols
                  0.120,  0.140,  0.130,  0.100, -9999., &  !  24 Xerosols
                  0.020,  0.050,  0.050, -9999., -9999., &  !  25 Yermosols
                  0.120,  0.140,  0.050, -9999., -9999., &  !  26 Solonchaks
                  0.020,  0.100, -9999., -9999., -9999., &  !  27 Dunes or shifting sands
                  0.050,  0.100, -9999., -9999., -9999., &  !  28 Salt flats
                  0.030,  0.100, -9999., -9999., -9999., &  !  29 Rock debris or desert detritus
                  0.000, -9999., -9999., -9999., -9999., &  !  30 Glaciers
                  0.000, -9999., -9999., -9999., -9999./    !  31 Water body

! Dry soil density (kg/m**3)
!                                Horizon
!                   1       2       3       4       5          Soil type
data rogdr/        900.,  1100.,  1400.,  1600.,  2300., &  !   1 Acrisols
!!!                   500.,   900.,  1300.,  1600.,  2300., &  !   2 Cambisols
                   500.,  1600.,  1700.,  1700.,  2300., &  !   2 Cambisols
                   500.,  1000.,  1300.,  1500.,  2300., &  !   3 Chernozems
                   500.,  1200.,  1500.,  1500.,  2300., &  !   4 Podzoluvisols
                   500.,  1500.,  2300., -9999., -9999., &  !   5 Rendzinans
                  1100.,  1300.,  1600.,  2300., -9999., &  !   6 Ferrasols
                  1100.,  1400.,  1600.,  2300., -9999., &  !   7 Gleysols
                  1000.,  1400.,  1600.,  2300., -9999., &  !   8 Phaeozems
                  1200.,  2300., -9999., -9999., -9999., &  !   9 Lithosols
                  1000.,  1400.,  2300., -9999., -9999., &  !  10 Fluvisols
                  1100.,  1300.,  1500.,  2300., -9999., &  !  11 Kastanozems
                   500.,  1100.,  1400.,  1500.,  2300., &  !  12 Luvisols
                  1000.,  1400.,  1600.,  2300., -9999., &  !  13 Greyzems
                  1100.,  1400.,  1600.,  2300., -9999., &  !  14 Nitosols
                   300.,  1400.,  1600.,  2300., -9999., &  !  15 Histosols
                   500.,  1300.,  1600.,  1500.,  2300., &  !  16 Podzols
                  1100.,  1500.,  2300., -9999., -9999., &  !  17 Arenosols
                  1200.,  1600.,  2300., -9999., -9999., &  !  18 Regosols
                  1200.,  1500.,  1600.,  2300., -9999., &  !  19 Solonetz
                   900.,  1000.,  2300., -9999., -9999., &  !  20 Andosols
                  1100.,  1500.,  2300., -9999., -9999., &  !  21 Rankers
                  1100.,  1300.,  1500.,  2300., -9999., &  !  22 Vertisols
                   600.,  1500.,  2300., -9999., -9999., &  !  23 Planosols
                  1100.,  1400.,  1500.,  2300., -9999., &  !  24 Xerosols
                  1100.,  1400.,  2300., -9999., -9999., &  !  25 Yermosols
                  1100.,  1500.,  2300., -9999., -9999., &  !  26 Solonchaks
                  1200.,  2300., -9999., -9999., -9999., &  !  27 Dunes or shifting sands
                  1400.,  2300., -9999., -9999., -9999., &  !  28 Salt flats
                  1400.,  2300., -9999., -9999., -9999., &  !  29 Rock debris or desert detritus
                   915., -9999., -9999., -9999., -9999., &  !  30 Glaciers
                  1000., -9999., -9999., -9999., -9999./    !  31 Water body

! Dry soil heat capacity (J/kg/K)
!                                Horizon
!                   1       2       3       4       5          Soil type
data cgdry/        900.,   900.,   900.,   800.,   800., &  !   1 Acrisols
!!!                   900.,   850.,   900.,   900.,   800., &  !   2 Cambisols
                  1000.,  1320.,  1320.,  1320.,  1300., &  !   2 Cambisols
                   900.,   900.,   850.,   850.,   800., &  !   3 Chernozems
                   900.,   850.,   800.,   800.,   800., &  !   4 Podzoluvisols
                   900.,   850.,   800., -9999., -9999., &  !   5 Rendzinans
                   850.,   900.,   900.,   800., -9999., &  !   6 Ferrasols
                   850.,   750.,   700.,   800., -9999., &  !   7 Gleysols
                   850.,   850.,   700.,   800., -9999., &  !   8 Phaeozems
                   850.,   800., -9999., -9999., -9999., &  !   9 Lithosols
                   850.,   800.,   800., -9999., -9999., &  !  10 Fluvisols
                   850.,   850.,   800.,   800., -9999., &  !  11 Kastanozems
                   900.,   850.,   800.,   800.,   800., &  !  12 Luvisols
                   850.,   850.,   700.,   800., -9999., &  !  13 Greyzems
                   900.,   800.,   800.,   800., -9999., &  !  14 Nitosols
                  1000.,   900.,   900.,   900., -9999., &  !  15 Histosols
                   900.,   850.,   800.,   800.,   800., &  !  16 Podzols
                   850.,   800.,   800., -9999., -9999., &  !  17 Arenosols
                   850.,   800.,   800., -9999., -9999., &  !  18 Regosols
                   850.,   850.,   850.,   800., -9999., &  !  19 Solonetz
                   800.,   800.,   800., -9999., -9999., &  !  20 Andosols
                   850.,   800.,   800., -9999., -9999., &  !  21 Rankers
                   850.,   800.,   800.,   800., -9999., &  !  22 Vertisols
                   900.,   850.,   800., -9999., -9999., &  !  23 Planosols
                   850.,   850.,   800.,   800., -9999., &  !  24 Xerosols
                   800.,   800.,   800., -9999., -9999., &  !  25 Yermosols
                   800.,   800.,   800., -9999., -9999., &  !  26 Solonchaks
                   800.,   800., -9999., -9999., -9999., &  !  27 Dunes or shifting sands
                   800.,   800., -9999., -9999., -9999., &  !  28 Salt flats
                   850.,   800., -9999., -9999., -9999., &  !  29 Rock debris or desert detritus
                 2093.4, -9999., -9999., -9999., -9999., &  !  30 Glaciers
                 4186.8, -9999., -9999., -9999., -9999./    !  31 Water body

! Saturated moisture potential of the soil (at QG=QGMAX) (m)
!                                Horizon
!                   1       2       3       4       5          Soil type
data psig0/      -0.350, -0.400, -0.450, -0.500, -0.300, &  !   1 Acrisols
                 -0.300, -0.400, -0.400, -0.450, -0.300, &  !   2 Cambisols
                 -0.300, -0.250, -0.200, -0.150, -0.200, &  !   3 Chernozems
                 -0.330, -0.130, -0.500, -0.350, -0.100, &  !   4 Podzoluvisols
                 -0.200, -0.150, -0.100, -9999., -9999., &  !   5 Rendzinans
                 -0.300, -0.500, -0.500, -0.300, -9999., &  !   6 Ferrasols
                 -0.200, -0.150, -0.400, -0.300, -9999., &  !   7 Gleysols
                 -0.500, -0.300, -0.150, -0.100, -9999., &  !   8 Phaeozems
                 -0.150, -0.100, -9999., -9999., -9999., &  !   9 Lithosols
                 -0.500, -0.500, -0.300, -9999., -9999., &  !  10 Fluvisols
                 -0.200, -0.150, -0.150, -0.200, -9999., &  !  11 Kastanozems
                 -0.300, -0.300, -0.650, -0.450, -0.300, &  !  12 Luvisols
                 -0.500, -0.300, -0.150, -0.100, -9999., &  !  13 Greyzems
                 -0.300, -0.600, -0.600, -0.300, -9999., &  !  14 Nitosols
                 -0.100, -0.600, -0.600, -0.300, -9999., &  !  15 Histosols
                 -0.350, -0.130, -0.600, -0.350, -0.100, &  !  16 Podzols
                 -0.100, -0.150, -0.300, -9999., -9999., &  !  17 Arenosols
                 -0.150, -0.200, -0.300, -9999., -9999., &  !  18 Regosols
                 -0.250, -0.350, -0.450, -0.300, -9999., &  !  19 Solonetz
                 -0.050, -0.020, -0.050, -9999., -9999., &  !  20 Andosols
                 -0.200, -0.150, -0.100, -9999., -9999., &  !  21 Rankers
                 -0.400, -0.500, -0.500, -0.300, -9999., &  !  22 Vertisols
                 -0.200, -0.600, -0.300, -9999., -9999., &  !  23 Planosols
                 -0.300, -0.400, -0.350, -0.300, -9999., &  !  24 Xerosols
                 -0.100, -0.080, -0.100, -9999., -9999., &  !  25 Yermosols
                 -0.300, -0.350, -0.300, -9999., -9999., &  !  26 Solonchaks
                 -0.100, -0.300, -9999., -9999., -9999., &  !  27 Dunes or shifting sands
                 -0.200, -0.300, -9999., -9999., -9999., &  !  28 Salt flats
                 -0.100, -0.300, -9999., -9999., -9999., &  !  29 Rock debris or desert detritus
                   2.20, -9999., -9999., -9999., -9999., &  !  30 Glaciers
                   0.57, -9999., -9999., -9999., -9999./    !  31 Water body

! Saturated hydraulic conductivity (m/s)
!                                Horizon
!                   1       2        3       4        5          Soil type
data ykghy/      2.0e-7, 2.5e-7,  3.0e-7, 4.5e-7, 0.5e-7, &  !   1 Acrisols
!!!                 5.0e-7, 7.5e-7,  8.5e-7,10.0e-7, 1.0e-7, &  !   2 Cambisols
                 5.0e-7,12.8e-7,  9.3e-7, 2.3e-7, 1.0e-7, &  !   2 Cambisols
                 5.0e-7, 9.0e-7, 10.0e-7, 2.0e-7, 0.5e-7, &  !   3 Chernozems
                 3.0e-7, 8.0e-7, 12.5e-7, 1.5e-7, 0.5e-7, &  !   4 Podzoluvisols
                 4.5e-7, 5.0e-7,  1.0e-7, -9999., -9999., &  !   5 Rendzinans
                 1.5e-7, 6.5e-7,  1.5e-7, 1.0e-7, -9999., &  !   6 Ferrasols
                 4.5e-7,11.0e-7,  3.5e-7, 1.0e-7, -9999., &  !   7 Gleysols
                 7.0e-7,10.0e-7,  5.0e-7, 1.0e-7, -9999., &  !   8 Phaeozems
                12.0e-7, 1.0e-7,  -9999., -9999., -9999., &  !   9 Lithosols
                 1.2e-7, 1.5e-7,  1.0e-7, -9999., -9999., &  !  10 Fluvisols
                 8.0e-7,11.0e-7,  4.0e-7, 0.5e-7, -9999., &  !  11 Kastanozems
                 5.0e-7, 6.0e-7,  8.0e-7,10.0e-7, 1.0e-7, &  !  12 Luvisols
                 7.0e-7,10.0e-7,  5.0e-7, 1.0e-7, -9999., &  !  13 Greyzems
                 1.0e-7, 2.0e-7,  2.5e-7, 0.5e-7, -9999., &  !  14 Nitosols
                 6.0e-7, 2.0e-7,  2.5e-7, 0.5e-7, -9999., &  !  15 Histosols
                 3.0e-7, 8.0e-7, 10.0e-7, 1.5e-7, 0.5e-7, &  !  16 Podzols
                12.0e-7, 3.5e-7,  1.0e-7, -9999., -9999., &  !  17 Arenosols
                 6.0e-7, 11.e-7,  1.0e-7, -9999., -9999., &  !  18 Regosols
                 2.0e-7, 4.5e-7,  6.0e-7, 1.0e-7, -9999., &  !  19 Solonetz
                12.0e-7,20.0e-7,  2.0e-7, -9999., -9999., &  !  20 Andosols
                 8.0e-7,10.0e-7,  1.0e-7, -9999., -9999., &  !  21 Rankers
                 3.0e-7, 6.0e-7, 10.0e-7, 1.0e-7, -9999., &  !  22 Vertisols
                 2.0e-7, 3.0e-7,  1.0e-7, -9999., -9999., &  !  23 Planosols
                 8.0e-7,12.5e-7,  1.5e-7, 0.5e-7, -9999., &  !  24 Xerosols
                 5.5e-7,11.0e-7,  0.5e-7, -9999., -9999., &  !  25 Yermosols
                 3.0e-7, 4.0e-7,  1.0e-7, -9999., -9999., &  !  26 Solonchaks
                15.0e-7, 1.0e-7,  -9999., -9999., -9999., &  !  27 Dunes or shifting sands
                 1.0e-7, 1.0e-7,  -9999., -9999., -9999., &  !  28 Salt flats
                12.0e-7, 1.0e-7,  -9999., -9999., -9999., &  !  29 Rock debris or desert detritus
                     0., -9999.,  -9999., -9999., -9999., &  !  30 Glaciers
                     0., -9999.,  -9999., -9999., -9999./    !  31 Water body

! Exponent B (dimensionless) for moisture potential equation:
! PSI=PSIG*(QGMAX/QG)**B
!                                Horizon
!                   1       2       3       4       5          Soil type
data bghyd/        8.00,   9.50,   8.00,  11.00,   3.00, &  !   1 Acrisols
                   5.00,   4.50,   7.50,   9.00,   3.00, &  !   2 Cambisols
                   5.00,   5.00,   4.50,   4.00,   3.00, &  !   3 Chernozems
                   5.00,   4.50,   8.50,   7.50,   3.00, &  !   4 Podzoluvisols
                   5.00,   4.10,   3.00, -9999., -9999., &  !   5 Rendzinans
                   5.00,   9.00,  10.00,   5.00, -9999., &  !   6 Ferrasols
                   4.90,   4.20,  11.00,   5.00, -9999., &  !   7 Gleysols
                   5.30,   4.90,   4.10,   3.00, -9999., &  !   8 Phaeozems
                   4.50,   3.00, -9999., -9999., -9999., &  !   9 Lithosols
                   5.00,  10.50,   3.00, -9999., -9999., &  !  10 Fluvisols
                   4.70,   4.50,   4.00,   3.00, -9999., &  !  11 Kastanozems
                   5.00,   6.00,   8.00,   4.00,   3.00, &  !  12 Luvisols
                   5.30,   4.90,   4.10,   3.00, -9999., &  !  13 Greyzems
                   8.50,  11.50,  11.00,   5.00, -9999., &  !  14 Nitosols
                   4.00,  11.50,  11.00,   5.00, -9999., &  !  15 Histosols
                   5.00,   4.50,  10.50,   7.50,    3.0, &  !  16 Podzols
                   4.20,  10.00,   3.00, -9999., -9999., &  !  17 Arenosols
                   6.00,   7.00,   3.00, -9999., -9999., &  !  18 Regosols
                   7.00,  10.00,  11.00,     3., -9999., &  !  19 Solonetz
                   4.10,   4.00,   3.00, -9999., -9999., &  !  20 Andosols
                   4.50,   4.00,   3.00, -9999., -9999., &  !  21 Rankers
                   7.50,   8.00,   8.00,   3.00, -9999., &  !  22 Vertisols
                   5.50,  11.00,   3.00, -9999., -9999., &  !  23 Planosols
                   8.00,   9.00,   8.50,   3.00, -9999., &  !  24 Xerosols
                   4.00,   4.50,   3.00, -9999., -9999., &  !  25 Yermosols
                   8.00,  10.00,   3.00, -9999., -9999., &  !  26 Solonchaks
                   4.00,   3.00, -9999., -9999., -9999., &  !  27 Dunes or shifting sands
                   6.00,   3.00, -9999., -9999., -9999., &  !  28 Salt flats
                   4.50,   3.00, -9999., -9999., -9999., &  !  29 Rock debris or desert detritus
                   0.00, -9999., -9999., -9999., -9999., &  !  30 Glaciers
                   0.00, -9999., -9999., -9999., -9999./    !  31 Water body

! Exponent C (dimensionless) for ... equation:
! 
!                                Horizon
!                   1       2       3       4       5          Soil type
data cghyd/        3.00,   3.00,   3.00,   3.00,   3.00, &  !   1 Acrisols
                   3.00,   3.00,   3.00,   3.00,   3.00, &  !   2 Cambisols
                   3.00,   3.00,   3.00,   3.00,   3.00, &  !   3 Chernozems
                   3.00,   3.00,   3.00,   3.00,   3.00, &  !   4 Podzoluvisols
                   3.00,   3.00,   3.00, -9999., -9999., &  !   5 Rendzinans
                   3.00,   3.00,   3.00,   3.00, -9999., &  !   6 Ferrasols
                   3.00,   3.00,   3.00,   3.00, -9999., &  !   7 Gleysols
                   3.00,   3.00,   3.00,   3.00, -9999., &  !   8 Phaeozems
                   3.00,   3.00, -9999., -9999., -9999., &  !   9 Lithosols
                   3.00,   3.00,   3.00, -9999., -9999., &  !  10 Fluvisols
                   3.00,   3.00,   3.00,   3.00, -9999., &  !  11 Kastanozems
                   3.00,   3.00,   3.00,   3.00,   3.00, &  !  12 Luvisols
                   3.00,   3.00,   3.00,   3.00, -9999., &  !  13 Greyzems
                   3.00,   3.00,   3.00,   3.00, -9999., &  !  14 Nitosols
                   3.00,   3.00,   3.00,   3.00, -9999., &  !  15 Histosols
                   3.00,   3.00,   3.00,   3.00,   3.00, &  !  16 Podzols
                   3.00,   3.00,   3.00, -9999., -9999., &  !  17 Arenosols
                   3.00,   3.00,   3.00, -9999., -9999., &  !  18 Regosols
                   3.00,   3.00,   3.00,   3.00, -9999., &  !  19 Solonetz
                   3.00,   3.00,   3.00, -9999., -9999., &  !  20 Andosols
                   3.00,   3.00,   3.00, -9999., -9999., &  !  21 Rankers
                   3.00,   3.00,   3.00,   3.00, -9999., &  !  22 Vertisols
                   3.00,   3.00,   3.00, -9999., -9999., &  !  23 Planosols
                   3.00,   3.00,   3.00,   3.00, -9999., &  !  24 Xerosols
                   3.00,   3.00,   3.00, -9999., -9999., &  !  25 Yermosols
                   3.00,   3.00,   3.00, -9999., -9999., &  !  26 Solonchaks
                   3.00,   3.00, -9999., -9999., -9999., &  !  27 Dunes or shifting sands
                   3.00,   3.00, -9999., -9999., -9999., &  !  28 Salt flats
                   3.00,   3.00, -9999., -9999., -9999., &  !  29 Rock debris or desert detritus
                   0.00, -9999., -9999., -9999., -9999., &  !  30 Glaciers
                   0.00, -9999., -9999., -9999., -9999./    !  31 Water body

! Relative (0-1) water content of the soil at wilting point (proportion)
!                                Horizon
!                   1       2       3       4       5          Soil type
data &
 qg_rel_wilt/      0.45,   0.42,   0.35,   0.40,   0.10, &  !   1 Acrisols
!!!                   0.40,   0.30,   0.35,   0.40,   0.10, &  !   2 Cambisols
                   0.08,   0.04,   0.03,   0.10,   0.10, &  !   2 Cambisols
                   0.40,   0.20,   0.15,   0.15,   0.10, &  !   3 Chernozems
                   0.30,   0.17,   0.30,   0.25,   0.10, &  !   4 Podzoluvisols
                   0.20,   0.10,   0.10, -9999., -9999., &  !   5 Rendzinans
                   0.25,   0.40,   0.50,   0.20, -9999., &  !   6 Ferrasols
                   0.18,   0.15,   0.40,   0.20, -9999., &  !   7 Gleysols
                   0.25,   0.20,   0.10,   0.10, -9999., &  !   8 Phaeozems
                   0.15,   0.10, -9999., -9999., -9999., &  !   9 Lithosols
                   0.22,   0.40,   0.10, -9999., -9999., &  !  10 Fluvisols
                   0.17,   0.30,   0.25,   0.10, -9999., &  !  11 Kastanozems
                   0.40,   0.20,   0.35,   0.30,   0.10, &  !  12 Luvisols
                   0.25,   0.20,   0.10,   0.10, -9999., &  !  13 Greyzems
                   0.35,   0.40,   0.40,   0.20, -9999., &  !  14 Nitosols
                   0.50,   0.40,   0.40,   0.20, -9999., &  !  15 Histosols
                   0.30,   0.17,   0.40,   0.25,   0.10, &  !  16 Podzols
                   0.15,   0.30,   0.10, -9999., -9999., &  !  17 Arenosols
                   0.20,   0.30,   0.10, -9999., -9999., &  !  18 Regosols
                   0.25,   0.30,   0.40,   0.10, -9999., &  !  19 Solonetz
                   0.10,   0.10,   0.10, -9999., -9999., &  !  20 Andosols
                   0.20,   0.15,   0.10, -9999., -9999., &  !  21 Rankers
                   0.35,   0.40,   0.40,   0.10, -9999., &  !  22 Vertisols
                   0.20,   0.40,   0.10, -9999., -9999., &  !  23 Planosols
                   0.30,   0.40,   0.35,   0.10, -9999., &  !  24 Xerosols
                   0.10,   0.15,   0.10, -9999., -9999., &  !  25 Yermosols
                   0.30,   0.40,   0.10, -9999., -9999., &  !  26 Solonchaks
                   0.10,   0.10, -9999., -9999., -9999., &  !  27 Dunes or shifting sands
                   0.60,   0.10, -9999., -9999., -9999., &  !  28 Salt flats
                   0.15,   0.10, -9999., -9999., -9999., &  !  29 Rock debris or desert detritus
                   0.00, -9999., -9999., -9999., -9999., &  !  30 Glaciers
                   0.00, -9999., -9999., -9999., -9999./    !  31 Water body

! Relative (0-1) water content of soil at reference point for evapotranspiration (proportion)
!                                Horizon
!                   1       2       3       4       5          Soil type
data &
 qg_rel_refw/      0.70,   0.70,   0.75,   0.68,   0.65, &  !   1 Acrisols
                   0.70,   0.70,   0.75,   0.68,   0.65, &  !   2 Cambisols
                   0.70,   0.80,   0.65,   0.68,   0.65, &  !   3 Chernozems
                   0.68,   0.67,   0.65,   0.68,   0.65, &  !   4 Podzoluvisols
                   0.70,   0.67,   0.65, -9999., -9999., &  !   5 Rendzinans
                   0.80,   0.70,   0.70,   0.65, -9999., &  !   6 Ferrasols
                   0.70,   0.68,   0.67,   0.65, -9999., &  !   7 Gleysols
                   0.70,   0.68,   0.67,   0.65, -9999., &  !   8 Phaeozems
                   0.67,   0.65, -9999., -9999., -9999., &  !   9 Lithosols
                   0.80,   0.68,   0.65, -9999., -9999., &  !  10 Fluvisols
                   0.70,   0.65,   0.68,   0.65, -9999., &  !  11 Kastanozems
                   0.70,   0.70,   0.70,   0.68,   0.65, &  !  12 Luvisols
                   0.70,   0.68,   0.67,   0.65, -9999., &  !  13 Greyzems
                   0.75,   0.68,   0.67,   0.65, -9999., &  !  14 Nitosols
                   0.65,   0.68,   0.67,   0.65, -9999., &  !  15 Histosols
                   0.68,   0.67,   0.65,   0.68,   0.65, &  !  16 Podzols
                   0.67,   0.70,   0.65, -9999., -9999., &  !  17 Arenosols
                   0.67,   0.66,   0.65, -9999., -9999., &  !  18 Regosols
                   0.68,   0.67,   0.66,   0.65, -9999., &  !  19 Solonetz
                   0.66,   0.65,   0.65, -9999., -9999., &  !  20 Andosols
                   0.68,   0.67,   0.65, -9999., -9999., &  !  21 Rankers
                   0.72,   0.70,   0.68,   0.65, -9999., &  !  22 Vertisols
                   0.68,   0.70,   0.65, -9999., -9999., &  !  23 Planosols
                   0.75,   0.68,   0.67,   0.65, -9999., &  !  24 Xerosols
                   0.67,   0.66,   0.65, -9999., -9999., &  !  25 Yermosols
                   0.68,   0.67,   0.65, -9999., -9999., &  !  26 Solonchaks
                   0.67,   0.65, -9999., -9999., -9999., &  !  27 Dunes or shifting sands
                   0.65,   0.65, -9999., -9999., -9999., &  !  28 Salt flats
                   0.67,   0.65, -9999., -9999., -9999., &  !  29 Rock debris or desert detritus
                   0.00, -9999., -9999., -9999., -9999., &  !  30 Glaciers
                   0.00, -9999., -9999., -9999., -9999./    !  31 Water body

! Dry soil (qg=qgmin) albedo
!                                Soil type
!                 Acrisols    Cambisols   Chernozems   Podzoluvisols Rendzinans
data soildalbed/    0.20,        0.20,       0.10,         0.20,       0.10, & 
!                Ferrasols    Gleysols    Phaeozems    Lithosols     Fluvisols
!!!                    0.15,        0.20,       0.10,         0.22,       0.15, & 
                    0.15,        0.20,       0.10,         0.15,       0.15, & 
!                Kastanozems  Luvisols    Greyzems     Nitosols      Histosols 
                    0.15,        0.20,       0.10,         0.18,       0.15, & 
!                Podzols      Arenosols   Regosols     Solonetz      Andosols
                    0.20,        0.25,       0.25,         0.17,       0.05, & 
!                Rankers      Vertisols   Planosols    Xerosols      Yermosols
!!!                    0.20,        0.17,       0.20,         0.24,       0.22, & 
                    0.20,        0.17,       0.20,         0.17,       0.16, & 
!                Solonchaks   Dunes       Salt flats   Rock debris   Glaciers
!!!                    0.22,        0.30,       0.35,         0.25,       0.30, &
                    0.20,        0.20,       0.25,         0.20,       0.30, &
!                Water body
                    0.07/

! Wet soil (qg=qgmax) albedo
!                                Soil type
!                 Acrisols    Cambisols   Chernozems   Podzoluvisols Rendzinans
data soilwalbed/    0.10,        0.10,       0.05,         0.10,       0.05, & 
!                Ferrasols    Gleysols    Phaeozems    Lithosols     Fluvisols
!!!                    0.08,        0.13,       0.05,         0.14,       0.08, & 
                    0.08,        0.13,       0.05,         0.08,       0.08, & 
!                Kastanozems  Luvisols    Greyzems     Nitosols      Histosols 
                    0.08,        0.10,       0.05,         0.10,       0.05, & 
!                Podzols      Arenosols   Regosols     Solonetz      Andosols
                    0.10,        0.15,       0.15,         0.14,       0.03, & 
!                Rankers      Vertisols   Planosols    Xerosols      Yermosols
!!!                    0.10,        0.13,       0.10,         0.14,       0.12, & 
                    0.10,        0.13,       0.10,         0.08,       0.07, & 
!                Solonchaks   Dunes       Salt flats   Rock debris   Glaciers
!!!                    0.18,        0.18,       0.20,         0.18,       0.30, &
                    0.18,        0.13,       0.20,         0.18,       0.30, &
!                Water body
                    0.07/

! Dry soil (qg=qgmin) emissivity in broadband window
!                                Soil type
!                 Acrisols    Cambisols   Chernozems   Podzoluvisols Rendzinans
data soildemis1/   0.960,       0.960,      0.975,        0.960,      0.970, & 
!                Ferrasols    Gleysols    Phaeozems    Lithosols     Fluvisols
!!!                   0.965,       0.960,      0.975,        0.940,      0.965, & 
                   0.965,       0.960,      0.975,        0.960,      0.965, & 
!                Kastanozems  Luvisols    Greyzems     Nitosols      Histosols 
                   0.970,       0.960,      0.975,        0.970,      0.970, & 
!                Podzols      Arenosols   Regosols     Solonetz      Andosols
                   0.960,       0.940,      0.940,        0.970,      0.980, & 
!                Rankers      Vertisols   Planosols    Xerosols      Yermosols
!!!                   0.960,       0.960,      0.960,        0.950,      0.960, & 
                   0.960,       0.960,      0.960,        0.960,      0.965, & 
!                Solonchaks   Dunes       Salt flats   Rock debris   Glaciers
!!!                   0.950,       0.930,      0.940,        0.950,      0.970, &
                   0.950,       0.945,      0.945,        0.960,      0.970, &
!                Water body
                   0.970/

! Wet soil (qg=qgmax) emissivity in broadband window
!                                Soil type
!                 Acrisols    Cambisols   Chernozems   Podzoluvisols Rendzinans
data soilwemis1/   0.985,       0.985,      0.995,        0.985,      0.990, & 
!                Ferrasols    Gleysols    Phaeozems    Lithosols     Fluvisols
!!!                   0.985,       0.985,      0.995,        0.980,      0.985, & 
                   0.985,       0.985,      0.995,        0.985,      0.985, & 
!                Kastanozems  Luvisols    Greyzems     Nitosols      Histosols 
                   0.990,       0.985,      0.995,        0.990,      0.990, & 
!                Podzols      Arenosols   Regosols     Solonetz      Andosols
                   0.985,       0.970,      0.970,        0.990,      0.995, & 
!                Rankers      Vertisols   Planosols    Xerosols      Yermosols
!!!                   0.985,       0.980,      0.985,        0.975,      0.980, & 
                   0.985,       0.980,      0.985,        0.980,      0.985, & 
!                Solonchaks   Dunes       Salt flats   Rock debris   Glaciers
!!!                   0.975,       0.960,      0.980,        0.960,      0.970, &
                   0.980,       0.975,      0.975,        0.970,      0.970, &
!                Water body
                   0.970/

! Dry soil (qg=qgmin) emissivity in 8-12 micron window
!                                Soil type
!                 Acrisols    Cambisols   Chernozems   Podzoluvisols Rendzinans
data soildemis2/   0.950,       0.950,      0.960,        0.950,      0.955, & 
!                Ferrasols    Gleysols    Phaeozems    Lithosols     Fluvisols
!!!                   0.950,       0.950,      0.960,        0.935,      0.950, & 
                   0.950,       0.950,      0.960,        0.945,      0.950, & 
!                Kastanozems  Luvisols    Greyzems     Nitosols      Histosols 
                   0.955,       0.950,      0.960,        0.955,      0.960, & 
!                Podzols      Arenosols   Regosols     Solonetz      Andosols
                   0.950,       0.920,      0.940,        0.960,      0.970, & 
!                Rankers      Vertisols   Planosols    Xerosols      Yermosols
!!!                   0.950,       0.940,      0.950,        0.930,      0.940, & 
                   0.950,       0.940,      0.950,        0.945,      0.950, & 
!                Solonchaks   Dunes       Salt flats   Rock debris   Glaciers
!!!                   0.940,       0.890,      0.920,        0.920,      0.980, &
                   0.940,       0.930,      0.925,        0.930,      0.980, &
!                Water body
                   0.970/

! Wet soil (qg=qgmax) emissivity in 8-12 micron window
!                                Soil type
!                 Acrisols    Cambisols   Chernozems   Podzoluvisols Rendzinans
data soilwemis2/   0.970,       0.970,      0.985,        0.955,      0.990, & 
!                Ferrasols    Gleysols    Phaeozems    Lithosols     Fluvisols
!!!                   0.970,       0.955,      0.985,        0.940,      0.970, & 
                   0.970,       0.955,      0.985,        0.965,      0.970, & 
!                Kastanozems  Luvisols    Greyzems     Nitosols      Histosols 
                   0.980,       0.970,      0.985,        0.975,      0.985, & 
!                Podzols      Arenosols   Regosols     Solonetz      Andosols
                   0.970,       0.940,      0.960,        0.980,      0.990, & 
!                Rankers      Vertisols   Planosols    Xerosols      Yermosols
!!!                   0.970,       0.950,      0.970,        0.950,      0.960, & 
                   0.970,       0.950,      0.970,        0.960,      0.970, & 
!                Solonchaks   Dunes       Salt flats   Rock debris   Glaciers
!!!                   0.960,       0.920,      0.950,        0.930,      0.980, &
                   0.960,       0.950,      0.955,        0.960,      0.980, &
!                Water body
                   0.970/

! Approximate depth (m) of water table top
!                                Soil type
!                 Acrisols    Cambisols   Chernozems   Podzoluvisols Rendzinans
data &
!!! water_table_depth/  10.,          4.,        20.,           2.,        50., & 
 water_table_depth/  10.,         20.,        25.,           2.,        50., & 
!                Ferrasols    Gleysols    Phaeozems    Lithosols     Fluvisols
!!!                     20.,          2.,        20.,         100.,         2., & 
                     20.,          5.,        20.,         100.,         5., & 
!                Kastanozems  Luvisols    Greyzems     Nitosols      Histosols 
!!!                     20.,          4.,        10.,          15.,        0.5, & 
                     20.,         20.,        10.,          15.,        0.5, & 
!                Podzols      Arenosols   Regosols     Solonetz      Andosols
                      3.,         50.,       100.,          50.,       100., & 
!                Rankers      Vertisols   Planosols    Xerosols      Yermosols
                     10.,         10.,        1.5,          50.,        50., & 
!                Solonchaks   Dunes       Salt flats   Rock debris   Glaciers
                     50.,        100.,        50.,          50.,         0., &
!                Water body
                      0./

 print *
 print *,"Definition of physical parameters of soil (FAO soil classification), parameters change depending on soil level"

 nsmooth = 0
! nsmooth = 1
 if (nsmooth > 0) print*, "nsmooth in presoil", nsmooth

! Definition of physical parameters of soil

! Definition of soil parameters at soil levels for various soil types using data
! about soil horizons location

! Correction of hydraulic conductivity following the results of experiment:
! real atmospheric condition in arid zone (Greece), no (!) presipitation
! during 3 summer month, it is necessary reduce hydraulic conductivity in 10 (!)
! times to obtain the decreasing of top soil water content

 do itype=1,nst
   do ihorizon=1,soil_horizon_number
     if (int(ykghy(ihorizon,itype)) /= -9999) ykghy(ihorizon,itype)=ykghy(ihorizon,itype)*0.1
   enddo
 enddo

 do itype=1,nst
   do ilev=1,nlevel_soil
     do ihorizon=1,soil_horizon_number
       if (level_soil(ilev)-horizon_bottom(ihorizon,itype) <= 0.) exit
     enddo
     qgmax_lev(ilev,itype) = qgmax(ihorizon,itype)
     qgmin_lev(ilev,itype) = qgmin(ihorizon,itype)
     rogdr_lev(ilev,itype) = rogdr(ihorizon,itype)
     cgdry_lev(ilev,itype) = cgdry(ihorizon,itype)
     psig0_lev(ilev,itype) = psig0(ihorizon,itype)
     ykghy_lev(ilev,itype) = ykghy(ihorizon,itype)
     bghyd_lev(ilev,itype) = bghyd(ihorizon,itype)
     cghyd_lev(ilev,itype) = cghyd(ihorizon,itype)
     qg_rel_wilt_lev(ilev,itype) = qg_rel_wilt(ihorizon,itype)
     qg_rel_refw_lev(ilev,itype) = qg_rel_refw(ihorizon,itype)
   enddo
 enddo

 soil_qmax(:,:,:) = 0.
 soil_qmin(:,:,:) = 0.
 soil_c(:,:,:) = 0.
 soil_rho(:,:,:) = 0.
 soil_psi(:,:,:) = 0.
 soil_k(:,:,:) = 0.
 soil_par_b(:,:,:) = 0.
 soil_par_c(:,:,:) = 0.
 soil_qrel_wilt(:,:,:) = 0.
 soil_qrel_ref(:,:,:) = 0.
 soil_water_table(:,:) = 0.
 soil_albedo_dry(:,:) = 0.
 soil_albedo_wet(:,:) = 0.
 soil_emiss_1_dry(:,:) = 0.
 soil_emiss_1_wet(:,:) = 0.
 soil_emiss_2_dry(:,:) = 0.
 soil_emiss_2_wet(:,:) = 0.

 do jp=1,nlat
 do ip=1,nlon

   do itype=1,nst
     n11 = itype
     zfrac = soil_fao(ip,jp,n11+1)

! Weighted (with soil fraction) average of the various soil parameters.
! The following fills glacier and water points by soil type "Podzols", for subsequent spatial smoothing only

      if (n11 >= nst-1.and.nsmooth /= 0) n11 = 16

      do ilev=1,nlevel_soil
        soil_qmax(ip,jp,ilev) = soil_qmax(ip,jp,ilev) + qgmax_lev(ilev,n11)*zfrac
        soil_qmin(ip,jp,ilev) = soil_qmin(ip,jp,ilev) + qgmin_lev(ilev,n11)*zfrac
        soil_c(ip,jp,ilev) = soil_c(ip,jp,ilev) + cgdry_lev(ilev,n11)*zfrac
        soil_rho(ip,jp,ilev) = soil_rho(ip,jp,ilev) + rogdr_lev(ilev,n11)*zfrac
        soil_psi(ip,jp,ilev) = soil_psi(ip,jp,ilev) + psig0_lev(ilev,n11)*zfrac
        soil_k(ip,jp,ilev) = soil_k(ip,jp,ilev) + ykghy_lev(ilev,n11)*zfrac
        soil_par_b(ip,jp,ilev) = soil_par_b(ip,jp,ilev) + bghyd_lev(ilev,n11)*zfrac
        soil_par_c(ip,jp,ilev) = soil_par_c(ip,jp,ilev) + cghyd_lev(ilev,n11)*zfrac
        soil_qrel_wilt(ip,jp,ilev) = soil_qrel_wilt(ip,jp,ilev) + qg_rel_wilt_lev(ilev,n11)*zfrac
        soil_qrel_ref(ip,jp,ilev) = soil_qrel_ref(ip,jp,ilev) + qg_rel_refw_lev(ilev,n11)*zfrac
      enddo

      soil_water_table(ip,jp) = soil_water_table(ip,jp) + water_table_depth(n11)*zfrac
      soil_albedo_dry(ip,jp) = soil_albedo_dry(ip,jp) + soildalbed(n11)*zfrac
      soil_albedo_wet(ip,jp) = soil_albedo_wet(ip,jp) + soilwalbed(n11)*zfrac
      soil_emiss_1_dry(ip,jp) = soil_emiss_1_dry(ip,jp) + soildemis1(n11)*zfrac
      soil_emiss_1_wet(ip,jp) = soil_emiss_1_wet(ip,jp) + soilwemis1(n11)*zfrac
      soil_emiss_2_dry(ip,jp) = soil_emiss_2_dry(ip,jp) + soildemis2(n11)*zfrac
      soil_emiss_2_wet(ip,jp) = soil_emiss_2_wet(ip,jp) + soilwemis2(n11)*zfrac

   enddo

 enddo
 enddo

! Horizontal smoothing of soil parameters at all soil levels

 if (nsmooth.ne.0) then

   do ilev=1,nlevel_soil
     zzz1(:,:) = soil_qmax(:,:,ilev)
     call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
     soil_qmax(:,:,ilev) = zzz2(:,:)
   enddo

   do ilev=1,nlevel_soil
     zzz1(:,:) = soil_qmin(:,:,ilev)
     call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
     soil_qmin(:,:,ilev) = zzz2(:,:)
   enddo

   do ilev=1,nlevel_soil
     zzz1(:,:) = soil_c(:,:,ilev)
     call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
     soil_c(:,:,ilev) = zzz2(:,:)
   enddo

   do ilev=1,nlevel_soil
     zzz1(:,:) = soil_rho(:,:,ilev)
     call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
     soil_rho(:,:,ilev) = zzz2(:,:)
   enddo

   do ilev=1,nlevel_soil
     zzz1(:,:) = soil_psi(:,:,ilev)
     call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
     soil_psi(:,:,ilev) = zzz2(:,:)
   enddo

   do ilev=1,nlevel_soil
     zzz1(:,:) = soil_k(:,:,ilev)
     call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
     soil_k(:,:,ilev) = zzz2(:,:)
   enddo

   do ilev=1,nlevel_soil
     zzz1(:,:) = soil_par_b(:,:,ilev)
     call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
     soil_par_b(:,:,ilev) = zzz2(:,:)
   enddo

   do ilev=1,nlevel_soil
     zzz1(:,:) = soil_par_c(:,:,ilev)
     call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
     soil_par_c(:,:,ilev) = zzz2(:,:)
   enddo

   do ilev=1,nlevel_soil
     zzz1(:,:) = soil_qrel_wilt(:,:,ilev)
     call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
     soil_qrel_wilt(:,:,ilev) = zzz2(:,:)
   enddo

   do ilev=1,nlevel_soil
     zzz1(:,:) = soil_qrel_ref(:,:,ilev)
     call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
     soil_qrel_ref(:,:,ilev) = zzz2(:,:)
   enddo

   zzz1(:,:) = soil_water_table(:,:)
   call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
   soil_water_table(:,:) = zzz2(:,:)

   zzz1(:,:) = soil_albedo_dry(:,:)
   call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
   soil_albedo_dry(:,:) = zzz2(:,:)

   zzz1(:,:) = soil_albedo_wet(:,:)
   call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
   soil_albedo_wet(:,:) = zzz2(:,:)

   zzz1(:,:) = soil_emiss_1_dry(:,:)
   call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
   soil_emiss_1_dry(:,:) = zzz2(:,:)

   zzz1(:,:) = soil_emiss_1_wet(:,:)
   call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
   soil_emiss_1_wet(:,:) = zzz2(:,:)

   zzz1(:,:) = soil_emiss_2_dry(:,:)
   call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
   soil_emiss_2_dry(:,:) = zzz2(:,:)

   zzz1(:,:) = soil_emiss_2_wet(:,:)
   call smooth_soil(zzz1,zzz2,nlon,nlat,0.2,nsmooth)
   soil_emiss_2_wet(:,:) = zzz2(:,:)

 endif

! Water body and glacier parameters

 do jp=1,nlat
 do ip=1,nlon

   n11 = soil_fao(ip,jp,1)
   if (n11 >= nst-1) then
     soil_qmax(ip,jp,:) = qgmax(1,n11)
     soil_qmin(ip,jp,:) = qgmin(1,n11)
     soil_c(ip,jp,:) = cgdry(1,n11)
     soil_rho(ip,jp,:) = rogdr(1,n11)
     soil_psi(ip,jp,:) = psig0(1,n11)
     soil_k(ip,jp,:) = ykghy(1,n11)
     soil_par_b(ip,jp,:) = bghyd(1,n11)
     soil_par_c(ip,jp,:) = cghyd(1,n11)
     soil_qrel_wilt(ip,jp,:) = qg_rel_wilt(1,n11)
     soil_qrel_ref(ip,jp,:) = qg_rel_refw(1,n11)
     soil_water_table(ip,jp) = water_table_depth(n11)
     soil_albedo_dry(ip,jp) = soildalbed(n11)
     soil_albedo_wet(ip,jp) = soilwalbed(n11)
     soil_emiss_1_dry(ip,jp) = soildemis1(n11)
     soil_emiss_1_wet(ip,jp) = soilwemis1(n11)
     soil_emiss_2_dry(ip,jp) = soildemis2(n11)
     soil_emiss_2_wet(ip,jp) = soilwemis2(n11)
   else
     soil_qmin(ip,jp,:) = min(soil_qmin(ip,jp,:), soil_qmax(ip,jp,:)-1.e-8)
     soil_qrel_wilt(ip,jp,:) = min(soil_qrel_wilt(ip,jp,:), soil_qrel_ref(ip,jp,:)-1.e-8)
   endif
 enddo
 enddo

return
end subroutine definition_soil_phys_param
!======================================================================
subroutine definition_soil_bottom(soil_type, veg_type, t_amplitude, t_climate, horog, &
 veg_water_table, soil_water_table, water_table, &
 qrel_surf, soil_qmax, soil_qmin, soil_c, soil_rho, soil_psi, soil_par_b, &
 lat_geo, level_soil, nlon, nlat, nst, nvt, nlevel_soil_tot, &
 nlevel_soil, nlevel_soil_water, depth, t_bottom, q_bottom, qrel_bottom)

! Procedure defines bottom condition of soil:
! index of bottom soil level at each grid point and bottom relative soil water content;
! at the bottom soil level climate temperature and bottom relative soil water content
! will be used as the bottom condition.

! Input variables:

!      soil_type:  soil types (classification FAO, 31 types) defined on the grid points, 3D array,
!                1-st and 2-nd array index are grid point index, 3-rd index is soil type number plus 1,
!                if 3-rd index has value 1, then the variable means dominate soil type,
!                values 2-32 of 3-rd index means fraction (of area, proportion) of the each of 31 soil types;

!      veg_type:   vegetation (landuse) types (classification GLC2000, 22 types) defined on the grid points, 3D array,
!               1-st and 2-nd array index are grid point index, 3-rd index is vegetation types number plus 1,
!               if 3-rd index has value 1, then the variable means dominate vegetation type,
!               values 2-23 of 3-rd index means fraction (of area, proportion) of the each of 14 vegetation types

!      t_amplitude: climate amplitude of temperature (K);
!      t_climate: climate temperature (K);
!      horog: orography height (m);
!      veg_water_table: depth of ground water table (m) following to vegetation type;
!      soil_water_table: depth of ground water table (m) following to soil type;
!      water_table: depth of ground water table (m);
!      qrel_surf: relative soil water content (proportion) in top soil layer;
!      soil_qmax: maximum specific volumetric soil water content (m^3/m^3) at all soil levels;
!      soil_qmin: minimum specific volumetric soil water content (m^3/m^3) at all soil levels;
!      soil_c: dry soil thermal capacity (J/K/m^3) at all soil levels;
!      soil_rho: dry soil density (kg/m^3) at all soil levels;
!      soil_psi: idraulic potential of saturated soil (m) at all soil levels;
!      soil_par_b: parameter in soil water transport equation (dimensionless) at all soil levels;

!      lat_geo:   geographical latitude of the grid points (degree, -90...90)

!      level_soil:   soil levels values (m) at total (maximum possible) soil levels;

!      nlon:  grid point number along axis of rotated longitude;
!      nlat:  grid point number along axis of rotated latidude;
!      nst:   total number of soil types;
!      nvt:   total number of vegetation types;
!      nlevel_soil_tot: number (maximum possible) of soil levels.

! Output variables:

!      nlevel_soil: index of bottom (climate) soil level for heat transport at each grid point;
!      nlevel_soil_water: index of bottom soil level for water transport at each grid point;
!      depth: depth (m) of bottom (climate) soil level at each grid point;
!      t_bottom: soil temperature (K) at bottom soil level;
!      q_bottom: soil water content (m^3/m^3) at bottom soil level;
!      qrel_bottom: relative soil water content (proportion) at bottom soil level..

implicit none

! Input:

integer :: nlon, nlat, nst, nvt, nlevel_soil_tot
real, dimension(nlon,nlat,nst+1) :: soil_type
real, dimension(nlon,nlat,nvt+1) :: veg_type
real, dimension(nlon,nlat) :: t_amplitude, t_climate, horog, veg_water_table, soil_water_table, qrel_surf, lat_geo
real, dimension(nlon,nlat,nlevel_soil_tot) :: soil_qmax, soil_qmin, soil_c, soil_rho, soil_psi, soil_par_b
real, dimension(nlevel_soil_tot) :: level_soil

! Output:

integer, dimension(nlon,nlat) :: nlevel_soil, nlevel_soil_water
real, dimension(nlon,nlat) :: depth, water_table, t_bottom, q_bottom, qrel_bottom

! Internal:

real, parameter :: water_table_veg_weight=0.7
integer, parameter :: nveg_water_table=3
integer, dimension(nveg_water_table) :: vegtype_water_table=(/16, 17, 18/)
real, dimension(nveg_water_table) :: weight_work
integer :: i, ii, j, lev, lev_ref, iflag, itype
real :: depth_ref=1.5, period=31536000., pi, amplitude_ref=1., depth_prec_season=4., depth_min=2., &
 qrel_ref, q, psi, lambda, kappa, cw=4186.8, rhow=1000., frac_bare, zdepth, z1, weight, weight_orog

 print *
 print *,"Definition of bottom condition for soil variables"

 pi=abs(acos(-1.))

! Reference soil level

 do lev=1,nlevel_soil_tot
   if (level_soil(lev) >= depth_ref) exit
 enddo
 lev_ref=max( min(lev, nlevel_soil_tot), 1)

 do j=1,nlat
 do i=1,nlon

   if (soil_psi(i,j,1) < 0.) then ! Not water body or glacier

! Reference relative soil water content

     qrel_ref=max(qrel_surf(i,j), 0.5)

!  Reference soil water content

     q=qrel_ref*(soil_qmax(i,j,lev_ref)-soil_qmin(i,j,lev_ref))+soil_qmin(i,j,lev_ref)

! Reference soil hydraulic potential

     psi=soil_psi(i,j,lev_ref)*((soil_qmax(i,j,lev_ref)/q)**soil_par_b(i,j,lev_ref))

! Thermal conductivity (dynamic) is the function of the soil hydraulic potential
! Soil thermal conductivity at reference level is used for the whole soil column

     lambda=min( 0.5+1400.*exp(-6.-alog10(abs(psi))), 2.5)

! Kinematic thermal conductivity of soil

     kappa=lambda/(soil_c(i,j,lev_ref)*soil_rho(i,j,lev_ref)+q*cw*rhow)

! Depth of attenuation of season temperature oscilation to refenence value

     depth(i,j)=(-alog((amplitude_ref/t_amplitude(i,j)**2)))/(2.*sqrt(pi/(kappa*period)))

     depth(i,j)=max(depth(i,j), depth_min)

! Control of presence of vegeation with is typical for expressed season precipitation regime

     iflag=0

     itype=2 ! Tree cover, broadleaved deciduous, closed (>40%)
     if (veg_type(i,j,itype+1) > 0.3) iflag=1

     itype=3 ! Tree cover, broadleaved deciduous, open (15-40%)
     if (veg_type(i,j,itype+1) > 0.3) iflag=1

     itype=11 ! Shrubcover, closed to open (>15%), evergreen(broadleaved or needleleaved)
     if (veg_type(i,j,itype+1) > 0.3) iflag=1

     if (iflag == 1) depth(i,j)=max(depth(i,j), depth_prec_season)

   else ! Water body or glacier

     depth(i,j)=100.

     iflag=0
     itype=nst-1 ! glacier
     if (soil_type(i,j,itype+1) > 0.5) iflag=1
     if (iflag == 1) depth(i,j)=10.

   endif

! Index of soil level with be used as bottom condition level for heat transport simulation

   do lev=1,nlevel_soil_tot
     if (level_soil(lev) >= depth(i,j)) exit
   enddo
   nlevel_soil(i,j)=min(lev, nlevel_soil_tot) 

! Index of soil level with be used as bottom condition level for water transport simulation

   if (soil_psi(i,j,1) < 0.) then ! Not water body or glacier
!     itype=19 ! Bare Areas 
!     frac_bare=veg_type(i,j,itype+1)
!     water_table(i,j)=veg_water_table(i,j)*(1.-frac_bare)+soil_water_table(i,j)*frac_bare 
!     water_table(i,j)=min(veg_water_table(i,j), soil_water_table(i,j)) 
     weight=water_table_veg_weight
     weight_work(1:nveg_water_table)=water_table_veg_weight
     do ii=1,nveg_water_table
       itype=vegtype_water_table(ii)+1
       weight_work(ii)=veg_type(i,j,itype)*(0.2-water_table_veg_weight)+water_table_veg_weight
     enddo
     weight=minval(weight_work(1:nveg_water_table))
     water_table(i,j)=veg_water_table(i,j)*weight+soil_water_table(i,j)*(1.-weight) 
     if (horog(i,j) <= 500.) then
       weight_orog=0.
     elseif (horog(i,j) >=1500.) then
       weight_orog=1.
     else
       weight_orog=1./1000.*horog(i,j)-500./1000.
     endif
     water_table(i,j)=water_table(i,j)*(1.-weight_orog)+50.*weight_orog
   else ! Water body or glacier
     water_table(i,j)=0.
   endif

   do lev=1,nlevel_soil_tot
     if (level_soil(lev) >= water_table(i,j)) exit
   enddo
   nlevel_soil_water(i,j)=min(lev, nlevel_soil_tot) 
   nlevel_soil_water(i,j)=min(nlevel_soil_water(i,j), nlevel_soil(i,j)) 

! Bottom condition of temperature and water content

   if (soil_psi(i,j,1) < 0.) then ! Not water body or glacier

     t_bottom(i,j)=t_climate(i,j)

     zdepth=level_soil(nlevel_soil(i,j))

     if (water_table(i,j) > zdepth) then
       z1=zdepth/water_table(i,j)
       qrel_bottom(i,j)=min(1.*z1+qrel_surf(i,j)*(1.-z1), 1.)
     else
       qrel_bottom(i,j)=1.
     endif

     lev=nlevel_soil_water(i,j)
     q_bottom(i,j)=qrel_bottom(i,j)*(soil_qmax(i,j,lev)-soil_qmin(i,j,lev))+soil_qmin(i,j,lev)

   else ! Water body or glacier

     t_bottom(i,j)=t_climate(i,j)

     iflag=0
     itype=nst-1 ! glacier
     if (soil_type(i,j,itype+1) > 0.5) iflag=1
     if (iflag == 1) t_bottom(i,j)=min(t_climate(i,j), 263.15)

     qrel_bottom(i,j)=1.
     nlevel_soil_water(i,j)=0

   endif

 enddo
 enddo

return
end subroutine definition_soil_bottom
!======================================================================
      subroutine smooth_soil(vt,vout,ni,nj,w,nsmooth)

!  Performs smoothing of matrix VT(NI,NJ) using a 5-point filter (laplacian).
!  It is assumed that the grid distances in x and y are similar
!  The filtering degree is prop. to param. W (between 0 and 1)
!  (for stronger filtering it is advised to increase the number of filtering iterations
!  NSMOOTH, keeping W not larger than about 0.5).
!  Output in VOUT(NI,NJ) and also in VT (caution: VT is redefined!)

      dimension vt(ni,nj), vout(ni,nj)

      do ns=1,nsmooth

      do j=2,nj-1
      do i=2,ni-1
      vout(i,j) = (1.-w)*vt(i,j)+.25*w*(vt(i+1,j)+vt(i-1,j)+vt(i,j+1)+vt(i,j-1))
      enddo
      enddo

!  Border points (1-D filter)

      do i=2, ni-1
      vout(i,1) = (1.-.7*w)*vt(i,1)+.35*w*(vt(i+1,1)+vt(i-1,1))
      vout(i,nj) = (1.-.7*w)*vt(i,nj)+.35*w*(vt(i+1,nj)+vt(i-1,nj))
      enddo
      do j=2, nj-1
      vout(1,j) = (1.-.7*w)*vt(1,j)+.35*w*(vt(1,j+1)+vt(1,j-1))
      vout(ni,j) = (1.-.7*w)*vt(ni,j)+.35*w*(vt(ni,j+1)+vt(ni,j-1))
      enddo

!  Corner points (unchanged)

      vout(1,1) = vt(1,1)
      vout(1,nj) = vt(1,nj)
      vout(ni,nj) = vt(ni,nj)
      vout(ni,1) = vt(ni,1)

      do j=1,nj
      do i=1,ni
      vt(i,j) = vout(i,j)
      enddo
      enddo

      enddo
      return
      end
!=======================================================================
subroutine interp(alfa, ex1, ex2, npi, xi, g, x, f, nval)

! (Similar to INTERP_SPLINE_1D, but with different input, and does not require
!  call of SUBR. NEAR - moreover, it checks that input coordinates are monotonic).

!  Interpolates with splines with tension in one dimension.
!  The spline is defined imposing that the second derivative is the average
!  of second derivatives computed at the two adjacent points.
!  At interval extremes the second derivative is assumed null.
!  This subroutine also extrapolates out of the interval where the input funtion G is defined

!  INPUT:  function G defined at coordinates XI (CAUTION: can be changed by this subroutine)
!          G(1:NPI) values at irregular but strictly growing coordinates XI(1:NPI)
!  OUTPUT: F(1:NVAL) interpolated values at arbitrary coordinates X(1:NVAL)

!  ALFA: spline tension parameter, comprised between 0 and 1
!  If ALFA=1, pure linear interpolation; if ALFA=0, pure spline

!  EX1: param. determining extrapolation for X < XI(1)
!  EX2: param. determining extrapolation for X > XI(NPI)
!  If EX1=0 OR EX2=0, constant value extrapolation is used at corresponding extreme
!  If EX1=1 OR EX2=1, linear extrapolation is used at corresponding extreme
!  Intermediate values of EX1 AND EX2 give intermediate extrapolation values

  real, dimension(npi )  :: xi, g
  real, dimension(nval)  :: x,  f

  if(alfa.lt..0.or.alfa.gt.1) then
  print*, 'CAUTION: in INTERP, ALFA out of interval 0-1'
  endif
  if(ex1.lt..0.or.ex1.gt.1) then
  print*, 'CAUTION: in INTERP, EX1 out of interval 0-1'
  endif
  if(ex2.lt..0.or.ex2.gt.1) then
  print*, 'CAUTION: in INTERP, EX2 out of interval 0-1'
  endif

! Fix for the case in which coordinates of the input function are not strictly increasing
! Note that this changes the input coordinates in the calling programme

  do k=2,npi
   if(xi(k).le.xi(k-1)) then
   print*, "CAUTION: in INTERP, coordinates of input function changed because not monotonic!"
   exit
   endif
  enddo

  zeps=(xi(npi)-xi(1))*1.e-6   ! Small deviation used to set apart interlaced coordinates
200 do k=2,npi
     if(xi(k).le.xi(k-1)) then
     ximed=0.5*(xi(k)+xi(k-1))
     xi(k-1)=ximed-zeps
     xi(k)=ximed+zeps
     gmed=0.5*(g(k)+g(k-1))
     g(k-1)=gmed
     g(k)=gmed
     endif
    enddo

 do k=2,npi
  if(xi(k).le.xi(k-1)) then
  goto 200
  endif
 enddo

 do 100 jval =1, nval

!  2 cases of extrapolation

 if(x(jval).lt.xi(1)) then
 f(jval) = g(1) + ex1*(g(1)-g(2))/(xi(1)-xi(2)) * (x(jval)-xi(1))
 go to 100
 elseif (x(jval).ge.xi(npi)) then
 f(jval) = g(npi) + ex2*(g(npi)-g(npi-1))/(xi(npi)-xi(npi-1)) * (x(jval)-xi(npi))
 go to 100
 endif

 ir = 0

!  IR is a reference index determining the interpolation interval
!  The interpolation expression is applied also if X=XI(J)

 do j = 1, npi
 if (x(jval).ge.xi(j)) ir = ir + 1
 enddo

 if (ir.eq.1) then
 fmm = 2*g(1) - g(2)
 xmm = 2*xi(1) - xi(2)
 fpp = g(ir+2)
 xpp = xi(ir+2)
 elseif (ir.eq.(npi-1)) then
 fpp = 2*g(npi) - g(npi-1)
 xpp = 2*xi(npi) - xi(npi-1)
 fmm = g(ir-1)
 xmm = xi(ir-1)
 else
 fmm = g(ir-1)
 xmm = xi(ir-1)
 fpp = g(ir+2)
 xpp = xi(ir+2)
 endif

 fm     = g(ir)
 xm     = xi(ir)
 fp     = g(ir+1)
 xp     = xi(ir+1)
 delx   = xp - xm
 delxp  = xpp - xp
 delxm  = xm - xmm
 delx1  = x(jval) - xm
 delx2  = xp - x(jval)
 delxs  = delx**2
 delx1s = delx1**2
 delx2s = delx2**2

!  Spline contribution to interpolation

 spl = fm*(delx2/delx + delx1*delx2s/(delxs*delxm) - delx1s*     &
       delx2/((delx+delxp)*delxs)) + fp*(delx1/delx +            &
       delx1s*delx2/(delxs*delxp) - delx1*delx2s/((delx+delxm)*  &
       delxs)) - fmm * delx1*delx2s/((delx+delxm)*delx*delxm) -  &
       fpp * delx1s*delx2/((delx+delxp)*delx*delxp)

!  Linear interpolation contribution

 clin = (fm*delx2 + fp*delx1)/delx

!  Final interpolation combined using ALFA

 f(jval) = alfa*clin + (1.-alfa)*spl

 100  continue

 return
end subroutine interp
!=======================================================================
