*********************************************************; * Update master DATA1 file for modifications *; * - log-e transform selected data *; * - drop selected data *; * - define new terms *; * - set Columbia River boundary to annual load * and shut off other catchment loadings (12/19/08, DRW)*; * - added code to calculate reciprocal hydraulic load values * new MRB_RF1 reservoirs (12/31/08). * - added code to make suggested corrections from review of residuals. * - corrected routine for replacing missing values for log-transformed variables (2/6/09) *********************************************************; *MODIFIED 07/06/09: Dan Wise. Inactivated code that removes reaches w/o hydseq number. *MODIFIED 07/08/09: Dan Wise. Corrected errors in setting boundary conditions, added east-west versions of some variables; *MODIFIED 07/14/09: Dan Wise. Converted permeability and natural P to SI units; *MODIFIED 07/20/09: Dan Wise. Broke out catchment characteristics by location and elevation; *MODIFIED 07/23/09: Dan Wise. Calculated additional variables; *MODIFIED xx/xx/xx: Dan Wise. Final RF1 version. Corrected errors in point-source loads. Converted avg_permeability to units of cm/hr; LIBNAME DIR_D1 'E:\SPARROW\RF1_report\nawqa_archiving\makedata_programs'; data DIR_D1.data1mod; set DIR_D1.data1; * Add other variables not included in first phase of data c *; LIBNAME newvar 'E:\SPARROW\ancillary_data\new_variables'; data newvar.data1mod2; merge DIR_D1.data1mod newvar.newvar newvar.east_west newvar.dem_slp newvar.nhd_perennial_total; /*File desctriptions newvar: variables representing red alder distribution, natural phosphorus in soil, and irrigated agriculture. east_west: indicator of whether catchment is located on the west of east side of the Cascade Range divide. dem_slp: estimates of minimum, mean,and maximum slope based on digitial elevation models for each RF1 catchment. nhd_perennial_total: summaries of National Hydrographic Dataset (NHD) stream reaches for each RF1 catchment. */ by mrb_id; *Add labels for new variables and provide data sources; Label alder_m2 = 'basal area of alder forest (square meters)' /*Data source: Red Alder Basal Area, by Stream Reach, for the Pacific Northwest. Available from: http://water.usgs.gov/GIS/metadata/usgswrd/XML/redAlderSources_MRB7.xml.*/ natural_p_ppm = 'background concentration of phosphorus (ppm)' /*Data source: Natural Phosphorus Sources for the Pacific Northwest. Available from: http://water.usgs.gov/GIS/metadata/usgswrd/XML/naturalPhosphorusSources_MRB7.xml.*/ east_west = 'indicator of catchment location: east(0) and west (1)' /*Data soure: catchments were assigned to east or west side of Cascade Range divide.*/ nlcd_forest_km2 = 'area of forested land (square kilometers)' /*Data source: Previously defined variable nlcd_forest converted from square meters to square kilometers*/ west_nlcd_forest_km2 = 'area of west-side forested land (square kilometers)' /*Data source: nlcd_forest_km2 located on the west side of the Cascades drainage divide.*/ east_nlcd_forest_km2 = 'area of east-side forested land (square kilometers)' /*Data source: nlcd_forest_km2 located on the east side of the Cascades drainage divide.*/ nlcd_devel_km2 = 'area of developed land (square kilometers)' /*Data source: Previously defined variable nlcd_devel converted from square meters to square kilometers.*/ nlcd_range_km2 = 'area of range land (square kilometers)' /*Data source: Previously defined variable nlcd_range converted from square meters to square kilometers.*/ nlcd_agric_km2 = 'area of agricultural land (square kilometers)' /*Data source: Previously defined variable nlcd_agric converted from square meters to square kilometers.*/ nlcd_other_km2 = 'area of other land (square kilometers)' /*Data source: Previously defined variable nlcd_other converted from square meters to square kilometers.*/ irrag_km2 = 'area of irrigated agriculture (squared kilometers)' /*Data source: Data set compiled by MRB7 and has not yet been reviewed and approved. Data set can be obtained by contacting: Daniel Wise, Hydrologist USGS Oregon Water Science Center 2130 SW 5th Avenue Portland, OR 97201 503-252-3213 dawise@usgs.gov*/ dryag_km2 = 'area of dryland agriculture (squared kilometers)' /*Data source: Area of nlcd_agric_km2 located east of Cascades drainage divide, but excluding area of irrag_km2.*/ westag_km2 = 'area of westside agriculture (squared kilometers)' /*Data source: Area of nlcd_agric_km2 located west of Cascades drainage divide.*/ dem_slp_min = 'minimum slope value from DEM (percent)' dem_slp_max = 'maximum slope value from DEM (percent)' dem_slp_mean = 'mean slope value from DEM (percent)' /*Data source: Values obtained using GRID functions in ARC Workstation on 30 meter elevation grids. The results were summarized for each RF1 catchment.*/ nhd_per_length = 'total length of perennial NHD streams (meters)' nhd_per_dens = 'drainage density for perennial NHD streams (meters/squared kilometer)'; /*Data source: NHD data set. Perennial NHD streams were defined as all reaches described as "perennial", "intermittent", or "artifical path" where the mean average Q was > 1.0 cfs, in addition to all reaches described as canal/ditch.*/ *Calculate additonal variables based on exisitng data; nhd_per_dens = nhd_per_length/demiarea; if nhd_per_length = . then nhd_per_length = 0; if nhd_per_dens = . then nhd_per_dens = 0; nlcd_forest_km2 = nlcd_forest/1000000; nlcd_devel_km2 = nlcd_devel/1000000; nlcd_range_km2 = nlcd_range/1000000; nlcd_agric_km2 = nlcd_agric/1000000; nlcd_other_km2 = nlcd_other/1000000; percent_forest = 100*nlcd_forest_km2/demiarea; perc_west_forest = east_west * percent_forest; if perc_west_forest = . then perc_west_forest = 0; perc_east_forest = -(east_west-1) * percent_forest; if perc_east_forest = . then perc_east_forest = 0; percent_devel = 100*nlcd_devel_km2/demiarea; percent_agric = 100*nlcd_agric_km2/demiarea; percent_range = 100*nlcd_range_km2/demiarea; percent_other = 100*nlcd_other_km2/demiarea; if(percent_forest<50 and percent_devel<50 and percent_range<50 and percent_agric<50 and percent_other<50) then nlcd_mixed_km2 = demiarea; if nlcd_mixed_km2 = . then nlcd_mixed_km2 = 0; if(huc>17100100 and huc<17100300) then coast_forest_km2 = nlcd_forest_km2; if coast_forest_km2 = . then coast_forest_km2 = 0; west_nlcd_forest_km2 = east_west * nlcd_forest_km2; east_nlcd_forest_km2 = -(east_west-1) * nlcd_forest_km2; if(huc<17100100 or huc>17100300) then non_coast_west_forest_km2 = west_nlcd_forest_km2; if non_coast_west_forest_km2 = . then non_coast_west_forest_km2 = 0; west_nlcd_devel_km2 = east_west * nlcd_devel_km2; east_nlcd_devel_km2 = -(east_west-1) * nlcd_devel_km2; west_nlcd_range_km2 = east_west * nlcd_range_km2; east_nlcd_range_km2 = -(east_west-1) * nlcd_range_km2; west_atdep_load = east_west * atdep_load; east_atdep_load = -(east_west-1) * atdep_load; west_ps_n_02_all = east_west * ps_n_02_all; east_ps_n_02_all = -(east_west-1) * ps_n_02_all; west_ps_p_02_all = east_west * ps_p_02_all; east_ps_p_02_all = -(east_west-1) * ps_p_02_all; west_nlcd_devel_km2 = east_west * nlcd_devel_km2; east_nlcd_devel_km2 = -(east_west-1) * nlcd_devel_km2; west_nlcd_agric_km2 = east_west * nlcd_agric_km2; east_nlcd_agric_km2 = -(east_west-1) * nlcd_agric_km2; west_dem_slp_mean = east_west * dem_slp_mean; east_dem_slp_mean = -(east_west-1) * dem_slp_mean; *specify mean slope for coastal drainages*; if (HUC > 17100100 and HUC < 17100300) then coast_slope = dem_slp_mean; if coast_slope = . then coast_slope = 0; *specify mean slope for noncoastal drainages*; if (HUC < 17100100 or HUC > 17100300) then noncoast_slope = dem_slp_mean; if noncoast_slope = . then noncoast_slope = 0; *create elevation variable for high-elevation catchments *; if elevcm_mean > 1000 then high_elev = elevcm_mean; if high_elev = . then high_elev = 0; *create slope variable for high-elevation catchments *; if elevcm_mean > 1000 then high_slope = dem_slp_mean; if high_slope = . then high_slope = 0; *create slope variable for high-elevation west-side catchments *; if (elevcm_mean > 1000 and east_west = 1) then whigh_slope = dem_slp_mean; if whigh_slope = . then whigh_slope = 0; *create slope variable for high-elevation east-side catchments *; if (elevcm_mean > 1000 and east_west = 0) then ehigh_slope = dem_slp_mean; if ehigh_slope = . then ehigh_slope = 0; *create slope variable for low-elevation catchments *; if elevcm_mean <= 1000 then low_slope = dem_slp_mean; if low_slope = . then low_slope = 0; *create variable indicating high-elevation catchments *; if elevcm_mean > 1000 then highcatch = 1; if highcatch = . then highcatch = 0; *create variable indicating high-elevation forested catchments *; if (elevcm_mean > 1000 and percent_forest > 90) then highforest = 1; if highforest = . then highforest = 0; *create variable indicating low-elevation forested catchments *; if elevcm_mean < 1000 then low_perc_forest = percent_forest ; if low_perc_forest = . then low_perc_forest = 0; *create variable indicating high-elevation forested catchments *; if elevcm_mean > 1000 then high_perc_forest = percent_forest ; if high_perc_forest = . then high_perc_forest = 0; *create variables indicating high-elevation range catchments *; if (elevcm_mean > 1000 and percent_range > 90) then highrange = 1; if highrange = . then highrange = 0; if (elevcm_mean > 1000) then highrange_km2 = nlcd_range_km2; if highrange_km2 = . then highrange_km2 = 0; if (elevcm_mean > 1000) then perc_highrange = percent_range; if percent_range = . then percent_range = 0; *convert natural_p_percent to natural_p_ppm; natural_p_ppm = natural_p_percent*10000; *convert average soil permeability from inches/hr to cm/sec; soil_PERMAVE_SI = PERMAVE/1417.3; *convert average soil permeability from inches/hr to cm/hr; soil_PERMAVE_cm_hr = PERMAVE*2.54; /* Convert natural_p_ppm to a loading (2008-10-30, HMJ) */ natural_p_load = natural_p_ppm * demiarea; west_natural_p_load = east_west * natural_p_load; east_natural_p_load = -(east_west-1) * natural_p_load; label natural_p_load = 'background phosphorus concentration (ppm) x incremental area (squared km)'; run; data newvar.data1mod3; set newvar.data1mod2; *delete data for reaches that do not transport load*; *if hydseq = . then delete; run; /* ACCOUNT FOR LOAD FROM UNMONITORED CANADIAN DRAINAGES */; data newvar.data1mod4; set newvar.data1mod3; /*Correct fadrain_percent values for catchments 56740 and 56741*/ if mrb_id = 56740 then fadrain_percent = 0.90; if mrb_id = 56741 then fadrain_percent = 0.99; can50 = .; canbasin = .; fadrain_scaler = 1.0; if (fadrain_percent > 0.0 and fadrain_percent <=0.50) then fadrain_scaler = 1/(1-fadrain_percent); /*Extrapolate area-dependent load variables into unmeasured Canadian areas of catchments (<= 50% of catchment)*/; /* Don't extrapolate for atmospheric N and natural P loads - already scaled by catchment area*/ /* Don't extrapolate delivery variables - should be independent of catchment area*/ ps_p_02_all = ps_p_02_all * fadrain_scaler; west_ps_p_02_all = west_ps_p_02_all * fadrain_scaler; east_ps_p_02_all = east_ps_p_02_all * fadrain_scaler; manc_p = manc_p * fadrain_scaler; manu_p = manu_p * fadrain_scaler; total_manure_p = total_manure_p * fadrain_scaler; farm_p = farm_p * fadrain_scaler; natural_p_load = natural_p_load * fadrain_scaler; nofarm_p = nofarm_p * fadrain_scaler; ps_n_02_all = ps_n_02_all * fadrain_scaler; west_ps_n_02_all = west_ps_n_02_all * fadrain_scaler; east_ps_n_02_all = east_ps_n_02_all * fadrain_scaler; manc_n = manc_n * fadrain_scaler; manu_n = manu_n * fadrain_scaler; total_manure_n = total_manure_n * fadrain_scaler; farm_n = farm_n * fadrain_scaler; nofarm_n = nofarm_n * fadrain_scaler; alder_m2 = alder_m2 * fadrain_scaler; alder_km2 = alder_km2 * fadrain_scaler; nlcd_forest_km2 = nlcd_forest_km2 * fadrain_scaler; west_nlcd_forest_km2 = west_nlcd_forest_km2 * fadrain_scaler; east_nlcd_forest_km2 = east_nlcd_forest_km2 * fadrain_scaler; nlcd_devel_km2 = nlcd_devel_km2 * fadrain_scaler; west_nlcd_devel_km2 = west_nlcd_devel_km2 * fadrain_scaler; east_nlcd_devel_km2 = east_nlcd_devel_km2 * fadrain_scaler; /*Set boundary load variable for unmonitored Canadian drainages(> 50% of catchment)*/; if fadrain_percent > 0.50 then can50 = 1.0; else can50 = 0; canbasin = can50 * demiarea; /*Set landscape loadings for Canadian drainages(> 50% of catchment) to zero*/; if fadrain_percent > 0.50 then ps_p_02_all = 0; if fadrain_percent > 0.50 then west_ps_p_02_all = 0; if fadrain_percent > 0.50 then east_ps_p_02_all = 0; if fadrain_percent > 0.50 then manc_p = 0; if fadrain_percent > 0.50 then manu_p = 0; if fadrain_percent > 0.50 then total_manure_p = 0; if fadrain_percent > 0.50 then pfarm_p = 0; if fadrain_percent > 0.50 then pnofarm_p = 0; if fadrain_percent > 0.50 then natural_p_load = 0; if fadrain_percent > 0.50 then west_natural_p_load = 0; if fadrain_percent > 0.50 then east_natural_p_load = 0; if fadrain_percent > 0.50 then atdep_load = 0; if fadrain_percent > 0.50 then west_atdep_load = 0; if fadrain_percent > 0.50 then east_atdep_load = 0; if fadrain_percent > 0.50 then ps_n_02_all = 0; if fadrain_percent > 0.50 then west_ps_n_02_all = 0; if fadrain_percent > 0.50 then east_ps_n_02_all = 0; if fadrain_percent > 0.50 then manc_n = 0; if fadrain_percent > 0.50 then manu_n = 0; if fadrain_percent > 0.50 then total_manure_n = 0; if fadrain_percent > 0.50 then farm_n = 0; if fadrain_percent > 0.50 then nofarm_n = 0; if fadrain_percent > 0.50 then alder_m2 = 0; if fadrain_percent > 0.50 then alder_km2 = 0; if fadrain_percent > 0.50 then nlcd_forest_km2 = 0; if fadrain_percent > 0.50 then west_nlcd_forest_km2 = 0; if fadrain_percent > 0.50 then east_nlcd_forest_km2 = 0; if fadrain_percent > 0.50 then nlcd_devel_km2 = 0; if fadrain_percent > 0.50 then west_nlcd_devel_km2 = 0; if fadrain_percent > 0.50 then east_nlcd_devel_km2 = 0; run; * COMPUTATIONS AND TRANSFORMATION OF NEW VARIABLES: ROUND 1 *; data newvar.data1mod5; set newvar.data1mod4; /* Rename variable that captures the NHD+ drainage density (2008-11-11, HMJ) */ drainden_nhd = .; drainden_nhd = offrf1_density_m_km2; /* Area-normalize: irrag_km2, offrf1_res_area_km2, offrf1_res_n (2008-10-28, HMJ)*/ irrag_pct = .; offrf1_res_sa_den = .; offrf1_res_n_den = .; irrag_pct = irrag_km2 / demiarea * 100; /* percent */ offrf1_res_sa_den = offrf1_res_area_km2 / demiarea; /* km2 per km2 */ offrf1_res_n_den = offrf1_res_n / demiarea; /* number reservoirs per km2 */ /*Calculate effective precipitation by accounting for irrigated agriculture (2009-06-30, DRW)*/ /* Assume all irrigated areas receive 36 inches (914.4 mm) of irrigation water each year*/ /* Multiply irrigation water depth by 100 to be consistent with precipitation scale*/ irrig_water = irrag_pct*91440*0.01; /*Estimate effective precipitation*/; precip_eff = ppt30mean + irrig_water; prec_eff_cm = precip_eff/1000; /* Initialize variables for farm nutrient loadings and calculate westside farm nutrient loadings */ if east_west = 1 then do; farm_n_west = .; farm_p_west = .; farm_n_west = farm_n + manc_n; farm_p_west = farm_p + manc_p; end; else if east_west = 0 then do; irrag_frac = .; farm_n_dry = .; farm_p_dry = .; farm_n_irr = .; farm_p_irr = .; irrag_frac = irrag_km2 / nlcd01ag_km2; dryag_frac = dryag_km2 / demiarea; end; /*Replace missing values for irrag_frac with 0*/ if irrag_frac = . then irrag_frac = 0; if dryag_frac = . then dryag_frac = 0; /* Disaggregate farm-applied N & P to different irrigation types */ if east_west = 1 then do; irrag_frac = .; end; else if east_west = 0 then do; farm_n_dry = (1 - irrag_frac) * ( farm_n + manc_n ); farm_p_dry = (1 - irrag_frac) * ( farm_p + manc_p ); farm_n_irr = irrag_frac * ( farm_n + manc_n ); farm_p_irr = irrag_frac * ( farm_p + manc_p ); end; if farm_n_west = . then farm_n_west = 0; if farm_n_irr = . then farm_n_irr = 0; if farm_n_dry = . then farm_n_dry = 0; if farm_p_west = . then farm_p_west = 0; if farm_p_irr = . then farm_p_irr = 0; if farm_p_dry = . then farm_p_dry = 0; /* Calculate total farm nutrient loading */ total_farm_n = farm_n + manc_n; total_farm_p = farm_p + manc_p; /* Calculate total agricultural nutrient loading */ total_ag_n = farm_n + manc_n + manu_n; total_ag_p = farm_p + manc_p + manu_p; west_total_ag_n = east_west * total_ag_n; east_total_ag_n = -(east_west-1) * total_ag_n; west_total_ag_p = east_west * total_ag_p; east_total_ag_p = -(east_west-1) * total_ag_p; /* Natural-log transformations (June/July 2008) */ ln_prec_eff_cm = log(prec_eff_cm); ln_soil_permave = log(soil_PERMAVE_cm_hr); run; *REPLACE MISSING VALUES FOR LOG-TRANSFORMED VARIABLES WITH MEAN VALUES (except for irrigation water)*; Proc means data=newvar.data1mod5 noprint; var ln_prec_eff_cm ln_soil_permave ; output out=mean mean= mean_ln_prec_eff_cm mean_ln_soil_permave ; run; data newvar.data1mod5; set newvar.data1mod5; do point=1 to nobs; set mean point=point nobs=nobs; output; end; data newvar.data1mod5a; set newvar.data1mod5; if ln_prec_eff_cm = . then ln_precip_eff_cm = mean_ln_prec_eff_cm; if ln_soil_permave = . then ln_soil_permave = mean_ln_soil_permave; drop mean_ln_soil_permave mean_ln_prec_eff_cm _TYPE_ _FREQ_; run; *OVERWRITE RF1 FRAC VALUES WITH VALUES BASED ON RF1 FLOW ACCUMULATION* /*Read in data flow accumulation ratios */; data newvar.data1mod5b; merge newvar.data1mod5a newvar.qratios; by mrb_id; *remove records for reaches not on modeling network; *if hydseq = . then delete; /*Reset frac values to flow accumulation ratios if ratio is less than 1.0 */ if ratio < 1.0 then frac = ratio; /*Retain original RF1 frac values for braided streams*/ if mrb_id = 55899 then frac = 0.001; if mrb_id = 55898 then frac = 0.999; if mrb_id = 56060 then frac = 0.002; if mrb_id = 56062 then frac = 0.998; if mrb_id = 54155 then frac = 0.005; if mrb_id = 54157 then frac = 0.995; if mrb_id = 55917 then frac = 0.030; if mrb_id = 55915 then frac = 0.070; if mrb_id = 51807 then frac = 0.037; if mrb_id = 51810 then frac = 0.963; if mrb_id = 53851 then frac = 0.436; if mrb_id = 53850 then frac = 0.564; /*Retain original RF1 frac values that represent stream diversions*/ if mrb_id = 56658 then frac = 0.637; if mrb_id = 55133 then frac = 0.693; if mrb_id = 56685 then frac = 0.719; if mrb_id = 53416 then frac = 0.807; if mrb_id = 92036 then frac = 0.807; if mrb_id = 52631 then frac = 0.809; if mrb_id = 52633 then frac = 0.834; if mrb_id = 53282 then frac = 0.882; if mrb_id = 91884 then frac = 0.916; if mrb_id = 48988 then frac = 0.942; if mrb_id = 52619 then frac = 0.948; if mrb_id = 65466 then frac = 0.963; if mrb_id = 56613 then frac = 0.984; *MAKE CORRECTIONS FROM REVIEWS OF MODELING RESULTS /*Correct point-source loadings /*Transfer loads for Portland and Gresham WWTP's from Columbia Slough to Columbia River*/; if mrb_id = 54156 then PS_N_02_all = 0; if mrb_id = 54156 then PS_P_02_all = 0; if mrb_id = 92085 then PS_N_02_all = PS_N_02_all + 2599820; if mrb_id = 92085 then PS_P_02_all = PS_P_02_all + 295819; */Add load from Yakima WWTP (not included in original point-source estimates, use average of 1992 and 1997); if mrb_id = 65465 then PS_N_02_all = PS_N_02_all + 224705; if mrb_id = 65465 then PS_P_02_all = PS_P_02_all + 25743; */Add load from Driggs WWTP; if mrb_id = 49048 then PS_N_02_all = PS_N_02_all + 6526; if mrb_id = 49048 then PS_P_02_all = PS_P_02_all + 1296; */Add load from St. Anthony WWTP; if mrb_id = 48977 then PS_N_02_all = PS_N_02_all + 8797; if mrb_id = 48977 then PS_P_02_all = PS_P_02_all + 1758; */Add load from Horseshoe Bend WWTP; if mrb_id = 50799 then PS_N_02_all = PS_N_02_all + 2084; if mrb_id = 50799 then PS_P_02_all = PS_P_02_all + 394; */Add load from Rexberg WWTP; if mrb_id = 65591 then PS_N_02_all = PS_N_02_all + 43635; if mrb_id = 65591 then PS_P_02_all = PS_P_02_all + 7931; */Add load from McCall fish hatchery; if mrb_id = 91952 then PS_N_02_all = PS_N_02_all + 12572; if mrb_id = 91952 then PS_P_02_all = PS_P_02_all + 1548; */Add load from Stanley fish hatchery; if mrb_id = 51575 then PS_N_02_all = PS_N_02_all + 12518; if mrb_id = 51575 then PS_P_02_all = PS_P_02_all + 1542; */Add load from Lake County/Charlo WWTP; if mrb_id = 47671 then PS_N_02_all = PS_N_02_all + 2487; if mrb_id = 47671 then PS_P_02_all = PS_P_02_all + 543; */Add load from Brookings WWTP; if mrb_id = 80232 then PS_N_02_all = PS_N_02_all + 13587; if mrb_id = 80232 then PS_P_02_all = PS_P_02_all + 2434; */Add load from Pasco WWTP; if mrb_id = 48793 then PS_N_02_all = PS_N_02_all + 34876; if mrb_id = 48793 then PS_P_02_all = PS_P_02_all + 6491; */Add load from Jackson WWTP; if mrb_id = 48888 then PS_N_02_all = PS_N_02_all + 23964; if mrb_id = 48888 then PS_P_02_all = PS_P_02_all + 5011; */Add load from GT Acquisitions LLC (Ski Resort) WWTP; if mrb_id = 49049 then PS_N_02_all = PS_N_02_all + 25279; if mrb_id = 49049 then PS_P_02_all = PS_P_02_all + 3203; */Add load from Thayne WWTP; if mrb_id = 48922 then PS_N_02_all = PS_N_02_all + 3734; if mrb_id = 48922 then PS_P_02_all = PS_P_02_all + 649; */Add load from Lynnwood WWTP; if mrb_id = 80330 then PS_N_02_all = PS_N_02_all + 81852; if mrb_id = 80330 then PS_P_02_all = PS_P_02_all + 13738; */Add load from Lebanon WWTP; if mrb_id = 53946 then PS_N_02_all = PS_N_02_all + 8448; if mrb_id = 53946 then PS_P_02_all = PS_P_02_all + 1940; */Add load from Astoria WWTP; if mrb_id = 53661 then PS_N_02_all = PS_N_02_all + 63532; if mrb_id = 53661 then PS_P_02_all = PS_P_02_all + 8100; */Add load from Grants Pass WWTP; if mrb_id = 91902 then PS_N_02_all = PS_N_02_all + 76857; if mrb_id = 91902 then PS_P_02_all = PS_P_02_all + 14546; */Add load from Chehalis WWTP; if mrb_id = 55523 then PS_N_02_all = PS_N_02_all + 28388; if mrb_id = 55523 then PS_P_02_all = PS_P_02_all + 4947; */Add load from Cowlitz WWTP; if mrb_id = 53516 then PS_N_02_all = PS_N_02_all + 28642; if mrb_id = 53516 then PS_P_02_all = PS_P_02_all + 2864; */Add load from Port Townsend Paper Company; if mrb_id = 80267 then PS_N_02_all = PS_N_02_all + 29953; if mrb_id = 80267 then PS_P_02_all = PS_P_02_all + 75804; */Add TN load from Snokist (Yakima) processing plant; if mrb_id = 65465 then PS_N_02_all = PS_N_02_all + 4324; */Add load from Dworshak fish hatchery; if mrb_id = 52527 then PS_N_02_all = PS_N_02_all + 40602; if mrb_id = 52527 then PS_P_02_all = PS_P_02_all + 5000; */Add load from Rapid River fish hatchery; if mrb_id = 52158 then PS_N_02_all = PS_N_02_all + 12756; if mrb_id = 52158 then PS_P_02_all = PS_P_02_all + 1571; */Add load from Ellis fish hatchery; if mrb_id = 51643 then PS_N_02_all = PS_N_02_all + 8434; if mrb_id = 51643 then PS_P_02_all = PS_P_02_all + 1039; */Add load from Wyoming Game and Fish Commission fish hatchery; if mrb_id = 48927 then PS_N_02_all = PS_N_02_all + 3694; if mrb_id = 48927 then PS_P_02_all = PS_P_02_all + 455; */Remove load from Transalta Centralia Mining Co. (no nutrient load in 2002); if mrb_id = 55550 then PS_N_02_all = 0; if mrb_id = 55550 then PS_P_02_all = 0; /*Change reach type flags*/; if mrb_id = 48111 then rchtype = 0; if mrb_id = 52248 then rchtype = 2; /*Calculate new reservoir properties*/ if mrb_id = 52248 then hload = 236; if mrb_id = 52248 then rhload = 1/hload; /*Change frac values based on USGS gage data*/ if mrb_id = 50797 then frac = 0.89; if mrb_id = 50492 then frac = 0.49; if mrb_id = 92078 then frac = 0.25; if mrb_id = 54248 then frac = 0.85; if mrb_id = 92168 then frac = 0.77; if mrb_id = 56115 then frac = 0.36; if mrb_id = 49110 then frac = 0.78; if mrb_id = 49006 then frac = 0.38; if mrb_id = 49014 then frac = 0.37; if mrb_id = 51777 then frac = 0.78; if mrb_id = 56159 then frac = 0.22; if mrb_id = 91869 then frac = 0.48; if mrb_id = 91874 then frac = 0.38; if mrb_id = 53134 then frac = 0.27; if mrb_id = 48947 then frac = 0.975; if mrb_id = 48946 then frac = 0.975; if mrb_id = 48945 then frac = 0.975; if mrb_id = 48943 then frac = 0.975; if mrb_id = 48941 then frac = 0.975; if mrb_id = 56685 then frac = 0.72; if mrb_id = 48977 then frac = 0.80; if mrb_id = 91953 then frac = 0.78; if mrb_id = 56598 then frac = 0.53; /* changed 1/26/09 based on comparison of USBR gage RDBW and USGS 12484500 */ if mrb_id = 57840 then frac = 0.77; /* added for 2/2/09 version (Malad River system)*/ if mrb_id = 49601 then frac = 0.40; /* added for 2/2/09 version (Malad River system)*/ if mrb_id = 49706 then frac = 0.40; /* added for 2/2/09 version (Malad River system)*/ /*Create new delivery variable that accounts for gain in reach load due to underestimated load from unmapped tributaries and ground water upwelling. Variable is non-zero for reaches where the ratio of the cumulative streamflow flowing into the reach is more 1.25 the streamflow for the reach (based on Dave Wolock's mean annual streamflow estimates), otherwise it is equal to ratio*/; if ratio > 1.25 then ln_qgain = log(ratio); else ln_qgain = -10; /*Prevent load from being transfered from canals: done in SPARROW control file*/ /*Change load estimates*/ *remove sites located too close to other sites*; if station_id = '3701054' then LOAD_A_00600 = .; if station_id = '3701054' then LOAD_A_00665 = .; if station_id = '3701002' then LOAD_A_00600 = .; if station_id = '3701002' then LOAD_A_00665 = .; if station_id = '3701116' then LOAD_A_00600 = .; if station_id = '3701116' then LOAD_A_00665 = .; if station_id = 'A438' then LOAD_A_00600 = .; if station_id = 'A438' then LOAD_A_00665 = .; if station_id = '12473520' then LOAD_A_00600 = .; if station_id = '12473520' then LOAD_A_00665 = .; if station_id = '402112' then LOAD_A_00600 = .; if station_id = '402112' then LOAD_A_00665 = .; if station_id = '404201' then LOAD_A_00600 = .; if station_id = '404201' then LOAD_A_00665 = .; if station_id = 'SKA_2' then LOAD_A_00600 = .; if station_id = 'SKA_2' then LOAD_A_00665 = .; if station_id = '12391950' then LOAD_A_00600 = .; if station_id = '12391950' then LOAD_A_00665 = .; if station_id = '12424500' then LOAD_A_00600 = .; if station_id = '12424500' then LOAD_A_00665 = .; if station_id = '12472600' then LOAD_A_00600 = .; if station_id = '12472600' then LOAD_A_00665 = .; if station_id = '402020' then LOAD_A_00600 = .; if station_id = '402020' then LOAD_A_00665 = .; if station_id = '402012' then LOAD_A_00600 = .; if station_id = '402012' then LOAD_A_00665 = .; if station_id = '402000' then LOAD_A_00600 = .; if station_id = '402000' then LOAD_A_00665 = .; if station_id = '10433' then LOAD_A_00600 = .; if station_id = '10433' then LOAD_A_00665 = .; if station_id = '402075' then LOAD_A_00600 = .; if station_id = '402075' then LOAD_A_00665 = .; if station_id = '404168' then LOAD_A_00600 = .; if station_id = '404168' then LOAD_A_00665 = .; *remove sites with inadequate WQ record*; if station_id = '26B100' then LOAD_A_00600 = .; if station_id = '26B100' then LOAD_A_00665 = .; *change loads based on updated streamflow or wq record*; if station_id = '13135000' then LOAD_A_00600 = 6765826; if station_id = '13135000' then LOAD_A_00665 = 567787; if station_id = '12027' then LOAD_A_00600 = 1510534; if station_id = '12027' then LOAD_A_00665 = 205836; /* Create new load variable for international boundary reaches /* Set value for reach MRB_ID = 65462 to mean annual load for Columbia River at Northport (12400520), /* MRB_ID 57366 to Okanogan River at Oroville to 49A190, and MRB_ID 48478 to Similkamen R. at Oroville (49B070) /* Set values for all other reaches to zero /* Set values for catchment loadings for reach MRB_ID_65462, MRB_ID_57366, and MRB_ID 92219 to zero /* Need to shut off load transfer from upstream reaches in SPARROW control file*/; boundary_n = 0; boundary_p = 0; if mrb_id = 65462 then boundary_n = load_A_00600; if mrb_id = 65462 then boundary_p = load_A_00665; if mrb_id = 65462 then atdep_load = 0; if mrb_id = 65462 then east_atdep_load = 0; if mrb_id = 65462 then ps_n_02_all = 0; if mrb_id = 65462 then east_ps_n_02_all = 0; if mrb_id = 65462 then septic_cnt = 0; if mrb_id = 65462 then manu_n = 0; if mrb_id = 65462 then total_ag_n = 0; if mrb_id = 65462 then east_total_ag_n = 0; if mrb_id = 65462 then farm_n_west = 0; if mrb_id = 65462 then farm_n_irr = 0; if mrb_id = 65462 then farm_n_dry = 0; if mrb_id = 65462 then nlcd_devel_km2 = 0; if mrb_id = 65462 then east_nlcd_devel_km2 = 0; if mrb_id = 65462 then nlcd_forest_km2 = 0; if mrb_id = 65462 then east_nlcd_forest_km2 = 0; if mrb_id = 65462 then alder_km2 = 0; if mrb_id = 65462 then ps_p_02_all = 0; if mrb_id = 65462 then east_ps_p_02_all = 0; if mrb_id = 65462 then manu_p = 0; if mrb_id = 65462 then total_ag_p = 0; if mrb_id = 65462 then east_total_ag_p = 0; if mrb_id = 65462 then farm_p_west = 0; if mrb_id = 65462 then farm_p_irr = 0; if mrb_id = 65462 then farm_p_dry = 0; if mrb_id = 65462 then natural_p_load = 0; if mrb_id = 65462 then east_natural_p_load = 0; if mrb_id = 65462 then west_natural_p_load = 0; if mrb_id = 57366 then boundary_n = load_A_00600; if mrb_id = 57366 then boundary_p = load_A_00665; if mrb_id = 57366 then atdep_load = 0; if mrb_id = 57366 then east_atdep_load = 0; if mrb_id = 57366 then ps_n_02_all = 0; if mrb_id = 57366 then east_ps_n_02_all = 0; if mrb_id = 57366 then septic_cnt = 0; if mrb_id = 57366 then manu_n = 0; if mrb_id = 57366 then total_ag_n = 0; if mrb_id = 57366 then east_total_ag_n = 0; if mrb_id = 57366 then farm_n_west = 0; if mrb_id = 57366 then farm_n_irr = 0; if mrb_id = 57366 then farm_n_dry = 0; if mrb_id = 57366 then nlcd_devel_km2 = 0; if mrb_id = 57366 then east_nlcd_devel_km2 = 0; if mrb_id = 57366 then nlcd_forest_km2 = 0; if mrb_id = 57366 then east_nlcd_forest_km2 = 0; if mrb_id = 57366 then alder_km2 = 0; if mrb_id = 57366 then ps_p_02_all = 0; if mrb_id = 57366 then east_ps_p_02_all = 0; if mrb_id = 57366 then manu_p = 0; if mrb_id = 57366 then total_ag_p = 0; if mrb_id = 57366 then east_total_ag_p = 0; if mrb_id = 57366 then farm_p_west = 0; if mrb_id = 57366 then farm_p_irr = 0; if mrb_id = 57366 then farm_p_dry = 0; if mrb_id = 57366 then natural_p_load = 0; if mrb_id = 57366 then east_natural_p_load = 0; if mrb_id = 57366 then west_natural_p_load = 0; if mrb_id = 92219 then boundary_n = load_A_00600; if mrb_id = 92219 then boundary_p = load_A_00665; if mrb_id = 92219 then atdep_load = 0; if mrb_id = 92219 then east_atdep_load = 0; if mrb_id = 92219 then ps_n_02_all = 0; if mrb_id = 92219 then east_ps_n_02_all = 0; if mrb_id = 92219 then septic_cnt = 0; if mrb_id = 92219 then manu_n = 0; if mrb_id = 92219 then total_ag_n = 0; if mrb_id = 92219 then east_total_ag_n = 0; if mrb_id = 92219 then farm_n_west = 0; if mrb_id = 92219 then farm_n_irr = 0; if mrb_id = 92219 then farm_n_dry = 0; if mrb_id = 92219 then nlcd_devel_km2 = 0; if mrb_id = 92219 then east_nlcd_devel_km2 = 0; if mrb_id = 92219 then nlcd_forest_km2 = 0; if mrb_id = 92219 then east_forest_km2 = 0; if mrb_id = 92219 then alder_km2 = 0; if mrb_id = 92219 then ps_p_02_all = 0; if mrb_id = 92219 then east_ps_p_02_all = 0; if mrb_id = 92219 then manu_p = 0; if mrb_id = 92219 then total_ag_p = 0; if mrb_id = 92219 then east_total_ag_p = 0; if mrb_id = 92219 then farm_p_west = 0; if mrb_id = 92219 then farm_p_irr = 0; if mrb_id = 92219 then farm_p_dry = 0; if mrb_id = 92219 then natural_p_load = 0; if mrb_id = 92219 then east_natural_p_load = 0; if mrb_id = 92219 then west_natural_p_load = 0; /* Create new load variable for load from ground water inflow*/ gwinflow_n = 0; gwinflow_p = 0; /* Create a modified canopy cover variable */ if east_west = 1 then do; CNPY_WEST = CNPY_MEAN; end; /* Area-normalize HLR16 */ HLR_HUMTPSI1_PCT = HLR_16 / (DEMIAREA * 1000000); /* Area-normalize HLR15 */ HLR_SAMTISPB_PCT = HLR_15 / (DEMIAREA * 1000000); /* Area-normalize HLR20 */ HLR_HUMTPSI2_PCT = HLR_20 / (DEMIAREA * 1000000); /* Combined HLR16 + HLR20 */ hlr_humtpsib_pct = HLR_HUMTPSI1_PCT + HLR_HUMTPSI2_PCT; /* Area-normalize HLR10*/ HLR_ARPTISPB_PCT = HLR_10 / (DEMIAREA * 1000000); run; * DROP UNUSED VARIABLES *; data newvar.data1mod6; set newvar.data1mod5b; *drop xxxx; run; *delete reaches outside of modeling domain and store final data in data1mod file*; data dir_d1.data1mod; set newvar.data1mod6; if mrb_id = 0 then delete; if mrb_id = 3 then delete; proc sort; by mrb_id; run;