subroutine ch_rtmusk !! ~ ~ ~ PURPOSE ~ ~ ~ !! this subroutine routes a daily flow through a reach using the !! Muskingum method !! code provided by Dr. Valentina Krysanova, Pottsdam Institute for !! Climate Impact Research, Germany !! Modified by Balaji Narasimhan !! Spatial Sciences Laboratory, Texas A&M University use basin_module use channel_data_module use channel_module use hydrograph_module !, only : ob, icmd, jrch, isdch, fp_stor, ch_stor, wet use time_module use channel_velocity_module use sd_channel_module use climate_module use reservoir_module use reservoir_data_module use water_body_module use hru_module, only : hru use conditional_module implicit none integer :: ii !none |current day of simulation integer :: ihru integer :: iihru integer :: icha integer :: irtstep integer :: isubstep integer :: ires integer :: ihyd integer :: irel real :: ch_stor_init !m3 |storage in channel at beginning of day real :: fp_stor_init !m3 |storage in flood plain above wetlands emergency spillway at beginning of day real :: wet_stor_init !m3 |storage in flood plain wetlands at beginning of day real :: tot_stor_init real :: inout !m3 |inflow - outflow for day real :: del_stor !m3 |change in storage of channel + flood plain + wetlands real :: topw !m |top width of main channel real :: qinday !units |description real :: qoutday !units |description real :: inflo !m^3 |inflow water volume real :: inflo_rate !m^3/s |inflow rate real :: dep_flo = 0. !m |depth of flow real :: ttime !hr |travel time through the reach real :: outflo !m^3 |outflow water volume real :: tl !m^3 |transmission losses during time step real :: trans_loss = 0. !m^3 |transmission losses during day real :: ev !m^3 |evaporation during time step real :: evap = 0. !m^3 |evaporation losses during day real :: rto real :: rto1 real :: rto_w real :: rto_emer real :: outflo_rate real :: dts !seconds |time step interval for substep real :: dthr real :: scoef real :: vol_ch real :: sum_inflo, sum_outflo real :: dep real :: evol_m3, pvol_m3 real :: wet_evol jrch = isdch jhyd = sd_dat(jrch)%hyd qinday = 0 qoutday = 0 ht2 = hz ob(icmd)%hyd_flo = 0. hyd_rad = 0. trav_time = 0. flo_dep = 0. trans_loss = 0. evap = 0. !***jga !ob(icmd)%tsin = (/0., 800., 2000., 4200., 5200., 4400., 3200., 2500., 2000., 1500., 1000., 700., 400., & ! 0., 0., 0., 0., 0., 1000000., 1000000., 1000000., 1000000., 1000000., 1000000./) sum_inflo = sum (ob(icmd)%tsin) !! total wetland volume at start of day wet_stor(jrch) = hz wet_evol = 0. do ihru = 1, sd_ch(jrch)%fp%hru_tot iihru = sd_ch(jrch)%fp%hru(ihru) wet_stor(jrch) = wet_stor(jrch) + wet(iihru) wet_evol = wet_evol + wet_ob(iihru)%evol end do wet_stor_init = wet_stor(jrch)%flo ch_stor_init = ch_stor(jrch)%flo fp_stor_init = fp_stor(jrch)%flo tot_stor_init = ch_stor_init + fp_stor_init irtstep = 1 isubstep = 0 dts = time%dtm / sd_ch(jrch)%msk%substeps * 60. dthr = dts / 3600. !! subdaily time step do ii = 1, sd_ch(jrch)%msk%nsteps !! water entering reach during time step - substeps for stability isubstep = isubstep + 1 if (isubstep > sd_ch(jrch)%msk%substeps) then irtstep = irtstep + 1 isubstep = 1 end if inflo = ob(icmd)%tsin(irtstep) / sd_ch(jrch)%msk%substeps inflo_rate = inflo / dts !! interpolate rating curve using inflow rates icha = jrch call rcurv_interp_flo (icha, inflo_rate) ch_rcurv(jrch)%in2 = rcurv !! save variables at each routing time step for sediment routing if (isubstep == 1 .and. rcurv%wet_perim > 1.e-6) then hyd_rad(irtstep) = rcurv%xsec_area / rcurv%wet_perim trav_time(irtstep) = rcurv%ttime flo_dep(irtstep) = rcurv%dep end if !! add inflow to total storage if (ht1%flo > 1.e-6) then rto = inflo / ht1%flo rto = Max(0., rto) rto = Min(1., rto) tot_stor(jrch) = tot_stor(jrch) + rto * ht1 end if !! if no water in channel - skip routing and set rating curves to zero if (tot_stor(jrch)%flo < 1.e-6) then ch_rcurv(jrch)%in1 = rcz ch_rcurv(jrch)%out1 = rcz sd_ch(jrch)%in1_vol = 0. sd_ch(jrch)%out1_vol = 0. else if (bsn_cc%rte == 1) then !! Muskingum flood routing method outflo = sd_ch(jrch)%msk%c1 * inflo + sd_ch(jrch)%msk%c2 * sd_ch(jrch)%in1_vol + & sd_ch(jrch)%msk%c3 * sd_ch(jrch)%out1_vol outflo = Min (outflo, tot_stor(jrch)%flo) outflo = Max (outflo, 0.) !! save inflow/outflow volumes for next time step (and day) for Muskingum sd_ch(jrch)%in1_vol = inflo sd_ch(jrch)%out1_vol = outflo else !! Variable Storage Coefficent method - sc=2*dt/(2*ttime+dt) - ttime=(in2+out1)/2 scoef = 2. * dthr / (ch_rcurv(jrch)%in2%ttime + ch_rcurv(jrch)%out1%ttime + dthr) scoef = Min (scoef, 1.) outflo = scoef * tot_stor(jrch)%flo end if !! compute outflow rating curve for next time step outflo_rate = outflo / dts !convert to cms call rcurv_interp_flo (jrch, outflo_rate) ch_rcurv(jrch)%out2 = rcurv !! add outflow to daily hydrograph and subdaily flow rto = outflo / tot_stor(jrch)%flo rto = Min (1., rto) ht2 = ht2 + rto * tot_stor(jrch) ob(icmd)%hyd_flo(1,irtstep) = ob(icmd)%hyd_flo(1,irtstep) + outflo !! subtract outflow from total storage tot_stor(jrch) = (1. - rto) * tot_stor(jrch) !! set rating curve for next time step ch_rcurv(jrch)%in1 = ch_rcurv(jrch)%in2 ch_rcurv(jrch)%out1 = ch_rcurv(jrch)%out2 !! partition channel and flood plain based on bankfull volume if (tot_stor(jrch)%flo > ch_rcurv(jrch)%elev(2)%vol_ch) then !! fill channel to bank full if below rto = (tot_stor(jrch)%flo - ch_rcurv(jrch)%elev(2)%vol_ch) / tot_stor(jrch)%flo fp_stor(jrch) = rto * tot_stor(jrch) ch_stor(jrch) = (1. - rto) * tot_stor(jrch) else ch_stor(jrch) = tot_stor(jrch) fp_stor(jrch) = hz end if !! if flood plain link - fill wetlands to emergency if flood plain storage available do ihru = 1, sd_ch(jrch)%fp%hru_tot iihru = sd_ch(jrch)%fp%hru(ihru) ires= hru(iihru)%dbs%surf_stor !! wetland storage can't go above emergency in release dtbl - it becomes flood plain storage if (ires > 0 .and. fp_stor(jrch)%flo > 0.) then if (fp_stor(jrch)%flo > (wet_ob(iihru)%evol - wet(iihru)%flo)) then rto = (wet_ob(iihru)%evol - wet(iihru)%flo) / fp_stor(jrch)%flo if (rto > 1.e-6) then wet(iihru) = wet(iihru) + rto * fp_stor(jrch) hru(iihru)%wet_obank_in = (rto * fp_stor(jrch)%flo) / (10. * hru(iihru)%area_ha) rto1 = 1. - rto fp_stor(jrch) = rto1 * fp_stor(jrch) end if else wet(iihru) = wet(iihru) + fp_stor(jrch) hru(iihru)%wet_obank_in = fp_stor(jrch)%flo / (10. * hru(iihru)%area_ha) fp_stor(jrch) = hz end if end if end do tot_stor(jrch) = ch_stor(jrch) + fp_stor(jrch) end if ! tot_stor(jrch)%flo < 1.e-6 end do ! end of sub-daily loop !! compute water balance - precip, evap and seep !! km * m * 1000 m/km * ha/10000 m2 = ha ch_wat_d(ich)%area_ha = sd_ch(ich)%chl * sd_ch(ich)%chw / 10. !! m3 = 10. * mm * ha ch_wat_d(ich)%precip = 10. * wst(iwst)%weat%precip * ch_wat_d(ich)%area_ha ch_wat_d(ich)%evap = 10. * bsn_prm%evrch * wst(iwst)%weat%pet * ch_wat_d(ich)%area_ha ch_wat_d(ich)%seep = 10. * sd_ch(ich)%chk * ch_wat_d(ich)%area_ha !k units to mm/d !! add precip ht2%flo = ht2%flo + ch_wat_d(ich)%precip !! calculate transmission losses if (ht2%flo > 1.e-6) then !! mm/hr * km * m * 24 / nsteps = m3 trans_loss = sd_ch(jrch)%chk * sd_ch(jrch)%chl * rcurv%wet_perim * 24. trans_loss = Min(trans_loss, tot_stor(jrch)%flo) !! subtract transmission loses from outflow if (ht2%flo > trans_loss) then rto = trans_loss / ht2%flo ht2 = (1. - rto) * ht2 ob(icmd)%hyd_flo(1,:) = (1. - rto) * ob(icmd)%hyd_flo(1,:) else ht2 = hz ob(icmd)%hyd_flo(1,:) = 0. end if end if !! calculate evaporation losses if (ht2%flo > 1.e-6) then !! calculate width of channel at water level - flood plain evap calculated in wetlands if (dep_flo <= sd_ch(jrch)%chd) then topw = ch_rcurv(jrch)%out2%surf_area else topw = 1000. * sd_ch(jrch)%chl * sd_ch(jrch)%chw end if iwst = ob(icmd)%wst !! mm/day * m2 / (1000. * sd_ch(jrch)%msk%nsteps) evap = bsn_prm%evrch * wst(iwst)%weat%pet * topw / 1000. if (evap < 0.) evap = 0. if (ht2%flo > evap) then if (outflo > 1.e-6) then rto = evap / outflo rto = Min (1., rto) else rto = 0. end if ht2 = (1. - rto) * ht2 ob(icmd)%hyd_flo(1,:) = (1. - rto) * ob(icmd)%hyd_flo(1,:) else ht2 = hz ob(icmd)%hyd_flo(1,:) = 0. end if end if !! check water balance at end of day sum_outflo = ht2%flo inout = sum_inflo - sum_outflo - trans_loss - evap !! total wetland volume at end of day wet_stor(jrch) = hz do ihru = 1, sd_ch(jrch)%fp%hru_tot iihru = sd_ch(jrch)%fp%hru(ihru) wet_stor(jrch) = wet_stor(jrch) + wet(iihru) end do del_stor = (ch_stor(jrch)%flo - ch_stor_init) + (fp_stor(jrch)%flo - fp_stor_init) + & (wet_stor(jrch)%flo - wet_stor_init) ch_fp_wb(jrch)%inflo = sum_inflo ch_fp_wb(jrch)%outflo = sum_outflo ch_fp_wb(jrch)%tl = trans_loss ch_fp_wb(jrch)%ev = evap ch_fp_wb(jrch)%ch_stor_init = ch_stor_init ch_fp_wb(jrch)%ch_stor = ch_stor(jrch)%flo ch_fp_wb(jrch)%fp_stor_init = fp_stor_init ch_fp_wb(jrch)%fp_stor = fp_stor(jrch)%flo ch_fp_wb(jrch)%tot_stor_init = tot_stor_init ch_fp_wb(jrch)%tot_stor = tot_stor(jrch)%flo ch_fp_wb(jrch)%wet_stor_init = wet_stor_init ch_fp_wb(jrch)%wet_stor = wet_stor(jrch)%flo return end subroutine ch_rtmusk