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 C///// Total Flow 20, Groundwater Stream Input 21 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 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 (p86=13,p85=14,p84=15,a86=16,a85=17,a84=18,fnum=19) parameter (outflow=20, grndflow=21) 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 GHMFac(7),GHYFac(7),GWTFac(7) real RHMFac(7),RHYFac(7),RWTFac(7) real RBODfc(7),GBODfc(7) C note - air temp determined directly at river mouth - so 8 factors. real airfac(8) 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), tLd(365), real temp integer soil integer S90(7), S85(7) real TeffDO(41), effOD(41) real RO(365) real GFac,RFac, growth, Otcng real pcpm86(365), pcpl85(365), pcph84(366), pcp(365), pcpmu integer year,OtA90,OtDev,Ot85,P1995,P1990 integer CPopG, SAcres(7),CtoG integer strt,which,stp, kind, BMP,order(6),cnt, mainBMP real Gdis(365),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) real SChange(7), ROperA,thecrv(41) real cover(7), initgrnd(7) integer CN(7,4), BMPCN(7,4) integer currentcn real p, s, gwstor(365), et(365) real alpha, beta(7), gamma, curalp, curgam, shade(7) C/// alpha is the percentage of groundwater storage that outflows every day C/// Nitrate Factors data GNFac/0.946,0.946,7.392,2.23,5.23,0.946,5.23/ data RNFac/0.544,0.544,9.326,3.086,5.108,0.5,5.108/ C/// Phosphate Factors data GPFac/0.007,0.007,0.007,0.197,0.057,0.057,0.057/ data RPFac/0.093,0.093,1.836,0.408,0.58,0.02,0.58/ C/// Sediment Factors data GSdFac/0,0,0,0,0,0,0/ data RSdFac/5,5,140,140,140,3,140/ C/// Heavy Metal Factors data GHMFac/0,0,0,0,0,0,0/ data RHMFac/0,0,0,0.103,0.892,0.1,0.892/ C/// Hydrocarbon Factors data GHYFac/0,0,0,0,0,0,0/ data RHYFac/0,0,0,2,10,0,10/ C/// Water Temperature Factors (not used) 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 data RBODfc/58,58,92,91,68,21,92/ data GBODfc/7*20/ C/// the crv is zerve / tld is 0-de / 'nd inflow low data thecrv/41*0.0/ data tLd/365*0.0/ data tinflow/365*0.0/ C/// Acreage in each terrain data S90/89139,89139,105387,15950,7672,2429,2357/ data S85/91434,91434,108945,10671,7192,2463,2261/ 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,.75,.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)/5e8,5e8,5e8,5e7/ data (initgrnd(i),i=5,7)/0,1.2e7,1e7/ 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/// helps drive surface evaporation, which will lower this factor. data beta/1,1,1,.5,1,1,.45/ 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******************************************************************************* C***** Reference Equations ***** C******************************************************************************* C G_Nitrogen(i) = Gdis[]*3.8{L/gal}*GNFac[i]{mg/L Total N} C Runoff_Nitrogen(i) = RunOff[]*3.8{L/gal}*RNFac[i]{mg/L N} C Nitrogen_conc(i) = (Runoff_Nitrogen(i)+ G_Nitrogen(i))/(inflow*3.8){L/gal} C Nitrogen_Load(i) = G_Nitrogen(i)+Runoff_Nitrogen(i) C** NitCon(8) = NitLd(8)/(inflow(8)*3.8) need to move to end C******************************************************************************* C G_Phosphorous_1 = Gdis_1*3.8{L/gal}*0.007{mg/L N} C Runoff_Phosphorous_1 = RunOff_1*3.8{L/gal}*0.093{mg/L P} C Phosphorous_conc_1 = (Runoff_Phosphorous_1+G_Phosphorous_1)/(inflow_1*3.8){L/gal} C** PhsCon(8) = PhsLd(8)/(inflow(8)*3.8) Need to move to end C******************************************************************************* C G_Sediment_1 = Gdis_1*3.8{L/gal}*0{mg/L TSS} C Runoff_Sediment_1 = RunOff_1*3.8{L/gal}*5{mg/L TSS} C Sediment_conc_1 = (Runoff_Sediment_1+G_Sediment_1)/(inflow_1*3.8){L/gal} C Sediment_Load_1 = G_Sediment_1+Runoff_Sediment_1 C** SedCon(8) = SedLd(8)/(inflow(8)*3.8) Need to move to end C******************************************************************************* C G_Heavy_Metals_1 = Gdis_1*3.8{L/gal}*0{mg/L Metals} C Runoff_Heavy_Metals_1 = RunOff_1*3.8{L/gal}*0{mg/L Heavy Metals} C Heavy_Metals_conc_1 = (Runoff_Heavy_Metals_1+G_Heavy_Metals_1)/(inflow_1*3.8){L/gal} C Heavy_Metals_Load_1 = G_Heavy_Metals_1+Runoff_Heavy_Metals_1 C** HMCon(8) = HMLd(8)/(inflow(8)*3.8) Need to move to end C******************************************************************************* C G_Hydrocarbons_1 = Gdis_1*3.8{L/gal}*0{mg/L hydrocarbons} C Runoff_Hydrocarbons_1 = RunOff_1*3.8{L/gal}*0{mg/L Hydrocarbons} C Hydrocarbons_conc_1 = (Runoff_Hydrocarbons_1+G_Hydrocarbons_1)/(inflow_1*3.8){L/gal} C Hydrocarbons_Load_1 = G_Hydrocarbons_1+Runoff_Hydrocarbons_1 C** HYCon(8) = HYLd(8)/(inflow(8)*3.8) Need to move to end C******************************************************************************* C G_Water_Temp_1 = Gdis_1*3.8{L/gal}*0.681*Air_Temp_1+2.28{'C} C Runoff_Water_Temp_1 = RunOff_1*3.8{L/gal}*Air_Temp_1{'C} C Water_Temperature_1 = (Runoff_Water_Temp_1+G_Water_Temp_1)/(inflow_1*3.8){L/gal} C Water_Temp_Load_1 = G_Water_Temp_1+Runoff_Water_Temp_1 C** WTemp(8) = WTLd(8)/(inflow(8)*3.8) Need to move to end C******************************************************************************* C--bodR_1 = RunOff_1*3.8{L/gal}*58{8+50 mg/L BOD+COD} C--bodG_1 = Gdis_1*3.8{L/gal}*20{mg/L BOD+COD} C--BOD_conc_1 = (bodR_1+bodG_1)/(inflow_1*3.8){L/gal} C******************************************************************************* C--moles/L H3O+ mean pH for precipitation in Bay area pH 4.5 C--Runoff_pH_1 = RunOff_1*3.8{L/gal}*10**(-4.5) C--{moles/L H2O+assume alkalinity of soils results in pH of 7.0 for groundwater} C--G_pH_1 = Gdis_1*3.8{L/gal}*10^(-7.0) C--pH_1 = -LOG10((Runoff_pH_1+G_pH_1)/(inflow_1*3.8)){L/gal} C******************************************************************************* C---------------------------------- OtA90 = 4898 year = 1986 C/// This is setting the year to use for the precipitation and acreage 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 C 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 44 i=1,7 SAcres(i)=S90(i) 44 continue C------ check where this goes---- C// This formulation is the change of acres per day. C// 1825 is number of days in 5 years. C// There is no trade off - everything increases - perhaps in C// "other developed"?? C// do 113 i=1,7 C// SChange(i)=(S90(i)-S85(i))/1825*SAcres(i) C// SAcres(i)=SAcres(i)+SChange(i)*dt C// 113 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 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 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---------------------------------------------------- C********For last station******* if (mainBMP.eq.1) then if (which.eq.8) then BMP = 0 do 556 i = 1, 6 if (order(i).eq.(strt-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 (*,*) "using BMP" currentcn=bmpcn(strt,soil) else C write (*,*) "using present" currentcn=cn(strt,soil) end if C write (*,*) currentcn 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 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 RO(i)=(((pcp(i)-0.2*s)**2)/(pcp(i)+(0.8*s))) RO(i) = RO(i)*SAcres(strt) RO(i) = RO(i)*CtoG end if else RO(i)=0 end if tcon(i) = RO(i)/1000 + tcon(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. p=.21 if (i.lt.334) p=.22 if (i.lt.303) p=.25 if (i.lt.273) p=.28 if (i.lt.243) p=.31 if (i.lt.212) p=.33 if (i.lt.181) p=.34 if (i.lt.151) p=.32 if (i.lt.120) p=.30 if (i.lt.90) p=.27 if (i.lt.59) p=.24 if (i.lt.31) p=.22 et(i)=(.46*(airfac(strt)*Airtmp(i)+2.28)+8) et(i)=SAcres(strt)*cover(strt)*et(i)*p/2.54*0.1*CtoG 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=.15*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". temp=beta(strt)*(pcp(i+1)-(RO(i+1)/SAcres(strt))/CtoG) if (temp.lt.0) temp=0 C* temp holds the pcp (inches) that ends up in storage after each i. if (gwstor(i).gt.et(i)) then gwstor(i+1)=gwstor(i)-et(i)+(temp)*SAcres(strt)*CtoG else gwstor(i+1)=(temp)*SAcres(strt)*CtoG end if curalp=alpha if (gwstor(i)*2.lt.initgrnd(strt)) curalp=.5*alpha if (gwstor(i)*5.lt.initgrnd(strt)) curalp=.2*alpha if (gwstor(i)*10.lt.initgrnd(strt)) curalp=.1*alpha if (gwstor(i)*20.lt.initgrnd(strt)) curalp=.05*alpha gdis(i+1)=0.1*curalp*gwstor(i-2) gdis(i+1)=gdis(i+1)+0.15*curalp*gwstor(i-1) gwstor(i-1)=gwstor(i-1)-0.15*curalp*gwstor(i-1) gdis(i+1)=gdis(i+1)+0.25*curalp*gwstor(i) gwstor(i)=gwstor(i)-0.25*curalp*gwstor(i) gdis(i+1)=gdis(i+1)+0.5*curalp*gwstor(i+1) gwstor(i+1)=gwstor(i+1)-0.5*curalp*gwstor(i+1) C********** C* Groundwater is basically a storage basin. basically (inexactly) C* alpha*fourdaytimeavgd(gwstor) is the total amount entering the stream from gw. C* This algorithm probably needs a bit more thought and/or justification. C********** C write (*,*) i,":", " gw:",gwstor(i)," gd:",gdis(i)," et:",et(i) C write (*,*) i, " gw:",gwstor(i)," gi:",temp*SAcres(strt)*CtoG C write (*,*) i, " ro:",ro(i), " gd:",gdis(i)," et:",et(i) 94 continue C******River section inflow is runoff + ground discharge do 121 i=1,365 inflow(i)=RO(i)+Gdis(i) tinflow(i)=tinflow(i)+inflow(i) tgrndflow(i)=tgrndflow(i)+gdis(i) 121 continue C----------------------------------------------- C - Calculate the quality indicator requested C----------------------------------------------- C////----------RunOff----------//// if (kind.eq.tRo) then do 667 i=1, 365 RO(i) = RO(i)/1000 667 continue C////----------Nitrogen----------//// else if (kind.eq.nit) then GFac=GNFac(strt) RFac=RNFac(strt) call clcprm(Gdis,GFac,RFac,RO,inflow,NitCon,NitLd,tLd) call concFor7(strt, which, tcon,tLd, tinflow) C////----------Phosphorous----------//// else if (kind.eq.phos) then GFac=GPFac(strt) RFac=RPFac(strt) call clcprm(Gdis,GFac,RFac,RO,inflow,PhsCon,PhsLd,tLd) call concFor7(strt, which, tcon,tLd, tinflow) C////----------Sediment----------//// else if (kind.eq.sed) then GFac=GSdFac(strt) RFac=RSdFac(strt) call clcprm(Gdis,GFac,RFac,RO,inflow,SedCon,SedLd,tLd) call concFor7(strt, which, tcon,tLd, tinflow) C////----------Heavy Metals---------//// else if (kind.eq.hm) then GFac=GHMFac(strt) RFac=RHMFac(strt) call clcprm(Gdis,GFac,RFac,RO,inflow,HMCon,HMLd, tLd) call concFor7(strt, which, tcon,tLd, tinflow) C////----------HydroCarbons---------//// else if (kind.eq.hy) then GFac=GHYFac(strt) RFac=RHYFac(strt) call clcprm(Gdis,GFac,RFac,RO,inflow,HYCon,HYLd,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) = 0 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) = 0 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**(-4.5) pHG(i) = Gdis(i)*3.8*10**(-7.0) 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 if (tcon(i-1).gt.0) then tcon(i) = tcon(i-1) else tcon(i) = 7 end if end if else if (inflow(i).gt.0) then pH(i) = -LOG10((pHR(i)+pHG(i))/(inflow(i)*3.8)) else pH(i) = pH(i-1) 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 GFac=GBODfc(strt) RFac=RBODfc(strt) call clcprm(Gdis,GFac,RFac,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) = 0 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) Wtemp(i)=tempLd(i)/(tinflow(i)*3.8) 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 GFac=GBODfc(strt) RFac=RBODfc(strt) call clcprm(Gdis,GFac,RFac,RO,tinflow,BODcon,BODLd,tLd) call ConcFor7(strt, which, tcon,tLd, tinflow) 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// if station 7, then write total if (which.eq.8) then do 707 i=1,365 if (kind.eq.outflow) then write(20,FMT='(F14.4)') tinflow(i)/1000 else if (kind.eq.grndflow) then write(20,FMT='(F14.4)') tgrndflow(i)/1000 else write(20,FMT='(F12.4)') tcon(i) end if 707 continue else if (kind.eq.nit) then do 708 i=1,365 write(20,FMT='(F12.4)') NitCon(i) 708 continue else if (kind.eq.phos) then do 709 i=1,365 write(20,FMT='(F12.4)') PhsCon(i) 709 continue else if (kind.eq.sed) then do 710 i=1,365 write(20,FMT='(F12.4)') SedCon(i) 710 continue else if (kind.eq.hm) then do 711 i=1,365 write (20,FMT='(F12.4)') HMCon(i) 711 continue else if (kind.eq.hy) then do 712 i=1,365 write(20,FMT='(F12.4)') HYCon(i) 712 continue else if (kind.eq.wt) then do 713 i=1,365 write(20,FMT='(F12.4)') WTemp(i) 713 continue else if (kind.eq.pPH) then do 714 i=1,365 write(20,FMT='(F12.4)') pH(i) 714 continue else if (kind.eq.disO) then do 715 i=1,365 write(20,FMT='(F12.4)') xDO(i) 715 continue else if (kind.eq.bod) then do 716 i=1,365 write(20,FMT='(F12.4)') BODCon(i) 716 continue else if (kind.eq.sdo) then do 717 i=1,365 write(20,FMT='(F12.4)') SatDO(i) 717 continue else if (kind.eq.atm) then do 718 i=1,365 write(20,FMT='(F12.4)') Atemp(i) 718 continue else if (kind.eq.tRo) then do 719 i=1,365 write(20,FMT='(F14.4)') RO(i) 719 continue else if (kind.eq.outflow) then do 720 i=1,365 write(20,FMT='(F14.4)') inflow(i)/1000 720 continue else if (kind.eq.grndflow) then do 721 i=1,365 write(20,FMT='(F14.4)') gdis(i)/1000 721 continue else error=.true. write (*,*) "error" close (20) stop 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 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) 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=(2.718281828)**(j-i-lag) cnt=cnt+wt tmavgt=tmavgt+wt*(afac*airtemp(j-lag)+2.28) else wt=(2.718281828)**(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 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) 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) parcon(i)=0 if (inflow(i).gt.0) parCon(i)=parLd(i)/(3.8* inflow(i)) tparLd(i)=tparLd(i)+parLd(i) 912 continue return end C*************************************************** C****************THE END*************************** C***************************************************