subroutine river_route(ihout, inum1, inum2, chc, chxk, conn, cpp, da9, dart, er, flu, ida, iy, iyr, revapst, runs, sbar, sub, susb, xysp, yd, yon, yph)
!**** PURPOSE: THIS SUBROUTINE CONTROLS THE CHANNEL ROUTING
! STEPS: 1) xxqd = varoute(2, )
! 2) xxqd recalc in rthyd
! 3) varoute(2, )=xxqd
!**** CALLED IN: MAIN
! ATTN! Look at the final line in the .fig file to define where
! the final basin output is written: in route() or in add()
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! PARAMETERS & VARIABLES
!
! >>>>> PARAMETERS IN TITLE
! icode = code to switch between routing subroutines (here: 2)
! ihout = hydrological storage location
! inum1 = reach number
! inum2 = inflow hydrograph(inum2 hydrograph is routed through inum1)
! >>>>>
! >>>>> COMMON PARAMETERS & VARIABLES
! accf(1) = accumulated runoff, m3/sec.
! accf(2) = accumulated sediment, t
! accf(3) = accumulated NO3-N, kg
! da = drainage area, km2
! dart(ih) = drainage area for subbasin, km2
! diver = diversion (not active)
! evp = evaporation from river surface, m3
! flu(j) = fraction of subbasin area in the basin
! ida = current day
! iy = current year as counter (1, ..., nbyr)
! iyr = current year
! revapst(j) = transmission losses in the reach, mm
! rflow = return flow (not acive)
! rl = seepage (not acive)
! runs(ida) = runoff simulated for the basin, m3/sec.
! sdti = inflow + storage, calc in rthyd, m3/sec.
! srch(18, j) = monthly reach outputs:
! srch(1, j) = water inflow, m3/s
! (2, j) = water outflow, m3/s
! (3, j) = sediment in, t
! (4, j) = sediment out, t
! (5, j) = sediment conc, mg/l
! (6, j) = org N in, kg
! (7, j) = org N out, kg
! (8, j) = sed P in, kg
! (9, j) = sed P out, kg
! (10, j) = evap, m3/s
! (11, j) = transmission losses, m3/s
! (12, j) = seepage, m3/s
! (13, j) = diversions, m3/s
! (14, j) = return flow, m3/s
! (15, j) = NO3-N in, kg
! (16, j) = NO3-N out, kg
! (17, j) = sol P in, kg
! (18, j) = sol P out, kg
! tlc = transmission losses, m3
! varoute(1:8, ih) = variables 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
! xxnit = NO3-N amount in the reach, kg
! xxqd = water amount in the reach, m3
! xxssf = subsurface flow, m3
! xysp = soluble P, kg
! yd = daily soil loss from subbasin caused by water erosion, t
! ydi = sediment, t
! yon = org N loss with erosion, kg
! yph = P org. loss with erosion, kg
! >>>>>
! >>>>> STATIC PARAMETERS
! cnit = local par
! ivar2neg = local variable: negative surface runoff
! ivar8neg = local variable: negative subsurface runoff
! j = local par, reach no.
! sedcon = local variable: sediment concentration, t/m3
! xx = local par
! >>>>>
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
real(dp), dimension(:), intent(in) :: chc
real(dp), dimension(:), intent(in) :: chxk
real(dp), intent(in) :: conn
real(dp), intent(in) :: cpp
real(dp), intent(in) :: da9
real(dp), dimension(:), intent(in) :: dart
real(dp), intent(inout) :: er
real(dp), dimension(:), intent(in) :: flu
integer, intent(in) :: ida
integer, intent(in) :: iy
integer, intent(in) :: iyr
real(dp), dimension(:), intent(inout) :: revapst
real(dp), dimension(366), intent(inout) :: runs
real(dp), dimension(:), intent(in) :: sbar
real(dp), dimension(30), intent(inout) :: sub
real(dp), dimension(:, :), intent(inout) :: susb
real(dp), intent(inout) :: xysp
real(dp), intent(inout) :: yd
real(dp), intent(inout) :: yon
real(dp), intent(inout) :: yph
integer, intent(in) :: ihout, inum1, inum2
integer :: ivar2neg, ivar8neg, j
real(dp) :: cnit, sedcon, xx
! transmission losses, mm
real(dp) :: tl_mm
! evaporative losses from river surface, mm
real(dp) :: evp_mm
!**** INITIALIZATION
j = inum1
tl_mm = 0.
evp_mm = 0.
tlc = 0.
evp = 0.
xxqd = varoute(2, inum2)
yd = varoute(3, inum2)
yon = varoute(4, inum2)
yph = varoute(5, inum2)
xxnit = varoute(6, inum2)
xysp = varoute(7, inum2)
xxssf = varoute(8, inum2)
if (ida .eq. 1 .and. j .eq. 1) ivar2neg = 0
if (ida .eq. 1 .and. j .eq. 1) ivar8neg = 0
if (ida .eq. 1 .and. j .eq. 1) then
accf(1) = 0.
accf(2) = 0.
accf(3) = 0.
end if
!#### CALL RTHYD to route hydrograph to basin outlet: recalc xxqd, xxssfn
call river_muskingum_routing(j, ihout, inum1, ida, iy, iyr)
!**** COMPUTE transmission losses in the reach
if (tlrch > 0. .OR. evrch > 0.) then
call river_transmission_loss (j)
tl_mm = tlc / sbar(j) * 1000. ! convert m3 to mm
evp_mm = evp / sbar(j) * 1000. ! convert m3 to mm
select case (tlgw)
case (0) ! add transmission losses to shallow ground water storage
revapst(j) = revapst(j) + tl_mm ! add to revap storage
susb(18, j) = susb(18, j) + tl_mm ! add to gw recharge
sub(18) = sub(18) + tl_mm * flu(j)
case (1) ! add transmission losses to deep ground water storage (lost from the system)
susb(17, j) = susb(17, j) + tl_mm ! add to deep gw (gwseep)
sub(17) = sub(17) + tl_mm * flu(j)
! add transmission losses to transmission losses variables
susb(10, inum1) = susb(10, inum1) + tl_mm
sub(10) = sub(10) + tl_mm * flu(j)
case (2) ! add transmission losses to both shallow and deep ground water storage to equal parts
revapst(j) = revapst(j) + tl_mm * .5 ! add to revap storage
susb(18, j) = susb(18, j) + tl_mm * .5 ! add to gw recharge
sub(18) = sub(18) + tl_mm * .5 * flu(j)
susb(17, j) = susb(17, j) + tl_mm ! add to deep gw (gwseep)
sub(17) = sub(17) + tl_mm * flu(j)
! add transmission losses to transmission losses variables
susb(10, inum1) = susb(10, inum1) + tl_mm * .5
sub(10) = sub(10) + tl_mm * .5 * flu(j)
end select
! add evaporative losses from river surface to ETact
susb(13, j) = susb(13, j) + evp_mm
sub(13) = sub(13) + evp_mm * flu(j)
end if ! ( tlrch > 0. .OR. evrch > 0. )
!#### CALL RTSED, ENRRT, RTORGN, RTPSED to route sediments, org N & org P
! rtsed - recalc yd
! enrrt - calc er - enrichment coefficient
! rtorgn - recalc yon
! rtpsed - recalc yph
if (sdti .gt. 0.) then
ydi = yd
call river_route_erosion(j, inum1, chc, chxk, yd)
call river_erosion_enritchment_ratio(inum1, da9, er, yd)
if (varoute(3, inum2) .gt. 0.) then
xx = yd / varoute(3, inum2)
else
xx = 0.
end if
call river_route_nitrate(conn, da9, er, yon)
yon = yon * xx
if (yon .gt. varoute(4, inum2)) yon = varoute(4, inum2)
call river_route_phosphorus(cpp, da9, yph)
yph = yph * xx
if (yph .gt. varoute(5, inum2)) yph = varoute(5, inum2)
end if
!**** Calculate monthly reach outputs
! output 5 is changed from deposition in tons to concentration in ppm
srch(1, inum1) = srch(1, inum1) + varoute(2, inum2) / 86400.
srch(2, inum1) = srch(2, inum1) + xxqd / 86400.
srch(3, inum1) = srch(3, inum1) + varoute(3, inum2)
srch(4, inum1) = srch(4, inum1) + yd
if (xxqd .gt. 0.1) then
sedcon = yd / xxqd * 1.e6
else
sedcon = 0.
end if
if (sedcon .gt. 200000.) sedcon = 200000.
srch(5, inum1) = srch(5, inum1) + sedcon
srch(6, inum1) = srch(6, inum1) + varoute(4, inum2)
srch(7, inum1) = srch(7, inum1) + yon * dart(ihout) * 100.
srch(8, inum1) = srch(8, inum1) + varoute(5, inum2)
srch(9, inum1) = srch(9, inum1) + yph * dart(ihout) * 100.
srch(10, inum1) = srch(10, inum1) + evp / 86400.
srch(11, inum1) = srch(11, inum1) + tlc / 86400.
srch(12, inum1) = srch(12, inum1) + rl / 86400.
srch(13, inum1) = srch(13, inum1) + diver / 86400.
srch(14, inum1) = srch(14, inum1) + rflow / 86400.
srch(15, inum1) = srch(15, inum1) + varoute(6, inum2)
srch(16, inum1) = srch(16, inum1) + xxnit
srch(17, inum1) = srch(17, inum1) + varoute(7, inum2) * 100.
srch(18, inum1) = srch(18, inum1) + xysp * dart(ihout) * 100.
! sl begin
srch(19, inum1) = srch(19, inum1) + (varoute(2, inum2) + varoute(8, inum2)) / 86400.
! sl end
!**** RECALCULATE varoute in Hydrological Storage Location ihout
varoute(2, ihout) = xxqd
varoute(3, ihout) = yd
varoute(4, ihout) = yon
varoute(5, ihout) = yph
varoute(6, ihout) = xxnit
varoute(7, ihout) = xysp
varoute(8, ihout) = xxssf
!**** WRITE DAILY OUTPUT in outlet (if j=1) to river output file unit=75
! WRITE is possible in variants - needed one should be opened
if (j .eq. 1) then
runs(ida) = (varoute(2, ihout) + varoute(8, ihout)) / 86400.
accf(1) = accf(1) + runs(ida)
accf(2) = accf(2) + varoute(3, ihout)
accf(3) = accf(3) + varoute(6, ihout)
!**** CALC N CONCENTRATION
cnit = varoute(6, ihout) / runs(ida) / 86.4
if (varoute(2, ihout) .lt. 0.) ivar2neg = ivar2neg + 1
if (varoute(8, ihout) .lt. 0.) ivar8neg = ivar8neg + 1
end if
return
end subroutine river_route