!!!! 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
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer | :: | j | ||||
| integer | :: | ihout | ||||
| integer | :: | inum1 | ||||
| integer, | intent(in) | :: | ida | |||
| integer, | intent(in) | :: | iy | |||
| integer, | intent(in) | :: | iyr |
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