C/////////////////////////////////////////////////////////////////////////// C///// RIVERWEB BACKEND C///// Last updated 1/14/06 C///// Whoever wrote this originally was on crack. I mean, seriously, C///// who would pass that the station is improved and then make the C///// program check to make sure it is improved by also passing an array? C///// Not cool, man, not cool. I commented out a bunch of the code that C///// rechecks to make sure BMP is on, because we're actually going C///// to trust our input. C///// - Kevin McGehee 1/14/06 C///// Having to debug someone else's five year old Fortan code is not fun. C///// Changed some variables to organize files a bit more. C///// -Kevin McGehee 6/9/05 C///// Newly commented 7/19/00 Old trw7 base C///// Kind Codes C///// C///// Nitrogen 1, Phosphorus 2, Sediment 3, Heavy Metals 4, Hydrocarbons 5 C///// Water Temp 6, pH 7, DO 8, BOD 9, Saturated DO 10, Air Temp 11, tro=13 C///// Total Flow 14, Groundwater Stream Input 15, 16 Pesticide C///// velocity 17, depth 18, width 19, Soil Saturation 20 C///// C///// BIO INDICATORS: C///// fishibi 21, Benthic IBI 22, # Fish Species 23, % Fish Species Expected 24 C///// Benthic Taxa num 25, % Benthic Species Expected 26, # EPT Species 27, C///// actual fish count (per 10 m^2) 28 C///// C///// loads: Nitrogen 51, Phosphorous 52, Sediment 53, HM 54, HY 55 C//////////////////////////////////////////////////////////////////////////// program riverweb C/////// Declare the variables and constants ////// C character*50 cid C character*75 iname,oname,aname logical error integer nit, phos, sed, hm, hy, wt, pPh, disO, bod, sdo, atm, tRO integer outflow, grndflow, pest, ck, bio, ioer, discr C//// parameters for kind codes/output parameter (nit=1, phos=2, sed=3, hm=4, hy=5,wt=6, pPh=7) parameter (disO=8, bod=9, SDO=10, atm=11, tRO=13) parameter (outflow=14, grndflow=15, pest=16, velo=17, dep=18) parameter (wid=19, soilsat=20, fishi=21, benthi=22) parameter (numfish=23, fishperc=24,numbenth=25, benthperc=26) parameter (numept=27, countfish=28) C/// Loads (need to be greater than all others. set ck mod at nitl-1) parameter (nitl=51, phosl=52, sedl=53, hml=54, hyl=55) C/// parameter (p86=13,p85=14,p84=15,a86=16,a85=17,a84=18,fnum=19) real airt85(365),airt86(365),airt84(366),AirTmp(365) C//// missing SDO: Uses std. function of water temp. (EffOD still is a curvefunction of C//// BOD as of 7/6) real GNFac(7), GPFac(7),GSdFac(7) real RNFac(7), RPFac(7), RSdFac(7) real RbNFac(5,7), RbPFac(5,7), RbSdFac(5,7) real RbHMFac(5,7),RbHYFac(5,7) real RbBodfc(5,7), rbph(5,7) real GHMFac(7),GHYFac(7) real RHMFac(7),RHYFac(7) real RBODfc(7),GBODfc(7), soilin(7) real rph(7), gph(7), perdev(7), bperdev(7), tperdev(7) C note - air temp determined directly at river mouth - so 8 factors. real airfac(8),yeah real NitCon(365),PhsCon(365),SedCon(365) real HMCon(365), HYCon(365), WTCon(365) real pHCon(365), BODCon(365), tCon(365) real NitLd(365),PhsLd(365), SedLd(365) real fishexpected, fishcount(365), fishibi(365) real fishper(365), fishnum(365) real benthibi(365), benthper(365), benthnum(365) real benthexpected, eptnum(365), holdyrph(365) real holdyrnit(365), wtdfishibi(365), wtdbenthibi(365) real wtdeptnum(365) real HMLd(365), HYLd(365), WTLd(365) real BODLd(365),PScon(365), PSld(365), tLd(365) real temp, totro(365), lag, curco, e, hold, tempron real vel(365), depth(365), width(365), distemp real holdcon integer soil, n, f integer S90(7) real effOD(41) real RO(365),cro(365), actrnfac(365) real GFac(365),RFac(365), gcfac, rcfac growth, Otcng real pcpm86(365), pcpl85(365), pcph84(366), pcp(365), pcpmu integer year integer SAcres(7), notnine real CtoM, CtoG integer strt,which,stp, kind, BMP,order(6),cnt, mainbmp real Gdis(366),inflow(365),tinflow(365),tgrndflow(365) real WTR(365),bodR(365),initcoeff(7) real WTemp(365),ATemp(365),SatDO(365) real bodG(365),pHR(365),pHG(365),pH(365),xDO(365) real templd(365),zyr(365), satur(365) real SChange(7), ROperA,thecrv(41) real cover(7), initgrnd(7), initsat(7,2) integer CN(7,4), BMPCN(7,4) integer currentcn, yes real p, s, gwstor(366), et(365), alpar(7) real alpha, beta(7), gamma, curalp, curgam, shade(7),tempacres real totacres real lagchar(7), slowlag(7), maxu, mini, avgn, summer real pestagri(12), pesturb(12), thepest(12), actipest(365) real pestzero(12), gwlagco(15), actlagco(15), falagco(5) C/// alpha is the percentage of groundwater storage that outflows every day C/// Nitrogen Factors [RO: CBP. GW is CBP for Stations 3,4,5,~6. 0,1, and 2 C/// were slightly altered from CBP to conform to the data from NCSU.] C/// I'm at about 60% confidence that all the factors are approriate. Part C/// of the problem is finding Nitrogen data in a world dominated by Nitrate C/// data. C older:data GNFac/0.4,0.45,7.392,2.23,4,0.4,4/ C old: data GNFac/0.4,0.45,8.4,2.23,4,0.4,4.4/ data GNFac/0.8,1.036,7.392,2.23,4,0.4,4.4/ data RNFac/0.544,0.544,9.326,3.086,5.108,0.5,5.108/ data (RbNFac(1,i),i=1,7)/0.54,0.544,9.326,2.35,5.108,0.5,3.89/ data (RbNFac(2,i),i=1,7)/0.54,0.544,9.326,3.086,5.108,0.5,5.108/ data (RbNFac(3,i),i=1,7)/0.54,0.544,9.326,3.086,5.108,0.5,5.108/ data (RbNFac(4,i),i=1,7)/0.54,0.544,9.326,3.086,5.108,0.5,5.108/ data (RbNFac(5,i),i=1,7)/0.54,0.544,9.326,3.086,5.108,0.5,5.108/ C/// Phosphorous Factors [RO Based on CC, TX, and CBP (avgs)] C/// "Phosphorous has low mobility in GW..."...GW=.01*RO C/// I have reasonably high confidence (>80%) in these numbers. C/// I would consider dropping the Phosphorous to a constant [why C/// variable if not contaminating effect?] in groundwater <=0.1 data GPFac/0.007,0.017,0.13,0.057,0.026,0.002,0.032/ data RPFac/0.07,0.17,1.3,0.57,0.26,0.02,0.32/ data (RbPFac(1,i),i=1,7)/0.07,0.17,1.3,0.24,0.26,0.02,0.16/ data (RbPFac(2,i),i=1,7)/0.07,0.17,1.3,0.57,0.26,0.02,0.32/ data (RbPFac(3,i),i=1,7)/0.07,0.17,1.3,0.57,0.26,0.02,0.32/ data (RbPFac(4,i),i=1,7)/0.07,0.17,1.3,0.57,0.26,0.02,0.32/ data (RbPFac(5,i),i=1,7)/0.07,0.17,1.3,0.57,0.26,0.02,0.32/ C/// Sediment Factors [Based on CBP factors] (not GSd. Untruthable? C/// Probably reasonable, though.) C/// I wonder if these factors should be combined with velocity. I'd like C/// to think so. data GSdFac/2.5,2.5,2.5,2.5,2.5,2.5,2.5/ data RSdFac/5,5,140,140,140,3,140/ data (RbSdFac(1,i),i=1,7)/5,5,140,70,140,3,140/ data (RbSdFac(2,i),i=1,7)/5,5,140,140,140,3,140/ data (RbSdFac(3,i),i=1,7)/5,5,140,140,140,3,140/ data (RbSdFac(4,i),i=1,7)/5,5,140,140,140,3,140/ data (RbSdFac(5,i),i=1,7)/5,5,140,140,140,3,140/ C/// Heavy Metal Factors [based on EMC MD (RO) for 2,3,4,6 (not 1,5)] C/// Reasonableness confidence ? 90% data GHMFac/0,0,0,0,0,0,0/ data RHMFac/0,0.03,0.202,0.111,.217,0.01,.226/ data (RbHMFac(1,i),i=1,7)/0,0.03,0.202,0.051,.217,0.01,.226/ data (RbHMFac(2,i),i=1,7)/0,0.03,0.202,0.111,.217,0.01,.226/ data (RbHMFac(3,i),i=1,7)/0,0.03,0.202,0.111,.217,0.01,.226/ data (RbHMFac(4,i),i=1,7)/0,0.03,0.202,0.111,.217,0.01,.226/ data (RbHMFac(5,i),i=1,7)/0,0.03,0.202,0.111,.217,0.01,.226/ C/// Hydrocarbon Factors (from EMC MD, oil and grease) C/// Reasonableness confidence ? 90% data GHYFac/0,0,0,0,0,0,0/ data RHYFac/0,0,0,1.7,9,0,9/ data (RbHYFac(1,i),i=1,7)/0,0,0,1,6,0,9/ data (RbHYFac(2,i),i=1,7)/0,0,0,1.7,9,0,9/ data (RbHYFac(3,i),i=1,7)/0,0,0,1.7,9,0,9/ data (RbHYFac(4,i),i=1,7)/0,0,0,1.7,9,0,9/ data (RbHYFac(5,i),i=1,7)/0,0,0,1.7,9,0,9/ C/// Air Factors [Totally made up. Then again, not unreasonable (or well C/// modelled)]...Should be calibrated with averages or something. data airfac/.9,.9,1.1,1,1.2,1.1,1.2,1/ C/// BOD factors {????} Currently CBP ...... I dunno. C/// Corpus Christi/20,20, 4,25.5,23,--,14/ C/// Maryland EMC /--,--,13.5,15,12,--,13/ C/// ???????? [notice big differences here???] Confidence <30% C/// Wonder if modelling DO this way will work at all. data RBODfc/58,58,92,91,68,21,92/ data (RbBODfc(1,i),i=1,7)/58,58,92,91,68,21,92/ data (RbBODfc(2,i),i=1,7)/58,58,92,91,68,21,92/ data (RbBODfc(3,i),i=1,7)/58,58,92,91,68,21,92/ data (RbBODfc(4,i),i=1,7)/58,58,92,91,68,21,92/ data (RbBODfc(5,i),i=1,7)/58,58,92,91,68,21,92/ C//// RbBod fac's not currently implemented. How change with BMP? C//// (* of course, this is a semi-universal question, but I've not even C//// bothered playing with this yet... *) data GBODfc/7*20/ C/// pH factors (These are totally made up. Makes sense that Forest>Agri) C/// Trying to make look like a typical watershed, which generally has a C/// fairly well-defined range of pH, caused by it's input, RO speed, and C/// dominant groundwater chemistry. I have 0% confidence in these C/// numbers and don't think this should be modelled this way... data Rph/6.8,6.8,6.3,6.6,6.5,6.2,6.4/ data (Rbph(1,i),i=1,7)/6.9,6.9,6.4,6.7,6.6,6.3,6.5/ data (Rbph(2,i),i=1,7)/6.5,6.5,5.35,6.1,5.9,5.7,5.9/ data (Rbph(3,i),i=1,7)/6.5,6.5,5.35,6.1,5.9,5.7,5.9/ data (Rbph(4,i),i=1,7)/6.5,6.5,5.35,6.1,5.9,5.7,5.9/ data (Rbph(5,i),i=1,7)/6.5,6.5,5.35,6.1,5.9,5.7,5.9/ data Gph/7.0,7.0,6.8,7.0,6.9,6.2,6.9/ C/// Set zero things. data thecrv/41*0.0/ data tLd/365*0.0/ data tinflow/365*0.0/ data zyr/365*0/ data mainbmp/0/ data fishcount/365*0.0/ data wtdfishibi/365*0.0/ data wtdbenthibi/365*0.0/ data wtdeptnum/365*0.0/ data totacres/0.0/ data discr/1/ C/// Acreage in each terrain [soon to be inputted.] C data S90/89139,89139,105387,15950,7672,2429,2357/ data S90/15000,5000,80000,25000,8000,1500,8000/ C/// Total Acreage (s90 2: 147500 acres = 230 square miles) C/// Total Acreage (s90 (original): 312073 acres = 490 square miles) C/// Ground cover which evapotranspires (evaporation and transpiration) data cover/1.35,1.1,1,.75,.55,.85,.55/ C/// Shade (groundcover @ streambed; wetlands I'm fudging as a special case, since the temperature regime there is far C/// different from a flowing stream system. I wonder if the "gw" storage temperature---having REALLY big tvavg pds needs to C/// be treated on its own, with the current method used only for RO. Incidentally, I think right now I'm thinking about this C/// strm temp calculation being some sort of combination of the two. data shade/1,.75,.25,.5,.2,.6,.15/ C/// Development percentages [-43 is artifactual trick to set ibi's of st.0 C/// to 5.0] data tperdev/-43.0,10.0,10.0,40.0,75.0,0.0,90.0/ data bperdev/-43.0,-10.0,0,30.0,60.0,-20.0,75.0/ data perdev/0.0,0.0,0.0,0.0,0.0,0.0,0.0/ C/// Curve Number Data (85%) data (CN(1,i),i=1,4)/30,55,70,77/ data (CN(2,i),i=1,4)/45,66,77,83/ data (CN(3,i),i=1,4)/67,78,85,89/ data (CN(4,i),i=1,4)/54,70,80,85/ data (CN(5,i),i=1,4)/81,88,90,93/ data (CN(6,i),i=1,4)/0,0,0,0/ data (CN(7,i),i=1,4)/74,82,87,89/ C/// (75 % Confidence) data (BMPcn(1,i),i=1,4)/30,55,70,77/ data (BMPcn(2,i),i=1,4)/30,55,70,77/ data (BMPcn(3,i),i=1,4)/62,71,78,81/ data (BMPcn(4,i),i=1,4)/51,68,74,84/ data (BMPcn(5,i),i=1,4)/78,83,86,87/ data (BMPcn(6,i),i=1,4)/0,0,0,0/ data (BMPcn(7,i),i=1,4)/71,79,84,89/ C//// The calibration of the initcoeff [for the GW Storage] depend on C//// a) landuse characteristics: cover(et), alpha(gw); and b) pcp data initcoeff/16,12,15,12,11,16,17/ C//// initsat's are pretty much arbitrary. I wouldn't advise relying C//// to much on the saturation calculation at any rate. data (initsat(i,1),i=1,7)/.2,.2,.2,.2,0,1,.2/ data (initsat(i,2),i=1,7)/.15,.15,.15,.15,0,1,.15/ c//// The "Alpha" array. These are calibrated to achieve the proper c//// RO/GW/ET, FLOW/ET balances for _this (1986)_ precipitation data C//// set. data alpar/.01,.025,.008,.03,.05,.01,.04/ C/// Initial Percentage of PCP-RO that ends up in Gwstorage (currently C/// appears to be behave as kind of as a fudge factor). Also, gws appears C/// to be very sensitive to small changes in beta (for 3, specifically, C/// but presumably for 6, too) while significantly less sensitive to changes C/// if beta is around 1/2. This bears some examination. I believe that C/// some combination of initgrnd and beta will have meaning and that there C/// is a way this can be determined. It is, however, sensible that residential C/// and urban have lower beta, since man-made interference with drainage systems C/// help drive surface evaporation, which will lower this factor. C/// If you buy that explanation, that is. In reality, it seems more C/// theoretically justifialbe to set this to 1 everywhere and let cover C/// RO alp etc drive the water flow through the system. C/// I dunno. I guess this isn't necesary with the arrays of alpha clearly C/// being the controlling factor after much experimentation. Thus, I'll leave C/// these here, set to unity, in case someone decides they want more of the C/// water to mysteriously jump up and leave the system. This does seem most C/// advantageous for the more urbanized watersheds where ET is small. (rel. C/// speaking, of course). C data beta/1,1,1,.5,.5,1,.5/ data beta/1,1,1,1,1,1,1/ C/// Pesticide Data (1-12 arrays for "urban" and "agriculture" both 75 C/// percentile) Good data: 95% Confidence that this is good enough. data pestagri/.4,.3,.2,.4,1.9,2,1.9,.9,.4,.3,.3,.4/ data pesturb/.15,.35,.45,.5,.8,.4,.7,.75,.3,.45,.25,.6/ data pestzero/0,0,0,0,0,0,0,0,0,0,0,0/ data e/2.718281828/ C/// Soil Type data [has not been thoroughly tested as of 7/26/00.] data soilin/3,3,3,3,3,3,3/ C/// Lag Characteristics (smaller numbers, longer lag) C/// lagchar is for RO, gwlagco is for GW. lagchar is station dep, and C/// gwlagco is just an array of coefficients [basically] that holds for C/// all stations baseflow. C/// The numbers are fit to resemble expected flows for the given landuse. C/// The way I see it, there are two dominant effects on RO lag. C/// 1. Size of area being drained. 2. Velocity in stream. Since these C/// are also roughly inverse effects (vel. being dependent on discharge C/// which closely tracks area), the range 0-6 should be relatively small. C/// The wetlands (station 5, strt6) value is completely irrelevant, since C/// there is no RO. Plus, wetlands are basically nonsensical anyway. C/// I'm pretty (70%) confident in these. data lagchar/8,8,8,8,9,8,9/ data slowlag/5,5,5,5,6,5,6/ C data lagchar/6,6,7.5,10,15,3,15/ C/// gwlagco is a slow exponential decay. Sigma i=1-15 = 1. C/// This looks backward in time... C/// y(day)=.108726*(.93590921758*e^(-.066236797x)) R^2=.99998 C/// Remember, these are coefficients. data (gwlagco(i),i=1,5)/.10203, .09545, .08929, .08353, .07814/ data (gwlagco(i),i=6,10)/.07311, .06839, .06398, .05905, .05599/ data (gwlagco(i),i=11,15)/.05238,.049,.04584,.04289,.040929908/ C/// How do factors (1=Nit, 2=Pho, etc) "drain"? If a is big the loss is less. C/// Some ideas about sources, solubility, etc. should inform these numbers. data falagco/80.0,1.0,10.0,10.0,1.0/ C--------------------------- C Old(er) code (still in use) to convert BOD--> Effective oxygen demand. C I (CIF) don't really know how well this is working...but believable values C come out of the model, and for real DO values I think we'd need a total C methodological overhaul anyway. See SWQ book by Chapra for about 10 C chapters on DO. C--------------------------- data effOD /0.00, 0.6, 1.13, 1.80, 2.48, 3.08, 3.90, 4.65, 5.33 $ , 6.23, 7.00,30*0/ C---------------------------------- year = 1986 C/// This is setting the year to use for the precipitation data C/// Changing this at present will reek havoc. C---------------------------- C error=.false. C---------------------------------------------------------- C///// get Precipitation Data //// C//// Again, changing this to get different pcp file //// C//// will cause a lot of problems. bummer, eh? //// C//// If you do change things, look for the alpar //// C//// and initcoeff to need changes. -YOU HAVE BEEN //// C//// WARNED- :) //// C---------------------------------------------------------- C read from med86 file open (13, FILE='fort.13', STATUS='OLD') read (p86,*) (pcpm86(i), i=1,365) C read from lo85 read (p85,*) (pcpl85(i), i=1,365) C read from hi84 read (p84,*) (pcph84(i), i=1,366) pcpmu=3.68 if (year.EQ.1986) then do 77 j=1,365 pcp(j)=pcpm86(j)*pcpmu/3.68 77 continue else if (year.EQ.1985) then do 78 j=1,365 pcp(j)=pcpl85(j)*pcpmu/3.68 78 continue else do 79 j=1,365 pcp(j)=pcph84(j) 79 continue endif C-------------------------------------- C///// get Air Temperature Data ////// C///// (this will not reek havoc.) //// C-------------------------------------- C read from airt86 file read (a86,*) (airt86(i), i=1,365) C read from airt85 read (a85,*) (airt85(i), i=1,365) C read from airt84 read (a84,*) (airt84(i), i=1,366) IF(year.eq.1986)THEN do 177 j=1,365 AirTmp(j)=airt86(j) 177 continue ELSE IF(year.eq.1985)THEN do 178 j=1,365 AirTmp(j)=airt85(j) 178 continue ELSE do 179 j=1, 365 AirTmp(j)=airt84(j) 179 continue endif C--- Read Acreage... ioer=0 open (88, FILE='/tmp/riveracres', STATUS='OLD', IOSTAT=ioer) if (ioer.eq.0) then read (88,*) (SAcres(i), i=1,7) else do 180 i=1,7 SAcres(i)=S90(i) 180 continue end if close (88) C----------------------------------------- C///// set Acres and initial groundwater ////// C----------------------------------------- do 40 i=1,7 initgrnd(i)=initcoeff(i)*SAcres(i) C write (*,*) "Acres: ",Sacres(i), " Gnd: ", initgrnd(i) 40 continue C// We read from STDIN now! -KRM 6/17/05 C-----------Read strt subwatershed--------- read (*,83) which C--------0-7 from web corresponds to 1-8 in fortran program which = which + 1 C-------------Read kind of quality indicator read (*,84) kind C---------add one to kind kind = kind + 1 C-------------Read BMP; BMP = 1 for Best Management Practices read (*,81) BMP C-----------Read list of improvements------- C read (*,82) (order(i),i=1,6) C-----------Read soil Type (IN FUTURE...) C read (19,83) soil C soil=3 C close (19) C-------------Create output file names, which depends on BMP C mainBMP = BMP C if (mainBMP.eq.0) then C oname='tmpgif/riverout' C oname(16:15+j)=cid(1:j) C else C oname='tmpgif/riverbout' C oname(17:16+j)=cid(1:j) C end if C open (20, FILE=oname, STATUS='UNKNOWN') 82 format (6i1) C// Format needs to be changed if array length changes. 83 format (i1) 84 format (i2) 81 format (i1) C-------------Set up conversion factors----------- C Convert_to_Gallons = 1/12{ft/"}*43560{sq ft/acre} C *7.4805{gal/cu. feet} C 1m^3=264(G) (1m^3=1000 L) CtoM=(1/12.0)*43560*7.48056*(1/264.172) CtoG=(1/12.0)*43560*7.48056 if (which.eq.8) then strt=1 stp=7 else strt=which stp=which endif C---------------------------------------------------- C***** If strt = 7, the last station is being C***** analyzed, thus need to add up all 7 C***** stations and go through the loop 7 times C***** Otherwise, for which<8, only execute loop once C***** for that particular station 333 if (strt.gt.stp) then go to 777 else C************************************************ C*** **** C*** Done for all Indicators **** C*** **** C************************************************ C--------------------------------------------------- C------ Calculate runoff for every day of the year C------ Runoff depends on start watershed, C------ whether improvements have been made C-------and the amount of precipitation - pcp(i) C---------------------------------------------------- soil=soilin(strt) gwstor(1)=initgrnd(strt) satur(1)=initsat(strt,bmp+1) temp=0 bio=0 ck=mod(kind,50) if ((kind.ge.21).and.(kind.le.28)) bio=1 C* Not a kludge >>>>> C* alpha is the percentage in GW storage that enters into the stream. C* alpha=alpar(strt) C* <<<<< C* Kludges and schtuff. if ((strt.eq.6).and.(bmp.eq.1)) alpha=0.00075 if ((strt.eq.4).and.(bmp.eq.1)) then gnfac(strt)=1.65 gpfac(strt)=.024 end if if ((strt.eq.7).and.(bmp.eq.1)) then gnfac(strt)=3.4 gpfac(strt)=.01 end if C* gdis(1)=alpha*gwstor(1) C********For last station******* C********Check if current station (strt) is improved (BMP) C if (mainBMP.eq.1) then C if (which.eq.8) then C BMP = 0 C do 556 i = 1, 6 C if (strt.eq.order(i)+1) then C if (order(i).gt.0) then C BMP = 1 C end if C end if C 556 continue C end if C end if do 94 i=1,365 C****Set Runoff using SCS curve number method. C****For more information on the equations see, e.g. C****"http://www.shodor.org/master/environmental/water/runoff/RunoffAlgorithm.html" C****** The Curve Number array (see top) comes from C****** "http://www.lmnoeng.com/Hydrology/hydrology.htm#Curve Numbers" C****** They got it from TR-55, which is available on the web...or, C****** at least at present, from the collection of articles that I saved. C**** Pick curve number (bmp or not) if (BMP.eq.1) then currentcn=bmpcn(strt,soil) else currentcn=cn(strt,soil) end if coeff=0 cro(i)=0 C**** Correct curve number if soil is saturated or dry. I should probably set C**** this up as an option, since it is unclear how well it works...especially C**** as relates to BMPs. Does it make sense that if you implement an impr. C**** at a site, that since the soil is saturated on occasions when it isn't C**** without improvement, the CN(bmp) > CN(normal)???? C call corcur(currentcn, satur(i)) C >>> Basic Curve Number stuff. if (currentcn.gt.0) then s=((1000-10*currentcn)/currentcn) if (i.eq.365) then end if if (pcp(i).gt.(0.2*s)) then cRO(i)=(((pcp(i)-0.2*s)**2)/(pcp(i)+(0.8*s))) cRO(i) = cRO(i)*SAcres(strt) cRO(i) = cRO(i)*CtoM end if else cRO(i)=0 end if C <<< C >>> Lagging Cro--> RO if (cRO(i).lt.1) cRO(i)=0 curco=1 coeff=1 j=0 if (cRO(i).gt.0) then if (which.eq.8) then lag=slowlag(strt) else lag=lagchar(strt) end if lag=log(cRO(i))/lag 67 if (curco.gt.0.01) then j=j+1 curco=e**(-(1/lag)*j) coeff=coeff+curco go to 67 end if do 68 k=i,i+j if (k.le.365) then RO(k)=RO(k)+cRO(i)*(1/coeff)*(e**((-1/lag)*(k-i))) else RO(k-365)=RO(k-365)+cRO(i)*(1/coeff)*(e**((-1/lag)*(k-i))) end if 68 continue end if if (ro(i).lt.10e-3) ro(i)=0 C <<< totro(i) = totro(i)+cRO(i) C********Groundwater Discharge, gwstor and Evapo/ration C* Initial conditions are dependent on alpha+beta+cover+size. C********Evapotranspiration gamma=1 C**** gamma is the percent of PET that actually is evapotranspired. As soil C**** is dried out, evapotranspiration becomes more difficult and the percent C**** of water in "storage" that flows into the stream also decreases as the C**** soil attempts to maintain its concentration. Thus, as et(i)/gwstorage(strt) C**** gets large, (% PET of GWS), gamma decreases. Here, it acts somewhat C**** discretely, which doesn't bother _me_, but might need future improvement. C*****This sets the "p" and pesticide factor values, which are C*****time-of-year dependent...Damn Fortran 77; no case statements. do 95 j=1,12 thepest(j)=pestzero(j) if (strt.eq.3) thepest(j)=0.001*pestagri(j) if (strt.eq.4) thepest(j)=0.001*pesturb(j) if (strt.eq.7) thepest(j)=0.001*pesturb(j) 95 continue actipest(i)=thepest(12) p=.21 if (i.lt.334) then p=.22 actipest(i)=thepest(11) end if if (i.lt.303) then p=.25 actipest(i)=thepest(10) end if if (i.lt.273) then p=.28 actipest(i)=thepest(9) end if if (i.lt.243) then p=.31 actipest(i)=thepest(8) end if if (i.lt.212) then p=.33 actipest(i)=thepest(7) end if if (i.lt.181) then p=.34 actipest(i)=thepest(6) end if if (i.lt.151) then p=.32 actipest(i)=thepest(5) end if if (i.lt.120) then p=.30 actipest(i)=thepest(4) end if if (i.lt.90) then p=.27 actipest(i)=thepest(3) end if if (i.lt.59) then p=.24 actipest(i)=thepest(2) end if if (i.lt.31) then p=.22 actipest(i)=thepest(1) end if C********(Potential) Ev.Tr. = (Acres)*(% plant land cover) C *p(seebelow)*(.46T+8)/2.54(in/cm)*(1/10 cm/mm) C********This is a slightly modified Blaney-Criddle Method. C********See "http://www.fao.org/docrep/S2022E/s2022e07.htm". et(i)=(.46*(airfac(strt)*Airtmp(i)+2.28)+8) et(i)=SAcres(strt)*cover(strt)*et(i)*p/2.54*0.1*CtoM curgam=gamma if (et(i)*5.gt.gwstor(i)) curgam=.75*gamma if (et(i)*3.gt.gwstor(i)) curgam=.5*gamma if (et(i)*2.gt.gwstor(i)) curgam=.2*gamma if (et(i)*1.5.gt.gwstor(i)) curgam=.05*gamma et(i)=curgam*et(i) C********** C* C* Groundwater is basically a storage basin. basically (inexactly) C* alpha*weighteddecayingtimeavg(gwstor) is the total amount entering the stream from gw. C* C********** if (Sacres(strt).ne.0) then temp=beta(strt)*(pcp(i)-(cRO(i)/SAcres(strt))/CtoM) else temp=0 end if if (temp.lt.0) then temp=0 end if C* temp holds the pcp (inches) that ends up in storage after each i. if (gwstor(i)+temp*Sacres(strt)*CtoM.gt.et(i)) then gwstor(i+1)=gwstor(i)-et(i)+temp*Sacres(strt)*CtoM else gwstor(i+1)=0 et(i)=gwstor(i)+temp*Sacres(strt)*CtoM end if C*** This is obviously superflous, since curalp used to hold the C*** dynamic alpha value, but now alpha isn't changing. Tooo lazy C*** to change: curalp=alpha n=0 do 73 k=1,15 actlagco(k)=gwlagco(k) 73 continue do 71 k=i-14,i hold=0 if (k.gt.0) then if ((k.eq.1).and.(n.gt.0)) then do 72 o=1,n hold=hold+gwlagco(16-o) 72 continue do 74 o=1,15 actlagco(o)=actlagco(o)/(1-hold) 74 continue end if gdis(i+1)=gdis(i+1)+actlagco(i-k+1)*curalp*gwstor(k) gwstor(i+1)=gwstor(i+1)-actlagco(i-k+1)*curalp*gwstor(k) else n=n+1 end if if (initgrnd(strt).ne.0) then satur(i+1)=0.5*satur(i)+0.5*initsat(strt, bmp+1)* $(gwstor(i+1)/initgrnd(strt)) if (satur(i+1).gt.1) satur(i+1)=1 else satur(i+1)=0 end if 71 continue if (gdis(i).lt.10e-3) gdis(i)=0 94 continue C* This is basically the end of the hydrology calculations... C******River section inflow is runoff + ground discharge do 121 i=1,365 inflow(i)=RO(i)+Gdis(i) tgrndflow(i)=tgrndflow(i)+gdis(i) tinflow(i)=tinflow(i)+inflow(i) 121 continue C if (which.ne.8) then C! write (*,*) C! write (*,*) "Station: ", strt-1 call wbudge(pcp, sacres(strt), ctom, et, gwstor, inflow, gdis) C end if C//// Total Velocity, Depth, and Width calculations _based_ on Leopold's regression C//// analysis. w/v calculuations are probably neither mouth nor at a station, really C/// [see pg. 176, A view of the River] if (which.ne.8) then do 122 i=1,365 distemp=inflow(i)/86400 vel(i)=((0.6)*(distemp**0.34)) depth(i)=(.5*(distemp**0.4)) width(i)=((10.0/3.0)*(distemp**0.26)) 122 continue else do 123 i=1,365 distemp=tinflow(i)/86400 vel(i)=(.5*distemp**0.15) depth(i)=(.4*(distemp**0.4)) width(i)=((5.0)*(distemp**0.45)) 123 continue end if C----------------------------------------------- C - Do Biology schtuff, if necessary (requires calculating some chemistry, C - {my apologies for the algorithm}. Easier to design this way than scrap C - old stuff and fix.) C_______________________________________________ C---- Fish. All regression analysis from MBSS, and Don Shaffer. C---- m=6, b=-10. m(survey piedmont)=5.5701 b=-8.1135 m(survey coastal)=65936 b=-13.0055 if (bio.eq.1) then do 224 i=1,7 if (bmp.eq.0) then perdev(i)=tperdev(i) else perdev(i)=bperdev(i) end if 224 continue tempacres=sacres(strt) fishexpected=int(6*log10(tempacres)-10) benthexpected=23 do 225 i=1,365 pHR(i) = RO(i)*3.8*10**(-rph(strt)) pHG(i) = Gdis(i)*3.8*10**(-gph(strt)) tLd(i) = pHG(i)+ pHR(i) + tLd(i) C**** 5-day Smoothed PH if (inflow(i).gt.0) then holdyrpH(i) = -LOG10((pHR(i)+pHG(i))/(inflow(i)*3.8)) else holdyrpH(i) = 9999.9999 end if 225 continue call dorf(RNfac(strt),RbNfac(1,strt),bmp,rfac,cro,falagco(1)) call dogf(GNfac(strt),gfac,gwstor,initgrnd(strt),strt,0.33,bmp,1) call vfacclc(Gdis,GFac,RFac,RO,inflow,holdyrnit,NitLd,tLd) C**** Careful at the condition where inflow=0 do 227 i=1,365 if (i.gt.4) then notnine=0 do 229 j=i-4,i if (holdyrph(j).eq.9999.9999) notnine=1 229 continue if (notnine.eq.0) then ph(i)=.2*(holdyrph(i-4)+holdyrph(i-3)+holdyrph(i-2) $ + holdyrph(i-1)+holdyrph(i)) nitcon(i)=0.2*(holdyrnit(i-4)+holdyrnit(i-3)+holdyrnit(i-2) $ + holdyrnit(i-1)+holdyrnit(i)) end if else notnine=0 do 230 j=1,i if (holdyrph(j).eq.9999.9999) notnine=1 230 continue if (notnine.eq.0) then hold=0 holdcon=0 do 231 j=1,i hold=hold+(1.0/i)*(holdyrph(j)) holdcon=holdcon+(1.0/i)*(holdyrnit(j)) 231 continue ph(i)=hold nitcon(i)=holdcon end if end if 227 continue C***** do 228 i=1,365 eptnum(i)=8 if (nitcon(i).gt.1) then eptnum(i)=int(-.4*nitcon(i)+8.0) end if if (nitcon(i).gt.10.0) eptnum(i)=4 228 continue do 226 i=1,365 benthibi(i)=5.0 fishibi(i)=5.0 if (ph(i).eq.9999.9999) then fishibi(i)=9999.9999 benthibi(i)=9999.9999 end if if ((ph(i).ge.8.0).and.(ph(i).ne.9999.9999)) then benthibi(i)=3.0 end if C Step Function or... if (discr.eq.1) then if (ph(i).le.6.5) then fishibi(i)=3.2 benthibi(i)=2.8 else if (ph(i).le.6.0) then fishibi(i)=2.7 benthibi(i)=2.48 else if (ph(i).le.5.5) then benthibi(i)=2.3 else if (ph(i).le.5.0) then benthibi(i)=2.20 else if (ph(i).le.4.5) then fishibi(i)=1.3 benthibi(i)=1.50 else if (ph(i).le.4.0) then fishibi(i)=1.0 benthibi(i)=1.0 end if else C Regression Analysis.... {take your pick...} if (ph(i).ne.9999.9999) then if (ph(i).le.6.5) then benthibi(i)=.695*(ph(i))-1.6 fishibi(i)=.89412*(ph(i))-2.644 if (ph(i).le.4.0) then benthibi(i)=1.0 fishibi(i)=1.0 end if end if end if end if if ((strt.le.3).or.(strt.eq.6)) then if (fishibi(i).ne.9999.9999) then if (fishibi(i).gt.(-0.0243*PerDev(strt)+3.957)) then fishibi(i)=(-0.0243*PerDev(strt)+3.957) end if end if else if (fishibi(i).ne.9999.9999) then if (fishibi(i).gt.(-0.0197*PerDev(strt)+3.5485)) then fishibi(i)=(-0.0197*PerDev(strt)+3.5485) end if end if end if if (benthibi(i).ne.9999.9999) then if (benthibi(i).gt.-0.0327*PerDev(strt)+3.3035) then benthibi(i)=(-0.0327*PerDev(strt)+3.3035) end if end if if ((benthibi(i).gt.eptnum(i)-3).and.(nitcon(i).gt.1)) then benthibi(i)=eptnum(i)-3 end if if (ph(i).ne.9999.9999) then fishper(i)=(.1278*fishibi(i)+.3798) fishnum(i)=int(fishper(i)*fishexpected) fishper(i)=fishnum(i)/fishexpected if (sacres(strt).gt.300) then fishcount(i)=1.0375*(fishibi(i))+2.8125 else fishcount(i)=0.0 end if benthper(i)=(.0861*benthibi(i)+.558399) benthnum(i)=int(benthper(i)*benthexpected) benthper(i)=benthnum(i)/benthexpected else fishcount(i)=0.0 fishper(i)=0.0 fishnum(i)=0.0 benthper(i)=0.0 benthnum(i)=0.0 end if 226 continue end if if ((bio.eq.1).and.(which.eq.8)) then totacres=totacres+sacres(strt) do 232 i=1,365 if (ph(i).ne.9999.9999) then wtdeptnum(i)=wtdeptnum(i)+sacres(strt)*eptnum(i) wtdfishibi(i)=wtdfishibi(i)+sacres(strt)*fishibi(i) wtdbenthibi(i)=wtdbenthibi(i)+sacres(strt)*benthibi(i) end if 232 continue if (strt.eq.7) then fishexpected=int(6*log10(totacres)-10) do 233 i=1,365 if (tinflow(i).ne.0.0) then eptnum(i)=wtdeptnum(i)/totacres fishibi(i)=wtdfishibi(i)/totacres benthibi(i)=wtdbenthibi(i)/totacres fishper(i)=(.1278*fishibi(i)+.3798) fishnum(i)=int(fishper(i)*fishexpected) fishper(i)=fishnum(i)/fishexpected if (totacres.gt.300) then fishcount(i)=1.0375*(fishibi(i))+2.8125 else fishcount(i)=0 end if benthper(i)=(.0861*benthibi(i)+.558399) benthnum(i)=int(benthper(i)*benthexpected) benthper(i)=benthnum(i)/benthexpected else eptnum(i)=0.0 fishibi(i)=9999.9999 benthibi(i)=9999.9999 fishper(i)=0.0 fishnum(i)=0.0 fishcount(i)=0.0 benthper(i)=0.0 benthnum(i)=0.0 end if 233 continue end if end if C----------------------------------------------- C - Calculate the quality indicator requested C----------------------------------------------- C////----------Nitrogen----------//// if (ck.eq.nit) then call dorf(RNfac(strt),RbNfac(1,strt),bmp,rfac,cro,falagco(ck)) call dogf(GNfac(strt),gfac,gwstor,initgrnd(strt),strt,0.33,bmp,1) call vfacclc(Gdis,GFac,RFac,RO,inflow,NitCon,NitLd,tLd) call concFor7(strt, which, tcon,tLd, tinflow) if (which.eq.8) then C do 344 i=1,365 C write (*,*) strt,": ",Nitld(i) C 344 continue end if C////----------Phosphorous----------//// else if (ck.eq.phos) then call dorf(RPfac(strt),RbPfac(1,strt),bmp,rfac,cro,falagco(ck)) call dogf(GPfac(strt),gfac,gwstor,initgrnd(strt),strt,0.0,bmp,2) call vfacclc(Gdis,GFac,RFac,RO,inflow,PhsCon,PhsLd,tLd) call concFor7(strt, which, tcon,tLd, tinflow) C////----------Sediment----------//// else if (ck.eq.sed) then call dorf(RSdfac(strt),RbSdfac(1,strt),bmp,rfac,cro,falagco(ck)) call dogf(GSdfac(strt),gfac,gwstor,initgrnd(strt),strt,2.0,bmp,3) call vfacclc(Gdis,GFac,RFac,RO,inflow,SedCon,SedLd,tLd) call concFor7(strt, which, tcon,tLd, tinflow) C////----------Heavy Metals---------//// else if (ck.eq.hm) then call dorf(RHmfac(strt),RbHmfac(1,strt),bmp,rfac,cro,falagco(ck)) call dogf(GHmfac(strt),gfac,gwstor,initgrnd(strt),strt,0.0,bmp,4) call vfacclc(Gdis,GFac,RFac,RO,inflow,HMCon,HMLd, tLd) call concFor7(strt, which, tcon,tLd, tinflow) C////----------HydroCarbons---------//// else if (ck.eq.hy) then call dorf(RHyfac(strt),RbHyfac(1,strt),bmp,rfac,cro,falagco(ck)) call dogf(GHyfac(strt),gfac,gwstor,initgrnd(strt),strt,-0.1,bmp,5) call vfacclc(Gdis,GFac,RFac,RO,inflow,HYCon,HYLd,tLd) call concFor7(strt, which, tcon,tLd, tinflow) C////----------Pesticide -------------/ else if (kind.eq.pest) then call vfacclc(Gdis,zyr,actipest,ro,inflow,psCon,psLd,tLd) call concFor7(strt, which, tcon,tLd, tinflow) C////----------Air Temperature---------//// C///Note - air temp determined directly at station8 C///not as a result of subwatershed air temperature else if (kind.eq.atm) then if (which.lt.8) then do 31 i=1,365 Atemp(i)=airfac(strt)*airTmp(i)+2.28 31 continue C///else if which = 8, working with last station else if (strt.eq.7) then do 831 i=1,365 tcon(i)=airfac(8)*airTmp(i)+2.28 831 continue end if C////----------Water Temperature---------//// else if (kind.eq.wt) then call WaterTemp(wtemp,airtmp,shade(strt),inflow,airfac(strt)) do 32 i=1,365 tLd(i)=tLd(i)+(wtemp(i)*inflow(i)*3.8) 32 continue call concFor7(strt, which, tcon,tLd, tinflow) C////----------Saturated Dissolved Oxygen---------//// C///Saturated_DO_1 = Temp_effect_DO_1 else if (kind.eq.SDO) then call WaterTemp(wtemp,airtmp,shade(strt),inflow,airfac(strt)) do 33 i=1,365 tLd(i)=tLd(i)+wtemp(i)*inflow(i) 33 continue do 30 i=1,365 if (inflow(i).gt.0) then SatDO(i) = sdot(WTemp(i),e) else SatDO(i) = 9999.9999 end if 30 continue C///If last pass through last station if ((strt.eq.7).and.(which.eq.8)) then do 577 i=1,365 if (tinflow(i).gt.0) then Wtemp(i)=(tLd(i))/(tinflow(i)*3.8) tcon(i) = sdot(WTemp(i),e) else tcon(i) = 9999.9999 end if 577 continue end if C////------------pH-------------//// else if (kind.eq.pPh) then do 25 i=1,365 pHR(i) = RO(i)*3.8*10**(-rph(strt)) pHG(i) = Gdis(i)*3.8*10**(-gph(strt)) tLd(i) = pHG(i)+ pHR(i) + tLd(i) if ((strt.eq.7).and.(which.eq.8)) then if (tinflow(i).gt.0) then tcon(i) = -LOG10((tLd(i))/(tinflow(i)*3.8)) else tcon(i) = 9999.9999 end if else if (inflow(i).gt.0) then pH(i) = -LOG10((pHR(i)+pHG(i))/(inflow(i)*3.8)) else pH(i) = 9999.9999 end if end if 25 continue C////------------Dissolved Oxygen-------------//// C////DO = Saturated_DO - effect_OD else if (kind.eq.disO) then call WaterTemp(wtemp,airtmp,shade(strt),inflow,airfac(strt)) do 133 i=1,365 C//For station 7 tempLd(i) = Wtemp(i)*inflow(i) + tempLd(i) 133 continue C///Calculate SatDO do 130 i = 1, 365 SatDO(i) = sdot(WTemp(i),e) 130 continue C//Calculate BOD GcFac=GBODfc(strt) RcFac=RBODfc(strt) call clcprm(Gdis,GcFac,RcFac,RO,inflow,BODcon,BODLd,tLd) call ConcFor7(strt, which, tCon, tLd, tinflow) do 126 i=1,365 if (inflow(i).gt.0) then xDO(i) = SatDO(i)-ROperA(BODcon(i),effOD,11,100) else xDO(i) = 9999.9999 end if 126 continue C///If last pass through last station C///If last pass through last station then C//replace Wtemp with culmulative values if ((strt.eq.7).and.(which.eq.8)) then do 127 i=1,365 if (tinflow(i).gt.0) then Wtemp(i)=(tempLd(i))/(tinflow(i)*3.8) else Wtemp(i)=9999.9999 end if SatDO(i) = sdot(WTemp(i),e) tcon(i) = SatDO(i)-ROperA(tcon(i),effOD,11,100) 127 continue end if C////------------BOD-------------//// else if (kind.eq.bod) then GcFac=GBODfc(strt) RcFac=RBODfc(strt) call clcprm(Gdis,GcFac,RcFac,RO,tinflow,BODcon,BODLd,tLd) call ConcFor7(strt, which, tcon,tLd, tinflow) end if if (which.eq.8) then do 234 l=1,365 ro(l)=0 gdis(l)=0 234 continue end if C----------------------------------------------- C************END indicators******************** C----------------------------------------------- C***********Population and OtDev change with time C CPopG = CPopG + (growth) * dt C OtDev = OtDev + (Otcng) * dt C----------------------------------------- strt = strt+1 go to 333 end if C**************LOOP TO TOP************** 777 continue C*** Find tld's if applicable. Two if statements to fit on lines. yes=0 if ((ck.eq.1).or.(ck.eq.2).or.(ck.eq.3)) then yes=1 call sumx(tld, summer) end if if ((ck.eq.4).or.(ck.eq.5)) then yes=1 call sumx(tld, summer) end if C*** C/// if station 7, find mins and maxs if (which.eq.8) then if (kind.eq.outflow) then call outputer(tinflow,0.0,0,1,mainbmp,kind,which) else if (kind.eq.grndflow) then call outputer(tgrndflow,0.0,0,1,mainbmp,kind,which) else if (kind.eq.tro) then call outputer(totro,0.0,0,1,mainbmp,kind,which) else if (kind.eq.velo) then call outputer(vel,0.0,0,1,mainbmp,kind,which) else if (kind.eq.dep) then call outputer(depth,0.0,0,1,mainbmp,kind,which) else if (kind.eq.wid) then call outputer(width,0.0,0,1,mainbmp,kind,which) else if (kind.eq.atm) then call outputer(tcon,0.0,0,1,mainbmp,kind,which) else if (kind.eq.wt) then call outputer(tcon,0.0,0,1,mainbmp,kind,which) else if (kind.eq.diso) then call outputer(tcon,0.0,0,0,mainbmp,kind,which) else if (kind.eq.sdo) then call outputer(tcon,0.0,0,0,mainbmp,kind,which) else if (kind.eq.pph) then call outputer(tcon,0.0,0,0,mainbmp,kind,which) else if (kind.eq.fishi) then call outputer(fishibi,0.0,0,0,mainbmp,kind,which) else if (kind.eq.numfish) then call outputer(fishnum,0.0,0,1,mainbmp,kind,which) else if (kind.eq.fishperc) then call outputer(fishper,0.0,0,1,mainbmp,kind,which) else if (kind.eq.countfish) then call outputer(fishcount,0.0,0,1,mainbmp,kind,which) else if (kind.eq.benthi) then call outputer(benthibi,0.0,0,0,mainbmp,kind,which) else if (kind.eq.numbenth) then call outputer(benthnum,0.0,0,1,mainbmp,kind,which) else if (kind.eq.benthperc) then call outputer(benthper,0.0,0,1,mainbmp,kind,which) else if (kind.eq.numept) then call outputer(eptnum,0.0,0,1,mainbmp,kind,which) else j=0 do 343 k=51,55 if (kind.eq.k) then call outputer(tld,summer/1e6,(10**6),1,mainbmp,kind,which) j=1 end if 343 continue if (j.eq.0) then call outputer(tcon,summer/1e6,0,1,mainbmp,kind,which) end if end if else if (kind.eq.nit) then call outputer(NitCon,summer/1e6,0,0,mainbmp,kind,which) else if (kind.eq.phos) then call outputer(PhsCon,summer/1e6,0,0,mainbmp,kind,which) else if (kind.eq.sed) then call outputer(SedCon,summer/1e6,0,0,mainbmp,kind,which) else if (kind.eq.hm) then call outputer(HmCon,summer/1e6,0,0,mainbmp,kind,which) else if (kind.eq.hy) then call outputer(HyCon,summer/1e6,0,0,mainbmp,kind,which) else if (kind.eq.pest) then call outputer(PsCon,summer/1e6,0,0,mainbmp,kind,which) else if (kind.eq.wt) then call outputer(Wtemp,0.0,0,1,mainbmp,kind,which) else if (kind.eq.pPH) then call outputer(ph,0.0,0,0,mainbmp,kind,which) else if (kind.eq.disO) then call outputer(xdo,0.0,0,0,mainbmp,kind,which) else if (kind.eq.bod) then call outputer(bodcon,0.0,0,0,mainbmp,kind,which) else if (kind.eq.sdo) then call outputer(satdo,0.0,0,0,mainbmp,kind,which) else if (kind.eq.atm) then call outputer(atemp,0.0,0,1,mainbmp,kind,which) else if (kind.eq.tRo) then call outputer(cro,0.0,0,1,mainbmp,kind,which) else if (kind.eq.outflow) then call outputer(inflow,0.0,0,1,mainbmp,kind,which) else if (kind.eq.grndflow) then call outputer(gdis,0.0,0,1,mainbmp,kind,which) else if (kind.eq.velo) then call outputer(vel,0.0,0,1,mainbmp,kind,which) else if (kind.eq.wid) then call outputer(width,0.0,0,1,mainbmp,kind,which) else if (kind.eq.dep) then call outputer(depth,0.0,0,1,mainbmp,kind,which) else if (kind.eq.soilsat) then call outputer(satur,0.0,0,1,mainbmp,kind,which) else if (kind.eq.fishi) then call outputer(fishibi,0.0,0,0,mainbmp,kind,which) else if (kind.eq.numfish) then call outputer(fishnum,0.0,0,1,mainbmp,kind,which) else if (kind.eq.fishperc) then call outputer(fishper,0.0,0,1,mainbmp,kind,which) else if (kind.eq.countfish) then call outputer(fishcount,0.0,0,1,mainbmp,kind,which) else if (kind.eq.benthi) then call outputer(benthibi,0.0,0,0,mainbmp,kind,which) else if (kind.eq.numbenth) then call outputer(benthnum,0.0,0,1,mainbmp,kind,which) else if (kind.eq.benthperc) then call outputer(benthper,0.0,0,1,mainbmp,kind,which) else if (kind.eq.numept) then call outputer(eptnum,0.0,0,1,mainbmp,kind,which) else if (kind.eq.nitl) then call outputer(nitld,summer/1e6,(10**6),1,mainbmp,kind,which) else if (kind.eq.phosl) then call outputer(phsld,summer/1e6,(10**6),1,mainbmp,kind,which) else if (kind.eq.sedl) then call outputer(sedld,summer/1e6,(10**6),1,mainbmp,kind,which) else if (kind.eq.hyl) then call outputer(hyld,summer/1e6,(10**6),1,mainbmp,kind,which) else if (kind.eq.hml) then call outputer(hmld,summer/1e6,(10**6),1,mainbmp,kind,which) else error=.true. write (*,*) "error" close (20) stop end if end if close (20) C! if (.not.error) write (*,*) "done" stop end C*************************************************** C** Subroutine - Find the Concentration of indicator C*************************************************** subroutine ConcFor7(start, inputStn, Totcon,TotLd, flow) integer start, inputStn real Totcon(365), TotLd(365), flow(365) if ((start.eq.7).and.(inputStn.eq.8)) then do 901 i=1,365 Totcon(i)=0 if (flow(i).gt.0) Totcon(i)= TotLd(i)/(flow(i)*3.8) 901 continue end if return end C*************************************************** C** Subroutine - Find the min of a data array (365) C*************************************************** subroutine mn(thing, mini) real thing(365), mini mini=1e15 do 910 i=1,365 if (thing(i).lt.mini) mini=thing(i) 910 continue return end C*************************************************** C** Subroutine - Find the max of a data array (365) C*************************************************** subroutine mx(thing, maxu) real thing(365), maxu maxu=-1e15 do 911 i=1,365 if (thing(i).ne.9999.9999) then if (thing(i).gt.maxu) maxu=thing(i) end if 911 continue return end C*************************************************** C** Subroutine - Find the avg of a data array (365) C*************************************************** subroutine avgx(thing, avgn) real thing(365), avgn avgn=0 do 913 i=1,365 avgn=thing(i)+avgn 913 continue avgn=avgn/365 return end C*************************************************** C** Subroutine - Find the sum of a data array [for lds] C*************************************************** subroutine sumx(thing, total) real thing(365), total total=0 do 914 i=1,365 total=thing(i)+total 914 continue return end C*************************************************** C** Subroutine - Find the Water Temperature C*************************************************** subroutine WaterTemp(wtemp,airtemp,shade,inflow,afac) integer i,j,sprd real shade, tmavgt, md, frac, wt, cnt,lag,afac real wtemp(365), airtemp(365), inflow(365) real e data e/2.718281828/ md=0 do 904 i=1,365 md=md+inflow(i) 904 continue md=md/365.0 do 905 i=1,365 sprd=0 frac=((inflow(i))/md) if (frac.lt.0.001) frac=0 if (frac.gt.1) then sprd=shade*(7.0+.333333*frac) else sprd=shade*(7.0-4.0*(1-frac)) end if lag=sprd/2 cnt=0 tmavgt=0 do 903 j=i-sprd,i+sprd if (j-lag.gt.0) then if (j.eq.i+lag) then cnt=cnt+1 tmavgt=tmavgt+(afac*airtemp(j-lag)+2.28) else if (j-lag.lt.i) then wt=e**(j-i-lag) cnt=cnt+wt tmavgt=tmavgt+wt*(afac*airtemp(j-lag)+2.28) else wt=e**(i-j+lag) cnt=cnt+wt tmavgt=tmavgt+wt*(afac*airtemp(j-lag)+2.28) end if end if 903 continue tmavgt=tmavgt/cnt wtemp(i)=tmavgt 905 continue return end C********************************************** C*** Linear Interpolation Function C********************************************** real function ROperA(rain,curve,crvct,range) C////note - range gives the difference between first x(0) and last x real rain, curve(41) integer crvct, range,i real civl, curvx,x1,x2,y1,y2,x,y civl=range/float(crvct-1) C// write (*,*) rain, crvct, range,civl i=1 curvx=0 111 if (rain.ge.curvx) then curvx=curvx+civl i=i+1 go to 111 endif y1=curve(i-1) y2=curve(i) x1=curvx-civl x2=curvx x=rain y=(y2-y1)*(x-x1)/(x2-x1)+y1 ROperA=y return end C********************************************** C*** SDO(t) Function C********************************************** real function sdot(watertemp, e) C////note - range gives the difference between first x(0) and last x real watertemp, e, t real o t=watertemp+273.15 o=(-139.3411)+(1.575701e5/t)-(6.642308e7/(t**2)) o=o+(1.243800e10/(t**3))-(8.621949e11/(t**4)) C write (*,*) o sdot=e**o return end C*************************************************** C** Subroutine for kind of quality indicator- C** Calculate parameters [Currently, Used for BOD] C*************************************************** subroutine clcprm(Gdis,GFac,RFac,RO,inflow,parCon,parLd,tparLd) real Gdis(365),GFac,RFac,inflow(365),tparLd(365) real parLd(365), parCon(365),parG(365),parR(365) real RO(365) C/// 3.8 is the conversion from gal --> L C/// Thus; Parld is the mg's of given factor C/// Parcon is mg/L of factor do 912 i=1,365 parG(i)=Gdis(i)*GFac parR(i)=RO(i)*RFac parLd(i)= 3.8*parG(i)+3.8*parR(i) if (inflow(i).gt.0) then parCon(i)=parLd(i)/(3.8* inflow(i)) else parcon(i)=9999.9999 end if tparLd(i)=tparLd(i)+parLd(i) 912 continue return end C*************************************************** C** Subroutine for kind of quality indicator- C** Calculate parameters [date-variable factors] C*************************************************** subroutine vfacclc(Gdis,GFac,RFac,RO,inflow,parCon,parLd,tparLd) real Gdis(365),GFac(365),RFac(365),inflow(365),tparLd(365) real parLd(365), parCon(365),parG(365),parR(365) real RO(365) C/// 3.8 is the conversion from gal --> L C/// Thus; Parld is the mg's of given factor C/// Parcon is mg/L of factor do 912 i=1,365 parG(i)=Gdis(i)*GFac(i) parR(i)=RO(i)*RFac(i) parLd(i)= (3.8*parG(i)+3.8*parR(i)) if (inflow(i).gt.0) then parCon(i)=parLd(i)/(3.8* inflow(i)) else parcon(i)=9999.9999 end if tparLd(i)=tparLd(i)+parLd(i) 912 continue return end C*************************************************** C** Output-er Subroutine C*************************************************** subroutine outputer(putter, sumld, diver, wantzero,bmp,kn,wh) real putter(365) real sumld, avgn, maxu, mini, old integer diver, wantzero, kn, wh,bmp, x character*18 space(4) if (diver.gt.0) then do 917 i=1,365 putter(i)=putter(i)/diver 917 continue end if call mn(putter, mini) call mx(putter, maxu) call avgx(putter, avgn) 110 format (' #', 5x, F12.4) write(*,110) mini write(*,110) maxu write(*,110) avgn if (sumld.gt.0) then write(*,110) sumld else write(*,*) "# n/a" end if do 918 i=1,365 112 format(2x,i3,2x,F12.4) if (putter(i).ne.9999.9999) then write(*,112) i, putter(i) else C// Printing hard returns at 0 flow if wantzero=0. if (wantzero.eq.1) then write(*,112) i, putter(i) else write(*,*) end if end if 918 continue return end C*************************************************** C** Water Budget Calculator C*************************************************** subroutine wbudge(pcp, acreage, ctom, et, gwstor, inflow, gdis) real pcp(365), ctom, et(365), gwstor(366), inflow(365), gdis(366) integer acreage, i real sumflow, sumet, sumin, sumro, sumdis, sumlost sumin=0 sumet=0 sumflow=0 sumdis=0 sumro=0 do 919 i=1,365 sumin=sumin+pcp(i)*acreage*ctom/1000 sumet=sumet+et(i)/1000 sumflow=sumflow+inflow(i)/1000 sumdis=sumdis+gdis(i)/1000 sumro=sumflow-sumdis 919 continue temp=gwstor(366)-gwstor(1) C write (*,*) "Totals: [in thousands of m^3]" C write (*,*) "Totalin:", sumin, " totET:", sumet, " totFLOW:" C $, sumflow," AddtoGwStor:",(gwstor(366)-gwstor(1))/1000 C write (*,*) C write (*,*) "Inflow Totals:" C if (sumflow.ne.0) then C write (*,*) "% GW dis:", sumdis/sumflow," % RO:", sumro/sumflow C else C write (*,*) "Sumflow is zero." C end if C write (*,*) C write (*,*) "percentages" C if (sumin.ne.0) then C write (*,*) "et (%):",sumet/sumin," flow:",sumflow/sumin C $, " Gw addition:", temp/(sumin*1000) C write (*,*) "Balance ?:", (sumet+sumflow+temp/1000)/sumin C else C write (*,*) "Sumin is zero." C end if return end C*************************************************** C** Correct Curve Number (for soil conditions) C*************************************************** subroutine corcur(cn, satu) integer cn real satu if (satu.gt.0.8) then cn=((23*cn)/(10.0+0.13*cn)) else if (satu.lt.0.2) then cn=((4.2*cn)/(10.0-0.058*cn)) end if return end C*************************************************** C** Assign Rfacs Subprocedure C*************************************************** subroutine dorf(norfac, bmpfac, bmp, rfac, cro, lagco) real norfac, bmpfac, rfac(365), cro(365), lagco, temp, e real agco integer bmp, j, k, l e=2.718281828 if (bmp.eq.1) then do 38 j=1,365 rfac(j)=bmpfac 38 continue else do 39 j=1,365 rfac(j)=norfac 39 continue end if do 41 j=1,365 if (cro(j).gt.0) then do 42 k=1,7 agco=lagco if (j+k.lt.365) then if (cro(j+k).gt.0) then do 43 l=1,15 if (bmp.eq.1) then if (rfac(j+k+l-1).lt.bmpfac) agco=agco*1.25 else if (rfac(j+k+l-1).lt.norfac) agco=agco*1.25 end if if ((j+k+l-1).le.365) then temp=1-(1/agco)*(e**((-(1.0/agco))*(k+(15-l)-1))) C if (temp.gt.0.99) temp=1 rfac(j+k+(15-l))=rfac(j+k+(15-l))*temp else temp=1-(1/agco)*(e**((-(1.0/agco))*(k+(15-l)-1))) C if (temp.gt.0.99) temp=1 rfac(j+k+(15-l)-366)=rfac(j+k+(15-l)-366)*temp end if 43 continue end if end if 42 continue end if 41 continue return end C*************************************************** C** Assign Gfacs Subprocedure C*************************************************** subroutine dogf(nogfac,gfac,gstr,ingws,strt,factin,bmp,code) real nogfac, gfac(365), gstr(366), lagger, y, factin real ingws, negfac, bfac, fact integer j, strt,k,bmp, code real wbfc(5,5) C/// Wetland BMP Groundwater factors data (wbfc(1,i),i=1,5)/0.28,0.28,0.28,0.28,0.28/ data (wbfc(2,i),i=1,5)/0.02,0.02,0.02,0.02,0.02/ data (wbfc(3,i),i=1,5)/2,2,2,2,2/ data (wbfc(4,i),i=1,5)/0.0,0.0,0.0,0.0,0.0/ data (wbfc(5,i),i=1,5)/0.0,0.0,0.0,0.0,0.0/ negfac=factin fact=factin if (strt.ne.6) then do 48 j=1,365 gfac(j)=nogfac 48 continue else if (bmp.eq.1) then nogfac=wbfc(code,1) fact=(1.0/3.0)*factin end if do 49 j=1,365 y=0 if (j-2.gt.0) then do 50 k=j-2,j if (fact.lt.0) negfac=-fact y=y+(1.0/3.0)*gstr(k)/(ingws*(2.0/negfac)) 50 continue else if (j-1.gt.0) then if (fact.lt.0) negfac=-fact y=(gstr(j-1)+gstr(j))/(ingws*(4.0/negfac)) else if (fact.lt.0) negfac=-fact y=gstr(j)/(ingws*(2.0/negfac)) end if gfac(j)=nogfac+(y-0.5)*(fact/5.0) if (gfac(j).lt.0) gfac(j)=0 49 continue end if return end C*************************************************** C****************THE END*************************** C***************************************************