Add supply from water user(s) to subbasins' hydrograph storage location (=j) If the subbasin is a reservoir, this step will be overwritten by the reservoir module
if external inputs from water user(s) are greater than 0 add supply of previous day to subbasins' hydrograph storage location (surface runoff component)
Withdraw water from subbasin outlet but only if subbasin is a headwater .and. not a reservoir. otherwise withdawals are taken during the routing process in subroutine 'add' Moreover, the water users %supplied and %tr_losses are calculated.
Check if current subbasin is a headwater...in this case, bWAMHydrograph is .true. at the array position of the subbasin number ATTENTION: The values of sda(2, j) and sda(8, j) may be changed!
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | ihout | |||
| integer, | intent(in) | :: | inum1 | |||
| logical, | intent(in) | :: | bSubcatch | |||
| real(kind=dp), | intent(in) | :: | da | |||
| real(kind=dp), | intent(in) | :: | da9 | |||
| integer, | intent(in) | :: | daycounter | |||
| integer, | intent(in) | :: | ida | |||
| integer, | intent(inout) | :: | ieap | |||
| integer, | intent(in) | :: | iy | |||
| integer, | intent(in) | :: | iyr | |||
| integer, | intent(in) | :: | mo | |||
| integer, | intent(in) | :: | nbyr | |||
| integer, | intent(in) | :: | nd |
subroutine runsubbasin(ihout, inum1, bSubcatch, da, da9, daycounter, ida, ieap, iy, iyr, mo, nbyr, nd)
!**** PURPOSE: THIS SUBROUTINE PERFORMS SUBBASIN OPERATIONS:
! 1. ESTABLISHES INITIAL CONDITIONS IN HYDROTOPES,
! 2. CALLS HYDROTOP TO CALCULATE HYDROTOP PROCESSES,
! 3. AGGREGATES HYDROTOP OUTPUTS,
! 4. CALCULATES EROSION IN SUBBASIN,
! 5. CALCULATES RETENTION OF NUTRIENTS,
! 6. CALCULATES VARIABLES FOR OUTPUT, and
! 7. SETS OUTPUTS FROM SUBBASIN FOR ROUTING
!**** CALLED IN: MAIN
! icode = 1, ihout = inum1 = inum2 = j
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! PARAMETERS & VARIABLES
!
! >>>>> INDICES
! i, ida = DAY
! j = SUBBASIN NR
! jea, je = HYDROTOPE NR
! k = SOIL NR
! n = LAND USE NR
! l = SOIL LAYER NR
! icr, iv = CROP or VEGETATION NR
! ih = HYDR. STORAGE LOCATON NR
! io = CROP OPERATION NR
! if = CROP FERTILISATION NR
! >>>>>
!
! >>>>> COMMON PARAMETERS & VARIABLES
! af = 1000. * da, km2
! aff = af * flu(j) = 1000 * da * flu(j), unit transf. coef.
! alai(j, je) = leaf area index
! alai0(n) = initial leaf area index
! ano3(j, je, l) = wno3(l, k) N-NO3 in sub j, HRU je, layer l, kg/ha
! anora(j, je, l) = wmn(l, k) active organic N in sub j, HRU je, layer l,
! kg/ha
! anors(j, je, l) = wn(l, k) stable organic N in sub j, HRU je, layer l,
! kg/ha
! ap(l, k) = labile (soluble) phosphorus, kg/ha (initialisation)
! bd(k) = bulk density of the upper (!) soil layer,
! calc from porosity, g/cm3
! canev = canopy evaporation, mm
! cklsp(j, je) = combined c, k, ls, p factor
! cn = current CN in hydrotope
! css(j) = sediment conc in return flow, ppm
! cv(j) = initial land cover, kg/ha
! cva(j, je) = vegetation cover, kg/ha
! da = basin area, km2
! dart(j) = drainage area for subbasin, km2
! degNgrw = N degradation in groundwater, (1/day)
! degNsub = N degradation in subsurface flow, (1/day)
! degNsur = N degradation in surface flow, (1/day)
! degPsur = P degradation in surface flow, (1/day)
! dm(j, je) = total biomass, kg/ha
! dur = flow duration, h
! eo = potential evapotranspiration, mm
! eopot = potential evapotranspiration, mm
! ep = plant transpirationm, mm
! es = soil evaporation, mm
! evasum(j, je) = SUM(ep + es + canev), mm
! fc(l, k) = field capacity, mm
! flu(j) = fraction of subbasin area in the basin
! frar(j, je) = fractional area of HYDROTOPE in SUBBASIN
! g(j, je) = fraction of heat unit accumulated
! gwchrg = groundwater recharge, mm
! gwq(j) = groundwaterw contribution to stream, mm
! gwseep = groundwater seepage, mm
! hsumfc(j, je) = sum of field capacity in soil, calc in subbasin, mm
! hsumul(j, je) = sum of upper limit water content in soil,
! calc in subbasin, mm
! hwss(2, j, je) = fun(field cap), calc i subbasin from hsumfc(j, je)
! ida = current day
! gis_m = switch code to write monthly results (water & crops) for GIS output
! gis_y = switch code to write annual results (water & crops) for GIS output
! gis_ave = switch code to write average annual results (water & crops) for GIS output
! ih1, ih2, ih3, ih4 = hydrotopes for HYDROTOPE PRINTOUT
! isb1, isb2, isb3, isb4 = subbasins for HYDROTOPE PRINTOUT
! isu1, isu2, isu3 = subbasins for SUBBASIN PRINTOUT
! iv = index for channel parameters chl(, ), chw(, ), chk(, )
! iy = current year as counter (1, ..., nbyr)
! iyr = current year
! mstruc(j, je, 2) = HRU structure vector to define k, n
! ndpri = day to write crop yield for GIS output
! neap(j) = number of HYDROTOPEs in subbasin, calc readbas
! nn = number of soil layers
! ns(k) = number of soil layers
! nsa(k) = number of layers for arable soils
! nucr(j, je) = crop number (database)
! op(l, k) = stable min. P, estimated from pmn(), kg/ha
! parsz(5, j) = particle size distribution, calc in subbasin
! percn = N leaching to g-w, kg/h
! plab(j, je, l) = ap(l, k) labile (soluable) P in sub j, HRU je,
! layer l, kg/ha
! pma(j, je, l) = pmn(l, k) active mineral P in sub j, HRU je,
! layer l, kg/ha
! pmn(l, k) = act. min. P, estimated in readsol, kg/ha
! pms(j, je, l) = op(l, k) stable mineral P in sub j, HRU je,
! layer l, kg/ha
! porg(j, je, l) = wpo(l, k) organic P in sub j, HRU je, layer l, kg/ha
! pr = peak runoff rate, m3/sec.
! gwrsum(j, je) = annual sum of percolation for GIS output, SUM(sep), mm
! precip = precipitation in subbasin, mm
! presum(j, je) = annual sum of precip for GIS output, mm
! psz(5, k) = particle size distribution
! q1 = daily surface runoff = xqd before call tran, mm
! qd = daily surface runoff in hydrotope, mm
! qi = daily surface runoff in hydrotope, mm
! qtl = transmission losses in subbasin, mm
! ra(j) = radiation in J/cm^2
! rd(j, je) = root depth, mm
! retNgrw = N retention in groundwater, days
! retNsub = N retention in subsurface flow, days
! retNsur = N retention in surface flow, days
! retPsur = P retention in surface flow, days
! revap = revaporation from groundwater, mm
! revapst(j) = transmission losses in the reach, calc in route, m3
! rsd(j, je, 2) = crop residue, kg/ha
! runsum(j, je) = annual sum of surface+subsurface runoff for GIS output,
! SUM(qd+ssf), mm
! sbp(j) = precipitation in subbasin, mm
! sda(1:10, j) = subbasin outputs for routing, analogue to varoute()
! varoute(1:10, ihout) = sda(1:10, j)
! sep = percolation to g-w, mm
! sml = snow melt, calc in snom, mm
! smq(j) = SUM(xqd) - surface runoff for subbasin, mm
! smsq(j) = SUM(xssf) - subsurface runoff for subbasin, mm
! sno(j) = accumulated value of snow cover, mm
! snoa(j, je) = snow water content in sub j, HRU je, mm
! snoev = snow evaporation, mm
! snow = precipitation as snow, mm
! ssf = subsurface flow in HYDROTOPE, mm
! ssfn = N loss with subsurface flow, kg/ha
! ste(j, je, l) = water content in sub j, HRU je, layer l, mm
! stin(l, k) = soil water content in soil layer, mm
! sub(1:28) = basin outputs (weighted sums) daily, see desc. below
! sub(1) = SUM(precip) precipitation, mm
! sub(2) = SUM(precip), if tx()<0 snowfall, mm
! sub(3) = SUM(sml) snowmelt(mm)
! sub(4) = SUM((tmx*flu) average "weighted" max temp., degree C
! sub(5) = SUM(tmn*flu) average "weighted" min temp., degree C
! sub(6) = SUM(ra*flu) radiation, J/cm^2
! sub(7) = SUM(sumcn*flu) curve number
! sub(8) = SUM(xqi*flu) surface runoff, mm
! sub(9) = SUM(xssf*flu) subsurface runoff, mm
! sub(10) = SUM(qtl*flu) channel transm. losses, mm
! sub(11) = SUM(xsep*flu) percolation, mm,
! sub(12) = SUM(xeo*flu) pot evapotr., mm
! sub(13) = SUM(xet*flu) evapotranspiration, mm
! sub(14) = SUM(snoev*flu) snow evaporation, mm
! sub(15) = SUM(gwq*flu) gr. water flow, mm
! sub(16) = SUM(revap*flu) revap from g-w to soil prof., mm
! sub(17) = SUM(gwseep*flu) g-w percolation, mm
! sub(18) = SUM(gwchrg*flu) g-w recharge, mm
! sub(19) = SUM(xswimd*flu) soil water index for basin
! sub(20) = SUM(wysb*flu) wysb=qi+ssf+gwq(j)-qtl, wat yld, mm
! sub(21) = SUM(yd/(100*da*flu) sed yield, t/ha
! sub(22) = SUM(yon*flu) org. N loss with sed, kg/ha
! sub(23) = SUM(yph*flu) org. P loss with sed, kg/ha
! sub(24) = SUM(ysp*flu) soluble P, kg/ha
! sub(25) = SUM(xyno3*flu) NO3-N in surface runoff, kg/ha
! sub(26) = SUM(xssfn*flu) NO3-N in subsurface runoff, kg/ha
! sub(27) = SUM(xpercn*flu) NO3-N leached to g-w, kg/ha
! sub(28) = SUM(xuno3*flu) NO3-N uptake by plants , kg/ha
! subp(j) = daily precipitation in subbasin, mm (from readcli)
! sumcn = weighted aver. Cur.Num. in subbasin (SUM(cn))
! susb(30, j) = subbasin outputs daily, analogue to sub() above
! swe(j, je) = water content in sub j, HRU je
! swin(k) = soil water content in soil profile, mm
! swind = soil water index for hydrotope:
! swind=swe(j, je)/sumfc(k), 1/1
! sym(j) = SUM(yd) = sum of sediment yield for subbasins, t
! tmn(j) = daily min temp in the subbasin, degree C
! tmx(j) = daily max temp in the subbasin, degree C
! tsnfall = threshold temperature for snow fall
! tx(j) = daily aver temp in the subbasin, degree C
! ul(l, k) = upper limit water content in layer, mm
! varoute(*, ih) = subbasin outputs for routing:
! Name Units Definition
! varoute(2, ih) |(m^3) |surface flow
! varoute(3, ih) |(tons) |sediment
! varoute(4, ih) |(kg) |organic N
! varoute(5, ih) |(kg) |organic P
! varoute(6, ih) |(kg) |nitrate N
! varoute(7, ih) |(kg) |soluble P
! varoute(8, ih) |(m^3) |subsurface + g-w flow
! vo = surface runoff, m3
! wmn(l, k) = active org. N content, kg/ha
! wn(l, k) = stable organic N content, kg/ha
! wno3(l, k) = NO3-N content, kg/ha
! wpo(l, k) = organic P content, kg/ha
! wysb = water yield in subbasin = xqi+xssf+gwq(j)-qtl, mm
! xcklsp = weighted aver. comb C*K*LS*P in subbasin, SUM(cklsp)
! xeo = weighted aver. pot. evapotr. in subbasin, SUM(eo), mm
! xet = weighted aver. act. evapotr. in subbasin,
! SUM(ep+es), mm
! xnorg = weighted aver. st.org.N in lay 1 in subb.,
! SUM(anors), g/t
! xnorgp = weighted aver. st.org.N in lay 1 in subb.,
! SUM(anors), kg/ha
! xpercn = weighted aver. N leaching to g-w in subb., kg/ha
! xporg = weighted aver. org.P in layer 1 in subb.,
! SUM(porg), g/t
! xporgp = weighted aver. org.P in layer 1 in subb.,
! SUM(porg), kg/ha
! xpsed = weighted aver. total P (porg+pms+pma) in subb., g/t
! xpsedp = weighted aver. total P (porg+pms+pma) in subb., kg/ha
! xqd = weighted aver. surface flow in subbasin, SUM(qd), mm
! xqi = weighted aver. surface flow in subbasin, SUM(qd), mm
! xsep = weighted aver. percolation in subbasin, SUM(sep), mm
! xsnoev = weighted aver. snow evap in subbasin, SUM(snoev), mm
! xssf = weighted aver. subsurface flow in subbasin,
! SUM(ssf), mm
! xssfn = weighted aver. N loss with ssf in subbasin,
! SUM(ssfn), mm
! xswind = weighted aver. soil wat. index in subbasin, SUM(swimd)
! xwysb = water yield for subbasin, m3
! xxswind = weighted aver. soil wat. index in the basin
! xyno3 = weighted aver. N loss with qd in subbasin,
! SUM(yno3), kg/ha
! xysp = weighted aver. soluble P leaching in subbasin, kg/ha
! yd = daily soil loss caused by water erosion, t
! yno3 = N loss with surface flow, kg/ha
! yon = org N loss with erosion, kg/ha
! yph = org P loss with erosion, kg/ha
! ysp = soluble P leaching, kg/ha
! >>>>>
! >>>>> STATIC PARAMETERS
! ii = local par
! im = local par
! j = subbasin nr
! jea = hydrotope nr
! jj = local par
! k = soil nr
! l = soil layer nr
! n = land use nr
! wabad = daily water balance
! >>>>>
!
! >>>>> UNITS TRANSFORMATION:
! g/t --> kg/ha *wt1=bden*dg/100., bden=bulk dens., dg=layer thickness
! % --> g/t *10000.
! % --> kg/ha *wt1*10000.
! kg/ha --> g/t 1/wt1 = 100./dg/bden()
! mg/kg = g/t
! mm --> m3 *(dart()*1000)
! kg/ha --> kg *(dart()*100)
! >>>>>
!~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
use crop, only : crop_yield_output, cva, ndgro, ndpri, nucr, ylda
use erosion, only : &
cklsp, &
css, &
erosion_enritchment_ratio, &
erosion_organic_nitrate_loss, &
erosion_phosphorus_loss, &
erosion_soil_loss, &
parsz, &
xcklsp, &
xnorg, &
xnorgp, &
xporg, &
xpsedp
use evapotranspiration, only : &
canev, &
eo, &
eopot, &
es, &
et, &
eta_output_id, &
etp_output_id, &
ra, &
snoev, &
soil_evaporation_output_id, &
tmn
use groundwater, only : &
additionalGwUptake, &
baseflow_output_id, &
bff, &
groundwater_process, &
groundwater_recharge_output_id, &
gwchrg, &
gwq, &
gwseep, &
revap, &
xet
use hydrotope, only : &
frar, &
hsumfc, &
hsumul, &
hydrotope_process, &
mstruc, &
precipitation_output_id, &
qi, &
soil_water_content_output_id, &
swe, &
swind, &
sym, &
tmax_output_id, &
tmean_output_id, &
tmin_output_id, &
tmpNgrw, &
tmpNsub, &
tmpNsur, &
tmpPsur, &
water_yield_output_id
use landuse, only : &
landuse_is_cropland, &
landuse_is_forest, &
landuse_is_natural_vegetation
use management, only : &
TSubbasin, &
bWAM_Module, &
bwamhydrograph, &
management_is_transfer_subbasin, &
management_subbasin_pointer, &
management_transfer_out
use nutrient, only : &
ano3, &
anora, &
anors, &
degNgrw, &
degNsub, &
degNsur, &
degPsur, &
percn, &
plab, &
pma, &
pms, &
porg, &
retNgrw, &
retNsub, &
retNsur, &
retPsur, &
ssf, &
ssfn, &
yno3, &
yon, &
yph, &
ysp
use output, only : &
area_tot_glacier, &
area_tot_snow, &
depth_ave_glacier, &
depth_ave_snow, &
output_store_hydrotope_value, &
output_store_subbasin_value, &
output_store_subbasin_values
use reservoir, only : xwysb, xxswind
use river, only : chk, chl, chw, pet_day, varoute
use snow, only : &
bSnowModule, &
gla, &
glac_acc_mm, &
glac_acc_mm0, &
precip_elev_cor, &
precipe, &
psnow, &
sml, &
smle, &
sno, &
snoa, &
snowVal, &
snow_acc_mm, &
snow_acc_mm0, &
snow_degree_day_melting, &
snow_depth_weq_output_id, &
snow_initialise_subbasin, &
snow_process, &
snowfall_weq_output_id, &
tmax, &
tmelt, &
tmin, &
tmit, &
tmx, &
tsnfall, &
vsn, &
xprecip, &
xprecip_elev_cor, &
xsml, &
xsnow, &
xtmax, &
xtmin, &
xtmit, &
xvsn
use soil, only : &
ap, &
bd, &
cn, &
cv, &
dur, &
fc, &
hrtc, &
hwss, &
iv, &
ms, &
nn, &
ns, &
nsa, &
op, &
percolation_output_id, &
pmn, &
pr, &
psz, &
q1, &
qd, &
rp, &
sc, &
sccor, &
sep, &
snum, &
soil_acc_mm, &
soil_curve_number_alpha, &
soil_curve_number_peak_runoff, &
soil_curve_transmission_losses, &
soil_water_index_output_id, &
ste, &
stin, &
subsurface_runoff_output_id, &
surface_runoff_output_id, &
swin, &
ul, &
vo, &
wmn, &
wn
use soil, only : &
wno3, &
wpo, &
xqd, &
z
use vegetation, only : &
alai, &
almn, &
dm, &
ep, &
g, &
gwqLastday, &
rd, &
transpiration_output_id, &
tsav, &
wsav
integer, intent(in) :: ihout, inum1
logical, intent(in) :: bSubcatch
real(dp), intent(in) :: da
real(dp), intent(in) :: da9
integer, intent(in) :: daycounter
integer, intent(in) :: ida
integer, intent(inout) :: ieap
integer, intent(in) :: iy
integer, intent(in) :: iyr
integer, intent(in) :: mo
integer, intent(in) :: nbyr
integer, intent(in) :: nd
integer :: ii, j, jea, jj, k, l, n, wet
real(dp) :: wabad
real(dp) :: pcp_tmp = 0. ! temporary precipitation values for output
real(dp) :: snow_tmp = 0. ! temporary snow values for output
real(dp) :: pcp_tot = 0. ! total preciptiation (rain + snow) for GIS output
! pointer on subbasin #### WATER MANAGEMENT MODULE ####
TYPE (TSubbasin), POINTER :: pS
!#### INITIALIZATION: call initsub
call subbasin_initialise_subbasin(canev, sml, xcklsp, xet, xnorg, xnorgp, xporg, xpsedp, xqd, yon, yph)
j = inum1
aff = da * 1000. * flu(j)
!**** COMPUTE TOTALS for precip:
precip = subp(j)
!###########################
!#### SNOW MODULE ####
!###########################
if (bSnowModule) then
call snow_initialise_subbasin
else
susb(1, j) = susb(1, j) + precip
sub(1) = sub(1) + precip * flu(j)
sbp(j) = sbp(j) + precip
!#### COMPUTE snow, recalc precip
! ATTN: Snow accumulation, melt, and snow evaporation
! are calculated at hydrotope level!
snowVal = 0.
if (tx(j) .le. tsnfall(j)) then
snowVal = precip
sub(2) = sub(2) + precip * flu(j)
call output_store_subbasin_value(snowfall_weq_output_id, j, precip)
susb(2, j) = susb(2, j) + precip
precip = 0.
end if
end if !( bSnowModule )
!###########################
!*********************************************************** START OF CYCLE 100
! MAIN CYCLE IN SUBBASIN: 100 - HYDROTOPE caclculation:
! 1) initial conditions,
! 2) call HYDROTOP,
! 3) compute weighted averages for subbasin.
do jea = 1, neap(j)
n = mstruc(j, jea, 1) ! land use type
k = mstruc(j, jea, 2) ! soil number
wet = mstruc(j, jea, 3) ! wetland 0 / 1
!###########################
!#### SNOW MODULE ####
!###########################
!#### CALL SNOM: COMPUTE snow & snowmelt, recalc precip
! ATTN: in snom() precip and sno() are recalculated
if (bSnowModule) then
call snow_process(j, jea, ida, mstruc, precip, tmn, tx)
else
if (snowVal > 0.) snoa(j, jea) = snoa(j, jea) + snowVal
if (tmx(j) > tmelt(j) .AND. snoa(j, jea) > 1E-4 ) call snow_degree_day_melting(j, jea)
end if
!###########################
!**** ESTABLISH number of soil layers: different for arable and non-arable soils
nn = ns(k)
if (landuse_is_cropland(n) ) nn = nsa(k)
!**** CALC hsumul(j, jea), hsumfc(j, jea), hwss(2, j, jea) and hrtc(j, jea, l)
call subbasin_assign_hsum(fc, hsumfc, hsumul, hwss, nn, snum, ul)
!#### CALL HRFLOWTT
call subbasin_flow_travel_time(j, jea, k, bSubcatch, bff, hrtc, nn, sc, sccor, z)
!**** INITIAL CONDITIONS for water & nutrients
if (iy == 1 .and. ida == 1) call subbasin_init_water_nutrients(n, nn, alai, almn, ano3, anora, anors, ap, cv, cva, dm, g, nucr, op, plab, pma, pmn, pms, porg, rd, sno, snoa, ste, stin, swe, swin, wmn, wn, wno3, wpo)
!**** INITIALIZATION of forest - in the beginning of every year
! ATTN! Forest is init. every year, as no forest management is included
if (ida .eq. 1) then
if (landuse_is_forest(n) ) then
dm(j, jea) = 200000.
rd(j, jea) = 2000.
end if
end if
!write(*,*) 'Inside subbasin.f95 ::: before olai(1, 1) = ', olai(1,1)
!#### CALL HYDROTOPE - calculate processes in HYDROTOPE
call hydrotope_process(j, jea, k, n, wet, bSubcatch, dart, daycounter, flu, frar, ida, iy, iyr, mo, nbyr, precip, sbar, tx)
!write(*,*) 'Inside subbasin.f95 ::: olai(1, 1) = ', olai(1,1)
! Subroutine wstress is called after ep was calculated in evap.f90
! This usually decreases ep but et is not updated thereafter.
! Important for hydrotope output!
et = ep + es
!**** COMPUTE WEIGHTED AVERAGES
! ATTN: xnorgp, xporgp - in kg/ha
! ATTN: xnorg, xporg - coef. (10./bd(k)) to tranf. kg/ha to g/t (I layer!)
! ATTN: xet - not from et, et is not recalc after cropmd
sumcn = sumcn + cn * frar(j, jea)
xqd = xqd + qd * frar(j, jea)
xqi = xqi + qi * frar(j, jea)
xssf = xssf + ssf * frar(j, jea)
xsep = xsep + sep * frar(j, jea)
xswind = xswind + swind * frar(j, jea)
xeo = xeo + eopot * frar(j, jea)
xet = xet + (ep + es + canev) * frar(j, jea)
xsnoev = xsnoev + snoev * frar(j, jea)
sno(j) = sno(j) + snoa(j, jea) * frar(j, jea)
! summarise total area covered by snow
if (snoa(j, jea) > 1E-4 ) then
area_tot_snow = area_tot_snow + frar(j, jea) * sbar(j)
depth_ave_snow = depth_ave_snow + snoa(j, jea) * frar(j, jea) * sbar(j) / (da * 10 ** 6)
call output_store_hydrotope_value(snow_depth_weq_output_id, j, jea, snoa(j, jea))
if (ida == nd) snow_acc_mm = snow_acc_mm + snoa(j, jea) * frar(j, jea) * sbar(j) / (da * 10 ** 6) - snow_acc_mm0
end if
if (ida == nd) soil_acc_mm = soil_acc_mm + swe(j, jea) * frar(j, jea) * sbar(j) / (da * 10 ** 6)
!###########################
!#### SNOW MODULE ####
!###########################
if (bSnowModule) then
if (gla(j, jea) > 1E-4 ) then
area_tot_glacier = area_tot_glacier + frar(j, jea) * sbar(j)
depth_ave_glacier = depth_ave_glacier + gla(j, jea) * frar(j, jea) * sbar(j) / (da * 10 ** 6)
if (ida == nd) glac_acc_mm = glac_acc_mm + gla(j, jea) * frar(j, jea) * sbar(j) / (da * 10 ** 6) - glac_acc_mm0
end if
xtmit = xtmit + tmit * frar(j, jea)
xtmax = xtmax + tmax * frar(j, jea)
xtmin = xtmin + tmin * frar(j, jea)
xprecip = xprecip + precipe * frar(j, jea) ! weighted subbasin rainfall or rainfall + snow melt (not total precipitation), mm
xprecip_elev_cor = xprecip_elev_cor + precip_elev_cor * frar(j, jea) ! weighted subbasin elevation - based corrected precipitation, mm
xsml = xsml + smle * frar(j, jea)
xsnow = xsnow + snowVal * frar(j, jea) ! weighted subbasin snow from snowfall, mm
xvsn = xvsn + vsn * frar(j, jea) ! weighted subbasin from snow pack, mm
end if
!###########################
xyno3 = xyno3 + yno3 * frar(j, jea)
xssfn = xssfn + ssfn * frar(j, jea)
xpercn = xpercn + percn * frar(j, jea)
xysp = xysp + ysp * frar(j, jea)
xcklsp = xcklsp + cklsp(j, jea) * frar(j, jea)
xnorg = xnorg + anors(j, jea, 1) * (10. / bd(k)) * frar(j, jea)
xnorgp = xnorgp + anors(j, jea, 1) * frar(j, jea)
xporg = xporg + porg(j, jea, 1) * (10. / bd(k)) * frar(j, jea)
xporgp = xporgp + porg(j, jea, 1) * frar(j, jea)
xpsed = xpsed + (porg(j, jea, 1) + pms(j, jea, 1) + pma(j, jea, 1)) * (10. / bd(k)) * frar(j, jea)
xpsedp = xpsedp + (porg(j, jea, 1) + pms(j, jea, 1) + pma(j, jea, 1)) * frar(j, jea)
!###########################
!#### SNOW MODULE ####
!###########################
if (bSnowModule) then
pcp_tmp = precip_elev_cor
pcp_tot = precipe + snowVal
else
pcp_tmp = precip
pcp_tot = subp(j)
end if
!###########################
call output_store_hydrotope_value(subsurface_runoff_output_id, j, jea, ssf)
call output_store_hydrotope_value(percolation_output_id, j, jea, sep)
call output_store_hydrotope_value(transpiration_output_id, j, jea, ep)
call output_store_hydrotope_value(soil_evaporation_output_id, j, jea, es)
call output_store_hydrotope_value(soil_water_index_output_id, j, jea, swind)
call output_store_hydrotope_value(soil_water_content_output_id, j, jea, swe(j, jea))
!**** CALC variables for GIS output
if (bSnowModule) then
call output_store_hydrotope_value(precipitation_output_id, j, jea, precipe + snowVal)
else
call output_store_hydrotope_value(precipitation_output_id, j, jea, subp(j))
end if
call output_store_hydrotope_value(surface_runoff_output_id, j, jea, qd)
call output_store_hydrotope_value(groundwater_recharge_output_id, j, jea, sep)
call output_store_hydrotope_value(eta_output_id, j, jea, ep + es + canev)
call output_store_hydrotope_value(etp_output_id, j, jea, eo)
!#### CALL CROP_GIS
if (iy > 1 .AND. ida == ndpri) call crop_yield_output(j, jea, k, ieap, ms, ndgro, tsav, wsav, ylda, landuse_is_cropland(n))
!**** CALC parsz() for subbasins
do jj = 1, 5
parsz(jj, j) = parsz(jj, j) + psz(jj, k) * frar(j, jea)
end do
end do
!*********************************************************** END OF CYCLE 100
pet_day(j) = xeo ! for transmission losses
!###########################
!#### SNOW MODULE ####
!###########################
if (bSnowModule) then
!#### recalculate precipitation and air temperature for the whole basin
precip = xprecip
tx(j) = xtmit
tmx(j) = xtmax
tmn(j) = xtmin
sml = xsml
psnow = xsnow
xprecip = 0.
xtmit = 0.
xtmax = 0.
xtmin = 0.
xsml = 0.
xsnow = 0.
end if
!###########################
!#### CALL alpha to calc r1 - alpha to be used in peakq
if (precip > 0.) then
call soil_curve_number_alpha(j, mo, precip, sml)
endif
!**** COMPUTE PEAK RUNOFF RATE, TRANSMISSION LOSSES & EROSION for subbasin
!#### CALL PEAKQ, TRAN, YSED
! xssf (mm), css (ppm) aff=1000*da*flu(j), yd (t)
if (xqd .ne. 0.) then
call soil_curve_number_peak_runoff(j)
if (xqd .gt. 0. .and. pr .gt. 0.) then
iv = 1
vo = xqd * da * flu(j) * 1000.
dur = vo / (pr * 3600.)
if (dur .gt. 24.) dur = 24.
! SL BEGIN: transmission losses are calculated in the routing process using subroutine tran_river
! SL only recalculate peak runoff rate (pr) in function tran(j)
q1 = xqd ! SL
call soil_curve_transmission_losses(j, chk, chl, chw)
!qtl = q1 - xqd
xqd = q1
! if ( qtl.lt.0. ) then
! xqd = q1
! qtl = 0.
! end if
! revapst(j) = revapst(j) + qtl
! SL END
end if
call erosion_soil_loss(aff, pr, xqd, yd)
end if
yd = yd + xssf * aff * css(j)
!**** COMPUTE ENRICHMENT RATIO, ORG N and P LOSS WITH EROSION
!#### CALL ENRSB, ORGNSED, PSED
if (qd > 0. .AND. pr > 0. .AND. precip > 0.) then
call erosion_enritchment_ratio(j, da, da9, flu, pr, precip, rp, xqd, yd)
call erosion_organic_nitrate_loss(j, da9, yd, yon)
call erosion_phosphorus_loss(j, da9, yd, yph)
end if
!**** COMPUTE GROUND WATER CONTRIBUTION TO FLOW & WATER YIELD
!#### CALL GWMOD
! NOTE that revap from shallow aquifer is add to xet
! call gwmod(j)
call groundwater_process(j, xsep, xet, xeo)
wysb = xqi + xssf + gwq(j) ! - qtl ! water yield, transmission losses are calculated in the routing process
runoff_mm(j) = wysb
call output_store_subbasin_value(river_runoff_output_id, j, wysb)
!**** COMPUTE RETENTION of NUTRIENTS in SUBBASIN
! Retention of N (xyno3, xssfn, xpercn) and soluble P (xysp)
! Method: from Fred Hattermann 2003
! Updated 11.03.2008 from Shaochun's file
xyno3 = (1. - exp(- 1. / retNsur - degNsur)) * xyno3 / (1. + retNsur * degNsur) + tmpNsur(j) * exp(- 1. / retNsur - degNsur)
xssfn = (1. - exp(- 1. / retNsub - degNsub)) * xssfn / (1. + retNsub * degNsub) + tmpNsub(j) * exp(- 1. / retNsub - degNsub)
xpercn = (1. - exp(- 1. / retNgrw - degNgrw)) * xpercn / (1. + retNgrw * degNgrw) + tmpNgrw(j) * exp(- 1. / retNgrw - degNgrw)
xysp = (1. - exp(- 1. / retPsur - degPsur)) * xysp / (1. + retPsur * degPsur) + tmpPsur(j) * exp(- 1. / retPsur - degPsur)
tmpNsur(j) = xyno3
tmpNsub(j) = xssfn
tmpNgrw(j) = xpercn
tmpPsur(j) = xysp
!**** COMPUTE TOTALS for water: sub(1...13)
!###########################
!#### SNOW MODULE ####
!###########################
if (bSnowModule) then
susb(1, j) = susb(1, j) + xprecip_elev_cor ! precip
sub(1) = sub(1) + xprecip_elev_cor * flu(j) ! precip * flu(j)
sbp(j) = sbp(j) + xprecip_elev_cor ! precip !precipe
sub(2) = sub(2) + psnow * flu(j)
call output_store_subbasin_value(snowfall_weq_output_id, j, psnow)
susb(2, j) = susb(2, j) + psnow
end if
!###########################
sub(3) = sub(3) + sml * flu(j)
susb(3, j) = susb(3, j) + sml
sub(4) = sub(4) + tmx(j) * flu(j)
call output_store_subbasin_value(tmax_output_id, j, tmx(j))
susb(4, j) = susb(4, j) + tmx(j)
sub(5) = sub(5) + tmn(j) * flu(j)
call output_store_subbasin_value(tmin_output_id, j, tmn(j))
susb(5, j) = susb(5, j) + tmn(j)
sub(29) = sub(29) + tx(j) * flu(j) ! / get_days_in_month(mo, iyr)
call output_store_subbasin_value(tmean_output_id, j, tx(j))
susb(29, j) = susb(29, j) + tx(j) ! / get_days_in_month(mo, iyr)
sub(6) = sub(6) + ra(j) * flu(j)
susb(6, j) = susb(6, j) + ra(j)
sub(7) = sub(7) + sumcn * flu(j)
susb(8, j) = susb(8, j) + xqd
sub(8) = sub(8) + xqd * flu(j)
smq(j) = smq(j) + xqd
susb(9, j) = susb(9, j) + xssf
sub(9) = sub(9) + xssf * flu(j)
smsq(j) = smsq(j) + xssf
! SL Transmission losses are calculated during the routing process...
! susb(10, j) = susb(10, j) + qtl
! sub(10) = sub(10) + qtl * flu(j)
! revapst(j) = revapst(j) + qtl
susb(11, j) = susb(11, j) + xsep
sub(11) = sub(11) + xsep * flu(j)
susb(12, j) = susb(12, j) + xeo
sub(12) = sub(12) + xeo * flu(j)
susb(13, j) = susb(13, j) + xet
sub(13) = sub(13) + xet * flu(j)
!**** COMPUTE TOTALS for ground water: sub(15...19)
susb(15, j) = susb(15, j) + gwq(j)
sub(15) = sub(15) + gwq(j) * flu(j)
call output_store_subbasin_value(baseflow_output_id, j, gwq(j))
susb(16, j) = susb(16, j) + revap
sub(16) = sub(16) + revap * flu(j)
susb(17, j) = susb(17, j) + gwseep
sub(17) = sub(17) + gwseep * flu(j)
susb(18, j) = susb(18, j) + gwchrg
sub(18) = sub(18) + gwchrg * flu(j)
susb(19, j) = susb(19, j) + xswind
sub(19) = sub(19) + xswind * flu(j)
xxswind = xxswind + xswind * flu(j)
!**** Calc daily water balance
wabad = sub(1) - sub(8) - sub(9) - sub(11) - sub(13)
!**** COMPUTE TOTALS for water yield (wysb) & sediments: sub(20...25)
susb(20, j) = susb(20, j) + wysb
sub(20) = sub(20) + wysb * flu(j)
call output_store_subbasin_value(water_yield_output_id, j, wysb)
sub(21) = sub(21) + yd / (100. * da * flu(j))
sym(j) = sym(j) + yd
susb(21, j) = susb(21, j) + yd / (100. * da * flu(j))
sub(22) = sub(22) + yon * flu(j)
susb(22, j) = susb(22, j) + yon
sub(23) = sub(23) + yph * flu(j)
susb(23, j) = susb(23, j) + yph
sub(24) = sub(24) + ysp * flu(j)
susb(24, j) = susb(24, j) + ysp
!**** COMPUTE SUBBASIN OUTPUTS FOR ROUTING: sda(), varoute()
! ATTN: sda(6, j) = sum of 3 fluxes after retention, new version
! ATTN: coef (dart()*1000) to transform from mm to m3 (42, 48)
! ATTN: coef (dart()*100) to transform kg/ha to kg (44-47)
! ATTN: wysb(mm)*dart(km2)*1000 = wysb*dart*1000 (m3)
! ATTN: sda(2, j), sda(8, j) in m3
! ATTN: xyno3 in kg/ha, sda(6) & varoute(6, .) in kg
sda(1, j) = precip
sda(2, j) = xqi * dart(ihout) * 1000.
sda(3, j) = yd
sda(4, j) = yon * dart(ihout) * 100
sda(5, j) = yph * dart(ihout) * 100
sda(6, j) = (xyno3 + xssfn + xpercn) * dart(ihout) * 100
sda(7, j) = xysp * dart(ihout) * 100
sda(8, j) = (wysb - xqi) * dart(ihout) * 1000.
xwysb = xwysb + wysb * dart(ihout) * 1000.
!cc riparian zone take water form GW flow as demanded
gwqLastday(j) = sda(8, j) * 0.9
sda(8, j) = sda(8, j) - additionalGwUptake(j)
!#################################
!#### WATER MANAGEMENT MODULE ####
!#################################
if (bWAM_Module) then
! check if water transfers take place in current subbasin
if (management_is_transfer_subbasin(j) ) then
pS => management_subbasin_pointer(j)
!-------------------------------------------------------------
! Add supply from water user(s) to subbasins' hydrograph storage location (=j)
! If the subbasin is a reservoir, this step will be overwritten by the reservoir module
!-------------------------------------------------------------
! if external inputs from water user(s) are greater than 0
! add supply of previous day to subbasins' hydrograph storage location (surface runoff component)
if (daycounter > 1) then
if (pS % inflow(daycounter - 1) > 0. ) &
sda(2, j) = sda(2, j) + pS % inflow(daycounter - 1) * 86400. ! convert input from m3 / s to m3 / d
end if
!-------------------------------------------------------------
!-------------------------------------------------------------
! Withdraw water from subbasin outlet
! but only if subbasin is a headwater .and. not a reservoir.
! otherwise withdawals are taken during the routing process in subroutine 'add'
! Moreover, the water users %supplied and %tr_losses are calculated.
!-------------------------------------------------------------
! Check if current subbasin is a headwater...in this case,
! bWAMHydrograph is .true. at the array position of the subbasin number
! ATTENTION: The values of sda(2, j) and sda(8, j) may be changed!
if (bWAMHydrograph(j) ) then
call management_transfer_out(j, sda(2, j), sda(8, j), daycounter, ida, iyr)
end if
!-------------------------------------------------------------
pS % Q_act(daycounter) = real(sda(2, j) / 86400. + sda(8, j) / 86400., 4)
end if
end if
!#################################
do ii = 1, 8
varoute(ii, ihout) = sda(ii, j)
end do
!###########################
!#### SNOW MODULE ####
!###########################
if (bSnowModule) then
pcp_tmp = xprecip_elev_cor
snow_tmp = psnow
else
pcp_tmp = precip
snow_tmp = snowVal
end if
!###########################
return
!******************************************************************************
contains
! Subroutines inside subroutine subbasis
subroutine subbasin_assign_hsum(fc, hsumfc, hsumul, hwss, nn, snum, ul)
real(dp), dimension(:, :), intent(in) :: fc
real(dp), dimension(:, :), intent(inout) :: hsumfc
real(dp), dimension(:, :), intent(inout) :: hsumul
real(dp), dimension(:, :, :), intent(inout) :: hwss
integer, intent(in) :: nn
integer, dimension(:), intent(in) :: snum
real(dp), dimension(:, :), intent(in) :: ul
! Correction from Fred Hattermann
! Reason: nn is defined for HRU, & sumul(k) and sumfc(k) are not appropriate
hsumul(j, jea) = 0.
hsumfc(j, jea) = 0.
do l = 1, nn
hsumul(j, jea) = hsumul(j, jea) + ul(l, k)
hsumfc(j, jea) = hsumfc(j, jea) + fc(l, k)
end do
if (hsumfc(j, jea) .lt. 0.001) then
call log_error("subbasin_assign_hsum", "hsumfc less than 0.001, is this soil in soil.cio? (subbasin, hru, landuse, soil)", &
ints=(/j, jea, n, snum(k)/))
end if
hwss(1, j, jea) = - 1.2677 +log(hsumfc(j, jea))
hwss(2, j, jea) = 1.6234 / hsumfc(j, jea)
end subroutine subbasin_assign_hsum
!******************************************************************************
!******************************************************************************
subroutine subbasin_init_water_nutrients(n, nn, alai, almn, ano3, anora, anors, ap, cv, cva, dm, g, nucr, op, plab, pma, pmn, pms, porg, rd, sno, snoa, ste, stin, swe, swin, wmn, wn, wno3, wpo)
use landuse, only : LULC, landuse_index, landuse_is_cropland, landuse_is_natural_vegetation
real(dp), dimension(:, :), intent(inout) :: alai
real(dp), dimension(:), intent(in) :: almn
real(dp), dimension(:, :, :), intent(inout) :: ano3
real(dp), dimension(:, :, :), intent(inout) :: anora
real(dp), dimension(:, :, :), intent(inout) :: anors
real(dp), dimension(:, :), intent(in) :: ap
real(dp), dimension(:), intent(in) :: cv
real(dp), dimension(:, :), intent(inout) :: cva
real(dp), dimension(:, :), intent(inout) :: dm
real(dp), dimension(:, :), intent(inout) :: g
integer, dimension(:, :), intent(inout) :: nucr
real(dp), dimension(:, :), intent(in) :: op
real(dp), dimension(:, :, :), intent(inout) :: plab
real(dp), dimension(:, :, :), intent(inout) :: pma
real(dp), dimension(:, :), intent(in) :: pmn
real(dp), dimension(:, :, :), intent(inout) :: pms
real(dp), dimension(:, :, :), intent(inout) :: porg
real(dp), dimension(:, :), intent(inout) :: rd
real(dp), dimension(:), intent(in) :: sno
real(dp), dimension(:, :), intent(inout) :: snoa
real(dp), dimension(:, :, :), intent(inout) :: ste
real(dp), dimension(:, :), intent(in) :: stin
real(dp), dimension(:, :), intent(inout) :: swe
real(dp), dimension(:), intent(in) :: swin
real(dp), dimension(:, :), intent(in) :: wmn
real(dp), dimension(:, :), intent(in) :: wn
real(dp), dimension(:, :), intent(in) :: wno3
real(dp), dimension(:, :), intent(in) :: wpo
! assign initial conditions for water & nutrients at the first
! day of simulation
! land use type
integer, intent(in) :: n
! number of soil arable layers
integer, intent(in) :: nn
integer :: l = 0
do l = 1, nn
ste(j, jea, l) = stin(l, k)
end do
swe(j, jea) = swin(k)
!###########################
!#### SNOW MODULE ####
!###########################
if (.NOT. bSnowModule) then
snoa(j, jea) = sno(j)
end if
!###########################
!**** Initial conditions for nutrients
do l = 1, nn
ano3(j, jea, l) = wno3(l, k)
anora(j, jea, l) = wmn(l, k)
anors(j, jea, l) = wn(l, k)
plab(j, jea, l) = ap(l, k)
porg(j, jea, l) = wpo(l, k)
pma(j, jea, l) = pmn(l, k)
pms(j, jea, l) = op(l, k)
end do
!**** Correction of initial values ano3(), plab() for noncropland
if (landuse_is_cropland(n) ) then
do l = 1, nn
ano3(j, jea, l) = ano3(j, jea, l) * .1
plab(j, jea, l) = plab(j, jea, l) * .1
end do
end if
!**** INITIALIZATION of crop parameters
g(j, jea) = 0.
!alai(j, jea)=alai0(get_lu_index(n))
! set LAI to minimum LAI as indicated in crop.dat
if (LULC % veg_code(landuse_index(n)) > 0 ) then
alai(j, jea) = almn(LULC % veg_code(landuse_index(n)))
else
alai(j, jea) = 0.
end if
!**** Correction of initial values ano3(), plab() for noncropland
if (landuse_is_cropland(n) ) then
cva(j, jea) = cv(j)
else
cva(j, jea) = 0.
nucr(j, jea) = 0 ! crop number (crop.dat)
dm(j, jea) = 0. ! total biomass
rd(j, jea) = 0. ! root depth
end if
if (landuse_is_natural_vegetation(n) ) then
dm(j, jea) = 3.
rd(j, jea) = 300.
end if
end subroutine subbasin_init_water_nutrients
!******************************************************************************
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
end subroutine runsubbasin