subroutine sd_channel_control3 use sd_channel_module use channel_velocity_module use basin_module use hydrograph_module use constituent_mass_module use conditional_module use channel_data_module use channel_module use ch_pesticide_module use climate_module use water_body_module use time_module use water_allocation_module, only: wallo use ch_salt_module !rtb salt use ch_cs_module !rtb cs use gwflow_module, only: flood_freq !rtb gwflow use ch_pesticide_module !!! nbs added 7-20-23 use channel_velocity_module implicit none !real :: rcharea !m^2 |cross-sectional area of flow real :: flo_rt !m^3/s |flow rate in reach for day integer :: isd_db ! | integer :: iob ! | integer :: idb !none |channel data pointer integer :: ihyd ! | integer :: ipest ! | integer :: isalt ! |salt ion counter (rtb salt) integer :: ihru ! | integer :: iru ! | integer :: ise ! | integer :: ielem ! | integer :: id integer :: iter real :: ebtm_m !m |erosion of bottom of channel real :: ebank_m !m |meander cut on one side real :: erode_bank_cut !cm |widening caused by downcutting (both sides) real :: ebtm_t !tons |bottom erosion real :: ebank_t !tons |bank erosion real :: sedout !mg |sediment out of waterway channel real :: washld !tons |wash load real :: bedld !tons |bed load real :: dep !tons |deposition real :: hc_sed !tons |headcut erosion real :: chside !none |change in horizontal distance per unit ! |change in vertical distance on channel side ! |slopes; always set to 2 (slope=1/2) real :: a !m^2 |cross-sectional area of channel real :: b !m |bottom width of channel real :: c !none |inverse of channel side slope real :: p !m |wetting perimeter real :: rh !m |hydraulic radius real :: qman !m^3/s or m/s |flow rate or flow velocity real :: frac !0-1 |fraction of hydrograph real :: valint ! | integer :: ivalint ! | real :: tbase !none |flow duration (fraction of 24 hr) real :: tb_pr ! | real :: tb ! | real :: vol_ovb ! | real :: const ! | integer :: ics !none |counter real :: ob_const ! | integer :: ii !none |counter real :: sum_vol ! | real :: xx ! | integer :: ic ! | real :: vol_overmx ! | real :: flood_dep ! | real :: dep_e ! | real :: rto !none |cloud cover factor real :: sumtime ! | real :: vc !m/s |flow velocity in reach real :: pr_ratio ! | real :: shear_btm_cr ! | real :: shear_btm ! | real :: hc !m/yr |head cut advance integer :: max ! | integer :: iaq integer :: iaq_ch real :: det !hr |time step real :: scoef !none |Storage coefficient real :: flo_ls real :: vel, cohes, vel_cr, b_coef, qcms, veg real :: rad_curv, cla, pk_rto, vel_bend, vel_rch real :: arc_len, hyd_radius, prot_len real :: gw_salt_in !kg |salt loading to channel from aquifer real :: gw_cs_in !kg |constituent loading to channel from aquifer real :: seep_mass !kg |salt mass in seepage water real :: salt_conc(8) !kg |salt concentration in channel water real :: cs_conc(8) !kg |constituent concentration in channel water real :: bf_flow !m3/s |bankfull flow rate * adjustment factor ich = isdch isd_db = sd_dat(ich)%hyd iwst = ob(icmd)%wst !rtb floodplain if(bsn_cc%gwflow.eq.1) flood_freq(ich) = 0 !! set ht1 to incoming hydrograph ht1 = ob(icmd)%hin !! set outgoing flow and sediment - ht2 ht2 = hz !! add water transfer if (ob(icmd)%trans%flo > 1.e-6) then ht1 = ht1 + ob(icmd)%trans ob(icmd)%trans = hz end if !set constituents to incoming loads (rtb salt; rtb cs) if (cs_db%num_tot > 0) then hcs1 = obcs(icmd)%hin(1) endif chsd_d(ich)%flo_in = ht1%flo / 86400. !flow for morphology output ch_in_d(ich) = ht1 !set inflow om hydrograph ch_in_d(ich)%flo = ht1%flo / 86400. !flow for om output !rtb hydrograph separation hdsep1%flo_surq = ob(icmd)%hdsep_in%flo_surq hdsep1%flo_latq = ob(icmd)%hdsep_in%flo_latq hdsep1%flo_gwsw = ob(icmd)%hdsep_in%flo_gwsw hdsep1%flo_swgw = ob(icmd)%hdsep_in%flo_swgw hdsep1%flo_satex = ob(icmd)%hdsep_in%flo_satex hdsep1%flo_satexsw = ob(icmd)%hdsep_in%flo_satexsw hdsep1%flo_tile = ob(icmd)%hdsep_in%flo_tile !rtb hydrograph separation !! adjust precip and temperature for elevation using lapse rates w = wst(iwst)%weat if (bsn_cc%lapse == 1) call cli_lapse wst(iwst)%weat = w ht1%temp = 5.0 + 0.75 * wst(iwst)%weat%tave wtemp = 5.0 + 0.75 * wst(iwst)%weat%tave !! if connected to aquifer - add flow if (sd_ch(ich)%aqu_link > 0) then iaq = sd_ch(ich)%aqu_link iaq_ch = sd_ch(ich)%aqu_link_ch if (aq_ch(iaq)%ch(iaq_ch)%flo_fr > 0.) then chsd_d(ich)%aqu_in = (aq_ch(iaq)%ch(iaq_ch)%flo_fr * aq_ch(iaq)%hd%flo) / 86400. chsd_d(ich)%aqu_in_mm = (aq_ch(iaq)%ch(iaq_ch)%flo_fr * aq_ch(iaq)%hd%flo) / (10. * ob(icmd)%area_ha) ht1 = ht1 + aq_ch(iaq)%ch(iaq_ch)%flo_fr * aq_ch(iaq)%hd !rtb salt do isalt=1,cs_db%num_salts gw_salt_in = aq_ch(iaq)%ch(iaq_ch)%flo_fr * aq_chcs(iaq)%hd(1)%salt(isalt) !kg chsalt_d(ich)%salt(isalt)%gw_in = gw_salt_in !kg hcs1%salt(isalt) = hcs1%salt(isalt) + gw_salt_in !kg enddo !rtb cs do ics=1,cs_db%num_cs gw_cs_in = aq_ch(iaq)%ch(iaq_ch)%flo_fr * aq_chcs(iaq)%hd(1)%cs(ics) !kg chcs_d(ich)%cs(ics)%gw_in = gw_cs_in !kg hcs1%cs(ics) = hcs1%cs(ics) + gw_cs_in !kg enddo aq_ch(iaq)%ch(iaq_ch)%flo_fr = 0. end if end if !if gwflow is active, calulate aquifer interactions (ht1 is updated) if(bsn_cc%gwflow.eq.1) then call gwflow_gwsw(ich) !channel <--> groundwater call gwflow_canl(ich) !channel --> canal seepage call gwflow_tile(ich) !groundwater --> channel call gwflow_satx(ich) !groundwater --> channel endif !! set inflow hyds for printing chsd_d(ich)%flo_in = ht1%flo / 86400. !flow for morphology output - m3/s chsd_d(ich)%flo_in_mm = ht1%flo / (10. * ob(icmd)%area_ha) !flow in mm ch_in_d(ich) = ht1 !set inflow om hydrograph ch_in_d(ich)%flo = ht1%flo / 86400. !flow for om output - m3/s !set constituents (rtb salt) to incoming loads if (cs_db%num_tot > 0) then hcs1 = obcs(icmd)%hin(1) end if !! zero outgoing flow and sediment - ht2 ht2 = hz !call Muskingum and variable storage coefficient flood routing method call ch_rtmusk call sd_channel_sediment3 !! route constituents call ch_rtpest !! call mike winchell's new routine for pesticide routing ! call ch_rtpest2 call ch_rtpath !salt and constituent concentrations (g/m3) for inflow water if(cs_db%num_salts > 0 .or. cs_db%num_cs > 0) then hcs2 = hcs1 !set outflow to inflow do isalt=1,cs_db%num_salts if(ht2%flo > 0) then salt_conc(isalt) = (hcs2%salt(isalt) * 1000.) / ht2%flo !g/m3 = mg/L else salt_conc(isalt) = 0. endif enddo do ics=1,cs_db%num_cs if(ht2%flo > 0) then cs_conc(ics) = (hcs2%cs(ics) * 1000.) / ht2%flo !g/m3 = mg/L else cs_conc(ics) = 0. endif enddo endif !salt mass in seepage do isalt=1,cs_db%num_salts seep_mass = salt_conc(isalt) * ch_wat_d(ich)%seep !g/m3 * m3 = g seep_mass = seep_mass / 1000. !kg if(seep_mass > hcs2%salt(isalt)) then seep_mass = hcs2%salt(isalt) endif hcs2%salt(isalt) = hcs2%salt(isalt) - seep_mass !kg chsalt_d(ich)%salt(isalt)%seep = seep_mass !kg (channel salt output) enddo !constituent mass in seepage do ics=1,cs_db%num_cs seep_mass = cs_conc(ics) * ch_wat_d(ich)%seep !g/m3 * m3 = g seep_mass = seep_mass / 1000. !kg if(seep_mass > hcs2%cs(ics)) then seep_mass = hcs2%cs(ics) endif hcs2%cs(ics) = hcs2%cs(ics) - seep_mass !kg chcs_d(ich)%cs(ics)%seep = seep_mass !kg (channel constituent output) enddo !! subtract evaporation if (ht2%flo < ch_wat_d(ich)%evap) then ch_wat_d(ich)%evap = ht2%flo ht2%flo = 0. else ht2%flo = ht2%flo - ch_wat_d(ich)%evap end if !! total outgoing to output to SWIFT ob(icmd)%hout_tot = ob(icmd)%hout_tot + ht2 !compute stream temperature ! Call Subroutune for Ficklin Model, Linear Equation Model, Energy Balance Model !ht2%temp = "output from subroutine" !salt, constituent mass if(cs_db%num_salts > 0 .or. cs_db%num_cs > 0) then hcs3 = hcs2 + ch_water(ich) !incoming + storage endif !rtb hydrograph separation if (rttime > det) then ! ht1 = incoming + storage !! travel time > timestep -- then all incoming is stored and frac of stored is routed hdsep2%flo_surq = scoef * ch_stor_hdsep(ich)%flo_surq hdsep2%flo_latq = scoef * ch_stor_hdsep(ich)%flo_latq hdsep2%flo_gwsw = scoef * ch_stor_hdsep(ich)%flo_gwsw hdsep2%flo_swgw = scoef * ch_stor_hdsep(ich)%flo_swgw hdsep2%flo_satex = scoef * ch_stor_hdsep(ich)%flo_satex hdsep2%flo_satexsw = scoef * ch_stor_hdsep(ich)%flo_satexsw hdsep2%flo_tile = scoef * ch_stor_hdsep(ich)%flo_tile ch_stor_hdsep(ich)%flo_surq = (frac*ch_stor_hdsep(ich)%flo_surq) + hdsep1%flo_surq ch_stor_hdsep(ich)%flo_latq = (frac*ch_stor_hdsep(ich)%flo_latq) + hdsep1%flo_latq ch_stor_hdsep(ich)%flo_gwsw = (frac*ch_stor_hdsep(ich)%flo_gwsw) + hdsep1%flo_gwsw ch_stor_hdsep(ich)%flo_swgw = (frac*ch_stor_hdsep(ich)%flo_swgw) + hdsep1%flo_swgw ch_stor_hdsep(ich)%flo_satex = (frac*ch_stor_hdsep(ich)%flo_satex) + hdsep1%flo_satex ch_stor_hdsep(ich)%flo_satexsw = (frac*ch_stor_hdsep(ich)%flo_satexsw) + hdsep1%flo_satexsw ch_stor_hdsep(ich)%flo_tile = (frac*ch_stor_hdsep(ich)%flo_tile) + hdsep1%flo_tile else !! travel time < timestep -- route all stored and frac of incoming hdsep2%flo_surq = scoef * hdsep1%flo_surq hdsep2%flo_latq = scoef * hdsep1%flo_latq hdsep2%flo_gwsw = scoef * hdsep1%flo_gwsw hdsep2%flo_swgw = scoef * hdsep1%flo_swgw hdsep2%flo_satex = scoef * hdsep1%flo_satex hdsep2%flo_satexsw = scoef * hdsep1%flo_satexsw hdsep2%flo_tile = scoef * hdsep1%flo_tile hdsep2%flo_surq = hdsep2%flo_surq + ch_stor_hdsep(ich)%flo_surq hdsep2%flo_latq = hdsep2%flo_latq + ch_stor_hdsep(ich)%flo_latq hdsep2%flo_gwsw = hdsep2%flo_gwsw + ch_stor_hdsep(ich)%flo_gwsw hdsep2%flo_swgw = hdsep2%flo_swgw + ch_stor_hdsep(ich)%flo_swgw hdsep2%flo_satex = hdsep2%flo_satex + ch_stor_hdsep(ich)%flo_satex hdsep2%flo_satexsw = hdsep2%flo_satexsw + ch_stor_hdsep(ich)%flo_satexsw hdsep2%flo_tile = hdsep2%flo_tile + ch_stor_hdsep(ich)%flo_tile ch_stor_hdsep(ich)%flo_surq = frac * hdsep1%flo_surq ch_stor_hdsep(ich)%flo_latq = frac * hdsep1%flo_latq ch_stor_hdsep(ich)%flo_gwsw = frac * hdsep1%flo_gwsw ch_stor_hdsep(ich)%flo_swgw = frac * hdsep1%flo_swgw ch_stor_hdsep(ich)%flo_satex = frac * hdsep1%flo_satex ch_stor_hdsep(ich)%flo_satexsw = frac * hdsep1%flo_satexsw ch_stor_hdsep(ich)%flo_tile = frac * hdsep1%flo_tile end if ob(icmd)%hdsep%flo_surq = hdsep2%flo_surq ob(icmd)%hdsep%flo_latq = hdsep2%flo_latq ob(icmd)%hdsep%flo_gwsw = hdsep2%flo_gwsw ob(icmd)%hdsep%flo_swgw = hdsep2%flo_swgw ob(icmd)%hdsep%flo_satex = hdsep2%flo_satex ob(icmd)%hdsep%flo_satexsw = hdsep2%flo_satexsw ob(icmd)%hdsep%flo_tile = hdsep2%flo_tile !store outflow components for writing (and convert from m3 --> m3/sec) hyd_sep_array(ich,1) = hdsep2%flo_surq / 86400. hyd_sep_array(ich,2) = hdsep2%flo_latq / 86400. hyd_sep_array(ich,3) = hdsep2%flo_gwsw / 86400. hyd_sep_array(ich,4) = hdsep2%flo_swgw / 86400. hyd_sep_array(ich,5) = hdsep2%flo_satex / 86400. hyd_sep_array(ich,6) = hdsep2%flo_satexsw / 86400. hyd_sep_array(ich,7) = 0. !hdsep2%flo_tile / 86400. !rtb hydrograph separation !end if ich = isdch !! check decision table for flow control - water diversion if (ob(icmd)%ruleset /= "null" .and. ob(icmd)%ruleset /= "0") then id = ob(icmd)%flo_dtbl d_tbl => dtbl_flo(id) call conditions (ich, id) call actions (ich, icmd, id) end if !! check decision table for water allocation if (sd_ch(isdch)%wallo > 0) then call wallo_control (sd_ch(isdch)%wallo) end if !! set outflow hyd to ht2 after diverting water ob(icmd)%hd(1) = ht2 !channel salt updates if(cs_db%num_salts > 0) then do isalt=1,cs_db%num_salts hcs2%salt(isalt) = scoef * hcs3%salt(isalt) ch_water(ich)%salt(isalt) = hcs3%salt(isalt) - hcs2%salt(isalt) enddo endif !channel constituent updates if(cs_db%num_cs > 0) then do ics=1,cs_db%num_cs hcs2%cs(ics) = scoef * hcs3%cs(ics) ch_water(ich)%cs(ics) = hcs3%cs(ics) - hcs2%cs(ics) enddo endif !! calculate stream temperature ob(icmd)%hd(1)%temp = 5. + .75 * wst(iwst)%weat%tave ht2%temp = 5. + .75 * wst(iwst)%weat%tave ch_stor(isdch)%temp = 5. + .75 * wst(iwst)%weat%tave !! set constituents for routing if (cs_db%num_pests > 0) then obcs(icmd)%hd(1)%pest = hcs2%pest end if if (cs_db%num_salts > 0) then !rtb salt obcs(icmd)%hd(1)%salt = hcs2%salt end if if (cs_db%num_cs > 0) then !rtb cs obcs(icmd)%hd(1)%cs = hcs2%cs end if !! output channel organic-mineral ch_out_d(isdch) = ob(icmd)%hd(1) !set outflow om hydrograph ch_out_d(isdch)%flo = ob(icmd)%hd(1)%flo / 86400. !m3 -> m3/s !! output channel morphology chsd_d(isdch)%flo = ob(icmd)%hd(1)%flo / 86400. !adjust if overbank flooding is moved to landscape chsd_d(isdch)%flo_mm = ob(icmd)%hd(1)%flo / (10. * ob(icmd)%area_ha) !flow out in mm chsd_d(isdch)%peakr = peakrate chsd_d(isdch)%sed_in = ob(icmd)%hin%sed chsd_d(isdch)%sed_out = sedout chsd_d(isdch)%sed_stor = ch_stor(isdch)%sed chsd_d(isdch)%washld = washld chsd_d(isdch)%bedld = bedld chsd_d(isdch)%dep = dep chsd_d(isdch)%deg_btm = ebtm_t chsd_d(isdch)%deg_bank = ebank_t chsd_d(isdch)%hc_sed = hc_sed chsd_d(isdch)%width = sd_ch(isdch)%chw chsd_d(isdch)%depth = sd_ch(isdch)%chd chsd_d(isdch)%slope = sd_ch(isdch)%chs chsd_d(isdch)%deg_btm_m = ebtm_m chsd_d(isdch)%deg_bank_m = ebank_m chsd_d(isdch)%n_tot = ob(icmd)%hd(1)%orgn + ob(icmd)%hd(1)%no3 + ob(icmd)%hd(1)%nh3 + ob(icmd)%hd(1)%no2 chsd_d(isdch)%p_tot = ob(icmd)%hd(1)%sedp + ob(icmd)%hd(1)%solp chsd_d(ich)%dep_bf = sd_ch_vel(ich)%dep_bf chsd_d(ich)%velav_bf = sd_ch_vel(ich)%vel_bf !! set pesticide output variables do ipest = 1, cs_db%num_pests chpst_d(isdch)%pest(ipest)%tot_in = obcs(icmd)%hin(1)%pest(ipest) chpst_d(isdch)%pest(ipest)%sol_out = frsol * obcs(icmd)%hd(1)%pest(ipest) chpst_d(isdch)%pest(ipest)%sor_out = frsrb * obcs(icmd)%hd(1)%pest(ipest) chpst_d(isdch)%pest(ipest)%react = chpst%pest(ipest)%react chpst_d(isdch)%pest(ipest)%volat = chpst%pest(ipest)%volat chpst_d(isdch)%pest(ipest)%settle = chpst%pest(ipest)%settle chpst_d(isdch)%pest(ipest)%resus = chpst%pest(ipest)%resus chpst_d(isdch)%pest(ipest)%difus = chpst%pest(ipest)%difus chpst_d(isdch)%pest(ipest)%react_bot = chpst%pest(ipest)%react_bot chpst_d(isdch)%pest(ipest)%bury = chpst%pest(ipest)%bury chpst_d(isdch)%pest(ipest)%water = ch_water(ich)%pest(ipest) chpst_d(isdch)%pest(ipest)%benthic = ch_benthic(ich)%pest(ipest) end do !rtb salt - set salt output variables do isalt = 1, cs_db%num_salts chsalt_d(ich)%salt(isalt)%tot_in = hcs1%salt(isalt) + chsalt_d(ich)%salt(isalt)%gw_in !mass entering the channel (upstream + groundwater) chsalt_d(ich)%salt(isalt)%tot_out = hcs2%salt(isalt) !mass leaving the channel (for current day) chsalt_d(ich)%salt(isalt)%water = ch_water(ich)%salt(isalt) !mass stored in channel (for current day) !concentration of channel water (= concentration of outflow water) if(ht2%flo > 0) then chsalt_d(ich)%salt(isalt)%conc = (hcs2%salt(isalt) * 1000.) / ht2%flo !g/m3 = mg/L else chsalt_d(ich)%salt(isalt)%conc = 0. endif enddo !rtb cs - set constituent output variables do ics = 1, cs_db%num_cs chcs_d(ich)%cs(ics)%tot_in = hcs1%cs(ics) + chcs_d(ich)%cs(ics)%gw_in !mass entering the channel (upstream + groundwater) chcs_d(ich)%cs(ics)%tot_out = hcs2%cs(ics) !kg mass leaving the channel (for current day) chcs_d(ich)%cs(ics)%water = ch_water(ich)%cs(ics) !kg mass stored in channel (for current day) if(ht2%flo > 0) then !concentration of outflow water chcs_d(ich)%cs(ics)%conc = (hcs2%cs(ics) * 1000.) / ht2%flo !g/m3 = mg/L else chcs_d(ich)%cs(ics)%conc = 0. endif enddo !! set values for recharge hydrograph - should be trans losses !ob(icmd)%hd(2)%flo = perc return end subroutine sd_channel_control3