subroutine soil_read_input(mb, mstruc, neap, cn2a, cn2b, cn2c, cn2d)
!**** PURPOSE: THIS SUBROUTINE READS SOIL DATA
!**** CALLED IN: MAIN
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! PARAMETERS & VARIABLES
!
! >>>>> COMMON PARAMETERS & VARIABLES
! abd(k) = sum of porosity for soil type, mm, used in solt
! alai0(n) = initial Leaf Area Index for land use type
! ap(l, k) = labile (soluble) phosphorus, g/t, then trans. to kg/ha
! asc(1:17) = empirical coef. (estimated here) to calc sc(l, k)
! awc(l, k) = available water capacity, % , /100. to transfer to m/m
! bd(k) = bulk density of the upper (!) soil layer, g/cm3,
! calc from porosity
! bden(l, k) = bulk density, input, g/cm3
! bsc(1:17) = empirical coef to calc sc(l, k)
! cbn(l, k) = organic carbon content, %
! cla(l, k) = clay content, %
! clay(1, k) = clay content, %
! cn2(k, n, j) = Curve Numbers for soil k, land use n, and subbasin
! cn2a(n) = Curve Numbers for soil group A
! cn2b(n) = Curve Numbers for soil group B
! cn2c(n) = Curve Numbers for soil group C
! cn2d(n) = Curve Numbers for soil group D
! cv(j) = initial land cover, kg/ha
! ek(k) = soil erodibility factor
! fc(l, k) = field capacity, %, then trans. to mm
! hk(l, k) = beta coef. to estimate hydr. cond., used in perc
! hum(l, k) = humus content, kg/ha
! isc = code for saturated cond.: 0 - read, 1 - calc
! nn = number of soil layers
! ns(k) = number of soil layers
! nsa(k) = number of layers for arable soils
!block nsolgr(k) = number of soils group for soil type, 1=A, 2=B, 3=C, 4=D
! op(l, k) = stable min. P, estimated from pmn(l, k), coef = 4., kg/ha
! pmn(l, k) = act. min. P, estimated from ap(l, k), psp=0.5, kg/ha
! por(l, k) = porosity, m/m (calc)
! poros(l, k) = porosity, % (input)
!block psp = parameter = 0.5
! psz(5, k) = particle size distribution
!block rtn = parameter = .15
! san(l, k) = sand content, %
! sand(1, k) = sand content, %
! sc(l, k) = saturated conductivity, mm/h, calc, if isc = 1
! sccor = correction factor for saturated conductivity (all soils)
! sil(l, k) = silt content, %
! silt(1, k) = silt content, %
! snam(k) = soil name
! stin(l, k) = water content in soil layer, mm
! stinco = init. water content coef, from readbas
! sumfc(k) = sum of field capacity in soil, mm
! sumul(k) = sum of upper limit water content in soil, mm
! swin(k) = soil water content in soil profile, mm
! ul(l, k) = upper limit water content in layer, mm
! up(l, k) = upper limit water content, m/m (used only here to calc fc)
! wmn(l, k) = active org. N, estim. from wn(l, k), rtn=.15, kg/ha
! wn(l, k) = stable organic N content, %,
! recalc from % to g/t OR
! estim. from C using C:N=10. and 10000. coef from % to g/t,
! then trans. to kg/ha
! wno3(l, k) = NO3-N content (kg/ha), estimated in g/t, trans to kg/ha
! wp(l, k) = wilting point, in m/m (used here and in perc)
! wpo(l, k) = org. P, estimated from wn(), coef=.125, kg/ha
! z(l, k) = depth to bottom of layer, mm
! >>>>>
! >>>>> STATIC PARAMETERS
! ccf = correction coef for up and por by silt content
! dg = layer thickness
! i = cycle par
! j = par
! jj = par
! k = soil type number
! ka = par
! l = par
! n = par
! numsb = actual number of soils in database
! sfile = soil file name
! tex = sum of clay, silt and sand
! titldum = text
! wt1 = coef. to transform from g/t to kg/ha
! xx = par
! xz = par
! zz = par
! >>>>>
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!**** Include common parameters
use utilities, only : open_file, check_range, to_string
use input, only : read_string_column
integer, intent(in) :: mb
integer, dimension(:, :, :), intent(in) :: mstruc
integer, dimension(:), intent(in) :: neap
real(dp), dimension(:), intent(in) :: cn2a, cn2b, cn2c, cn2d
! soil counter
integer :: k
! layer counter
integer :: l
integer :: jj, n
integer :: j, jea
character(len = 100) :: titldum, snam, sn
character(len=path_max_length) :: soildat
real(dp) :: dg, wt1, xx, xz, zz
real(dp), dimension(17):: asc
integer :: fid
call log_debug('soil_read_input', 'SOIL PARAMETERS:')
call log_debug('soil_read_input', 'Number of soils in soil.csv = ', int=msdb)
!**** Reading Soil parameters from soilXX.dat files
call read_string_column(soil_input_file_id, "soil_file", soil_filenames, "")
do k = 1, msdb
if (k .gt. msdb .or. k .eq. 0) then
call log_error("soil_read_input", 'Soil type (k > msdb) = max soil number, or 0', &
i1=k, i2=msdb)
end if
soildat = soil_filenames(k)
! open soil parameter files
fid = open_file(trim(soil_input_dir)//'/'//soildat)
!**** READ SOIL HYDRAULIC PROPERTIES
! THE SOIL IS DIVIDED VERTICALLY INTO UP TO 10 LAYERS
! LAYER THICKNESS CAN BE VARIABLE
! ns - number of layers, nsa - number of layers for arable soils
read(fid, *, ERR=55) snum(k), ns(k), nsa(k), nsolgr(k)
nn = ns(k)
if (ns(k) .gt. ml .or. ns(k) .eq. 0 ) then
call log_error("soil_read_input", 'Number of layers for (soil > ml) = max number, or 0', i1=k, i2=ml)
end if
if (nsa(k) .gt. ns(k) .or. nsa(k) .eq. 0 ) then
call log_error("soil_read_input", 'Number of arrable layers grater than number of layers (soil) or 0', i1=k, ints=(/nsa(k), ns(k)/))
end if
read(fid, "(a)") snam
read(fid, "(a)") titldum
read(fid, "(a)") titldum
read(fid, *, ERR=55) (z(l, k), l = 1, nn)
read(fid, *, ERR=55) (cla(l, k), l = 1, nn)
read(fid, *, ERR=55) (sil(l, k), l = 1, nn)
read(fid, *, ERR=55) (san(l, k), l = 1, nn)
read(fid, *, ERR=55) (bden(l, k), l = 1, nn)
read(fid, *, ERR=55) (poros(l, k), l = 1, nn)
read(fid, *, ERR=55) (awc(l, k), l = 1, nn)
read(fid, *, ERR=55) (fc(l, k), l = 1, nn)
read(fid, *, ERR=55) (cbn(l, k), l = 1, nn)
read(fid, *, ERR=55) (wn(l, k), l = 1, nn)
read(fid, *, ERR=55) (sc(l, k), l = 1, nn)
read(fid, *, ERR=55) ek(k)
close(fid)
if (k .gt. ms .or. k .eq. 0) call log_error("soil_read_input", &
'Soil type number > ms = max soil number, or 0', int=k)
sn = to_string(int=k)
call check_range(z(:, k), 'depth of soil '//trim(sn), (/0.001, 1e5/))
call check_range(cla(:, k) + sil(:, k) + san(:, k), 'texture of soil '//trim(sn), &
(/99., 101./))
call check_range(bden(:, k), 'bulk density of soil '//trim(sn), (/0.001, 100./))
call check_range(poros(:, k), 'porosity of soil '//trim(sn), (/0.01, 100./))
call check_range(awc(:, k), 'water capacity of soil '//trim(sn), (/0.01, 100./))
call check_range(fc(:, k), 'field capacity of soil '//trim(sn), (/0.01, 100./))
call check_range(sc(:, k), 'saturated conductivity of soil '//trim(sn), (/0.01, 1e6/))
call check_range((/ek(k)/), "erodibility of soil "//trim(sn), (/0.01, 1./), index=k)
!**** ESTIMATE CARBON FOR LOWER LAYERS IF AVAILABLE ONLY FOR UPPER LAYER
if (cbn(3, k) .le. 0) then
xx = z(2, k)
do l = 3, nn
dg = (z(l, k) - xx)
if (cbn(l, k) .eq. 0.) cbn(l, k) = cbn(l - 1, k) * exp(- .001 * dg)
xx = z(l, k)
end do
end if
!**** CALCULATE INITIAL NUTRIENT CONTENT: wno3, wn, wmn, wpo, ap, pmn, op, hum
! INITIALIZE sc() if isc = 0
! Units transformation:
! g/t --> kg/ha *wt1=bden()*dg/100., bden() - bulk density
! % --> g/t *10000.
! % --> kg/ha *wt1*10000.
! kg/ha --> g/t 1/wt1=100./dg/bden()
! mg/kg = g/t
xx = 0.
do l = 1, nn
if (isc .eq. 1) sc(l, k) = 0.
dg = (z(l, k) - xx)
wt1 = bden(l, k) * dg / 100.
if (wno3(l, k) .le. 0.) then
wno3(l, k) = 10. * exp(- z(l, k) / 1000.)
wno3(l, k) = wno3(l, k) * wt1
endif
wn(l, k) = 10000. * wn(l, k)
if (wn(l, k) .le. 0.) wn(l, k) = 10000. * cbn(l, k) / 10.
wn(l, k) = wn(l, k) * wt1
wmn(l, k) = wn(l, k) * rtn
wn(l, k) = wn(l, k) * (1. - rtn)
wpo(l, k) = .125 * wn(l, k)
if (ap(1, k) .le. 0.) ap(1, k) = 5.
ap(l, k) = ap(1, k) * wt1
pmn(l, k) = ap(l, k) * (1. - psp) / psp
op(l, k) = 4. * pmn(l, k)
hum(l, k) = cbn(l, k) * wt1 * 172.
xx = z(l, k)
end do
!**** CALCULATE PARTICLE SIZE DISTRIBUTION psz
silt(1, k) = sil(1, k) / 100.
clay(1, k) = cla(1, k) / 100.
sand(1, k) = san(1, k) / 100.
psz(1, k) = (1. - clay(1, k)) ** 2.49 * sand(1, k)
psz(2, k) = .13 * silt(1, k)
psz(3, k) = .20 * clay(1, k)
if (clay(1, k) .le. .5) then
if (clay(1, k) .lt. .25) then
psz(4, k) = 2. * clay(1, k)
else
psz(4, k) = .28 * (clay(1, k) - .25) + .5
end if
else
psz(4, k) = .57
end if
psz(5, k) = 1. - psz(1, k) - psz(2, k) - psz(3, k) - psz(4, k)
!**** CALCULATE/RECALCULATE HYDRAULIC PARAMETERS wp, up, por, awc &
! CALCULATE SATURATED CONDUCTIVITY sc IN CASE isc=1
! CORRECTION of up & por by silt content: Fred Hattermann 2002, ccf=.05
! Removed 2011 from Tobias Vetter
do l = 1, nn
! wp(l, k) = 0.4 * cla(l, k) * bden(l, k) / 100.
wp(l, k) = fc(l, k) / 100 - awc(l, k) / 100
! ccf = 0.05
! up(l, k) = wp(l, k) + awc(l, k)/100. + ccf * sil(1, k) / 100.
up(l, k) = fc(l, k) / 100
! por(l, k) = 1. - bden(l, k) / 2.65 + ccf * sil(1, k) / 100.
por(l, k) = poros(l, k) / 100
if (up(l, k) .ge. por(l, k)) then
up(l, k) = por(l, k) - .05
wp(l, k) = up(l, k) - awc(l, k) / 100.
if (wp(l, k) .le. 0.) then
up(l, k) = por(l, k) * .75
wp(l, k) = por(l, k) * .25
end if
end if
awc(l, k) = up(l, k) - wp(l, k)
if (isc .eq. 1) then
asc(1) = 1.
asc(2) = san(l, k)
asc(3) = cla(l, k)
asc(4) = poros(l, k) / 100.
asc(5) = san(l, k) * san(l, k)
asc(6) = cla(l, k) * cla(l, k)
asc(7) = poros(l, k) * poros(l, k) / 10000.
asc(8) = san(l, k) * cla(l, k)
asc(9) = san(l, k) * poros(l, k) / 100.
asc(10) = cla(l, k) * poros(l, k) / 100.
asc(11) = asc(5) * cla(l, k)
asc(12) = asc(6) * poros(l, k) / 100.
asc(13) = asc(5) * poros(l, k) / 100.
asc(14) = asc(6) * san(l, k)
asc(15) = asc(7) * cla(l, k)
asc(16) = asc(5) * asc(7)
asc(17) = asc(6) * asc(7)
do jj = 1, 17
sc(l, k) = sc(l, k) + asc(jj) * bsc(jj)
end do
sc(l, k) = exp(sc(l, k)) * 10.
end if
end do
ek(k) = ek(k) * ekc0
end do ! loop over number of soil files
!*********************************************************** END OF SOIL CYCLE
!**** DETERMINE DRAINAGE COEFS ul, sumul, fc, sumfc, hk, wss &
! WRITE SELECTED SOIL PARAMETERS
do k = 1, ms
if (z(1, k) .gt. 0.) then
nn = ns(k)
bd(k) = 2.65 * (1. - por(1, k))
abd(k) = 0.
zz = 0.
xx = 0.
swin(k) = 0.
do l = 1, nn
dg = z(l, k) - xx
xz = por(l, k) * dg * 10.
abd(k) = abd(k) + xz
ul(l, k) = (por(l, k) - wp(l, k)) * dg
sumul(k) = sumul(k) + ul(l, k)
fc(l, k) = dg * (up(l, k) - wp(l, k))
sumfc(k) = sumfc(k) + fc(l, k)
stin(l, k) = fc(l, k) * stinco
hk(l, k) = - 2.655 / log10(fc(l, k) / ul(l, k))
xx = z(l, k)
swin(k) = swin(k) + stin(l, k)
end do
abd(k) = abd(k) / (10. * z(nn, k))
end if
end do
!**** DEFINE CURVE NUMBERS
! cn2(k, n)- CN, conditions 2
do j = 1, mb
do jea = 1, neap(j)
n = mstruc(j, jea, 1)
k = mstruc(j, jea, 2)
select case(nsolgr(k))
case(1)
cn2(j, jea) = cn2a(n)
case(2)
cn2(j, jea) = cn2b(n)
case(3)
cn2(j, jea) = cn2c(n)
case(4)
cn2(j, jea) = cn2d(n)
end select
end do
end do
cv = 0.01
call log_info('soil_read_input', 'Soil parameters read for N soils: ', int=k)
return
55 continue
call log_error('soil_read_input', "Failed to read soil file for soil", int=k)
end subroutine soil_read_input