! 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) 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