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//////////////////////////////////////////////////////////////////////////// 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) 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, DO: Uses Curves (TeffDO, EffOD) real GNFac(7), GPFac(7),GSdFac(7) real RNFac(7), RPFac(7), RSdFac(7) real RbNFac(7), RbPFac(7), RbSdFac(7) real GHMFac(7),GHYFac(7),GWTFac(7) real RHMFac(7),RHYFac(7),RWTFac(7) real RbHMFac(7),RbHYFac(7),RbWTFac(7) real RBODfc(7),GBODfc(7), rbbodfc(7), soilin(7) real rph(7), rbph(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 integer soil, n, f integer S90(7) real TeffDO(41), 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),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), 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/0.544,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/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/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/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/0,0,0,1.7,9,0,9/ C/// Water Temperature Factors (not used [see shade]) data GWTFac/0.681,0.681,0.681,0.681,0.681,0.681,0.681/ data RWTFac/0,0,0,0,0,0,0/ C/// Air Factors data airfac/0.681,0.681,0.881,0.781,0.981,0.781,0.981,0.881/ 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/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/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,25000,80000,25000,15000,2500,30000/ C/// Ground cover which evapotranspires data cover/1.4,1,1,.75,.15,1.5,.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,9.76e5,2.675e6,3.75e5/ data (initgrnd(i),i=5,7)/0,7e4,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/// 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. C data lagchar/5,5.5,7,8,9,5,9/ data lagchar/2,2,2.5,3.5,5,1,5/ 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 TeffDO /14.6, 14.2, 13.8, 13.4, 13.1, 12.8, 12.4, 12.1 $ , 11.8, 11.6, 11.3, 11.0, 10.8, 10.5 $ , 10.3, 10.1, 9.86, 9.65, 9.45, 9.26, 9.08, 8.90, 8.73, 8.56 $ , 8.40, 8.24, 8.09 $ , 7.95, 7.81, 7.67, 7.54, 7.41, 7.29, 7.17, 7.05, 6.94, 6.82 $ , 6.72, 6.61, 6.51, 6.41/ 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 lag=log(cRO(i))/lagchar(strt) 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----------------------------------------------- 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(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(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(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(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(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) = ROperA(WTemp(i),TeffDO,41,40) 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) = ROperA(WTemp(i),TeffDO,41,40) 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) = ROperA(WTemp(i),TeffDO,41,40) 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) = ROperA(Wtemp(i),TeffDO,41,40) 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.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 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** 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****************THE END*************************** C***************************************************