!     FILE route.f
!
!     SUBROUTINES IN THIS FILE 	 CALLED FROM
!     subroutine route(icode, ihout, inum1, inum2) main
!     subroutine add(icode, ihout, inum1, inum2) main
!     subroutine transfer(icode, ihout, inum1, inum2) (main, not active)
module river

  use utilities, only : dp, path_max_length, log_debug, log_warn

  implicit none

  ! coef. for routing: correction coef to calculate xkm
  real(dp), save, dimension(:), allocatable :: bsn_roc2
  ! coef. for routing: correction coef to calculate xkm
  real(dp), save, dimension(:), allocatable :: bsn_roc4
  ! coefficient to correct channel width (from .bsn)
  real(dp), save :: chwc0 = 1.
  ! storage correction coef (from .bsn)
  real(dp), save :: storc1 = 0.5
  ! correction coef. for channel N (Manning) factor
  real(dp), save :: chnnc0 = 1.
  ! main channel slope, m / m
  real(dp), save, dimension(:), allocatable :: chs
  ! channel N value
  real(dp), save, dimension(:), allocatable :: chn
  ! overland flow N value
  real(dp), save, dimension(:), allocatable :: ovn
  ! return flow travel time, days. If 0, then estimated
  real(dp), save, dimension(:), allocatable :: rt
  ! daily routed and added discharge per subbasin
  real(dp), save, dimension(:), allocatable :: runsub_m3s
  real(dp), save, allocatable, dimension(:) :: xkm_qd
  real(dp), save, allocatable, dimension(:) :: c1_qd
  real(dp), save, allocatable, dimension(:) :: c2_qd
  real(dp), save, allocatable, dimension(:) :: c3_qd
  real(dp), save, allocatable, dimension(:) :: c4_qd
  real(dp), save, allocatable, dimension(:) :: xkm_ssf
  real(dp), save, allocatable, dimension(:) :: c1_ssf
  real(dp), save, allocatable, dimension(:) :: c2_ssf
  real(dp), save, allocatable, dimension(:) :: c3_ssf
  real(dp), save, allocatable, dimension(:) :: c4_ssf
  real(dp), save :: tlrch = 1.
  real(dp), save :: evrch = 1.
  real(dp), save :: tlc
  real(dp), save :: evp
  ! 0 = add transmission losses to shallow ground water; 1 = add to deep ground water
  integer, save :: tlgw = 0
  real(dp), save, dimension(:), allocatable :: pet_day
  !
  real(dp), save :: xxqd
  !
  real(dp), save :: sdti
  !
  real(dp), save :: ydi
  !
  real(dp), save :: diver
  !
  real(dp), save :: rflow
  !
  real(dp), save :: rl
  ! coef. to estimate peak runoff in stream
  real(dp), save :: prf = 1.
  !
  real(dp), save :: spcon
  ! exponent for estimation of sediment transport
  real(dp), save :: spexp = 1.
  !
  real(dp), save :: xxssf
  !
  real(dp), save :: xxnit
  !
  real(dp), save :: roc1
  !
  real(dp), save :: roc3
  !
  real(dp), save, dimension(10) :: accf = 0.
  ! coefs for each subbasin
  real(dp), save, dimension(:), allocatable :: roc2
  real(dp) :: roc2_0 = 9.
  ! coefs for each subbasin
  real(dp), save, dimension(:), allocatable :: roc4
  real(dp) :: roc4_0 = 9.
  ! main channel length, km
  real(dp), save, dimension(:, :), allocatable :: chl
  ! average width of main channel, m
  real(dp), save, dimension(:, :), allocatable :: chw
  ! effective hydraulic conductivity of main channel, mm / h
  real(dp), save, dimension(:, :), allocatable :: chk
  ! channel and flow parameters
  real(dp), save, dimension(:, :), allocatable :: phi
  ! initial water storage in subbasins, m3
  real(dp), save, dimension(:), allocatable :: sdtsav
  ! channel depth, m
  real(dp), save, dimension(:), allocatable :: chd
  ! channel slope, m / m
  real(dp), save, dimension(:), allocatable :: chss
  ! channel N value
  real(dp), save, dimension(:), allocatable :: chnn
  ! monthly REACH outputs (dif. components)
  real(dp), save, dimension(:, :), allocatable :: srch
  ! = qdinp(, ) in the last day of the year, m3
  real(dp), save, dimension(:), allocatable :: qdilast
  ! = qdout(, ) in the last day pof the year, m3
  real(dp), save, dimension(:), allocatable :: qdolast
  ! = qssinp(, ) in the last day of the year, m3
  real(dp), save, dimension(:), allocatable :: qsilast
  ! = qssout(, ) in the last day of the year, m3
  real(dp), save, dimension(:), allocatable :: qsolast
  ! varoute(1:19, ih) = variables for routing
  real(dp), save, dimension(:, :), allocatable :: varoute
  ! surface flow - daily input in reaches, m3
  real(dp), save, dimension(:, :), allocatable :: qdinp
  ! surface flow - daily output in reaches, m3
  real(dp), save, dimension(:, :), allocatable :: qdout
  ! subsurface + g - w flow - daily input in reaches, m3
  real(dp), save, dimension(:, :), allocatable :: qssinp
  ! subsurface + g - w flow - daily output in reaches, m3
  real(dp), save, dimension(:, :), allocatable :: qssout

  ! Output variable ids
  integer :: &
    river_discharge_output_id = 0

  namelist / RIVER_PARAMETERS / &
    roc2_0, &
    roc4_0, &
    chnnc0, &
    chwc0, &
    prf, &
    roc1, &
    roc3, &
    spcon, &
    spexp, &
    storc1, &
    tlgw, &
    tlrch, &
    evrch

