this is the time to empty the volume of water at the bankfull Q perform flood plain simulation increase Q in flood plain until all the volume can be emptied in one day 1 cm interval find the cross sectional area and depth for volrt 1 cm interval depth calculate width of channel at water level
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer | :: | j |
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