subroutine hydrotope_process(j, jea, k, n, wet, bSubcatch, dart, daycounter, flu, frar, ida, iy, iyr, mo, nbyr, precip, sbar, tx)
!**** PURPOSE: THIS SUBROUTINE CALCULATES ALL PROCESSES in HYDROTOPES
!**** CALLED IN: SUBBASIN
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! PARAMETERS & VARIABLES
!
! >>>>> COMMON PARAMETERS & VARIABLES
! bcv(j, je) = lag factor (residue and snow effect on temp)
! canstor(j, je) = canopy water storage, mm
! cn = current CN
! cn2(k, n, j) = Curve Numbers for soil k, land use n, and subbasin
! cva(j, je) = vegetation cover, kg/ha
! dflow(j, je, 20) = monthly flows for water and N (see writhru.f)
! eo = potential evapotranspiration, mm
! ep = plant transpiration, mm
! es = soil evaporation, mm
! icc = index for cover crop corr. number in crop database
! ida = current day
! idfe(n, if) = day of fertilization, if - number of fertilisation
! FORD = index for deciduous forest corr. number in crop database
! FORE = index for evergreen forest corr. number in crop database
! FORM = index for mixed forest corr. number in crop database
! RNGB = index for heather corr. number in crop database
! imea = index for meadows corr. number in crop database
! ipas = index for pasture corr. number in crop database
! iwet = index for wetland corr. number in crop database
! iwetf = index for forested wetland corr. number in crop database
! nn = number of soil layers
! nveg(j, je) = number of vegetation from crop database
! percn = N-NO3 leaching to g-w, kg/ha
! precip = precipitation, mm
! preinf(j, je) = precipitation adjusted for canopy storage, mm
! prk = lateral subsurface flow, mm, calc in perc for 4mm layers
! qd = surface flow in HYDROTOPE, mm
! qi = surface flow in HYDROTOPE, mm
! rain = preinf(j, je) - qd, mm
! sep = percolation, mm
! sml = snow melt, calc in snom(), mm
! snoa(j, jea) = snow water content in HYDROTOPE, mm
! snoev = snow evap. in HYDROTOPE, mm
! ssf = subsurface flow in HYDROTOPE, mm
! ssfn = N-NO3 in subsurface flow, kg/ha
! ste(j, je, l) = water storage, recalc here, mm
! strsn = N stress for plant growth
! strsp = P stress for plant growth
! sumfc(k) = sum of field capacity in soil, mm
! swe(j, je) = soil water content, mm
! swind = soil water index for hydrotope: swind=swe()/sumfc()
! ts = temp. stress
! tx(j) = daily average temperature in the subbasin, degree C
! uno3 = N uptake by plants, kg/ha
! vb = CN max - for output in main
! vl = CN min - for output in main
! ws(j, je) = water stress factor
! yno3 = N-NO3 loss with surface flow, kg/ha
! ysp = soluble P leaching, kg/ha
! >>>>>
! >>>>> STATI! PARAMETERS
! ii = local par
! l = local par
! xx = local par
! zz = local par
! >>>>>
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
use crop, only : &
be, &
crop_initialise_hydrotope, &
crop_process, &
cva, &
fen, &
feno, &
fep, &
fon, &
fop, &
huharv, &
hun, &
icc, &
idfe, &
igro, &
mfe, &
nucr, &
rdmx, &
rwt, &
snup, &
spup
use erosion, only : erosion_cklsp_factor, yone, yphe
use evapotranspiration, only : &
canstor, &
eo, &
es, &
evapotranspiration_process, &
humi, &
pit, &
ra, &
snoev, &
tmn, &
ylc, &
yls
use groundwater, only : additionalGwUptake
use landuse, only : &
LULC, &
landuse_index, &
landuse_is_cropland, &
landuse_is_forest, &
landuse_is_natural_vegetation
use management, only : &
TSubbasin, &
TWaterUser, &
bWAM_Module, &
management_is_active_period, &
management_is_transfer_subbasin, &
management_subbasin_pointer, &
management_user_pointer
use nutrient, only : &
dflow, &
nutrient_fertilisation, &
nutrient_leaching, &
percn, &
ssf, &
ssfn, &
strsn, &
strsp, &
uap, &
uno3, &
yno3, &
ysp
use snow, only : &
bSnowModule, &
precipe, &
rnew, &
rsn, &
sml, &
snoa, &
tmax, &
tmax_tmp, &
tmin, &
tmin_tmp, &
tmit, &
tmx, &
tx_tmp
use soil, only : &
avt, &
bd, &
cbn, &
cn, &
cn2, &
cncor, &
ek, &
fc, &
flate, &
nn, &
poe, &
prk, &
psp, &
qd, &
rtn, &
sep, &
soil_curve_number, &
soil_curve_number_runoff, &
soil_process, &
soil_temperature, &
ste, &
sumfc, &
te, &
ul, &
z
use vegetation, only : &
alai, &
blai, &
ep, &
olai, &
rsd, &
to, &
ts, &
vegetation_process, &
vegetation_store_output, &
ws
logical, intent(in) :: bSubcatch
real(dp), dimension(:), intent(in) :: dart
integer, intent(in) :: daycounter
real(dp), dimension(:), intent(in) :: flu
real(dp), dimension(:, :), intent(in) :: frar
integer, intent(in) :: ida
integer, intent(in) :: iy
integer, intent(in) :: iyr
integer, intent(in) :: mo
integer, intent(in) :: nbyr
real(dp), intent(in) :: precip
real(dp), dimension(:), intent(in) :: sbar
real(dp), dimension(:), intent(in) :: tx
! subbasin number
integer :: j
! number of hydrotope in subbasin
integer :: jea
! land use number from LULC map and * .lut file
integer :: n
! soil number
integer :: k
! wetland code 0 / 1
integer :: wet
integer :: ii
real(dp) :: xx, zz
! the curve number values passed to function curno
real(dp) :: cn2curno
!#### WATER MANAGEMENT MODULE ####
real(dp) :: area ! hydrotop area, m2
!#### WATER MANAGEMENT MODULE ####
TYPE (TSubbasin), POINTER :: pS ! pointer on subbasin
!#### WATER MANAGEMENT MODULE ####
TYPE (TWaterUser), POINTER :: pWU ! pointer on actual TWU
!###########################
!#### SNOW MODULE ####
!###########################
if (bSnowModule) then
tx_tmp = tmit
else
tx_tmp = tx(j)
end if
!###########################
!**** INITIALIZATION
qd = 0.
qi = 0.
ssf = 0.
sep = 0.
ep = 0.
es = 0.
snoev = 0.
ysp = 0.
yno3 = 0.
uno3 = 0.
ssfn = 0.
percn = 0.
ws(j, jea) = 1.
ts = 1.
strsn = 1.
strsp = 1.
!**** INITIZALIZE MANAGEMENT OPERATIONS
if (landuse_is_cropland(n) ) then
call crop_initialise_hydrotope(j, jea, iy, mstruc)
end if
!#### CALL CURNO: to set Curve Number parameters
if (bSubcatch) then
cn2curno = 0.
cn2curno = cn2(j, jea) * cncor(j)
if (cn2curno * cncor(j) > 100.) cn2curno = 100.
if (cn2curno * cncor(j) < 1.) cn2curno = 1.
call soil_curve_number(cn2curno, j, jea, hsumfc, hsumul)
else
call soil_curve_number(cn2(j, jea), j, k, hsumfc, hsumul)
end if
!**** CALC bcv() - lag factor for soil temperature
bcv(j, jea) = cva(j, jea) / (cva(j, jea) + exp(7.563-1.297e-4 * cva(j, jea)))
if (snoa(j, jea) .gt. 0.) then
if (snoa(j, jea) .le. 120.) then
xx = snoa(j, jea) / (snoa(j, jea) + exp(6.055 - .3002 * snoa(j, jea)))
else
xx = 1.
end if
bcv(j, jea) = amax1(real(xx, 4), real(bcv(j, jea), 4))
end if
if (bSnowModule) then
tx_tmp = tmit
tmax_tmp = tmax
tmin_tmp = tmin
else
tx_tmp = tx(j)
tmax_tmp = tmx(j)
tmin_tmp = tmn(j)
end if
!#### CALL SOLT - COMPUTE TEMP. of SOIL LAYERS
call soil_temperature(zz, j, jea, k, bcv, ida, mo, pit, preinf, swe, tmax_tmp, tmin_tmp, tx_tmp)
!#### CALL VOLQ TO CALC RUNOFF VOLUME,
!#### CALL ECKLSP TO CALC COMBINED CKLSP factor for hydrotope
! preinf = precipitation adjusted for canopy storage
! = precipitation + snow melt + canopy water storage
!###########################
!#### SNOW MODULE ####
!###########################
if (bSnowModule) then
!preinf(j, jea) = precipe + canstor(j, jea)
!preinf(j, jea) = precipe + canstor(j, jea) + vsn
preinf(j, jea) = precipe + canstor(j, jea)
!if ( vsn > 0. ) preinf(j, jea) = vsn + canstor(j, jea)
else
preinf(j, jea) = precip + sml + canstor(j, jea)
end if
!preinf(j, jea) = precip + sml + canstor(j, jea)
!###########################
!#################################
!#### WATER MANAGEMENT MODULE ####
!#################################
if (bWAM_Module .AND. daycounter > 1) then
if (management_is_transfer_subbasin(j) ) then
if (landuse_is_cropland(n) .AND. mstruc(j, jea, 7) >= 1 ) then ! if hydrotope is irrigated cropland
area = frar(j, jea) * sbar(j)
pS => management_subbasin_pointer(j)
!pWU => TWU(pS%pos_irr)
pWU => management_user_pointer(pS % pos_irr) ! pointer on current water user
! if sprinkler irrigation, add available irrigation volume to precipitation
! supply in m3/s, precipitation in mm
! pWU%supplied is the total amount of water for all irrigation hydrotopes within this subbasin.
! By converting to mm, the area weighted volume is maintained.
if (pWU % wu_opt == 4 .AND. pWU % irr_practice == 1 .AND. management_is_active_period(iyr, ida, pWU) ) &
preinf(j, jea) = preinf(j, jea) + pWU % supplied(daycounter - 1) * 1000. * 86400. / area !pWU % area
end if
end if
end if
!#################################
if (tx_tmp > 0.) then
if (preinf(j, jea) > 0. ) then
call soil_curve_number_runoff(j, jea, alai, blai, canstor, igro, nucr, preinf, LULC%canmx(landuse_index(n)))
qi = qd
call erosion_cklsp_factor(j, jea, k, cva, ek, igro, nucr, landuse_is_cropland(n), landuse_is_natural_vegetation(n), landuse_is_forest(n))
if (cn .ge. vl) then
if (cn .ge. vb) &
vb = cn
else
vl = cn
end if
end if
end if
!#### CALL PURK TO PERFORM SOIL WATER ROUTING
rain = preinf(j, jea) - qd
!#################################
!#### WATER MANAGEMENT MODULE ####
!#################################
if (bWAM_Module .AND. daycounter > 1) then
if (management_is_transfer_subbasin(j) ) then
if (landuse_is_cropland(n) .AND. mstruc(j, jea, 7) >= 1 ) then ! if hydrotope is irrigated cropland
area = frar(j, jea) * sbar(j)
pS => management_subbasin_pointer(j)
pWU => management_user_pointer(pS % pos_irr) ! pointer on current water user
! if drip irrigation, add available irrigation volume to precipitation, mm
! supply in m3/s, precipitation in mm
! pWU%supplied is the total amount of water for all irrigation hydrotopes within this subbasin.
! By converting to mm, the area weighted volume is maintained.
if (pWU % wu_opt == 4 .AND. pWU % irr_practice == 2 .AND. management_is_active_period(iyr, ida, pWU) ) &
rain = rain + pWU % supplied(daycounter - 1) * 1000. * 86400. / area !pWU % area
end if
end if
end if
!#################################
call soil_process(j, jea, k, rain, swe)
ssf = prk
!#### CALC swe = soil water content [mm]
! swe is re-calculated after vegetation growth below
! ste = water storage per layer [mm], calculated in purk
swe(j, jea) = 0.
swe(j, jea) = sum(ste(j, jea, :))
call evapotranspiration_process(j, jea, k, alai, cva, ep, ida, mo, nn, preinf, qd, snoa, ste, tx, z, bSnowModule, rnew, tmit, tx_tmp, rsn, LULC%ETcor(landuse_index(n)))
! re-calc soil water content
swe(j, jea) = 0.
swe(j, jea) = sum(ste(j, jea, :))
!#### CALL FERT, CRPMD & VEGMD to calculate vegetation growth & fertilisation
select case(LULC % lutype(landuse_index(n)))
case (0) !#### ! NO VEGETATION - - > plant transpiration = 0.
ep = 0.
case (1) !#### CROPLAND, managed
! Fertilizer application
do ii = 1, mfe
if (ida == idfe(ii) ) call nutrient_fertilisation(j, jea, ii, fen, feno, fep)
end do
call crop_process(j, jea, k, n, wet, additionalGwUptake, avt, bWAM_Module, dart, daycounter, es, fc, flu, frar, humi, ida, iy, iyr, mstruc, nbyr, nn, nveg, pit, ra, sbar, sep, ste, tmn, tx, uap, ylc, yls, z, bSnowModule, tmit)
call nutrient_leaching(j, jea, k, bd, cbn, flate, flu, fon, fop, frar, nn, poe, preinf, psp, qd, rsd, rtn, ste, te, ul, yone, yphe)
case(2) !#### NATURAL VEGETATION
nveg(j, jea) = LULC % veg_code(landuse_index(n))
call vegetation_process(j, jea, k, n, wet, additionalGwUptake, avt, bWAM_Module, be, cva, dart, daycounter, fc, flu, frar, huharv, humi, hun, icc, ida, iy, iyr, mstruc, nbyr, nn, nucr, nveg, ra, rdmx, rwt, sbar, sep, snup, spup, ste, sumfc, swe, tmn, to, tx, uap, z, alai, olai)
call nutrient_leaching(j, jea, k, bd, cbn, flate, flu, fon, fop, frar, nn, poe, preinf, psp, qd, rsd, rtn, ste, te, ul, yone, yphe)
case(3) !#### WATER
es = eo
case(4) !#### FORESTS
nveg(j, jea) = LULC % veg_code(landuse_index(n))
call vegetation_process(j, jea, k, n, wet, additionalGwUptake, avt, bWAM_Module, be, cva, dart, daycounter, fc, flu, frar, huharv, humi, hun, icc, ida, iy, iyr, mstruc, nbyr, nn, nucr, nveg, ra, rdmx, rwt, sbar, sep, snup, spup, ste, sumfc, swe, tmn, to, tx, uap, z, alai, olai)
call nutrient_leaching(j, jea, k, bd, cbn, flate, flu, fon, fop, frar, nn, poe, preinf, psp, qd, rsd, rtn, ste, te, ul, yone, yphe)
end select
call vegetation_store_output(j, jea)
!**** RE-CALC SOIL WATER CONTENT swe() & SOIL WATER INDEX swind
swe(j, jea) = 0.
swe(j, jea) = sum(ste(j, jea, :))
swind = swe(j, jea) / sumfc(k)
! dflow(j, je, 20) = monthly flows for water and N (see writhru.f)
dflow(j, jea, 1) = dflow(j, jea, 1) + qd
dflow(j, jea, 2) = dflow(j, jea, 2) + ssf
dflow(j, jea, 3) = dflow(j, jea, 3) + sep
dflow(j, jea, 4) = dflow(j, jea, 4) + es + ep
end subroutine hydrotope_process