contains

  subroutine river_initialise(mb, mch, mhyd, subbasin_input_file_id)
    use input, only : get_config_fid
    use output, only : output_register_subbasin_var
    integer, intent(in) :: mb, mch, mhyd, subbasin_input_file_id

    read(get_config_fid(), RIVER_PARAMETERS)

    river_discharge_output_id = output_register_subbasin_var("discharge")

    call river_allocate(mb, mch, mhyd)
    call river_read_input(subbasin_input_file_id)

    ! In the SWIM manual it is suggested to set roc1 and roc3 values to 0.
    roc1 = 0.
    roc3 = 0.

  end subroutine river_initialise

  subroutine river_read_input(subbasin_input_file_id)
    use input, only : read_real_column
    integer, intent(in) :: subbasin_input_file_id

    call read_real_column(subbasin_input_file_id, "chl", chl(1, :), range=(/0., 1e4/), closed="n")
    chl(2, :) = chl(1, :)
    call read_real_column(subbasin_input_file_id, "chs", chs, range=(/0., 3./), closed="u")
    chss = chs
    call read_real_column(subbasin_input_file_id, "chw", chw(1, :), range=(/0., 1e4/), closed="n")
    chw(2, :) = chw(1, :)
    chw(2, :) = chw(2, :) * chwc0
    call read_real_column(subbasin_input_file_id, "chd", chd, range=(/0., 100./), closed="n")
    call read_real_column(subbasin_input_file_id, "chk", chk(1, :), 0.37_dp)
    call read_real_column(subbasin_input_file_id, "chk", chk(2, :), 0.37_dp)
    call read_real_column(subbasin_input_file_id, "chn", chn, 0.075_dp)
    call read_real_column(subbasin_input_file_id, "ovn", ovn, 0.15_dp)
    call read_real_column(subbasin_input_file_id, "rt", rt, 0.0_dp)
    call read_real_column(subbasin_input_file_id, "sdtsav", sdtsav, 0.0_dp)

  end subroutine river_read_input


  subroutine river_allocate(mb, mch, mhyd)
    integer, intent(in) :: mb, mch, mhyd
    allocate(c1_qd(mb))
    allocate(c1_ssf(mb))
    allocate(c2_qd(mb))
    allocate(c2_ssf(mb))
    allocate(c3_qd(mb))
    allocate(c3_ssf(mb))
    allocate(c4_qd(mb))
    allocate(c4_ssf(mb))
    allocate(chd(mch))
    allocate(chk(2, mb))
    allocate(chl(2, mb))
    allocate(chn(mb))
    allocate(chnn(mch))
    allocate(chs(mb))
    allocate(chss(mch))
    allocate(chw(2, mb))
    allocate(ovn(mb))
    allocate(pet_day(mb)) ! for transmission losses
    allocate(phi(20, mb))
    allocate(qdilast(mhyd))
    allocate(qdinp(mhyd, 366))
    allocate(qdolast(mhyd))
    allocate(qdout(mhyd, 366))
    allocate(qsilast(mhyd))
    allocate(qsolast(mhyd))
    allocate(qssinp(mhyd, 366))
    allocate(qssout(mhyd, 366))
    allocate(roc2(mb))
    allocate(roc4(mb))
    allocate(rt(mb))
    allocate(sdtsav(mb))
    allocate(srch(20, mch))
    allocate(varoute(19, mhyd))
    allocate(xkm_qd(mb))
    allocate(xkm_ssf(mb))
    allocate(runsub_m3s(mb))

    runsub_m3s = 0.
    chd = 0.
    chk = 0.37
    chl = 0.
    chn = 0.075
    chnn = 0.075 * chnnc0
    chs = 0.
    chss = 0.
    chw = 0.
    ovn = 0.15
    pet_day = 0.
    phi = 0.
    qdilast = 0.
    qdinp = 0.
    qdolast = 0.
    qdout = 0.
    qsilast = 0.
    qsolast = 0.
    qssinp = 0.
    qssout = 0.
    roc2 = roc2_0
    roc4 = roc4_0
    rt = 0.
    sdtsav = 0.
    srch = 0.
    varoute = 0.
  end subroutine river_allocate

  subroutine river_allocate_subcatch(n_subcatch)
    integer, intent(in) :: n_subcatch
    allocate(bsn_roc2(n_subcatch))
    bsn_roc2 = 0.
    allocate(bsn_roc4(n_subcatch))
    bsn_roc4 = 0.
  end subroutine river_allocate_subcatch

  subroutine river_initialise_travel_time(mb)
    integer, intent(in) :: mb
    integer j

    do j = 1, mb
      call river_muskingum_travel_time(j)
    end do
  end subroutine river_initialise_travel_time

  subroutine dealloc_river
   deallocate(chd)
   deallocate(chk)
   deallocate(chl)
   deallocate(chn)
   deallocate(chnn)
   deallocate(chs)
   deallocate(chss)
   deallocate(chw)
   deallocate(ovn)
   deallocate(phi)
   deallocate(qdilast)
   deallocate(qdinp)
   deallocate(qdolast)
   deallocate(qdout)
   deallocate(qsilast)
   deallocate(qsolast)
   deallocate(qssinp)
   deallocate(qssout)
   deallocate(rt)
   deallocate(sdtsav)
   deallocate(srch)
   deallocate(varoute)

  end subroutine dealloc_river

  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



  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&



  subroutine river_route_add(bRunoffdat, ihout, inum1, inum2, additionalGwUptake, bWAM_Module, daycounter, ida, iyr, mb, nqobs, obssb, obs_discharge, runs, subouthyd, inum1s, bRsvModule, rsv_is_operational, rsvSubbasin)
    !**** PURPOSE: THIS SUBROUTINE ADDS OUTPUTS FOR ROUTING
    !**** CALLED IN: MAIN
    !     ATTN! Look at the final line in the .fig file to define
    !           where the output has to be written: in route() or in add()
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> PARAMETERS IN TITLE
    !tit icode = code to switch between routing subroutines (here: 2)
    !tit ihout = hydrological storage location
    !tit inum1 = reach number
    !tit inum2 = inflow hydrograph(inum2 hydrograph is routed through inum1)
    !      >>>>>

    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      accf(4) = accumulated org. N, m3/sec.
    !      accf(5) = accumulated org. P, t
    !      accf(6) = accumulated NO3-N, kg
    !      ida = current day
    !      iy = current year as counter (1, ..., nbyr)
    !      iyr = current year
    !      runs(ida) = runoff simulated for the basin, m3/sec.
    !      varoute(1:8, ih) = vector 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
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      cnit = NO3-N concentration, mg/l
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    use management, only : &
        TSubbasin, bWAMHydrograph, management_is_transfer_subbasin, &
        management_subbasin_pointer, management_transfer_out, management_transfer_out
    use output, only : output_store_subbasin_values

    logical, intent(in) :: bRunoffdat
    real(dp), dimension(:), intent(out) :: additionalGwUptake
    logical, intent(in) :: bWAM_Module
    integer, intent(in) :: daycounter
    integer, intent(in) :: ida
    integer, intent(in) :: iyr
    integer, intent(in) :: mb
    integer, intent(in) :: nqobs
    integer, dimension(100), intent(in) :: obssb
    real(dp), dimension(:, :), intent(in) :: obs_discharge
    real(dp), dimension(366), intent(inout) :: runs
    integer, dimension(:), intent(in) :: subouthyd
    integer, intent(in) :: ihout, inum1, inum2
    integer, dimension(:), intent(in) :: inum1s
    logical, intent(in) :: bRsvModule, rsv_is_operational
    integer, dimension(:), intent(in) :: rsvSubbasin

    integer :: i
    real(dp) :: cnit
    ! ### RESERVOIR MODULE
    logical :: bAdd
    !#### WATER MANAGEMENT MODULE ####
    TYPE (TSubbasin), POINTER :: pS ! pointer on actual TSub

    if (ida .eq. 1 .and. inum2 .eq. 1) then
      accf(4) = 0.
      accf(5) = 0.
      accf(6) = 0.
    end if

    !#######################
    ! ### RESERVOIR MODULE
    !#######################
    bAdd = .true.

    if (bRsvModule) then
      if (inum2 <= mb) then
        if (rsvSubbasin(inum2) /= 0 .AND. rsv_is_operational) then
          varoute(2, ihout) = varoute(2, inum1)
          varoute(3, ihout) = varoute(3, inum1)
          varoute(4, ihout) = varoute(4, inum1)
          varoute(5, ihout) = varoute(5, inum1)
          varoute(6, ihout) = varoute(6, inum1)
          varoute(7, ihout) = varoute(7, inum1)
          varoute(8, ihout) = varoute(8, inum1)
          bAdd = .false. !!! DO NOT CALL add_varoute !!!
        else
          !#################################
          !#### WATER MANAGEMENT MODULE ####
          !#################################
          ! check if water management module is switched on
          if (bWAM_Module) then
            ! check if water transfers take place in current subbasin
            if (bWAMHydrograph(ihout) ) then
              if (management_is_transfer_subbasin(inum2) ) then
                pS => management_subbasin_pointer(inum2)
                call river_route_add_varoute
                !-------------------------------------------------------------
                ! Withdraw water from subbasin outlet
                ! but only if subbasin is not a headwater .and. not a reservoir.
                !-------------------------------------------------------------
                ! ATTENTION: The values of varoute(2, j) and varoute(8, j) may be changed!
                call management_transfer_out(inum2, varoute(2, ihout), varoute(8, ihout), daycounter, ida, iyr)

                pS % Q_act(daycounter) = real(varoute(2, ihout) / 86400. + varoute(8, ihout) / 86400.)
                !-------------------------------------------------------------
                bAdd = .false. !!! DO NOT CALL add_varoute AGAIN !!!
              end if ! ( wam_is_transfer_subbasin(inum2) )
            end if
          end if ! ( bWAM_Module )
          !#################################

        end if ! ( rsvSubbasin(inum2) /= 0 )
      end if ! if inum2 > mb

    else ! if (.NOT.bRsvModule)
      !#################################
      !#### WATER MANAGEMENT MODULE ####
      !#################################
      ! check if water management module is switched on
      if (bWAM_Module) then
        ! check if water transfers take place in current subbasin
        ! check if water transfers take place in current subbasin
        if (bWAMHydrograph(ihout) ) then
          if (management_is_transfer_subbasin(inum2) ) then
            pS => management_subbasin_pointer(inum2)
            call river_route_add_varoute
            !-------------------------------------------------------------
            ! Withdraw water from subbasin outlet
            ! but only if subbasin is not a headwater .and. not a reservoir.
            !-------------------------------------------------------------
            ! ATTENTION: The values of varoute(2, j) and varoute(8, j) may be changed!
            call management_transfer_out(inum2, varoute(2, ihout), varoute(8, ihout), daycounter, ida, iyr)

            pS % Q_act(daycounter) = real(varoute(2, ihout) / 86400. + varoute(8, ihout) / 86400.)
            !-------------------------------------------------------------
            bAdd = .false. !!! DO NOT CALL add_varoute AGAIN !!!
          end if ! ( wam_is_transfer_subbasin(inum2) )
        end if
      end if ! ( bWAM_Module )
      !#################################

    end if ! ( bRsvModule )
    !#######################! END IF ( bRsvModule )

    if (bAdd) call river_route_add_varoute !**** ADD FLOWS: 2 - surface flow, 8 - subsurface flow

    !*** Overwrite simulated routed discharge with observed, eg. for calibration
    if (bRunoffdat .and. nqobs > 1) then
      do i = 2, nqobs
        if (inum2 .eq. obssb(i) .and. obs_discharge(ida, i) >= 0. ) then
          varoute(2, ihout) = obs_discharge(ida, i) * 86400.
          varoute(8, ihout) = 0.
        endif
      enddo
    endif

    !**** WRITE DAILY OUTPUT in outlet (if inum2=1) to river output file unit=76
    !     WRITE is possible in variants - needed one should be opened
    if (inum2 == 1) then ! if inum2 is outlet (usually the last add command in the .fig file)

      if (bRsvModule) then ! if reservoir module is active
        if (rsvSubbasin(inum2) /= 0 ) then ! if subbasin inum2 is a reservoir
          runs(ida) = (varoute(2, inum1) + varoute(8, inum1)) / 86400.
        else
          runs(ida) = (varoute(2, ihout) + varoute(8, ihout)) / 86400.
        end if
      else
        runs(ida) = (varoute(2, ihout) + varoute(8, ihout)) / 86400.
      end if
      !runs(ida)=(varoute(2, ihout)+varoute(8, ihout))/86400.
      accf(4) = accf(4) + runs(ida)
      accf(5) = accf(5) + varoute(3, ihout)
      accf(6) = accf(6) + varoute(6, ihout)

      !**** CALC N CONCENTRATION
      !        varoute(6, ) - kg, runs - m3/sec, cnit - mg/l
      cnit = varoute(6, ihout) / runs(ida) / 86.4

      if (river_discharge_output_id > 0) then
        do i = 1, mb
          if (bRsvModule) then
            if (rsvSubbasin(i) /= 0 ) then ! if subbasin is a reservoir
              runsub_m3s(i) = (varoute(2, inum1s(subouthyd(i))) + varoute(8, inum1s(subouthyd(i)))) / 86400.
            else
              runsub_m3s(i) = (varoute(2, subouthyd(i)) + varoute(8, subouthyd(i))) / 86400.
            end if
          else
            runsub_m3s(i) = (varoute(2, subouthyd(i)) + varoute(8, subouthyd(i))) / 86400.
          end if
        end do

        call output_store_subbasin_values(river_discharge_output_id, runsub_m3s)
      end if

      additionalGwUptake(1:mb) = 0

    end if ! (inum2 .eq. 1)

  contains
    ! Subroutines inside the `add` subroutine

    !--------------------------------------------------------
    subroutine river_route_add_varoute
      varoute(2, ihout) = varoute(2, inum1) + varoute(2, inum2)
      varoute(3, ihout) = varoute(3, inum1) + varoute(3, inum2)
      varoute(4, ihout) = varoute(4, inum1) + varoute(4, inum2)
      varoute(5, ihout) = varoute(5, inum1) + varoute(5, inum2)
      varoute(6, ihout) = varoute(6, inum1) + varoute(6, inum2)
      varoute(7, ihout) = varoute(7, inum1) + varoute(7, inum2)
      varoute(8, ihout) = varoute(8, inum1) + varoute(8, inum2)
    end subroutine river_route_add_varoute
    !--------------------------------------------------------

  end subroutine river_route_add



  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  subroutine river_transfer()
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !**** PURPOSE: THIS SUBROUTINE CONTROLS THE CHANNEL ROUTING, NOT USED!
    !**** CALLED IN: MAIN, "sleeping subroutine"
  end subroutine river_transfer
  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

  subroutine river_transmission_loss(j)
    ! *** PURPOSE: this subroutine calculate the transmssion and evaporation losses in river reaches (by Huang, 2013, based on SWAT code)
    ! *** CALLED IN: ROUTE
    !     METHOD: SWAT
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      chl(2, j) = channel length, km
    !      evrch |none |Reach evaporation adjustment factor.
    !                               |Evaporation from the reach is multiplied by
    !                               |EVRCH. This variable was created to limit the
    !                               |evaporation predicted in arid regions.
    !      phi(5, j) = "normal" flow
    !      sdti = water inflow + storage in the reach, m3/sec.
    !      sdtsav(ir) = water storage in the reach, m3
    !      xxqd = surface flow - daily input in reaches, m3
    !      xxssf = subsurface + g-w flow - daily input in reaches, m3
    !      tlc = transmission losses, m3
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      nreach = reach No.
    !      yy = local par
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    ! *** Include common parameters

    integer :: j, nreach
    real(dp) :: d, b, volrt, rchdep, p, rh, rcharea, rttime1
    real(dp) :: adddep, addp, addarea, tlc1, tlc2, evp1, evp2, topw
    real(dp) :: qdandssf, qd_ssf
    real(dp) :: x, w, Vconvert, k_convert, intercept, xx, decay, slope
    real(dp) :: slope2, intercept2, chside

    nreach = j

    chside = 2.
    d = chd(j)
    b = chw(2, j) - 2. * d * chside
    volrt = sdti
    sdti = 0.
    rchdep = 0.
    p = 0.
    rh = 0.
    qd_ssf = 0.
    tlc = 0.
    evp = 0.
    addp = 0.
    adddep = 0.
    rcharea = 0.
    rttime1 = 0.
    addarea = 0.
    tlc1 = 0.
    tlc2 = 0.
    evp1 = 0.
    evp2 = 0.
    topw = 0.
    qdandssf = 0.
    x = 0.
    w = 0.
    Vconvert = 0.
    k_convert = 0.
    intercept = 0.
    xx = 0.
    decay = 0.
    slope = 0.
    slope2 = 0.
    intercept2 = 0.

    !**** CALC bottom width b & channel side chside
    if (b <= 0.) then
      b = .5 * chw(2, j)
      chside = (chw(2, j) - b) / (2. * d)
    end if
    !**** calculate the current water level
    if (volrt >= phi(5, j) ) then
      rcharea = b * d + chside * d * d
      rchdep = d
      p = b + 2. * d * sqrt(chside * chside + 1)
      rh = rcharea / p
      sdti = phi(5, j)
      !! this is the time to empty the volume of water at the bankfull Q
      rttime1 = volrt * 86400. / (3600. * sdti)
      if (rttime1 > 24.) then
        !! perform flood plain simulation
        !! increase Q in flood plain until all the volume can be emptied in one day
        adddep = 0.
        do while ( rttime1 > 24.)
          !! 1 cm interval
          adddep = adddep + 0.01
          addarea = rcharea + ((chw(2, j) * 5.) + chside * adddep) * adddep
          addp = p + chw(2, j) * 4. + 2. * adddep * sqrt(17.)
          rh = addarea / addp
          sdti = river_mannings_discharge(addarea, rh, chnn(j), chss(j))
          rttime1 = volrt * 86400. / (3600. * sdti)
        end do
        rcharea = addarea
        rchdep = rchdep + adddep
        p = addp
      end if

    else
      do while (sdti .lt. volrt)
        !! find the cross sectional area and depth for volrt
        !! 1 cm interval depth
        rchdep = rchdep + 0.01
        rcharea = (b + chside * rchdep) * rchdep
        p = b + 2. * rchdep * sqrt(chside * chside + 1)
        rh = rcharea / p
        sdti = river_mannings_discharge(rcharea, rh, chnn(j), chss(j))
      end do
    end if ! ( volrt >= phi(5, j) )

    !**** sdti
    sdti = volrt

    if (sdti .gt. 0.) then

    !**** calculate the transmission losses in rivers
    qd_ssf = xxqd / (xxqd + xxssf)

    if (tlrch > 0.) then
      qdandssf = xxqd + xxssf
      if (qdandssf .gt. 0.) then
        if (volrt .ge. phi(5, j)) then ! SWAT method only for flood conditions
          tlc = tlrch * 24. * chk(2, j) * chl(2, j) * p
        else ! WASA method for normal conditions, water is within the channel
          x = chl(2, j) * 1000 * 0.00062 !river lenth, conversion kilo meter to miles 1000 * 0.00062
          w = 10. ** (log10(sdti) * 0.494 + 1.031) * 3.281 ! river width, convertion meter to feet, 3.281
          Vconvert = sdti * 0.00081 * 86400 ! conversion m3 / s to acre - feet
          k_convert = tlrch * 0.0394 * chk(2, j) ! mm / hr to inch / hr

          intercept = - 0.00465 * k_convert * 24.
          xx = (1 - 0.00545 * (k_convert * 24 / Vconvert))

          if (xx .gt. 0.) then
            decay = - 1.09 * log(1. - 0.00545 * (k_convert * 24. / Vconvert))
            slope = exp(- decay)
            slope2 = exp(- decay * w * x)

            if ((1. - slope) .gt. .01) then
              intercept2 = (intercept / (1. - slope)) * (1. - slope2)
            else
              intercept2 = (intercept / 0.01) * (1. - slope2)
            end if

            tlc = (Vconvert - (intercept2 + slope2 * Vconvert)) * 1234.6 !acre - feet to m3
          else
            tlc = 0.
          end if
        end if ! (volrt .ge. phi(5, j))

        tlc2 = tlc * sdtsav(nreach) / (xxqd + xxssf + sdtsav(nreach))

        if (sdtsav(nreach) .le. tlc2) then
          tlc2 = sdtsav(nreach)
          sdtsav(nreach) = sdtsav(nreach) - tlc2
          tlc1 = tlc - tlc2
          if (qdandssf .le. tlc1) then
            tlc1 = xxqd + xxssf
            xxqd = 0.
            xxssf = 0.
          else
            xxqd = xxqd - tlc1 * qd_ssf
            xxssf = xxssf - tlc1 * (1. - qd_ssf)
          end if
        else
          sdtsav(nreach) = sdtsav(nreach) - tlc2
          tlc1 = tlc - tlc2
          if (qdandssf .le. tlc1) then
            tlc1 = xxqd + xxssf
            xxqd = 0.
            xxssf = 0.
          else
            xxqd = xxqd - tlc1 * qd_ssf
            xxssf = xxssf - tlc1 * (1. - qd_ssf)
          end if
        end if
        tlc = tlc1 + tlc2
      end if ! (qdandssf .gt. 0.)
    end if ! ( tlrch > 0. )

    !**** calculate the evaporation losses in rivers (by Huang, 2013, based on SWAT code)
    if (evrch > 0.) then
        qdandssf = xxqd + xxssf
        if (qdandssf .gt. 0.) then
          !! calculate width of channel at water level
          topw = 0.
          if (rchdep .le. d) then
            !topw = b + 2. *rchdep*chside !SWAT version
            topw = 10. ** (log10(sdti) * 0.494 + 1.031) ! river width, WASA version
          else
            topw = 5. * chw(2, j) + 2. * (rchdep - chd(j)) * 4.
          end if

          evp = evrch * pet_day(j) * chl(2, j) * topw
          evp2 = evp * sdtsav(nreach) / (xxqd + xxssf + sdtsav(nreach))

          if (sdtsav(nreach) .le. evp2) then
            evp2 = sdtsav(nreach)
            sdtsav(nreach) = sdtsav(nreach) - evp2
            evp1 = evp - evp2
            if (qdandssf .le. evp1) then
              evp1 = xxqd + xxssf
              xxqd = 0.
              xxssf = 0.
            else
              xxqd = xxqd - evp1 * qd_ssf
              xxssf = xxssf - evp1 * (1 - qd_ssf)
            end if
          else
            sdtsav(nreach) = sdtsav(nreach) - evp2
            evp1 = evp - evp2
            if (qdandssf .le. evp1) then
              evp1 = xxqd + xxssf
              xxqd = 0.
              xxssf = 0.
            else
              xxqd = xxqd - evp1 * qd_ssf
              xxssf = xxssf - evp1 * (1. - qd_ssf)
            end if
          end if

          evp = evp1 + evp2
        end if ! (qdandssf .gt. 0.)
      end if ! evrch > 0.

    else
      xxqd = 0.
      xxssf = 0.
      sdti = 0.
    end if ! (sdti .gt. 0.)

    if (xxqd .lt. 0.) xxqd = 0.
    if (xxssf .lt. 0.) xxssf = 0.
    if (sdti .lt. 0.) sdti = 0.
    if (sdtsav(nreach) .lt. 0.) sdtsav(nreach) = 0.

  end subroutine river_transmission_loss


  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&


  subroutine river_muskingum_routing(j, ihout, inum1, ida, iy, iyr)
    !**** PURPOSE: this subroutine routes hydrograph to basin outlet
    !**** CALLED IN: ROUTE
    !     METHOD: MUSKINGUM
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !tit j=inum1 = Reach No.
    !tit inum2 = Inflow Hyd.(inum2 hydrograph is routed through
    !                 inum1 hydrograph)

    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      chl(2, j) = channel length, km
    !      ida = current day
    !      iyr = current year
    !      phi(5, j) = "normal" flow
    !      phi(10, j) = prelim. estimation of xkm for bankfull depth
    !      phi(13, j) = prelim. estimation of xkm for 0.1 bankfull depth
    !      qdilast(ih) = qdinp(, ) in the last day of the year, m3
    !      qdinp(ih, ida) = surface flow - daily input in reaches, m3
    !      qdolast(ih) = qdout(, ) in the last day pof the year, m3
    !      qdout(ih, ida) = surface flow - daily output in reaches, m3
    !      qsilast(ih) = qssinp(, ) in the last day of the year, m3
    !      qsolast(ih) = qssout(, ) in the last day of the year, m3
    !      qssinp(ih, ida) = subsurface + g-w flow - daily input in reaches, m3
    !      qssout(ih, ida) = subsurface + g-w flow - daily output in reaches, m3
    !      roc1 = 0
    !      roc2 = calibration parameter for routing
    !      roc3 = 0
    !      roc4 = calibration parameter for routing
    !      sdti = water inflow + storage in the reach, m3/sec.
    !      sdtsav(ir) = water storage in the reach, m3
    !      xxqd = surface flow - daily input in reaches, m3
    !      xxssf = subsurface + g-w flow - daily input in reaches, m3
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      c1 = local par
    !      c2 = local par
    !      c3 = local par
    !      c4 = local par
    !      det = local par
    !      nreach = reach No.
    !      qdin = xxqd
    !      ssfin = xxssf
    !      xat = dimentionless weighting factor for routing
    !      xkm = storage time constant for the reach
    !      yy = local par
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters

    integer, intent(in) :: ida
    integer, intent(in) :: iy
    integer, intent(in) :: iyr
    integer :: j, ihout, inum1, nreach
    !, det, xat, yy
    real(dp) :: qdin, ssfin

    ! SL only a preliminary solition
    if (xxqd <= 0.) then
      call log_warn("river_muskingum_routing", "1. xxqd in (day, subbasin) is negative! Will set to almost zero", &
        i1=ida, i2=j, real=xxqd)
      xxqd = max(0.00006_dp, xxqd)
    end if
    if (xxssf <= 0.) then
      call log_warn("river_muskingum_routing", "1. xxssf in (day, subbasin) is negative! Will set to almost zero", &
        i1=ida, i2=j, real=xxssf)
      xxssf = max(0.00006_dp, xxssf)
    end if

    nreach = inum1
    qdinp(ihout, ida) = xxqd
    qssinp(ihout, ida) = xxssf
    qdin = xxqd
    ssfin = xxssf

    !**** COMPUTE sdti = INFLOW + STORAGE, transfer m3 --> m3/sec
    sdti = (qdin + ssfin + sdtsav(nreach)) / 86400.

    !*********************************************************** STEP 1
    !**** COMPUTE SURFACE flow routing
    ! SL The variables xkm, c1, c2, c3, and c4 were calculated for surface and subsurface routing
    !     every time step this function is called.
    !     This is actually not necessary because they need to be computed based on subbasins'
    !     channell parameters only once during the pre-processing.
    !     This is now done in subroutine ttcoefi()

    if (iy .eq. 1 .and. ida .eq. 1) then
      xxqd = c1_qd(j) * xxqd + c2_qd(j) * sdtsav(inum1) + c3_qd(j) * sdtsav(inum1) + c4_qd(j)
    else
      if (ida .gt. 1) then
        xxqd = c1_qd(j) * xxqd + c2_qd(j) * qdinp(ihout, ida - 1) + c3_qd(j) * qdout(ihout, ida - 1)
      else
        xxqd = c1_qd(j) * xxqd + c2_qd(j) * qdilast(ihout) + c3_qd(j) * qdolast(ihout)
      end if
    end if

    !*********************************************************** STEP 2
    !**** COMPUTE SUBSURFACE flow routing
    ! SL The variables xkm, c1, c2, c3, and c4 were calculated for surface and subsurface routing
    !     every time step this function is called.
    !     This is actually not necessary because they need to be computed based on subbasins'
    !     channell parameters only once during the pre-processing.
    !     This is now done in subroutine ttcoefi()

    if (iy .eq. 1 .and. ida .eq. 1) then
      xxssf = c1_ssf(j) * xxssf + c2_ssf(j) * sdtsav(inum1) + c3_ssf(j) * sdtsav(inum1) + c4_ssf(j)
    else
      if (ida .gt. 1) then
        xxssf = c1_ssf(j) * xxssf + c2_ssf(j) * qssinp(ihout, ida - 1) + c3_ssf(j) * qssout(ihout, ida - 1)
      else
        xxssf = c1_ssf(j) * xxssf + c2_ssf(j) * qsilast(ihout) + c3_ssf(j) * qsolast(ihout)
      end if
    end if


    !sl begin
    !!!!!! DO NOT ALLOW NEGATIVE FLOWS
    !!!!!! This shouldn't be a long-term solution but an attempt to make the user aware of
    !!!!!! obviously wrong parameter settings
    if (xxqd < 0.) then
      call log_warn("river_muskingum_routing", "2. xxqd in (day, subbasin) is negative! Will set to almost zero", &
        i1=ida, i2=j, real=xxqd)
      xxqd = 0.000006
    end if
    if (xxssf < 0.) then
      call log_warn("river_muskingum_routing", "2. xxssf in (day, subbasin) is negative! Will set to almost zero", &
        i1=ida, i2=j, real=xxssf)
      xxssf = 0.000006
    end if
    !sl end

    !**** CALC outputs
    qdout(ihout, ida) = xxqd
    qssout(ihout, ida) = xxssf

    !**** CALC water storage in the reach
    sdtsav(nreach) = sdtsav(nreach) + qdin + ssfin - xxqd - xxssf
    if (sdtsav(nreach) .lt. 0.1) then
      sdtsav(nreach) = 0.
    endif

    !**** STORE input-output in the last day of the year:
    if (mod(iyr, 4) .eq. 0 .and. ida .eq. 366) then
      qdilast(ihout) = qdinp(ihout, ida)
      qdolast(ihout) = qdout(ihout, ida)
      qsilast(ihout) = qssinp(ihout, ida)
      qsolast(ihout) = qssout(ihout, ida)
    endif
    if (mod(iyr, 4) .ne. 0 .and. ida .eq. 365) then
      qdilast(ihout) = qdinp(ihout, ida)
      qdolast(ihout) = qdout(ihout, ida)
      qsilast(ihout) = qssinp(ihout, ida)
      qsolast(ihout) = qssout(ihout, ida)
    endif

  end subroutine river_muskingum_routing

  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&


  subroutine river_route_erosion(j, inum1, chc, chxk, yd)
    !**** PURPOSE: THIS SUBROUTINE ROUTES SEDIMENT FROM SUB-BASIN TO BASIN OUTLET
    !**** CALLED IN: ROUTE
    !      METHOD: from J. WILLIAMS (SWAT98)
    !      DEPOSITION IS BASED ON FALL VELOCITY AND DEGRADATION - ON STREAM POWER
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !      PARAMETERS & VARIABLES
    !
    !       >>>>> COMMON PARAMETERS & VARIABLES
    !       chc(j) = channel USLE C factor
    !       chnn(j) = channel N value
    !       chss(j) = channel slope, m/m
    !       chw(2, j) = average width of main channel, m
    !       chxk(j) = channel USLE K factor
    !       ida = current day
    !       irout = switch code to print from routfun()
    !       prf = coef. to estimate peak runoff in stream
    !       sdti = water inflow + storage in the reach, m3/sec.
    !       sdtsav(ir) = water storage in the reach, m3
    !       spcon = rate parameter for estimation of sediment transport, readbas
    !       spexp = exponent for estimation of sediment transport, readbas
    !       xxqd = surface flow - daily input in reaches, m3
    !       yd = daily soil loss from subbasin caused by water erosion, t
    !       >>>>>

    !       >>>>> STATIC PARAMETERS
    !       ach = local par
    !       bwd = local par
    !       cych = local par
    !       cyin = local par
    !       d = f(sdti, chnn, chw, chss)
    !       deg = degradation
    !       dep = deposition
    !       depnet = local par
    !       nsz = local par
    !       prst = peak runoff in stream
    !       qdin = surface inflow + water storage
    !       vc = f(sdti, chw, d)
    !       vs = local par
    !       ydin = yd
    !       >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !**** Include common parameters

    real(dp), dimension(:), intent(in) :: chc
    real(dp), dimension(:), intent(in) :: chxk
    real(dp), intent(inout) :: yd
    integer j, inum1
    real(dp) ach, cych, cyin, d, deg, dep, depnet, prst, qdin
    real(dp) vc, ydin

    prst = prf * sdti
    qdin = xxqd + sdtsav(inum1)
    ydin = yd
    if (qdin .le. 0.01) return

    !**** COMPUTE FLOW DEPTH
    d = ((prst * chnn(j)) / (chw(2, j) * chss(j) ** .5)) ** .6
    ach = d * chw(2, j)
    if (d .lt. .010) then
      vc = 0.01
    else
      vc = prst / ach
    endif
    if (vc .gt. 3.) vc = 3.

    !**** CALCULATE DEPOSITION & DEGRADATION (either - or)
    cyin = ydin / qdin
    cych = spcon * vc ** spexp
    depnet = qdin * (cych - cyin)
    if (depnet .gt. 0.) then
      deg = depnet * chxk(j) * chc(j)
      dep = 0.
    else
      dep = - depnet
      deg = 0.
    endif

    yd = yd + deg - dep
    if (yd .lt. 0.) yd = 0.

    return
  end subroutine river_route_erosion


  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&


  subroutine river_erosion_enritchment_ratio(inum1, da9, er, yd)
    !**** THIS SUBROUTINE CALCULATES ENRICHMENT RATIO
    !**** CALLED IN: ROUTE
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !      PARAMETERS & VARIABLES
    !
    !       >>>>> COMMON PARAMETERS & VARIABLES
    !       da9 = drainage area, ha
    !       er = enrichment coefficient
    !       sdtsav(ir) = water storage in the reach, m3
    !       xxqd = surface flow - daily input in reaches, m3
    !       ydi = yd, daily soil loss caused by water erosion, t,
    !                    in route before call river_route_erosion
    !       yd = daily soil loss caused by water erosion, t, recalc in river_route_erosion
    !       >>>>>

    !       >>>>> STATIC PARAMETERS
    !       cy = local par
    !       dr = local par
    !       qdin = inflow + watre storage
    !       x1 = local par
    !       x2 = local par
    !       >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters

    real(dp), intent(in) :: da9
    real(dp), intent(inout) :: er
    real(dp), intent(in) :: yd
    integer inum1
    real(dp) cy, dr, qdin, x1, x2

    qdin = xxqd + sdtsav(inum1)
    if (qdin .le. 0.1E-20) return

    !**** CALC enrichment coefficient
    dr = yd / (ydi + 1.e-10)
    if (dr .le. 0.01) dr = .01
    if (dr .ge. 10.) dr = 10.
    cy = 100000. * ydi / (da9 * qdin)
    x2 = - log10(dr) / 2.699
    x1 = 1. / .25 ** x2
    er = x1 * (cy + 1.e-10) ** x2
    if (er .lt. 1.) er = 1.
    if (er .gt. 3.) er = 3.

    return
  end subroutine river_erosion_enritchment_ratio



  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&



  subroutine river_route_nitrate(conn, da9, er, yon)
    !**** PURPOSE: THIS SUBROUTINE CALCULATES ORGANIC N ROUTING
    !**** CALLED IN: ROUTE
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !      PARAMETERS & VARIABLES
    !
    !       >>>>> COMMON PARAMETERS & VARIABLES
    !       conn = xnorg * er, g/t: N org. in I layer for subbasin,
    !               corrected for enrichment
    !       da9 = drainage area, ha
    !       er = enrichment coefficient
    !       ydi = yd, daily soil loss caused by water erosion, t,
    !               in route before call rtsed
    !       yon = org N loss with erosion, kg/ha, routed
    !       >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters

    real(dp), intent(in) :: conn
    real(dp), intent(in) :: da9
    real(dp), intent(in) :: er
    real(dp), intent(out) :: yon

    yon = .001 * conn * er * ydi / da9
    return
  end subroutine river_route_nitrate



  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&



  subroutine river_route_phosphorus(cpp, da9, yph)
    !**** PURPOSE: THIS SUBROUTINE CALCULATES PHOSPHORUS ROUTING
    !**** CALLED IN: ROUTE
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !      PARAMETERS & VARIABLES
    !
    !       >>>>> COMMON PARAMETERS & VARIABLES
    !       cpp = xporg * er, g/t: P org. in I layer in subbasin, g/t
    !                  corrected for enrichment
    !       da9 = drainage area, ha
    !       ydi = yd, daily soil loss caused by water erosion, t,
    !                 in route before call rtsed
    !       yph = P org. loss with erosion, kg/ha
    !       >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters

    real(dp), intent(in) :: cpp
    real(dp), intent(in) :: da9
    real(dp), intent(out) :: yph

    yph = .001 * cpp * ydi / da9
    return
  end subroutine river_route_phosphorus



  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&



  subroutine river_muskingum_travel_time(j)
    !**** THIS SUBROUTINE COMPUTES TRAVEL TIME COEFFICIENTS FOR ROTING,
    !      CALLS COEFS and QMAN
    !**** CALLED IN: MAIN
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !      PARAMETERS & VARIABLES
    !
    !       >>>>> COMMON PARAMETERS & VARIABLES
    !       chd(j) = channel depth, m
    !       chl(2, j) = main channel length, km
    !       chnn(j) = channel N value
    !       chss(j) = channel slope, m/m
    !       chw(2, j) = average width of main channel, m
    !       phi(1, j), phi(2, j) = coefs for travel time estimation, if higher than
    !                           normal flow = 1.2 depth
    !       phi(3, j), phi(4, j) = coefs for travel time estimation, if lower than
    !                           normal flow = 0.1 depth
    !       phi(5, j) = flow rate when reach is at bankfull depth, m3/sec
    !       phi(6, j) = bottom width of main channel, m
    !       phi(7, j) = depth of water when reach is at bankfull depth, m
    !       phi(8, j) = average velocity when reach is at bankfull depth, m/s
    !       phi(9, j) = wave celerity when reach is at bankfull depth, m/s
    !       phi(10, j) = prelim. estimation of xkm for bankfull depth:
    !                         storage time constant for reach at bankfull depth
    !                         (ratio of storage to discharge), h
    !       phi(11, j) = average velocity when reach is at 0.1 bankfull depth
    !                         (low flow), m/s
    !       phi(12, j) = wave celerity when reach is at 0.1 bankfull depth
    !                         (low flow), m/s
    !       phi(13, j) = prelim. estimation of xkm for 0.1 bankfull depth:
    !                         storage time constant for reach at 0.1 bankfull depth
    !                         (low flow) (ratio of storage to discharge), h
    !       >>>>>

    !       >>>>> STATIC PARAMETERS
    !       a = local par
    !       a1 = local par
    !       aa = local par
    !       b = local par
    !       celer = celerity
    !       chside = local par
    !       d = local par
    !       fpn = local par
    !       fps = local par
    !       p = local par
    !       p1 = from coefs
    !       pp2 = from coefs
    !       q2 = average flow
    !       qq1 = local par
    !       rh = hydraulic radius (m)
    !       tmne = local par
    !       tt1 = local par
    !       tt2 = local par
    !       velos = velocity
    !       ykm = xkm in hours (I estimation of xkm)
    !       >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters

    integer j
    real(dp) a, a1, aa, b, celer, chside, d, fpn, fps, p, p1
    real(dp) pp2, q2, qq1, rh, tmne, tt1, tt2, velos, ykm
    ! real(dp) qman
    real(dp) :: det, xat, yy

    chside = 2.
    fps = 4.
    fpn = chnn(j)
    d = chd(j)
    b = chw(2, j) - 2. * d * chside

    !**** CALC bottom width b & channel side chside
    if (b .le. 0.) then
      b = .5 * chw(2, j)
      chside = (chw(2, j) - b) / (2. * d)
    endif

    !**** CALC average flow q2
    phi(6, j) = b
    phi(7, j) = d
    p = b + 2. * d * sqrt(chside * chside + 1)
    a = b * d + chside * d * d
    rh = a / p
    q2 = river_mannings_discharge(a, rh, chnn(j), chss(j))

    !**** CALC velocity, celerity and storage time constant for bankfull depth
    aa = 1.
    velos = river_mannings_discharge(aa, rh, chnn(j), chss(j))
    celer = velos * 5. / 3.
    !      UNITS: chl() -in km!  km/ --> m/sec /3600 *1000
    ykm = chl(2, j) / celer / 3.6
    phi(5, j) = q2
    phi(8, j) = velos
    phi(9, j) = celer
    phi(10, j) = ykm
    tt2 = chl(2, j) * a / q2

    !**** CALC velocity, celerity and storage time constant for 1.2 bankfull depth
    !      depth = 1.2 * chd
    !#### CALL COEFS
    d = 1.2 * chd(j)
    a1 = chw(2, j) * chd(j) + fps * (d - chd(j)) ** 2
    p1 = 2. * chd(j) * sqrt(fps * fps + 1)
    rh = (a + a1) / (p + p1)
    tmne = (chnn(j) * a + fpn * a1) / (a + a1)
    a = a + a1
    qq1 = river_mannings_discharge(a, rh, tmne, chss(j))
    tt1 = chl(2, j) * a / qq1
    call river_travel_time_coefficients(qq1, q2, tt1, tt2, p1, pp2)
    phi(1, j) = p1
    phi(2, j) = pp2

    !**** CALC velocity, celerity and storage time constant for 0.1 bankfull depth
    !      depth = 0.1 * chd
    !#### CALL COEFS
    d = 0.1 * chd(j)
    p = b + 2. * d * sqrt(chside * chside + 1)
    a = b * d + chside * d * d
    rh = a / p
    qq1 = river_mannings_discharge(a, rh, chnn(j), chss(j))
    tt1 = chl(2, j) * a / qq1
    call river_travel_time_coefficients(qq1, q2, tt1, tt2, p1, pp2)
    phi(3, j) = p1
    phi(4, j) = pp2
    aa = 1.
    velos = river_mannings_discharge(aa, rh, chnn(j), chss(j))
    celer = velos * 5. / 3.
    !      UNITS: chl() -in km!  km/ --> m/sec /3600 *1000
    ykm = chl(2, j) / celer / 3.6
    phi(11, j) = velos
    phi(12, j) = celer
    phi(13, j) = ykm


    ! SL BEGIN ------------------------------------------------------------------------
    xat = 0.05 !0.1 ! 0.2 ! dimensionless weighting factor for routing
    det = 24.

    !#### Compute MUSKINGUM ROUTING PARAMETERS for surface flow component
    xkm_qd(j) = phi(10, j) * roc1 + phi(13, j) * roc2(j)
    yy = 2. * xkm_qd(j) * (1. - xat) + det
    c1_qd(j) = (det - 2. * xkm_qd(j) * xat) / yy
    c2_qd(j) = (det + 2. * xkm_qd(j) * xat) / yy
    c3_qd(j) = (2. * xkm_qd(j) * (1. - xat) - det) / yy
    c4_qd(j) = phi(5, j) * chl(2, j) * det / yy
    !*** Correction to avoid oszillation and negative discharge
    if (c3_qd(j) < 0. ) then
      c2_qd(j) = c2_qd(j) + c3_qd(j)
      c3_qd(j) = 0.
    end if
    if (c1_qd(j) < 0.) then
      c2_qd(j) = c2_qd(j) + c1_qd(j)
      c1_qd(j) = 0.
    end if

    !#### Compute MUSKINGUM ROUTING PARAMETERS for sub-surface flow component
    xkm_ssf(j) = phi(10, j) * roc3 + phi(13, j) * roc4(j)
    yy = 2. * xkm_ssf(j) * (1. - xat) + det
    c1_ssf(j) = (det - 2. * xkm_ssf(j) * xat) / yy
    c2_ssf(j) = (det + 2. * xkm_ssf(j) * xat) / yy
    c3_ssf(j) = (2. * xkm_ssf(j) * (1. - xat) - det) / yy
    c4_ssf(j) = phi(5, j) * chl(2, j) * det / yy
    !*** Correction to avoid oszillation and negative discharge
    if (c3_ssf(j) < 0. ) then
      c2_ssf(j) = c2_ssf(j) + c3_ssf(j)
      c3_ssf(j) = 0.
    end if
    if (c1_ssf(j) < 0.) then
      c2_ssf(j) = c2_ssf(j) + c1_ssf(j)
      c1_ssf(j) = 0.
    end if

    return
  end subroutine river_muskingum_travel_time


  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&


  subroutine river_travel_time_coefficients(qq1, q2, tt1, tt2, p1, pp2)
    !**** PURPOSE: THIS SUBROUTINE CALCULATES coefficients p1 and pp2
    !**** CALLED IN: TTCOEFI
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters

    real(dp) qq1, q2, tt1, tt2, p1, pp2

    pp2 = log(tt1 / tt2) / log(qq1 / q2)
    if (pp2 .lt. - 1.) pp2 = - 1.
    if (pp2 .gt. 1.5) pp2 = 1.5
    p1 = tt1 / (qq1 ** pp2)
    return
  end subroutine river_travel_time_coefficients



  !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&



  function river_mannings_discharge(a, rh, xn, chslope)
    !**** PURPOSE: THIS SUBROUTINE CALCULATES FLOW USING MANNINGS EQUATION
    !**** CALLED IN: TTCOEFI
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters

    real(dp) a, rh, xn, chslope, river_mannings_discharge

    river_mannings_discharge = a * rh ** .6666 * sqrt(chslope) / xn
    return
  end function river_mannings_discharge

end module river
