From 3bb22ed72cb2b1bb83233b12f1db0b49ef0ec499 Mon Sep 17 00:00:00 2001 From: Jeff Arnold Date: Tue, 31 Dec 2024 13:11:26 -0600 Subject: [PATCH] souirce code updates --- src/basin_module.f90 | 2 +- src/basin_prm_default.f90 | 2 +- src/cal_conditions.f90 | 5 + src/cal_parm_select.f90 | 13 +- src/cal_parmchg_read.f90 | 2 + src/ero_cfactor.f90 | 13 +- src/mgt_harvbiomass.f90 | 5 +- src/nut_nminrl.f90 | 59 +++++---- src/pl_dormant.f90 | 2 +- src/pl_fert.f90 | 3 +- src/plant_data_module.f90 | 2 +- src/plant_parm_read.f90 | 20 ++- src/smp_grass_wway.f90 | 259 +++++++++++++++++++------------------- 13 files changed, 210 insertions(+), 177 deletions(-) diff --git a/src/basin_module.f90 b/src/basin_module.f90 index 798eb56..cccfd8d 100644 --- a/src/basin_module.f90 +++ b/src/basin_module.f90 @@ -111,7 +111,7 @@ module basin_module real :: tlaps = 6.5 !! deg C/km |temperature lapse rate: deg C per km of elevation difference real :: nfixmx = 20.0 !! max daily n-fixation (kg/ha) real :: decr_min = 0.01 !! minimum daily residue decay - real :: rsd_covco = 0.30 !! residue cover factor for computing frac of cover + real :: rsd_covco = 0.75 !! residue cover factor for computing frac of cover real :: urb_init_abst = 1. !! maximum initial abstraction for urban areas when using Green and Ampt real :: petco_pmpt = 100.0 !! PET adjustment (%) for Penman-Montieth and Preiestly-Taylor methods real :: uhalpha = 1.0 !! alpha coeff for est unit hydrograph using gamma func diff --git a/src/basin_prm_default.f90 b/src/basin_prm_default.f90 index ea04a97..9e6432f 100644 --- a/src/basin_prm_default.f90 +++ b/src/basin_prm_default.f90 @@ -37,7 +37,7 @@ subroutine basin_prm_default !if (bsn_prm%dorm_hr < 1.e-6) bsn_prm%dorm_hr = -1. !! time threshold used to define dormant (hrs) if (bsn_prm%nfixmx < 1.e-6) bsn_prm%nfixmx = 20.0 !! max daily n-fixation (kg/ha) if (bsn_prm%decr_min < 1.e-6) bsn_prm%decr_min = 0.01 !! - if (bsn_prm%rsd_covco < 1.e-6) bsn_prm%rsd_covco = 0.30 !! residue cover factor for computing frac of cover + if (bsn_prm%rsd_covco < 1.e-6) bsn_prm%rsd_covco = 0.75 !! residue cover factor for computing C factor equation if (bsn_prm%urb_init_abst < 1.e-6) bsn_prm%urb_init_abst = 0. !! PET adjustment (%) for Penman-Montieth and Preiestly-Taylor methods if (bsn_prm%petco_pmpt < 0.5 .and. bsn_prm%petco_pmpt > 0.) bsn_prm%petco_pmpt = 0.0 !! reservoir sediment settling coeff bsn_prm%petco_pmpt = (100. + bsn_prm%petco_pmpt) / 100. !! convert to fraction diff --git a/src/cal_conditions.f90 b/src/cal_conditions.f90 index 0e1a9e4..b282f47 100644 --- a/src/cal_conditions.f90 +++ b/src/cal_conditions.f90 @@ -7,6 +7,7 @@ subroutine cal_conditions use hru_module, only : hru use soil_module use plant_module + use plant_data_module use time_module use climate_module, only : pcp, tmp @@ -74,6 +75,10 @@ subroutine cal_conditions if (pl_find == "n") cond_met = "n" exit end do + case ("pl_class") + if (cal_upd(ichg_par)%cond(ic)%targc /= pl_class(ielem)) then + cond_met = "n" + end if case ("landuse") !for hru if (cal_upd(ichg_par)%cond(ic)%targc /= hru(ielem)%land_use_mgt_c) then cond_met = "n" diff --git a/src/cal_parm_select.f90 b/src/cal_parm_select.f90 index 43ac894..5fd6573 100644 --- a/src/cal_parm_select.f90 +++ b/src/cal_parm_select.f90 @@ -34,7 +34,8 @@ subroutine cal_parm_select (ielem, ly, chg_parm, chg_typ, chg_val, absmin, absma use hydrograph_module use pesticide_data_module use plant_module - use gwflow_module + use plant_data_module + use gwflow_module implicit none @@ -80,12 +81,16 @@ subroutine cal_parm_select (ielem, ly, chg_parm, chg_typ, chg_val, absmin, absma hru(ielem)%lumv%usle_p = chg_par (hru(ielem)%lumv%usle_p, & chg_typ, chg_val, absmin, absmax) + case ("usle_c") + pldb(ielem)%usle_c = chg_par (pldb(ielem)%usle_c, & + chg_typ, chg_val, absmin, absmax) + case ("ovn") hru(ielem)%luse%ovn = chg_par (hru(ielem)%luse%ovn, & chg_typ, chg_val, absmin, absmax) case ("elev") - hru(ielem)%topo%elev = chg_par (hru(ielem)%topo%elev, & + hru(ielem)%topo%elev = chg_par (hru(ielem)%topo%elev, & chg_typ, chg_val, absmin, absmax) case ("slope") @@ -429,8 +434,8 @@ subroutine cal_parm_select (ielem, ly, chg_parm, chg_typ, chg_val, absmin, absma bsn_prm%n_updis = chg_par(bsn_prm%n_updis, & chg_typ, chg_val, absmin, absmax) - case ("p_updis") - bsn_prm%p_updis = chg_par(bsn_prm%p_updis, & + case ("rsd_covco") + bsn_prm%rsd_covco = chg_par(bsn_prm%rsd_covco, & chg_typ, chg_val, absmin, absmax) case ("dorm_hr") diff --git a/src/cal_parmchg_read.f90 b/src/cal_parmchg_read.f90 index 4ce5b29..d498b2e 100644 --- a/src/cal_parmchg_read.f90 +++ b/src/cal_parmchg_read.f90 @@ -102,6 +102,8 @@ subroutine cal_parmchg_read cal_upd(i)%num_elem = db_mx%dtbl_res case ("plt") cal_upd(i)%num_elem = sp_ob%hru + case ("pl_class") + cal_upd(i)%num_elem = db_mx%plantparm case ("lyr") cal_upd(i)%num_elem = sp_ob%hru case ("sol") diff --git a/src/ero_cfactor.f90 b/src/ero_cfactor.f90 index 2afc41a..3e81a68 100644 --- a/src/ero_cfactor.f90 +++ b/src/ero_cfactor.f90 @@ -103,23 +103,24 @@ subroutine ero_cfactor !! newer method using residue and biomass cover rsd_sumfac = (soil1(j)%rsd(1)%m +1.) / 1000. - grnd_covfact = 0. + grnd_sumfac = 0. can_covfact = 10000. do ipl = 1, pcom(j)%npl idp = pcom(j)%plcur(ipl)%idplt ab_gr_t = pl_mass(j)%ab_gr(ipl)%m / 1000. - !grnd_sumfac = grnd_sumfac + ab_gr_t - grnd_covfact = grnd_covfact + pldb(idp)%usle_c * ab_gr_t / (ab_gr_t + exp(1.175 - 1.748 * ab_gr_t)) + grnd_sumfac = grnd_sumfac + ab_gr_t + !! grnd_covfact = grnd_covfact + pldb(idp)%usle_c * ab_gr_t / (ab_gr_t + exp(1.175 - 1.748 * ab_gr_t)) can_covfact = amin1 (can_covfact, pcom(j)%plg(ipl)%cht) end do !grnd_covfact = grnd_sumfac / (grnd_sumfac + exp(1.175 - 1.748 * grnd_sumfac)) - rsd_covfact = exp(-0.75 * rsd_sumfac) + rsd_covfact = exp(-bsn_prm%rsd_covco * rsd_sumfac) can_frcov = amin1 (1., pcom(j)%lai_sum / 3.) can_covfact = 1. - can_frcov * Exp(-.328 * pcom(j)%cht_mx) - bio_covfac = 1. - grnd_covfact * exp(-0.1 * can_covfact) - c = Max(1.e-10, rsd_covfact * grnd_covfact * bio_covfac) + grnd_covfact = exp(-pldb(idp)%usle_c * grnd_sumfac) + !! bio_covfac = 1. - grnd_covfact * exp(-0.1 * can_covfact) + c = Max(1.e-10, rsd_covfact * grnd_covfact) ! * can_covfact) !! erosion output variables ero_output(j)%ero_d%c = c diff --git a/src/mgt_harvbiomass.f90 b/src/mgt_harvbiomass.f90 index e08dfc9..f05cfcc 100644 --- a/src/mgt_harvbiomass.f90 +++ b/src/mgt_harvbiomass.f90 @@ -40,11 +40,10 @@ subroutine mgt_harvbiomass (jj, iplant, iharvop) harv_seed = hi_tot * pl_mass(j)%seed(ipl) harv_leaf = hi_tot * pl_mass(j)%leaf(ipl) harv_stem = hi_tot * pl_mass(j)%stem(ipl) - pl_yield = harv_seed + harv_leaf - pl_yield = pl_yield + harv_stem + pl_yield = harv_seed + harv_leaf + harv_stem !! apply pest stress to harvest index - mass lost due to pests - don't add to residue - pl_yield = (1. - pcom(j)%plcur(ipl)%pest_stress) * (1. - harveff) * pl_yield + pl_yield = (1. - pcom(j)%plcur(ipl)%pest_stress) * pl_yield !! add plant carbon for printing hrc_d(j)%plant_surf_c = hrc_d(j)%plant_surf_c + pl_yield%c hpc_d(j)%harv_abgr_c = hpc_d(j)%harv_abgr_c + pl_yield%c diff --git a/src/nut_nminrl.f90 b/src/nut_nminrl.f90 index b56ac90..f57f993 100644 --- a/src/nut_nminrl.f90 +++ b/src/nut_nminrl.f90 @@ -32,7 +32,7 @@ subroutine nut_nminrl use septic_data_module use basin_module use organic_mineral_mass_module - use hru_module, only : rsdco_plcom, i_sep, ihru, ipl, isep + use hru_module, only : rsdco_plcom, i_sep, ihru, isep use soil_module use plant_module use output_landscape_module, only : hnb_d @@ -43,7 +43,7 @@ subroutine nut_nminrl integer :: k = 0 !none |counter (soil layer) integer :: kk = 0 !none |soil layer used to compute soil water and ! |soil temperature factors - integer :: idp = 0 + !integer :: idp = 0 real :: rmn1 = 0. !kg N/ha |amount of nitrogen moving from fresh organic ! |to nitrate(80%) and active organic(20%) ! |pools in layer @@ -62,7 +62,7 @@ subroutine nut_nminrl real :: cprf = 0. ! |carbon phosphorus ratio factor real :: ca = 0. ! | real :: decr = 0. ! | - real :: rdc = 0. ! | + !real :: rdc = 0. ! | real :: wdn = 0. !kg N/ha |amount of nitrogen lost from nitrate pool in ! |layer due to denitrification real :: cdg = 0. !none |soil temperature factor @@ -77,6 +77,8 @@ subroutine nut_nminrl hnb_d(j)%org_lab_p = 0. hnb_d(j)%act_sta_n = 0. hnb_d(j)%denit = 0. + hnb_d(j)%rsd_nitorg_n = 0. + hnb_d(j)%rsd_laborg_p = 0. !! compute humus mineralization of organic soil pools do k = 1, soil(j)%nly @@ -147,33 +149,36 @@ subroutine nut_nminrl cnrf = 1. end if - if (soil1(j)%rsd(k)%p > 1.e-4) then - cpr = soil1(j)%rsd(k)%c / soil1(j)%rsd(k)%p - if (cpr > 5000.) cpr = 5000. - cprf = Exp(-.693 * (cpr - 200.) / 200.) - else - cprf = 1. - end if + if (soil1(j)%rsd(k)%p > 1.e-4) then + cpr = soil1(j)%rsd(k)%c / soil1(j)%rsd(k)%p + if (cpr > 5000.) cpr = 5000. + cprf = Exp(-.693 * (cpr - 200.) / 200.) + else + cprf = 1. + end if - ca = Min(cnrf, cprf, 1.) + ca = Min(cnrf, cprf, 1.) - !! compute root and incorporated residue decomposition - !! all plant residue in soil is mixed - don't track individual plant residue in soil + !! compute root and incorporated residue decomposition + !! all plant residue in soil is mixed - don't track individual plant residue in soil - if (pcom(j)%npl > 0) then - decr = rsdco_plcom(j) / pcom(j)%npl * ca * csf - else - decr = 0.05 - end if - decr = Max(bsn_prm%decr_min, decr) - decr = Min(decr, 1.) - decomp = decr * soil1(j)%rsd(k) - soil1(j)%rsd(k) = soil1(j)%rsd(k) - decomp - soil1(j)%mn(k)%no3 = soil1(j)%mn(k)%no3 + .8 * decomp%n - soil1(j)%hact(k)%n = soil1(j)%hact(k)%n + .2 * decomp%n - soil1(j)%mp(k)%lab = soil1(j)%mp(k)%lab + .8 * decomp%p - soil1(j)%hsta(k)%p = soil1(j)%hsta(k)%p + .2 * decomp%p - + if (pcom(j)%npl > 0) then + decr = rsdco_plcom(j) / pcom(j)%npl * ca * csf + else + decr = 0.05 + end if + decr = Max(bsn_prm%decr_min, decr) + decr = Min(decr, 1.) + decomp = decr * soil1(j)%rsd(k) + soil1(j)%rsd(k) = soil1(j)%rsd(k) - decomp + soil1(j)%mn(k)%no3 = soil1(j)%mn(k)%no3 + .8 * decomp%n + soil1(j)%hact(k)%n = soil1(j)%hact(k)%n + .2 * decomp%n + soil1(j)%mp(k)%lab = soil1(j)%mp(k)%lab + .8 * decomp%p + soil1(j)%hsta(k)%p = soil1(j)%hsta(k)%p + .2 * decomp%p + + hnb_d(j)%rsd_nitorg_n = hnb_d(j)%rsd_nitorg_n + .8 * decomp%n + hnb_d(j)%rsd_laborg_p = hnb_d(j)%rsd_laborg_p + .8 * decomp%p + !! compute denitrification wdn = 0. if (i_sep(j) /= k .or. sep(isep)%opt /= 1) then diff --git a/src/pl_dormant.f90 b/src/pl_dormant.f90 index a887660..d5683ac 100644 --- a/src/pl_dormant.f90 +++ b/src/pl_dormant.f90 @@ -47,7 +47,7 @@ subroutine pl_dormant end if lai_drop = max (0., lai_drop) lai_drop = amin1 (1., lai_drop) - leaf_drop%m = lai_drop * pl_mass(j)%leaf(ipl)%m + leaf_drop%m = rto * lai_drop * pl_mass(j)%leaf(ipl)%m leaf_drop%n = leaf_drop%m * pcom(j)%plm(ipl)%n_fr leaf_drop%n = max (0., leaf_drop%n) leaf_drop%p = leaf_drop%m * pcom(j)%plm(ipl)%p_fr diff --git a/src/pl_fert.f90 b/src/pl_fert.f90 index 38b4838..0d2dbbd 100644 --- a/src/pl_fert.f90 +++ b/src/pl_fert.f90 @@ -19,7 +19,7 @@ subroutine pl_fert (ifrt, frt_kg, fertop) implicit none - real, parameter :: rtof=0.5 !none |weighting factor used to partition the + real :: rtof !none |weighting factor used to partition the ! |organic N & P concentration of septic effluent ! |between the fresh organic and the stable organic pools integer :: j = 0 !none |hru counter @@ -36,6 +36,7 @@ subroutine pl_fert (ifrt, frt_kg, fertop) j = ihru + rtof = 0.5 !! calculate c:n ratio for manure applications for SWAT-C if (bsn_cc%cswat == 2) then org_frt%m = frt_kg diff --git a/src/plant_data_module.f90 b/src/plant_data_module.f90 index 08b8d40..03e26d0 100644 --- a/src/plant_data_module.f90 +++ b/src/plant_data_module.f90 @@ -3,7 +3,7 @@ module plant_data_module implicit none character(len=40), dimension (:), allocatable :: plts_bsn !none |plant names simulated in current run - + character(len=25), dimension(:), allocatable :: pl_class !none |plant class - row crop, tree, grass, etc type plant_db character(len=40) :: plantnm = "" !none |crop name character(len=18) :: typ = "" !none |plant category diff --git a/src/plant_parm_read.f90 b/src/plant_parm_read.f90 index 8d5aa66..0bde64b 100644 --- a/src/plant_parm_read.f90 +++ b/src/plant_parm_read.f90 @@ -6,9 +6,10 @@ subroutine plant_parm_read implicit none - integer :: ic = 0 !none |counter - character (len=80) :: titldum = ""! |title of file - character (len=80) :: header = "" ! |header of file + integer :: ic = 0 !none |plant counter + character (len=80) :: titldum = "" ! |title of file + character (len=80) :: header = "" ! |header of file + character (len=80) :: plclass = "" ! |plant class - row crop, close grown, grass, tree, etc integer :: eof = 0 ! |end of file integer :: imax = 0 !none |determine max number for array (imax) and total number in file integer :: mpl = 0 ! | @@ -23,6 +24,7 @@ subroutine plant_parm_read if (.not. i_exist .or. in_parmdb%plants_plt == " null") then allocate (pldb(0:0)) allocate (plcp(0:0)) + allocate (pl_class(0:0)) else do open (104,file=in_parmdb%plants_plt) @@ -37,6 +39,7 @@ subroutine plant_parm_read end do allocate (pldb(0:imax)) allocate (plcp(0:imax)) + allocate (pl_class(0:imax)) rewind (104) read (104,*,iostat=eof) titldum @@ -44,8 +47,17 @@ subroutine plant_parm_read read (104,*,iostat=eof) header if (eof < 0) exit + read (104,*,iostat=eof) plclass + if (plclass /= "nam1") then + backspace (105) + end if + do ic = 1, imax - read (104,*,iostat=eof) pldb(ic) + if (plclass /= "nam1") then + read (104,*,iostat=eof) pldb(ic) + else + read (104,*,iostat=eof) pldb(ic), pl_class(ic) + end if if (eof < 0) exit pldb(ic)%mat_yrs = Max (1, pldb(ic)%mat_yrs) end do diff --git a/src/smp_grass_wway.f90 b/src/smp_grass_wway.f90 index 930ef0b..4e384a1 100644 --- a/src/smp_grass_wway.f90 +++ b/src/smp_grass_wway.f90 @@ -65,164 +65,167 @@ subroutine smp_grass_wway real :: xrem = 0. ! | integer :: k = 0 !m^3/s |Total number of HRUs plus this HRU number -!! set variables + !! set variables j = ihru - -!! do this only if there is surface runoff this day + !! do this only if there is surface runoff this day if (surfq(j) > 0.001) then -!! compute channel peak rate using SCS triangular unit hydrograph -!! Calculate average flow based on 3 hours of runoff + !! compute channel peak rate using SCS triangular unit hydrograph + !! calculate average flow based on 3 hours of runoff chflow_day = 1000. * surfq(j) * hru(ihru)%km - chflow_m3 = chflow_day/10800 - qp_cms = 2. * chflow_m3 / (1.5 * tc_gwat(j)) + chflow_m3 = chflow_day / 10800 + qp_cms = 2. * chflow_m3 / (1.5 * tc_gwat(j)) -!! if peak rate is greater than bankfull discharge - if (qp_cms > grwway_vel(j)%vel_bf) then - rcharea = grwway_vel(j)%area - rchdep = hru(j)%lumv%grwat_d - else -!! find the crossectional area and depth for todays flow -!! by iteration method at 1cm interval depth -!! find the depth until the discharge rate is equal to volrt - sdti = 0. - rchdep = 0. - - Do While (sdti < qp_cms) - rchdep = rchdep + 0.01 - rcharea = (grwway_vel(j)%wid_btm + 8 * rchdep) * rchdep - p = grwway_vel(j)%wid_btm + 2. * rchdep * Sqrt(1. + 8 * 8) - rh = rcharea / p - sdti = Qman(rcharea, rh, hru(j)%lumv%grwat_n, hru(j)%lumv%grwat_s) - end do - end if - -!! Sediment yield (kg) from fraction of area drained by waterway - - sedin = sedyld(ihru) * hru(ihru)%km -!! Calculate sediment losses in sheetflow at waterway sides - -!! calculate area of sheeflow in m^2 assumne *:1 side slope 8.06 = (8^2+1^2)^.5 - sf_area = (hru(j)%lumv%grwat_d - rchdep) * 8.06 * hru(j)%lumv%grwat_l * 1000 -!! Adjust Area to account for flow nonuniformities White and Arnold 2009 found half of flow in VFS -!!handled by 10% of VFS area. Waterways likely even more concentrated Assume only 20% of sideslope acts as filters - if (sf_area > 1.e-6) then - sf_area = sf_area * 0.20 -!! calculate runoff depth over sheetflow area in mm - sf_depth=surfq(j) * hru(ihru)%km * 1000000/sf_area -!! Calculate sediment load on sheetflow area kg/ha - sf_sed = sedin * 1000 / sf_area -!! Calculate runoff and sediment losses taken from mostly from filter.f - end if - - if (sf_area > 0.) then -!! surq_remove = 75.8 - 10.8 * Log(sf_depth) + 25.9 -!! & * Log(sol_k(1,j)) - !! Simpler form derived from vfsmod simulations. r2 = 0.57 Publication pending white and arnold 2008 - - surq_remove = 95.6 - 10.79 * Log(sf_depth) - if (surq_remove > 100.) surq_remove = 100. - if (surq_remove < 0.) surq_remove = 0. - - sed_remove = 79.0 - 1.04 * sf_sed + 0.213 * surq_remove - if (sed_remove > 100.) sed_remove = 100. - if (sed_remove < 0.) sed_remove = 0. - - Else - sed_remove = 0 - surq_remove = 0 - endif - sedint = sedin * (1. - sed_remove / 100.) - -!! calculate flow velocity - vc = 0.001 - if (rcharea > 1.e-4) then - vc = qp_cms / rcharea - if (vc > grwway_vel(j)%celerity_bf) vc = grwway_vel(j)%celerity_bf - end if - -!! compute deposition in the waterway - cyin = 0. - cych = 0. - depnet = 0. - deg = 0. - dep = 0. -!! if there is significant flow calculate - if (chflow_m3 > 1.e-4) then -!! Calculate sediment concentration in inflow mg/m^3 - cyin = sedint / chflow_day -!! Calculate sediment transport capacity mg/m^3 - cych = hru(j)%lumv%grwat_spcon * vc ** 1.5 -!! Calculate deposition in mg - depnet = chflow_day * (cyin - cych) - if (depnet < 0.) depnet = 0 - if (depnet > sedint) depnet = sedint - endif -!! Calculate sediment out of waterway channel - sedout = sedint - depnet + !! if peak rate is greater than bankfull discharge + if (qp_cms > grwway_vel(j)%vel_bf) then + rcharea = grwway_vel(j)%area + rchdep = hru(j)%lumv%grwat_d + else + !! find the crossectional area and depth for todays flow + !! by iteration method at 1cm interval depth + !! find the depth until the discharge rate is equal to volrt + sdti = 0. + rchdep = 0. + + Do While (sdti < qp_cms) + rchdep = rchdep + 0.01 + rcharea = (grwway_vel(j)%wid_btm + 8 * rchdep) * rchdep + p = grwway_vel(j)%wid_btm + 2. * rchdep * Sqrt(1. + 8 * 8) + rh = rcharea / p + sdti = Qman(rcharea, rh, hru(j)%lumv%grwat_n, hru(j)%lumv%grwat_s) + end do + end if -!! Calculate total fraction of sediment and surface runoff transported - if (sedyld(j) < .0001) sedyld(j) = .0001 - sed_frac = sedout/sedyld(j) + !! Sediment yield (t) from fraction of area drained by waterway + sedin = sedyld(j) + + !! Calculate sediment losses in sheetflow at waterway sides + !! calculate area of sheeflow in m^2 assumne *:1 side slope 8.06 = (8^2+1^2)^.5 + sf_area = (hru(j)%lumv%grwat_d - rchdep) * 8.06 * hru(j)%lumv%grwat_l * 1000. + !! limit area of sheet flow to 10% of hru area + sf_area = Min(0.1 * hru(j)%area_ha * 10000., sf_area) + + !! Adjust Area to account for flow nonuniformities White and Arnold 2009 found half of flow in VFS + !!handled by 10% of VFS area. Waterways likely even more concentrated Assume only 20% of sideslope acts as filters + if (sf_area > 1.e-6) then + sf_area = sf_area * 0.20 + !! calculate runoff depth over sheetflow area in mm + sf_depth = surfq(j) * hru(ihru)%area_ha * 10000. / sf_area + !! Calculate sediment load on sheetflow area kg/m2 + sf_sed = sedin * 1000. / sf_area + !! Calculate runoff and sediment losses taken from mostly from filter.f + end if + if (sf_area > 0.) then + !! surq_remove = 75.8 - 10.8 * Log(sf_depth) + 25.9 * Log(sol_k(1,j)) + !! Simpler form derived from vfsmod simulations. r2 = 0.57 Publication pending White and Arnold 2008 + surq_remove = 95.6 - 10.79 * Log(sf_depth) + if (surq_remove > 95.) surq_remove = 95. + if (surq_remove < 0.) surq_remove = 0. - surq_frac = 1 - surq_remove/100 + sed_remove = 79.0 - 1.04 * sf_sed + 0.213 * surq_remove + if (sed_remove > 95.) sed_remove = 95. + if (sed_remove < 0.) sed_remove = 0. + else + sed_remove = 0 + surq_remove = 0 + endif + sedint = sedin * (1. - sed_remove / 100.) + + !! calculate flow velocity + vc = 0.001 + if (rcharea > 1.e-4) then + vc = qp_cms / rcharea + if (vc > grwway_vel(j)%celerity_bf) vc = grwway_vel(j)%celerity_bf + end if -!! Subtract reductions from sediment, nutrients, bacteria, and pesticides NOT SURFACE RUNOFF to protect water balance - sedtrap = sedyld(j) * (1. - sed_frac) - sedyld(j) = sedyld(j) * sed_frac - sedminpa(j) = sedminpa(j) * sed_frac - sedminps(j) = sedminps(j) * sed_frac - sedorgp(j) = sedorgp(j) * sed_frac - surqsolp(j) = surqsolp(j) * surq_frac - sedorgn(j) = sedorgn(j) * sed_frac - surqno3(j) = surqno3(j) * surq_frac + !! compute deposition in the waterway + cyin = 0. + cych = 0. + depnet = 0. + deg = 0. + dep = 0. + !! if there is significant flow calculate + if (chflow_m3 > 1.e-4) then + !! Calculate sediment concentration in inflow mg/m^3 + cyin = sedint / chflow_day + !! Calculate sediment transport capacity mg/m^3 + cych = hru(j)%lumv%grwat_spcon * vc ** 1.5 + !! Calculate deposition in mg + depnet = chflow_day * (cyin - cych) + if (depnet < 0.) depnet = 0 + if (depnet > sedint) depnet = sedint + endif + + !! Calculate sediment out of waterway channel + sedout = sedint - depnet + + !! Calculate total fraction of sediment and surface runoff transported + if (sedyld(j) < .0001) sedyld(j) = .0001 + sed_frac = sedout / sedyld(j) + if (sed_frac > 75.) sed_frac = 75. + if (sed_frac < 0.) sed_frac = 5. + + surq_frac = 1 - surq_remove / 100. + if (sed_frac > 30.) sed_frac = 30. + if (sed_frac < 0.) sed_frac = 0. + + !! Subtract reductions from sediment, nutrients, bacteria, and pesticides NOT SURFACE RUNOFF to protect water balance + sedtrap = sedyld(j) * (1. - sed_frac) + sedyld(j) = sedyld(j) * sed_frac + sedminpa(j) = sedminpa(j) * sed_frac + sedminps(j) = sedminps(j) * sed_frac + sedorgp(j) = sedorgp(j) * sed_frac + surqsolp(j) = surqsolp(j) * surq_frac + sedorgn(j) = sedorgn(j) * sed_frac + surqno3(j) = surqno3(j) * surq_frac xrem = 0. - if (sedtrap <= lagyld(j)) then - lagyld(j) = lagyld(j) - sedtrap - else - xrem = sedtrap - lagyld(j) - lagyld(j) = 0. - if (xrem <= sanyld(j)) then - sanyld(j) = sanyld(j) - xrem + if (sedtrap <= lagyld(j)) then + lagyld(j) = lagyld(j) - sedtrap else - xrem = xrem - sanyld(j) - sanyld(j) = 0. - if (xrem <= sagyld(j)) then - sagyld(j) = sagyld(j) - xrem + xrem = sedtrap - lagyld(j) + lagyld(j) = 0. + if (xrem <= sanyld(j)) then + sanyld(j) = sanyld(j) - xrem else - xrem = xrem - sagyld(j) - sagyld(j) = 0. - if (xrem <= silyld(j)) then - silyld(j) = silyld(j) - xrem + xrem = xrem - sanyld(j) + sanyld(j) = 0. + if (xrem <= sagyld(j)) then + sagyld(j) = sagyld(j) - xrem else - xrem = xrem - silyld(j) - silyld(j) = 0. - if (xrem <= clayld(j)) then - clayld(j) = clayld(j) - xrem + xrem = xrem - sagyld(j) + sagyld(j) = 0. + if (xrem <= silyld(j)) then + silyld(j) = silyld(j) - xrem else - xrem = xrem - clayld(j) - clayld(j) = 0. + xrem = xrem - silyld(j) + silyld(j) = 0. + if (xrem <= clayld(j)) then + clayld(j) = clayld(j) - xrem + else + xrem = xrem - clayld(j) + clayld(j) = 0. + end if end if end if end if end if - end if sanyld(j) = Max(0., sanyld(j)) silyld(j) = Max(0., silyld(j)) clayld(j) = Max(0., clayld(j)) sagyld(j) = Max(0., sagyld(j)) lagyld(j) = Max(0., lagyld(j)) -!! Calculate pesticide removal -!! based on the sediment and runoff removal only + !! Calculate pesticide removal + !! based on the sediment and runoff removal only do k = 1, cs_db%num_pests hpestb_d(j)%pest(k)%surq = hpestb_d(j)%pest(k)%surq * surq_frac hpestb_d(j)%pest(k)%sed = hpestb_d(j)%pest(k)%sed * (1. - sed_remove / 100.) end do - end if + end if !! surfq(j) > 0.001 + return end subroutine smp_grass_wway \ No newline at end of file