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