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 C//////////////////////////////////////////////////////////////////////////// program riverweb C/////// Declare the variables and constants ////// character*6 cid character*20 iname,oname logical error integer nit, phos, sed, hm, hy, wt, pPh, disO, bod, sdo, atm, tRO integer outflow, grndflow, pest 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) 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) 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 HtLd(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 integer soil, n, f integer S90(7) real effOD(41) real RO(365),cro(365) real GFac(365),RFac(365), gcfac, rcfac growth, Otcng real pcpm86(365), pcpl85(365), pcph84(366), pcp(365), pcpmu integer year,OtA90,OtDev,Ot85,P1995,P1990 integer CPopG, SAcres(7) 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) real WTemp(365),ATemp(365),SatDO(365) real bodG(365),pHR(365),pHG(365),pH(365),xDO(365) real templd(365),zyr(365) real SChange(7), ROperA,thecrv(41) real cover(7), initgrnd(7) integer CN(7,4), BMPCN(7,4) integer currentcn, yes real p, s, gwstor(366), et(365) real alpha, beta(7), gamma, curalp, curgam, shade(7) 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. GN:basically, trying to create reasonable ranges.] C/// Thus, GW is CBP for Stations 3,4,5,~6. 0,1, and 2 were altered from CBP to C/// decrease the shortterm _range_ of concentrations.] C old:data GNFac/0.4,0.45,7.392,2.23,4,0.4,4/ data GNFac/0.4,0.45,8.4,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,3.086,5.108,0.5,5.108/ 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 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.57,0.26,0.02,0.32/ 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/// This is an even more problematic set of numbers as a whole than the others. C/// I don't even really know what is meant by sediments. (Total Suspended Solids? Dissolved?) C/// I dunno.... 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,140,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)] 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.111,.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) 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.7,9,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 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???] 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'm too lazy to C//// bother even 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. data Rph/6.5,6.5,5.35,6.1,5.9,5.7,5.9/ data (Rbph(1,i),i=1,7)/6.5,6.5,5.35,6.1,5.9,5.7,5.9/ 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,7,5.9,6.6,6.4,5.7,6.4/ 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/ C/// Acreage in each terrain C data S90/89139,89139,105387,15950,7672,2429,2357/ data S90/40000,15000,80000,25000,15000,2500,30000/ C/// Total Acreage (s90 2: 207500 acres = 324 square miles) C/// Total Acreage (s90 (original): 312073 acres = 490 square miles) C/// Ground cover which evapotranspires data cover/1.4,1,1,.75,.15,1,.35/ 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/// Curve Number Data 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)/72,81,88,91/ data (CN(4,i),i=1,4)/54,78,80,85/ data (CN(5,i),i=1,4)/89,92,94,95/ data (CN(6,i),i=1,4)/0,0,0,0/ data (CN(7,i),i=1,4)/77,85,90,92/ 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,79,84/ data (BMPcn(5,i),i=1,4)/85,89,91,92/ data (BMPcn(6,i),i=1,4)/0,0,0,0/ data (BMPcn(7,i),i=1,4)/77,85,90,92/ C/// Initial Groundwater Storage (These are all calibrated to the final C/// day values for the year...[These are soil type C] data (initgrnd(i),i=1,4)/1.204e6,5.9e5,2.675e6,3.75e5/ data (initgrnd(i),i=5,7)/0,2e5,7.5e5/ C/// Old [non m3. G] C data (initgrnd(i),i=1,4)/7.5e8,8.8e8,9e8,6.25e7/ C data (initgrnd(i),i=5,7)/0,1.2e7,1.5e8/ C/// C/// Initial Percentage of PCP-RO that ends up in Gwstorage (currently C/// appears to be behaving 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. data beta/1,1,1,.5,.5,1,.5/ C/// On second or third thought, it seems to me that beta should all be 1. C/// Of course, one way of instantly accomplishing this is to raise initgd C/// for 3,4, and 6; but it ends up having to get ludicrously large to C/// balance out. Another possible solution is to increase cover. This C/// works, but than sigmaET seems larger than you'd expect. Maybe, it C/// makes sense to increase the CNs. Or, something. It doesn't make C/// sense to have extremely large GWSTOR when the RO should be >> GW. C/// Pesticide Data (1-12 arrays for "urban" and "agriculture" both 75 percentile) 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 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 should be relatively small. C/// The wetlands (station 5, strt6) value is completely irrelevant, since C/// there is no RO. data lagchar/8,8,8,8,9,8,9/ data slowlag/1,1,2,2,2,1,2/ C data lagchar/6,6,7.5,10,15,3,15/ C/// gwlagco is a slow exponential decay. Sigma i=1-15 = 1. 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. data falagco/30,1,1,2,4/ 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---------------------------------- OtA90 = 4898 year = 1986 C/// This is setting the year to use for the precipitation data C---------------------------- P1995 = 78400 CPopG = P1995 P1990 = 71347 C growth = (P1995-P1990)/P1990/5/365{monthly growth rate}*CPopG growth = (P1995-P1990)/P1990/5/365*CPopG Ot85 = 4776 OtDev = OtA90 Otcng = (OtA90-Ot85)/Ot85/1825*OtDev TotAs = 222937 C error=.false. C------------------------------------- C///// get Precipitation Data ////// 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 C note check last else with DS 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-------------------------------------- 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----------------------------------------- C///// set Acres ////// C----------------------------------------- do 40 i=1,7 SAcres(i)=S90(i) 40 continue C//read the process id read (*,*) cid iname = '/tmp/riverdata' j=len(cid) iname(15:14+j) cid(1:j) open (19, FILE=iname, STATUS='OLD') C-----------Read strt subwatershed--------- read (19,83)which C--------0-7 from web corresponds to 1-8 in fortran program which = which + 1 C-------------Read kind of quality indicator read (19,84)kind C-------------Read BMP; BMP = 1 for Best Management Practices read (19,81) BMP C-----------Read list of improvements------- read (19,82)(order(i),i=1,6) C-----------Read soil Type (IN FUTURE...) C read (19,83) soil C soil=3 close (19) C-------------Create output file name, which depends on BMP mainBMP = BMP if (mainBMP.eq.0) then oname='/tmp/riverout' oname(14:13+j)=cid(1:j) else oname='/tmp/riverbout' oname(15:14+j)=cid(1:j) end if 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) C********For last station******* if (mainBMP.eq.1) then if (which.eq.8) then BMP = 0 do 556 i = 1, 6 c write (*,*) strt," ", order(i) if (strt.eq.order(i)+1) then if (order(i).gt.0) then BMP = 1 end if end if 556 continue end if end if C*****Pick curve number (bmp or not) if (BMP.eq.1) then C write (*,*) strt," using BMP" currentcn=bmpcn(strt,soil) else C write (*,*) strt," using present" currentcn=cn(strt,soil) end if 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" do 66 i=1,365 coeff=0 cro(i)=0 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 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 totro(i) = totro(i)+cRO(i) 66 continue C********Groundwater, gwheight and Evapo/ration alpha=.002 C* alpha is the percentage in GW storage that enters into the stream. C* gwstor(1)=initgrnd(strt) C* Initial groundwater conditions (calibrated with end of the current year C* "empirically") gdis(1)=alpha*gwstor(1) temp=0 do 94 i=1,365 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. (see the changing C*****alpha factors below. Thus, as et(i)/gwstorage(strt) gets large, C*****(% PET of GWS), gamma decreases. Here, it acts somewhat discretely, which C*****likely needs future improvement. C*****This sets the "p" and pesticide factor values, which are C*****time-of-year dependent 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*0.75*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 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********(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". C temp=beta(strt)*(pcp(i+1)-(cRO(i+1)/SAcres(strt))/CtoM) temp=beta(strt)*(pcp(i)-(cRO(i)/SAcres(strt))/CtoM) 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 curalp=alpha if (gwstor(i)*2.lt.initgrnd(strt)) curalp=.5*alpha if (gwstor(i)*4.lt.initgrnd(strt)) curalp=.25*alpha if (gwstor(i)*8.lt.initgrnd(strt)) curalp=.1*alpha if (gwstor(i)*12.lt.initgrnd(strt)) curalp=.075*alpha if (gwstor(i)*20.lt.initgrnd(strt)) curalp=.05*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 71 continue C********** C* Groundwater is basically a storage basin. basically (inexactly) C* alpha*weighteddecayingtimeavg(gwstor) is the total amount entering the stream from gw. C* This algorithm probably needs a bit more thought and/or justification. C********** if (gdis(i).lt.10e-3) gdis(i)=0 94 continue 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 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)=((2.0/3.0)*(distemp**0.3)) depth(i)=(.5*(distemp**0.3)) width(i)=(3*(distemp**0.4)) 122 continue else do 123 i=1,365 distemp=tinflow(i)/86400 vel(i)=((2.0/3.0)*(distemp**0.1)) depth(i)=(0.4*(distemp**0.4)) width(i)=(3.75*(distemp**0.5)) 123 continue end if C----------------------------------------------- C - Calculate the quality indicator requested C----------------------------------------------- C////----------Nitrogen----------//// if (kind.eq.nit) then do 48 j=1,365 if (bmp.eq.1) then rfac(j)=RbNfac(1,strt) else rfac(j)=RNfac(strt) end if 48 continue do 41 j=1,365 if (cro(j).gt.0) then do 46 k=1,20 rfac(j+k)=(e**(-1/(falagco(kind)*k)))*rfac(j+k) 46 continue end if gfac(j)=GnFac(strt) 41 continue call vfacclc(Gdis,GFac,RFac,RO,inflow,NitCon,NitLd,tLd) call concFor7(strt, which, tcon,tLd, tinflow) C////----------Phosphorous----------//// else if (kind.eq.phos) then do 49 j=1,365 if (bmp.eq.1) then rfac(j)=RbPfac(1,strt) else rfac(j)=RPfac(strt) end if 49 continue do 50 j=1,365 if (cro(j).gt.0) then do 51 k=1,20 rfac(j+k)=(e**(-1/(falagco(kind)*k)))*rfac(j+k) 51 continue end if gfac(j)=GPFac(strt) 50 continue call vfacclc(Gdis,GFac,RFac,RO,inflow,PhsCon,PhsLd,tLd) call concFor7(strt, which, tcon,tLd, tinflow) C////----------Sediment----------//// else if (kind.eq.sed) then do 52 j=1,365 if (bmp.eq.1) then rfac(j)=RbSdFac(1,strt) else rfac(j)=RSdFac(strt) end if 52 continue do 53 j=1,365 if (cro(j).gt.0) then do 54 k=1,20 rfac(j+k)=(e**(-1/(falagco(kind)*k)))*rfac(j+k) 54 continue end if gfac(j)=GSdFac(strt) 53 continue call vfacclc(Gdis,GFac,RFac,RO,inflow,SedCon,SedLd,tLd) call concFor7(strt, which, tcon,tLd, tinflow) C////----------Heavy Metals---------//// else if (kind.eq.hm) then do 55 j=1,365 if (bmp.eq.1) then rfac(j)=RbHmFac(1,strt) else rfac(j)=RHmFac(strt) end if 55 continue do 56 j=1,365 if (cro(j).gt.0) then do 57 k=1,20 rfac(j+k)=(e**(-1/(falagco(kind)*k)))*rfac(j+k) 57 continue end if gfac(j)=GHmFac(strt) 56 continue call vfacclc(Gdis,GFac,RFac,RO,inflow,HMCon,HMLd, tLd) call concFor7(strt, which, tcon,tLd, tinflow) C////----------HydroCarbons---------//// else if (kind.eq.hy) then do 58 j=1,365 if (bmp.eq.1) then rfac(j)=rbHYFac(1,strt) else rfac(j)=rHYFac(strt) end if 58 continue do 59 j=1,365 if (cro(j).gt.0) then do 60 k=1,20 rfac(j+k)=(e**(-1/(falagco(kind)*k)))*rfac(j+k) 60 continue end if gfac(j)=ghyfac(strt) 59 continue 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 ((kind.eq.1).or.(kind.eq.2).or.(kind.eq.3)) then yes=1 call sumx(tld, summer) end if if ((kind.eq.4).or.(kind.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) else if (kind.eq.grndflow) then call outputer(tgrndflow, 0.0, 0, 1) else if (kind.eq.tro) then call outputer(totro, 0.0, 0, 1) else if (kind.eq.velo) then call outputer(vel, 0.0, 0, 1) else if (kind.eq.dep) then call outputer(depth, 0.0, 0, 1) else if (kind.eq.wid) then call outputer(width, 0.0, 0, 1) else if (kind.eq.atm) then call outputer(tcon, 0.0, 0, 1) else if (kind.eq.wt) then call outputer(tcon, 0.0, 0, 1) else if (kind.eq.diso) then call outputer(tcon, 0.0, 0, 0) else if (kind.eq.sdo) then call outputer(tcon, 0.0, 0, 0) else if (kind.eq.pph) then call outputer(tcon, 0.0, 0, 0) else call outputer(tcon, summer/1e6, 0, 0) end if else if (kind.eq.nit) then call outputer(NitCon, summer/1e6, 0, 0) else if (kind.eq.phos) then call outputer(PhsCon, summer/1e6, 0, 0) else if (kind.eq.sed) then call outputer(SedCon, summer/1e6, 0, 0) else if (kind.eq.hm) then call outputer(HmCon, summer/1e6, 0, 0) else if (kind.eq.hy) then call outputer(HyCon, summer/1e6, 0, 0) else if (kind.eq.pest) then call outputer(PsCon, summer/1e6, 0, 0) else if (kind.eq.wt) then call outputer(Wtemp, 0.0, 0, 1) else if (kind.eq.pPH) then call outputer(ph, 0.0, 0, 0) else if (kind.eq.disO) then call outputer(xdo, 0.0, 0, 0) else if (kind.eq.bod) then call outputer(bodcon, 0.0, 0, 0) else if (kind.eq.sdo) then call outputer(satdo, 0.0, 0, 0) else if (kind.eq.atm) then call outputer(atemp, 0.0, 0, 1) else if (kind.eq.tRo) then call outputer(cro, 0.0, 0, 1) else if (kind.eq.outflow) then call outputer(inflow, 0.0, 0, 1) else if (kind.eq.grndflow) then call outputer(gdis, 0.0, 0, 1) else if (kind.eq.velo) then call outputer(vel, 0.0, 0, 1) else if (kind.eq.wid) then call outputer(width, 0.0, 0, 1) else if (kind.eq.dep) then call outputer(depth, 0.0, 0, 1) else error=.true. write (*,*) "error" close (20) stop end if end if close (20) 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/ do 904 i=1,365 md=md+inflow(i) 904 continue md=md/365 do 905 i=1,365 frac=((inflow(i))/md) if (frac.gt.1) then sprd=shade*(7+.333333*frac) else if (1/frac.lt.12) then sprd=shade*(7-.333333*(1/frac)) else sprd=shade*3 end if 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 C write (*,*) wtemp(i)," ",airtemp(i)," ",cnt," ",lag 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)) 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) real putter(365) real sumld, avgn, maxu, mini integer diver, wantzero 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(20,110) mini write(20,110) maxu write(20,110) avgn if (sumld.gt.0) then write(20,110) sumld else write(20,*) "# n/a" end if do 918 i=1,365 112 format(2x,i3,2x,F12.4) if (putter(i).ne.9999.9999) then write(20,112) i, putter(i) else C// Printing hard returns at 0 flow if wantzero=0. if (wantzero.eq.1) then write(20,112) i, putter(i) else write(20,*) 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 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 write (*,*) "Totals: [in thousands of m^3]" write (*,*) "Totalin:", sumin, " totET:", sumet, " totFLOW:" $, sumflow," AddtoGwStor:",(gwstor(366)-gwstor(1))/1000 write (*,*) write (*,*) "Inflow Totals:" write (*,*) "% GW dis:", sumdis/sumflow," % RO:", sumro/sumflow write (*,*) write (*,*) "percentages" write (*,*) "et (%):",sumet/sumin," flow:",sumflow/sumin $, " Gw addition:", (gwstor(366)-gwstor(1))/(sumin*1000) temp=gwstor(366)-gwstor(1) write (*,*) "Balance ?:", (sumet+sumflow+temp/1000)/sumin return end C*************************************************** C****************THE END*************************** C***************************************************