!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutines for: *!
!* Soil and Water - Programs *!
!* *!
!* contains: *!
!* SOIL *!
!* SOIL_INI *!
!* SOIL_WAT *!
!* UPT_WAT *!
!* FRED1 - ...11 *!
!* TAKE_WAT *!
!* BUCKET *!
!* SNOWPACK *!
!* HUM_ADD *!
!* BC_APPL: application of biochar *!
!* *!
!* Copyright (C) 1996-2018 *!
!* Potsdam Institute for Climate Impact Reserach (PIK) *!
!* Authors and contributors see AUTHOR file *!
!* This file is part of 4C and is licensed under BSD-2-Clause *!
!* See LICENSE file or under: *!
!* http://www.https://opensource.org/licenses/BSD-2-Clause *!
!* Contact: *!
!* https://gitlab.pik-potsdam.de/foresee/4C *!
!* *!
!*****************************************************************!
SUBROUTINE soil
! Soil processes (frame)
use data_climate
use data_depo
use data_out
use data_simul
use data_soil
use data_soil_cn
implicit none
call evapo
call intercep
call soil_wat
call soil_temp
if (flag_wurz .eq. 4 .or. flag_wurz .eq. 6) then
call soil_stress !calculate ground stress factors
call cr_depth !define root depth
endif
call soil_cn
call root_eff
END subroutine soil
!**************************************************************
SUBROUTINE soil_ini
! Initialisation of soil data and parameters
use data_inter
use data_evapo
use data_out
use data_par
use data_simul
use data_soil
use data_soil_cn
use data_species
use data_stand
implicit none
integer i,j,k
real d_0, t_0
! Table of quarz and clay content (mass%) versus wlam
real, dimension(17) :: xwlam = (/1.5, 1.15, 0.9, 0.67, 0.6, 0.5, 0.38, 0.37, 0.3, 0.29, 0.27, 0.26, 0.25, 0.24, 0.23, 0.22, 0.15/), &
yquarz = (/93.0,85.0,80.0, 82.0,76.0, 64.0, 65.0, 51.0,60.0, 30.0, 14.0, 10.0, 12.0, 20.0, 30.0, 43.0, 23.0/), &
yclay = (/3.0, 3.0, 3.0, 12.0, 6.0, 6.0, 10.0, 4.0,21.0, 12.0, 10.0, 37.0, 15.0, 40.0, 30.0, 35.0, 55.0/)
real value
real, dimension(nlay):: humush(nlay)
real, dimension(nlay):: xfcap, xwiltp, xpv ! output of addition for water capacities
! estimation of soil layer values
d_0 = 0.
do j = 1, nlay
t_0 = thick(j)
mid(j) = d_0 + 0.5*t_0
d_0 = d_0 + t_0
depth(j) = d_0
enddo
perc = 0.
wupt_r = 0.
wupt_ev = 0.
thick_1 = thick(1)
select case (flag_soilin)
case (0,2)
do i=1,nlay
if (i .gt. 1) then
call tab_int(xwlam, yquarz, 17, wlam(i), value)
sandv(i) = value / 100. ! Mass% of mineral fraction
call tab_int(xwlam, yclay, 17, wlam(i), value)
clayv(i) = value / 100.
siltv(i) = 1. - clayv(i) - sandv(i)
else
sandv(1) = 0.0
clayv(1) = 0.0
siltv(1) = 0.0
endif
enddo
case (1,3,4)
clayv = clayv / 100.
sandv = sandv / 100.
siltv = 1. - clayv - sandv
if ((sandv(1) .le. zero) .and. (clayv(1) .le. zero)) siltv(1) = 0.
skelv = skelv / 100.
humusv = humusv / 100.
end select
! Settings for subroutine take_wat
skelfact = 1.
pv = skelfact * pv_v * thick * 0.1 ! mm
wilt_p = skelfact * wilt_p_v * thick * 0.1 ! mm
field_cap = skelfact * f_cap_v * thick * 0.1 ! mm
wats = field_cap ! mm
watvol = f_cap_v ! vol%
n_ev_d = nlay
nlgrw = nlay+1
do i=1,nlay
if (w_ev_d .gt. depth(i)) n_ev_d = i
if (grwlev .gt. depth(i)) nlgrw = i+1
vol(i) = thick(i) * 10000.
enddo
! dry mass of first layer
dmass = vol * dens
rmass1 = dmass(1) - (C_hum(1) + C_opm(1)) / cpart ! corection term of first layer
humush = humusv
if (2.*C_hum(1) .lt. humusv(1)*dmass(1)) then
humusv(1) = C_hum(1) / (dmass(1) * cpart)
endif
do i=2, nlay
humusv(i) = C_hum(i) / (dmass(i) * cpart)
enddo
if (flag_bc .gt. 0) y_bc_n = 1 ! actual number of biochar application
! calculation of additions for water capacities
call hum_add(xfcap, xwiltp, xpv)
fcaph = f_cap_v - xfcap
wiltph = wilt_p_v - xwiltp
pvh = pv_v - xpv
! ground water
do i = nlgrw, nlay
wats(i) = pv(i)
enddo
interc_can = 0.
int_st_can = 0.
interc_sveg = 0.
int_st_sveg = 0.
aev_i = 0.
wat_tot = SUM(wats)
END subroutine soil_ini
!**************************************************************
SUBROUTINE soil_wat
use data_out
! soil water balance
use data_climate
use data_evapo
use data_inter
use data_out
use data_par
use data_simul
use data_soil
use data_species
use data_stand
implicit none
real :: eva_dem ! evaporation demand of soil
real :: p_inf = 0. ! infiltrated water
real :: pev ! local: soil evaporation
real :: watot, wetot ! total water content at start and end
real :: wutot ! total water uptake from the soil
real :: wutot_ev ! total water uptake by soil evaporation
real :: wutot_r ! total water uptake by roots
real, allocatable, dimension(:) :: upt ! local array for uptake
real, external :: wat_new
real enr, wa, we, percolnl, snow_sm, buckdepth
integer j
allocate (upt(nlay))
wupt_ev = 0.
aev_s = 0.
prec_stand = MAX(0., prec - interc_can - interc_sveg) ! stand precipitation
if (flag_int .gt. 1000) then
prec_stand = prec_stand * prec_stand_red / 100.
endif
call snowpack(snow_sm, p_inf, pev)
if (anz_coh .le. 0) pev = pet
eva_dem = MAX(0., p_inf) - pev ! evaporation demand of soil
! if all stand precipitation is evaporated and there is still a demand
! there is an uptake from soil layers (up to an certain depth)
if (eva_dem .lt. 0.) then
if (snow .le. 0) call take_wat(eva_dem, cover)
aev_s = aev_s + p_inf + SUM(wupt_ev)
else
aev_s = aev_s + pev
endif ! eva_dem
aet = aev_s + aev_i
p_inf = MAX(eva_dem, 0.)
upt = wupt_ev
watot = SUM(wats) ! total initial water content
do j = 1, nlgrw-1
enr = p_inf - upt(j)
wa = wats(j) - field_cap(j)
we = wat_new(wa, enr, j)
p_inf = enr + wa - we
perc(j) = MAX(p_inf, 0.)
wats(j) = MAX(we+field_cap(j), wilt_p(j))
enddo
do j = nlgrw, nlay
enr = p_inf - upt(j)
wa = wats(j) - field_cap(j)
we = pv(j) - field_cap(j) ! ground water level is constant!
p_inf = enr + wa - we
perc(j) = MAX(p_inf, 0.)
wats(j) = MAX(we+field_cap(j), wilt_p(j))
enddo
if (flag_wred .le. 10) then
call upt_wat
else
call upt_wat1
endif
! root uptake balanced imediate after calculation
upt = upt + wupt_r
wats = wats - wupt_r
watvol = 10.*wats/(thick * skelfact) ! estimation for complete layer in Vol% without skeleton (only soil substrate)
! total water quantities
wetot = SUM(wats) ! total final water content
wutot_ev = SUM(wupt_ev) ! total water uptake by soil evaporation
wutot_r = SUM(wupt_r) ! total water uptake by roots
wutot = wutot_ev + wutot_r ! total water uptake
aet = aet + wutot_r ! daily total aet
percolnl = perc(nlay)
trans_tree = 0.
trans_sveg = 0.
zeig => pt%first
do while (associated(zeig))
if (zeig%coh%species .le. nspec_tree) then
trans_tree = trans_tree + zeig%coh%supply
else
trans_sveg = trans_sveg + zeig%coh%supply
endif
zeig => zeig%next
enddo ! zeig (cohorts)
! cumulative water quantities
perc_cum = perc_cum + perc(nlay)
wupt_r_c = wupt_r_c + wutot_r
wupt_e_c = wupt_e_c + wutot_ev
wupt_cum = wupt_cum + wutot
aet_cum = aet_cum + aet
dew_cum = dew_cum + dew_rime
tra_tr_cum = tra_tr_cum + trans_tree
tra_sv_cum = tra_sv_cum + trans_sveg
call bucket(bucks_100, bucks_root, buckdepth)
! number of drought days per layer
where ((wats-0.2) .le. wilt_p) s_drought = s_drought+1
if (flag_dayout .ge. 2) then
write (unit_wat, '(2I5, 7F7.2, 24F8.2)') time_cur, iday, airtemp, prec, interc_can, int_st_can, &
interc_sveg, int_st_sveg, snow, snow_sm, pet, trans_dem, &
pev, aev_s, aev_i, perc(nlay), watot, wetot, wutot, wutot_ev, wutot_r,&
trans_tree,trans_sveg, eva_dem, gp_can, aet, ceppot_can, ceppot_sveg
endif
deallocate (upt)
END subroutine soil_wat
!**************************************************************
SUBROUTINE upt_wat
! Water uptake by roots
use data_simul
use data_evapo
use data_soil
use data_stand
use data_par
use data_species
use data_climate
implicit none
real, dimension(1:anz_coh) :: tr_dem ! auxiliary arrays for cohorts
real wat_ava, hdem, frtrel, frtrel_1, hupt, hupt_c, totfrt_2, hv, demand_mistletoe_canopy
real wat_at ! total available water per layer
real wat_ar ! total available water per layer with uptake resistance
real hred ! resistance coefficient
real, external :: fred1, fred2, fred3, fred4, fred5, fred6, fred7, fred11
integer i, ianz, j, nroot3
! Calculation of Water Demand
ianz = anz_coh
tr_dem = 0
hdem = 0
hv = pet-aev_i-pev_s
if (hv .lt. 0.) hv = 0.
select case (flag_eva)
case (0,1,3)
if((pet .gt. 0.) .and. (hv .gt. 0.)) then
trans_dem = hv * alfm * (1. - exp(-gp_tot/gpmax)) ! pet (potential evapotranspiration) is reduced by intereption evaporation and potential ground evaporation
else
trans_dem = 0.0
endif
if (gp_tot .gt. zero) then
hdem = trans_dem / gp_tot
else
hdem= 0.
endif
case (8)
! potential transpiration demand of each cohort
if ((gp_tot .gt. zero) .and. (hv .gt. 0.)) then
hdem = (pet-aev_i-aev_s) / gp_tot
else
hdem= 0.
endif
case (2,4)
trans_dem = 0.
case (6,16,36)
! Eucalyptus
hv = pet
if((pet .gt. 0.) .and. (hv .gt. 0.)) then
trans_dem = hv * alfm * (1. - exp(-gp_tot/gpmax))
else
trans_dem = 0.
endif
! preparation: potential transpiration demand of each cohort
if (gp_tot .gt. zero) then
hdem = trans_dem / gp_tot
else
hdem= 0.
endif
case (7,17,37)
trans_dem = hv
! potential transpiration demand of each cohort
if (gp_tot .gt. zero) then
hdem = trans_dem / gp_tot
else
hdem= 0.
endif
end select
hdem = max(0., hdem)
! Distribution of total Demand into Demands of Cohorts
!extraction of demand of mistletoe cohort (for case flag eva = 1,3,6,7...)
zeig => pt%first
do while (associated(zeig))
if (zeig%coh%species.eq.nspec_tree+2) then
demand_mistletoe_canopy=zeig%coh%gp * zeig%coh%ntreea * hdem
end if
zeig => zeig%next
enddo
zeig => pt%first
i = 1
do while (associated(zeig))
select case (flag_eva)
case (0, 1, 3, 6, 7, 16, 17, 36, 37)
!uppermost tree cohort (with flag mistletoe) gets additinal demand of mistletoe
if (zeig%coh%mistletoe.eq.1) then
zeig%coh%demand = zeig%coh%gp * zeig%coh%ntreea * hdem + demand_mistletoe_canopy
elseif (zeig%coh%species.eq.nspec_tree+2) then ! set demand of mistletoe to zero as it will be fullfilled by the tree
zeig%coh%demand=0. ! set to zero because demand has been added to the infested tree cohort
else
zeig%coh%demand = zeig%coh%gp * zeig%coh%ntreea * hdem ! all other cohorts get their demand
end if
case (2,4)
!uppermost tree cohort (with flag mistletoe) gets additinal demand, that of mistletoe
if (zeig%coh%mistletoe.eq.1) then
zeig%coh%demand = (max(0., zeig%coh%demand - zeig%coh%aev_i) + demand_mistletoe_cohort)
endif
if (zeig%coh%species.eq.nspec_tree+2) then ! set demand of mistletoe to zero as it will be fullfilled by the tree
zeig%coh%demand=0.
endif
if (zeig%coh%mistletoe.ne.1 .AND. zeig%coh%species.ne.nspec_tree+2) then
zeig%coh%demand = max(0., zeig%coh%demand - zeig%coh%aev_i)
end if
trans_dem = trans_dem + zeig%coh%demand
end select
tr_dem(i) = zeig%coh%demand ! demand of transpiration per cohort
i = i + 1
zeig => zeig%next
enddo ! zeig (cohorts)
! Calculation of Water Supply
frtrel_1 = 1.
select case (flag_wurz)
case (0)
if (nroot_max .gt. 5) then
nroot3 = 5
else
nroot3 = nroot_max
endif
case default
nroot3=nroot_max
end select
! layers with seedlings
do j = 1,nroot3
! determination of resisctance coefficient
select case (flag_wred)
case(1)
hred = fred1(j)
case(2)
hred = fred2(j)
case(3)
hred = fred3(j)
case(4)
hred = fred4(j)
case(5)
hred = 1.
case(6)
hred = 0.5
case(7)
hred = 0.25
case(8)
hred = fred6(j)
case(10)
hred = fred7(j)
case(11)
hred = fred11(j)
end select
wat_res(j) = hred
if (temps(j) .gt. -0.3) then
wat_at = max(wats(j) - wilt_p(j), 0.) ! total available water per layer
wat_ar = hred * wat_at ! total available water per layer with uptake resistance
hupt = 0.
else
wat_ar = 0. ! frost
wat_at = 0.
hupt = 0.
endif
! Distribution of Water Supply into the Cohorts
! Distribution of Fine Roots
zeig => pt%first
i = 1 ! cohort index
do while (associated(zeig))
if (zeig%coh%species .ne. nspec_tree+2) then ! not for mistletoe
frtrel = zeig%coh%frtrelc(j)
wat_ava = frtrel * wat_ar ! available water per tree cohort and layer
if (wat_ava .ge. tr_dem(i)) then
hupt_c = tr_dem(i)
tr_dem(i) = 0.
else
hupt_c = wat_ava
tr_dem(i) = tr_dem(i) - wat_ava
endif
select case (flag_dis) !xylem clogger disturbance
case (1,2)
hupt_c = hupt_c * xylem_dis
end select
xwatupt(i,j) = hupt_c ! water uptake per cohorte and layer
zeig%coh%supply = zeig%coh%supply + hupt_c
if (zeig%coh%supply .lt.0.) then
continue
endif
hupt = hupt + hupt_c
i = i + 1
end if ! exclusion of mistletoe
zeig => zeig%next
enddo ! zeig (cohorts)
wupt_r(j) = hupt
enddo ! j
! layers without seedlings
if (totfrt_p.gt.(seedlfrt+zero)) then
totfrt_2 = 1./(totfrt_p-seedlfrt)
do j = nroot3+1, nroot_max
! determination of resisctance coefficient
select case (flag_wred)
case(1)
hred = fred1(j)
case(2)
hred = fred2(j)
case(3)
hred = fred3(j)
case(4)
hred = fred4(j)
case(5)
hred = 1.
case(6)
hred = 0.5
case(7)
hred = 0.25
end select
wat_res(j) = hred
if (temps(j) .gt. -0.3) then
wat_at = max(wats(j) - wilt_p(j), 0.) ! total available water per layer
wat_ar = hred * wat_at ! total available water per layer with uptake resistance
hupt = 0.
else
wat_ar = 0.
endif
zeig => pt%first
i = 1 ! cohort index
do while (associated(zeig))
frtrel = zeig%coh%frtrelc(j)
wat_ava = frtrel * wat_ar ! available water per tree cohort and layer
if (wat_ava .ge. tr_dem(i)) then
hupt_c = tr_dem(i)
tr_dem(i) = 0.
else
hupt_c = wat_ava
tr_dem(i) = tr_dem(i) - wat_ava
endif
select case (flag_dis) !xylem clogger disturbance
case (1,2)
hupt_c = hupt_c * xylem_dis
end select
xwatupt(i,j) = hupt_c
zeig%coh%supply = zeig%coh%supply + hupt_c
hupt = hupt + hupt_c
i = i + 1
zeig => zeig%next
enddo ! zeig (cohorts)
wupt_r(j) = hupt
enddo ! j
endif
END subroutine upt_wat
!**************************************************************
SUBROUTINE upt_wat1
! Water uptake by roots
! 2. Version
use data_simul
use data_evapo
use data_soil
use data_stand
use data_par
implicit none
real, dimension(1:anz_coh) :: tr_dem,frt_rel ! help arrays for cohorts
real wat_ava, hdem, frtrel, frtrel_1, hupt, hupt_c, totfrt_2
real wat_at ! total available water per layer
real wat_ar ! total available water per layer with uptake resistance
real hred ! resistance coefficient
real, external :: fred1, fred2, fred3, fred4, fred5, fred6, fred7, fred11
integer i, ianz, j, nroot3
ianz = anz_coh
tr_dem=0
trans_dem = (pet-aev_i) * alfm * (1. - exp(-gp_can/gpmax)) ! pet NOT reduced by ground evaporation
if (trans_dem .lt. 0.) trans_dem = 0.
! potential transpiration demand of each cohort
if (gp_can .gt. zero) then
hdem = trans_dem / gp_can
else
hdem= 0.
endif
! Estimation of transpiration demand of tree cohorts and total fine root mass
! in layers with and without seedlings
zeig => pt%first
i = 1
do while (associated(zeig))
select case (flag_eva)
case (0, 1, 3)
zeig%coh%demand = zeig%coh%gp * zeig%coh%ntreea * hdem
case (2)
zeig%coh%demand = zeig%coh%demand - zeig%coh%aev_i
end select
tr_dem(i) = zeig%coh%demand
i = i + 1
zeig => zeig%next
enddo ! zeig (cohorts)
! uptake controlled by share of roots
frtrel_1 = 1.
! layers with seedlings
do j = 1,nroot_max
! determination of resisctance coefficient
select case (flag_wred)
case(1)
hred = fred1(j)
case(2)
hred = fred2(j)
case(3)
hred = fred3(j)
case(4)
hred = fred4(j)
case(5)
hred = 1.
case(6)
hred = 0.5
case(7)
hred = 0.25
case(8) ! BKL, ArcEGMO
hred = fred6(j)
case(10)
hred = fred7(j)
end select
wat_at = max(wats(j) - wilt_p(j), 0.) ! total available water per layer
wat_ar = hred * wat_at ! total available water per layer with uptake resistance
hupt = 0.
zeig => pt%first
i = 1 ! cohort index
do while (associated(zeig))
frtrel = zeig%coh%frtrel(j) * zeig%coh%x_frt * zeig%coh%ntreea * totfrt_1
wat_ava = frtrel * wat_ar ! available water per tree cohort and layer
if (wat_ava .ge. tr_dem(i)) then
hupt_c = tr_dem(i)
tr_dem(i) = 0.
else
hupt_c = wat_ava
tr_dem(i) = tr_dem(i) - wat_ava
endif
select case (flag_dis) !xylem clogger disturbance
case (1,2)
hupt_c = hupt_c * xylem_dis
end select
xwatupt(i,j) = hupt_c
zeig%coh%supply = zeig%coh%supply + hupt_c
hupt = hupt + hupt_c
i = i + 1
zeig => zeig%next
enddo ! zeig (cohorts)
wupt_r(j) = hupt
enddo ! j
END subroutine upt_wat1
!**************************************************************
real FUNCTION fred1(j)
! Function for calculating uptake resistance
! from CHEN (1993)
! empirical relation between soil water content and resistance
! fred1=1 if (field_cap - 10%*field_cap) <= wats <= (field_cap + 10%*field_cap)
use data_par
use data_soil
implicit none
real hf, f09, f11, wc, diff
integer j
f09 = 0.9 * field_cap(j)
f11 = 1.1 * field_cap(j)
wc = wats(j)
if (wc .lt. wilt_p(j)) then
hf = 0.
else if (wc .lt. f09) then
diff = f09-wilt_p(j)
if (diff .lt. zero) diff = 0.001
hf = 1. - (f09-wc) / diff
else if (wc .gt. f11) then
diff = pv(j)-f11
if (diff .lt. zero) diff = 0.001
hf = 0.3 + 0.7 * (pv(j)-wc) / diff
if (hf .lt. zero) hf = 0.001
else
hf = 1.
endif
fred1 = hf
END function fred1
!**************************************************************
real FUNCTION fred2(j)
! Function for calculating uptake resistance
! from Aber and Federer (f=0.04 fuer Wasser in cm)
! only 40% of total available water are plant available per day
implicit none
integer j
fred2 = 0.05
END function fred2
!**************************************************************
real FUNCTION fred3(j)
! Function for calculating uptake resistance
! from CHEN (1993)
! modified to a profile defined in fred
! fred3 may be described by a function (old version):
! fred3 = 0.0004*j*j - 0.0107*j + 0.0735
! this case: set from a root profile, defined by input of root_fr
use data_par
use data_soil
implicit none
real hf, f09, f11, wc, diff
! hf is a reduction factor in dependence on water content
real fred(15)
integer j
! uptake reduction depending on water content
f09 = 0.9 * field_cap(j)
f11 = 1.1 * field_cap(j)
wc = wats(j)
if (wc .lt. wilt_p(j)) then
hf = 0.
else if (wc .lt. f09) then
diff = f09-wilt_p(j)
if (diff .lt. zero) diff = 0.001
hf = 1. - (f09-wc) / diff
else if (wc .gt. f11) then
diff = pv(j)-f11
if (diff .lt. zero) diff = 0.001
hf = 0.3 + 0.7 * (pv(j)-wc) / diff
if (hf .lt. zero) hf = 0.001
else
hf = 1.
endif
fred3 = root_fr(j) * hf
END function fred3
!**************************************************************
real FUNCTION fred4(j)
! Function for calculating uptake resistance
! modified to a profile defined in fred
! profile at Beerenbusch
use data_soil
implicit none
real fred(15)
integer j
fred = (/ 0.0, 0.03, 0.03, 0.02, 0.02, 0.02, 0.02, 0.01, 0.01, 0.01, &
0.01, 0.01, 0.01, 0.01, 0.01 /) ! fred fuer Beerenbusch
fred4 = fred(j)
END function fred4
!**************************************************************
real FUNCTION fred6(j)
! Function for calculating uptake resistance
! from Kloecking (2006) simular to fred1
! empirical relation between soil water content and resistance
! fred6=1 if field_cap <= wats <= (field_cap + 10%*field_cap)
use data_soil
implicit none
real hf, f09, f11, wc
integer j
f09 = field_cap(j)
f11 = 1.1 * field_cap(j)
wc = wats(j)
if (wc .le. wilt_p(j)) then
hf = 0.
else if (wc .lt. f09) then
hf = 0.1 + (0.9 *(wc-wilt_p(j)) / (f09-wilt_p(j)))
else if (wc .gt. f11) then
hf = 0.3 + 0.7 * (pv(j)-wc) / (pv(j)-f11)
if (hf .lt. 0.) hf = 0.001
else
hf = 1.
endif
fred6 = hf
END function fred6
!**************************************************************
real FUNCTION fred7(j)
! Function for calculating uptake resistance
! from CHEN (1993)
! empirical relation between soil water content and resistance
! fred1=1 if (field_cap - 10%*field_cap) <= wats <= (field_cap + 10%*field_cap)
use data_par
use data_soil
implicit none
real hf, f09, f11, wc, diff
integer j
f09 = 0.9 * field_cap(j)
f11 = 1.1 * field_cap(j)
wc = wats(j)
if (wc .lt. wilt_p(j)) then
hf = 0.
else if (wc .lt. f09) then
diff = f09-wilt_p(j)
if (diff .lt. zero) diff = 0.001
hf = exp(-5.*(f09-wc) / diff)
else if (wc .gt. f11) then
diff = pv(j)-f11
if (diff .lt. zero) diff = 0.001
hf = 0.3 + 0.7 * (pv(j)-wc) / diff
if (hf .lt. zero) hf = 0.001
else
hf = 1.
endif
fred7 = hf
END function fred7
!**************************************************************
real FUNCTION fred11(j)
! Function for calculating uptake resistance, especially adapted for Mistletoe disturbance
! function after van Wijk, 2000
use data_par
use data_soil
implicit none
real hf, S, f11, wc, diff
integer j
f11 = 1.1 * field_cap(j)
wc = wats(j)
if (wc .lt. wilt_p(j)) then
hf = 0.
else if (wc .lt. field_cap(j)) then
S=(field_cap(j)-wc)/(field_cap(j)-wilt_p(j))
hf = exp(-30*S) !30 = strong reduction in water avail.
else if (wc .gt. f11) then
diff = pv(j)-f11
if (diff .lt. zero) diff = 0.001
hf = 0.3 + 0.7 * (pv(j)-wc) / diff
if (hf .lt. zero) hf = 0.001
else
hf = 1.
endif
fred11 = hf
END function fred11
!**************************************************************
SUBROUTINE take_wat(eva_dem, psi)
! Estimation of water taking out for uncovered soil
use data_soil
use data_simul
implicit none
!input:
real :: eva_dem ! evaporation demand
real :: psi ! covering
integer i, ii, j, ntag ! max. layer of taking out
real, allocatable, dimension(:) :: gj
real, external :: b_r, funcov
real diff, gj_j, depth_j, depth_n, rij, rmax, rr, rs, sr
allocate (gj(nlay))
do i=1,nlay
wupt_ev(i)=0.0
gj(i)=0.0
enddo
ntag = 0
rmax = 0.0
depth_n = depth(n_ev_d)
do i=1,n_ev_d
rij = 0.0
rr = depth_n/depth(i)
sr = 0.0
rs = 0.0
do j=1,i
! depth for uncovered take out
depth_j = depth(j)
gj(j) = FUNCOV(w_ev_d, rs*rr, rr*depth_j)
rs = depth_j
sr = sr + gj(j)
enddo ! i
if (sr.gt.1.E-7) then
sr = 1.0/sr
do j=1,i
! water take out
! (psi = 1.-psi) no soil evaporation in case of total covering
! and maximal evaporation for uncovered soil
gj_j = -B_R(wats(j), field_cap(j), wilt_p(j)) &
* eva_dem * (1.-psi) * gj(j) * sr
gj_j = max(gj_j,0.0)
gj(j)= gj_j
rij = rij + gj_j
enddo ! i
if (rij .gt. rmax) then
rmax = rij
ntag = i
do ii=1,ntag
wupt_ev(ii) = gj(ii)
enddo
endif ! rij
endif ! sr
enddo ! n_ev_d
! balance
do i=1,nlay
diff = wats(i) - wilt_p(i)
if (wupt_ev(i) .gt. diff) then
wupt_ev(i) = diff
endif
enddo ! nlay
deallocate (gj)
END subroutine take_wat
!*******************************************************************************
real FUNCTION B_R(water, f_cap, wilting)
! Reduction function for water taking out (uncovered soil)
implicit none
!input:
real :: water ! water storage
real :: f_cap ! field capacity
real :: wilting ! wilting point
b_r = 1.0
if (water .lt. f_cap) B_R = max((water-wilting)/(f_cap-wilting), 0.0)
END function B_R
!******************************************************************************
real FUNCTION funcov(wt_d, a, bb)
! take out density function for uncovered soil
implicit none
!input:
real :: wt_d ! depth of water taking out by evaporation (cm)
real :: a, bb ! relative upper and lower depth of actual layer
real fk, wt_5, b
fk = .455218234
wt_5 = 0.05 * wt_d
b = min(bb,wt_d)
funcov = (- b + a + 1.05*wt_d*log((b+wt_5)/(a+wt_5)))*fk/wt_d
END function funcov
!******************************************************************************
real FUNCTION wat_new(wat_us, wat_in, ilayer)
! FUNCTION WIEN(WIA,NIST,ALAM,DTI,TT,DICK)
! Estimation of additional water after infiltration and percolation
use data_par
use data_soil
implicit none
! input:
real :: wat_us ! water content in relation to field capacity
real :: wat_in ! water infiltration into actual layer
integer :: ilayer ! number of actual layer
real dti !time step
real awi, b1, b2, la, hsqr, exphelp
dti = 1.
fakt = 0.4
if (fakt .ge. 0.0) then ! percolation?
la = 100.0 * fakt * dti * wlam(ilayer)/thick(ilayer)**2
if (wat_us .le. zero) then ! water near zero?
if (wat_in .le. zero) then ! infiltrated water near zero?
wat_new = wat_us + wat_in
else
if (wat_us+wat_in .gt. zero) then
exphelp = sqrt(la*wat_in) * (1 + wat_us/wat_in)*1
if (exphelp .le.10.) then ! avoid underflow
b1 = -exp(-2. * exphelp)
else
b1 = 0.
endif
wat_new = sqrt(wat_in/la) * (1+b1)/(1-b1)
else
wat_new = wat_us + wat_in
endif
endif ! wat_in
else
if (wat_in .lt. 0.) then
awi = abs(wat_in)
b1 = atan(wat_us/sqrt(awi/la)) / sqrt(la * awi)
if (b1 .gt. 1) then
b2 = sqrt (awi * la)
b1 = sin(b2) / cos(b2)
b2 = sqrt(awi / la)
wat_new = b2 * (wat_us - b2*b1) / (b2 + wat_us*b1)
else
wat_new = wat_in * (1-b1)
endif ! b1
else
if (wat_in .gt. 0.) then
b1 = sqrt(wat_in / la)
hsqr = sqrt(la*wat_in)
if (hsqr .lt. 10.) then
b2 = (wat_us - b1) * exp(-2.* hsqr) / (wat_us + b1)
if (b2 .ge. 1.0) then
b2 = 0.99999
endif
else
b2 = 0.
endif
wat_new = b1 * (1.+b2) / (1.-b2)
else
wat_new = wat_us / (1. + la*wat_us)
endif
endif ! wat_in
endif ! wat_us
else
wat_new = wat_us
endif ! fakt
END function wat_new
!******************************************************************************
SUBROUTINE bucket(bucksize1, bucksize2, buckdepth)
! calculation of bucket size (1m; without humus layer)
use data_soil
implicit none
real bucksize1, & ! bucket size of 1 m depth (nFK)
bucksize2, & ! bucket size of rooting zone
buckdepth, diff
integer j
bucksize1 = 0.
bucksize2 = 0.
buckdepth = 0.
do j=2,nlay
if ((depth(j)-depth(1)) .lt. 100.) then
bucksize1 = bucksize1 + wats(j) - wilt_p(j)
buckdepth = depth(j) - depth(1)
else
diff = 100. - buckdepth
bucksize1 = bucksize1 + (wats(j) - wilt_p(j))*diff/thick(j)
buckdepth = 100.
exit
endif
enddo
do j=2,nroot_max
bucksize2 = bucksize2 + wats(j) - wilt_p(j)
enddo
END subroutine bucket
!******************************************************************************
SUBROUTINE snowpack(snow_sm, p_inf, pev)
! properties of snow
! calculation of soil surface temperature under snow pack
use data_climate
use data_evapo
use data_inter
use data_par
use data_simul
use data_soil
use data_soil_t
implicit none
real p_inf ! infiltrated water
real snow_sm
real pev
real airtemp_sm ! melting temperature
real snow_old ! old snow pack
real tc_snow ! thermal conductivity of snow J/cm/s/K
real thick_snow ! thickness of snow
real dens_snow ! density of snow
real:: dens_sn_new = 0.1 ! density of fresh snow
real fakta
snow_old = snow
!substract evaporation of snowcover from snow in both cases
if (airtemp .lt. temp_snow) then ! frost conditions
snow = snow + prec_stand ! precipitation as snow
snow_sm = 0.0 ! no snow melting
p_inf = 0.0 ! no infiltrated precipitation
pev = max((pev_s - aev_i), 0.) ! interc. evapor. reduces soil evapor.
else
airtemp_sm = max(airtemp, 0.)
snow_sm = airtemp*(0.45+0.2*airtemp) ! snow melting
snow_sm = MIN(snow_sm, snow)
snow = snow - snow_sm
p_inf = prec_stand + snow_sm ! infiltrated precipitation
pev = max((pev_s - aev_i), 0.) ! interc. evapor. reduces soil evapor.
end if ! airtemp
if (snow .ge. zero) then
snow_day = snow_day + 1
days_snow = days_snow + 1
if (pev .le. zero) then
pev = 0.
else
! snow sublimation
aev_s = max(min(snow, pev), 0.)
snow = snow - aev_s
pev = pev - aev_s
endif
! soil surface temperature under snow pack
! snow hight = 0.2598 * water equivalent + 8.6851; adjustment from measurement values (see Bodentemperatur.xls)
if (snow .ge. 0.05) then
dens_snow = dens_sn_new + snow_day*0.025
dens_snow = MIN(dens_snow, 1.)
dens_snow = 0.5*(dens_sn_new*prec_stand + dens_snow*snow_old)/snow
dens_snow = MIN(dens_snow, 1.)
tc_snow = 0.7938*EXP(3.808*dens_snow)*0.001 ! thermal conductivity of snow J/cm/s/K
thick_snow = snow / dens_snow
fakta = tc_snow * 86400. * (thick(1)/2.) / (t_cond(1) * thick_snow) ! s --> day
temps_surf = (0.5*temps(1) + fakta*airtemp) / (1. + fakta) ! CoupModel (Jansson, 2001)
endif
else
snow_day = 0
endif
END subroutine snowpack
!******************************************************************************
SUBROUTINE soil_stress
! Calculation of the stress factors
use data_soil
use data_species
use data_stand
use data_par
implicit none
integer :: i, k
real :: m_1, m_2, n_1, n_2
real :: wratio, wafpo
real, dimension (1:4) :: allstress, xvar, yvar
!temperature stress
do i=1,nlay
do k=1,nspecies
if (temps(i) .ge. spar(k)%tbase) then
svar(k)%tstress(i) = sin((pi/2)*(temps(i)-spar(k)%tbase)/(spar(k)%topt-spar(k)%tbase))
else
svar(k)%tstress(i) = 0.
endif
!soil strength
wratio=0.
if (dens(i) .le. BDopt(i)) then
svar(k)%BDstr(i) = 1
svar(k)%BDstr(i) = 1
elseif (dens(i) .ge. svar(k)%BDmax(i)) then
svar(k)%BDstr(i) = 0
else
svar(k)%BDstr(i) = (svar(k)%BDmax(i)-dens(i))/(svar(k)%BDmax(i)-BDopt(i))
endif
if (watvol(i) .lt. wilt_p_v(i)) then
wratio = 0.
elseif (watvol(i) .gt. f_cap_v(i)) then
wratio = 1.
else
wratio = (watvol(i)-wilt_p_v(i))/(f_cap_v(i)-wilt_p_v(i))
endif
svar(k)%sstr(i)=svar(k)%BDstr(i)*sin(1.57*wratio)
!aeration
wafpo=watvol(i)/pv_v(i)
if (wafpo .ge. svar(k)%porcrit(i)) then
svar(k)%airstr(i) = (1.-wafpo)/(1.-svar(k)%porcrit(i))
else
svar(k)%airstr(i) = 1.
endif
if (svar(k)%airstr(i) .lt. 0.) svar(k)%airstr(i) = 0.
!soil acidity
xvar=(/spar(k)%ph_min, spar(k)%ph_opt_min, spar(k)%ph_opt_max, spar(k)%ph_max/)
yvar=(/0,1,1,0/)
m_1=(yvar(1)-yvar(2))/(xvar(1)-xvar(2))
n_1=yvar(2)-m_1*xvar(2)
m_2=(yvar(3)-yvar(4))/(xvar(3)-xvar(4))
n_2=yvar(4)-m_2*xvar(4)
if (phv(i) .gt. spar(k)%ph_opt_max .and. phv(i) .le. spar(k)%ph_max ) then
svar(k)%phstr(i)=m_2*phv(i)+n_2
elseif (phv(i) .lt. spar(k)%ph_opt_min .and. phv(i) .ge. spar(k)%ph_min ) then
svar(k)%phstr(i)=m_1*phv(i)+n_1
elseif (phv(i) .gt. spar(k)%ph_max .or. phv(i) .lt. spar(k)%ph_min) then
svar(k)%phstr(i)=0.
else
svar(k)%phstr(i)=1.
endif
! total stress (Rstress) is taken as the largest of the four
allstress(1)=svar(k)%tstress(i)
allstress(2)=svar(k)%sstr(i)
allstress(3)=svar(k)%airstr(i)
allstress(4)=svar(k)%phstr(i)
svar(k)%Rstress(i)= minval(allstress)
svar(k)%Smean(i)=svar(k)%Rstress(i)+svar(k)%Smean(i)
enddo
enddo
END subroutine soil_stress
!*******************************************************************************
SUBROUTINE hum_add(xfcap, xwiltp, xpv)
! Soil parameter according to [Kuntze et al., Bodenkunde, 1994], S. 172
use data_simul
use data_soil
use data_soil_cn
implicit none
integer :: i, k
real :: fcapi, clayvi, siltvi, humvi, humvi2, wiltpi, pvi, nfki, hcbc
real, dimension(nlay):: xfcap, xwiltp, xpv ! output of addition mm/dm
xfcap(1) = 0.0
xwiltp(1) = 0.0
xpv(1) = 0.0
do i = 1, nlay
fcapi = 0.
wiltpi = 0.
pvi = 0.
clayvi = clayv(i)
humvi = humusv(i)*100.
humvi2 = humusv(i)*humusv(i)
if (humvi .lt. 15.) then
if (clayvi .le. 0.05) then
wiltpi = 0.0609 * humvi2 + 0.33 * humvi
pvi = 0.0436 * humvi2 + 0.631 * humvi
nfki = -0.0009 * humvi2 + 1.171 * humvi
fcapi = nfki + wiltpi
else if (clayvi .le. 0.12) then
wiltpi = 0.0357 * humvi2 + 0.0762 * humvi
pvi = 0.0441 * humvi2 + 0.5455 * humvi
nfki = 0.0252 * humvi2 + 0.7462 * humvi
fcapi = nfki + wiltpi
else if (clayvi .le. 0.17) then
wiltpi = 0.0374 * humvi2 - 0.1777 * humvi
pvi = 0.0552 * humvi2 + 0.2936 * humvi
nfki = 0.0324 * humvi2 + 0.6243 * humvi
fcapi = nfki + wiltpi
else if (clayvi .le. 0.35) then
wiltpi = 0.0179 * humvi2 - 0.0385 * humvi
pvi = 0.0681 * humvi2 + 0.0768 * humvi
nfki = 0.0373 * humvi2 + 0.3617 * humvi
fcapi = nfki + wiltpi
else if (clayvi .le. 0.65) then
wiltpi = 0.0039 * humvi2 + 0.0254 * humvi
pvi = 0.0613 * humvi2 + 0.0947 * humvi
nfki = 0.0338 * humvi2 + 0.0904 * humvi
fcapi = nfki + wiltpi
else
wiltpi = 0.0
pvi = 0.0613 * humvi2 + 0.0947 * humvi
nfki = 0.0104 * humvi2 + 0.2853 * humvi
fcapi = nfki + wiltpi
endif
else ! humvi > 15
! organic soils
continue
endif ! humvi
xfcap(i) = fcapi
xwiltp(i) = wiltpi
xpv(i) = pvi
enddo
if (flag_bc .gt. 0) then
do i = 1, nlay
if (C_bc(i) .gt. 0.) then
fcapi = f_cap_v(i)
clayvi = clayv(i)
siltvi = siltv(i)
humvi = humusv(i)*100.
hcbc = C_bc(i)*100.*100. / (cpart_bc(y_bc_n) * dmass(i))
if ((clayvi .le. 0.17) .and. (siltvi .le. 0.5)) then ! sand
fcapi = 0.0619 * hcbc
wiltpi = 0.0375 * hcbc
nfki = 7.0
elseif ((clayvi .le. 0.45) .and. (siltvi .gt. 0.17)) then ! loam
fcapi = 0.015 * hcbc
wiltpi = 0.0157 * hcbc
nfki = 10.
else ! clay
fcapi = -0.0109 * hcbc
wiltpi = -0.0318 * hcbc
nfki = 16.
endif
xfcap(i) = xfcap(i) + fcapi
xwiltp(i) = xwiltp(i) + wiltpi
endif
enddo
endif
END subroutine hum_add
!*******************************************************************************
SUBROUTINE bc_appl
! application of biochar
use data_out
use data_simul
use data_soil
use data_soil_cn
implicit none
character :: text
integer :: ios, inunit, j
logical :: ex
real :: hcbc
call testfile(valfile(ip),ex)
IF (ex .eqv. .true.) then
inunit = getunit()
ios=0
open(inunit,file=valfile(ip),iostat=ios,status='old',action='read')
if (.not.flag_mult8910) then
print *,'***** Reading application values of biochar from file ',valfile(ip),'...'
write (unit_err, *) 'Application values of biochar from file ',trim(valfile(ip))
endif
do
read(inunit,*) text
IF(text .ne. '!')then
backspace(inunit)
exit
endif
enddo
read (inunit,*,iostat=ios) n_appl_bc
allocate (C_bc_appl(n_appl_bc))
allocate (N_bc_appl(n_appl_bc))
allocate (bc_appl_lay(n_appl_bc))
allocate (cnv_bc(n_appl_bc))
allocate (dens_bc(n_appl_bc))
allocate (cpart_bc(n_appl_bc))
allocate (y_bc(0 : n_appl_bc + 1))
y_bc = 0
C_bc_appl = 0.
N_bc_appl = 0.
do j = 1, n_appl_bc
read (inunit,*,iostat=ios) y_bc(j), cpart_bc(j), cnv_bc(j), dens_bc(j)
read (inunit,*,iostat=ios) bc_appl_lay(j), C_bc_appl(j)
enddo
endif ! ex
END subroutine bc_appl
!*******************************************************************************