C/////// Page 550 red book C/////// gallons per day C////// Kind Codes 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???? need to change most xx(i) i=1,7 to xx(strt, i) i=1,365 C???? also need to carry out specific code depending on kind C////////////////////////////////////////////////////////////////////////// program riverweb C/////// Declare the variables and constants ////// character*6 cid character*20 iname,oname logical error real ROCrv(7,25), ROBMP(7,25) integer CrvCt(7),BMPCt(7) 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) parameter (pPh=7, 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) real airt85(365),airt86(365),airt84(366),AirTmp(365) C //// missing SDO, DO 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) C// Initialize the Rsta (river station stocks) integer Gsta(7), Rsta(8), GSInit(7), RSInit(8) 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 real Gdis(365),Grecrg(365),inflow(365),tinflow(365) real WTR(365),bodR(365) real WTemp(365),outflow(365),ATemp(365),SatDO(365) real bodG(365),pHR(365),pHG(365),pH(365),xDO(365) C/// Not used: GODLd(365) real SChange(7), ROperA,thecrv(41) 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/ 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/ data GSdFac/0,0,0,0,0,0,0/ data RSdFac/5,5,140,140,140,3,140/ data GHMFac/0,0,0,0,0,0,0/ data RHMFac/0,0,0,0.103,0.892,0.1,0.892/ data GHYFac/0,0,0,0,0,0,0/ data RHYFac/0,0,0,2,10,0,10/ 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/ data airfac/0.681,0.681,0.881,0.781,0.981,0.781,0.981,0.881/ data RBODfc/58,58,92,91,68,21,92/ data GBODfc/7*20/ data thecrv/41*0.0/ data RSInit/7*1000000,7000000/ data GSInit/7*1000000/ data tLd/365*0.0/ data tinflow/365*0.0/ data S90/89139,89139,105387,15950,7672,2429,2357/ data S85/91434,91434,108945,10671,7192,2463,2261/ data CrvCt/25,25,21,25,25,25,25/ data BMPCt/25,25,25,25,21,25,21/ data (ROCrv(1,i),i=1,12)/0.00, 0.00, 0.00, 0.09, 0.24, 0.46, 0.71 $ , 1.01, 1.33, 1.67, 2.04, 2.46/ data (ROCrv(1,i),i=13,25)/2.81, 3.24, 3.62, 4.08, 4.46, 4.92 $ , 5.33, 5.82, 6.22, 6.72, 7.13, 7.62, 8.05/ data (ROCrv(2,i),i=1,12)/0.00, 0.03, 0.08, 0.3, 0.56, 0.89 $ , 1.25, 1.64, 2.04, 2.46, 2.89, 3.36/ data (ROCrv(2,i),i=13,25)/3.78, 4.26, 4.69, 5.22, 5.63, 6.12 $ , 6.57, 7.02, 7.52, 7.98, 8.48, 9.00, 9.45/ data (ROBMP(2,i),i=1,12)/0.00, 0.00, 0.00, 0.09, 0.24, 0.46 $ , 0.71, 1.01, 1.33, 1.67, 2.04, 2.46/ data (ROBMP(2,i),i=13,25)/2.81, 3.24, 3.62, 4.08, 4.46, 4.92 $ , 5.33, 5.82, 6.22, 6.72, 7.13, 7.62, 8.05/ data (ROCrv(3,i),i=1,12)/0.00, 0.12, 0.27, 0.65, 1.18, 1.59 $ , 2.10, 2.70, 3.10, 3.60, 4.30, 4.86/ data (ROCrv(3,i),i=13,21)/5.40, 6.00, 6.54, 7.18, 7.68, 8.34 $ , 8.88, 9.54, 10.1/ data (ROBMP(3,i),i=1,12)/0.00, 0.02, 0.08, 0.29, 0.56, 0.89 $ , 1.25, 1.64, 2.04, 2.46, 2.89, 3.30/ data (ROBMP(3,i),i=13,25)/3.78, 4.26, 4.69, 5.16, 5.63, 6.12 $ , 6.57, 7.02, 7.52, 8.04, 8.48, 9.00, 9.45/ data (ROCrv(4,i),i=1,12)/0.00, 0.03, 0.08, 0.29, 0.56, 0.89 $ , 1.25, 1.64, 2.04, 2.46, 2.89, 3.30/ data (ROCrv(4,i),i=13,25)/3.78, 4.20, 4.69, 5.16, 5.63, 6.06 $ , 6.57, 6.96, 7.52, 7.92, 8.48, 9.00, 9.45/ data (ROBMP(4,i),i=1,12)/0.00, 0.01, 0.03, 0.17, 0.38, 0.65 $ , 0.96, 1.30, 1.67, 2.05, 2.45, 2.88/ data (ROBMP(4,i),i=13,25)/3.28, 3.66, 4.15, 4.62, 5.04, 5.52 $ , 5.95, 6.42, 6.88, 7.38, 7.81, 8.34, 8.76/ data (ROCrv(5,i),i=1,12)/0.00, 0.24, 0.56, 0.97, 1.48, 1.96 $ , 2.45, 2.94, 3.43, 3.92, 4.42, 4.98/ data (ROCrv(5,i),i=13,25)/5.41, 5.94, 6.41, 6.90, 7.40, 7.98 $ , 8.40, 8.88, 9.40, 9.90, 10.4, 10.9, 11.4/ data (ROBMP(5,i),i=1,12)/0.00, 0.12, 0.27, 0.65, 1.18, 1.59 $ , 2.10, 2.70, 3.10, 3.60, 4.30, 4.86/ data (ROBMP(5,i),i=13,21)/5.40, 6.00, 6.54, 7.18, 7.68, 8.34 $ , 8.88, 9.54, 10.1/ data (ROCrv(6,i),i=1,12)/0.00, 0.12, 0.32, 0.7, 1.09, 1.53 $ , 1.98, 2.45, 2.92, 3.40, 3.88, 4.38/ data (ROCrv(6,i),i=13,25)/4.85, 5.34, 5.82, 6.42, 6.81, 7.32 $ , 7.79, 8.34, 8.78, 9.30, 9.77, 10.3, 10.8/ data (ROCrv(7,i),i=1,12)/0.00, 0.12, 0.32, 0.7, 1.09, 1.53 $ , 1.98, 2.45, 2.92, 3.40, 3.88, 4.38/ data (ROCrv(7,i),i=13,25)/ 4.85, 5.34, 5.82, 6.42, 6.81, 7.32 $ , 7.79, 8.34, 8.78, 9.30, 9.77, 10.3, 10.8/ data (ROBMP(7,i),i=1,12)/0.00, 0.12, 0.27, 0.65, 1.18, 1.59 $ , 2.10, 2.70, 3.10, 3.60, 4.30, 4.86/ data (ROBMP(7,i),i=13,21)/ 5.40, 6.00, 6.54, 7.18, 7.68, 8.34 $ , 8.88, 9.54, 10.1/ 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/// Initialize River Station Stocks (Rsta) Do 10 i=1,8 Rsta(i)=RSInit(i) C TRANSIT TIME = 1 C INFLOW LIMIT = INF C CAPACITY = INF 10 continue Do 111 i=1,7 Gsta(i) = GSInit(i) 111 continue C---------------------------------- OtA90 = 4898 year = 1986 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 units of mg/L C ?????????where are these declared??????????? WetN = 0 subN = 2.25 urbN = 2.25 IndN = 1.75 AgN = 2.25 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 - so number of C// acres a day. 1825 is number of days in 5 years. C// There is no trade off - everything increases - perhaps in C// "other developed"?? do 113 i=1,7 SChange(i)=(S90(i)-S85(i))/1825*SAcres(i) SAcres(i)=SAcres(i)+SChange(i)*dt 113 continue C//read the process id C//Commented out ONLY for TESTING!!! C-- read (*,*) cid iname = '/tmp/riverdata1234' C-- j=len(cid) C-- 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,7) close (19) do 780 i = 1, 6 write (*,*) order(i) 780 continue C-------------Create output file name, which depends on BMP C////MODIFIED FILE NAMES if (which.eq.8) then BMP = 0 do 557 i = 1, 6 if (order(i).gt.0) then BMP = 1 end if 557 continue end if if (BMP.eq.0) then oname='/tmp/riverout1' oname(14:13+j)=cid(1:j) C/////Process ID's not used right now. else oname='/tmp/riverout2' oname(15:14+j)=cid(1:j) C/////Process Id's not used. end if open (20, FILE=oname, STATUS='UNKNOWN') 82 format (6i1) 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//DIVIDE BY A MILLION CtoG=(1/12.0)*43560*7.48056 if (which.eq.8) then strt=1 stp=7 else strt=which stp=which endif C------------------------------------------ C-?????? figure out where this goes??????----- C------------------------------------------ C-- outflow_1 = CONVEYOR OUTFLOW do 65 i=1,8 outflow(i)=rsta(i) 65 continue C------------------------------------------ C*****If strt = 8, the last station is being C********* analyzed, thus need to add up all 6 C*********stations and go through the loop 8 times C*****Otherwise, for strt<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 (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 C********For best management practice installed******* if (BMP.eq.1) then write (*,*) "using BMP" cnt=BMPCt(strt) do 55 i=1,cnt thecrv(i)=ROBMP(strt,i) 55 continue C********For best management practice NOT installed***** else cnt = CrvCt(strt) write (*,*) "using Present" do 56 i=1,cnt thecrv(i)=ROCrv(strt,i) 56 continue end if C********Set RO, Gdis, Grecrg do 66 i=1,365 RO(i)=SAcres(strt)*ROperA(pcp(i),thecrv,cnt,12)*CtoG Grecrg(i)=pcp(i)*SAcres(strt)*CtoG-RO(i) Gdis(i)=1000000 tcon(i) = RO(i)+ tcon(i) 66 continue C------------------------- 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) 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.at) then if (which.lt.8) then do 31 i=1,365 Atemp(i)=airfac(strt)*airTmp(i)+2.28 31 continue 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 do 32 i=1,365 Atemp(i)=airfac(strt)*airTmp(i)+2.28 GWT=Gdis(i)*3.8*Atemp(i) RWT=RO(i)*3.8*Atemp(i) WTLd(i)=GWT+RWT Wtemp(i)=WTLd(i)/(inflow(i)*3.8) tLd(i)=tLd(i)+WTLd(i) 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 do 33 i=1,365 Atemp(i)=airfac(strt)*airTmp(i)+2.28 GWT=Gdis(i)*3.8*Atemp(i) RWT=RO(i)*3.8*Atemp(i) WTLd(i)=GWT+RWT Wtemp(i)=WTLd(i)/(inflow(i)*3.8) tLd(i)=tLd(i)+WTLd(i) 33 continue do 30 i=1,365 SatDO(i) = ROperA(WTemp(i),TeffDO,41,40) 30 continue C///If last pass through last station if (strt.eq.7) then Wtemp(i)=tLd(i)/(inflow(i)*3.8) do 577 i=1,365 tcon(i) = ROperA(WTemp(i),TeffDO,41,40) 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) then tcon(i) = -LOG10((tLd(i))/(tinflow(i)*3.8)) else pH(i) = -LOG10((pHR(i)+pHG(i))/(inflow(i)*3.8)) end if 25 continue C////------------Dissolved Oxygen-------------//// C////DO = Saturated_DO - effect_OD else if (kind.eq.disO) then do 133 i=1,365 Atemp(i)=airfac(strt)*airTmp(i)+2.28 GWT=Gdis(i)*3.8*Atemp(i) RWT=RO(i)*3.8*Atemp(i) WTLd(i)=GWT+RWT Wtemp(i)=WTLd(i)/(inflow(i)*3.8) C//For station 7 tLd(i) = WTLd(i) + tLd(i) 133 continue do 130 i=1,365 SatDO(i) = ROperA(WTemp(i),TeffDO,41,40) 130 continue C///If last pass through last station if (strt.eq.7) then Wtemp(i)=tLd(i)/(inflow(i)*3.8) do 567 i=1,365 SatDO(i) = ROperA(WTemp(i),TeffDO,41,40) 567 continue end if GFac=GBODfc(strt) RFac=RBODfc(strt) call clcprm(Gdis,GFac,RFac,RO,inflow,BODcon,BODLd,tLd) call ConcFor7(strt, which, BODCon, BODLd, tinflow) do 126 i=1,365 xDO(i) = SatDO(i)-ROperA(BODcon(i),effOD,11,100) 126 continue C///If last pass through last station if (strt.eq.7) then do 127 i=1,365 tcon(i) = SatDO(i)-ROperA(BODcon(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 CPopG = CPopG + (growth) * dt OtDev = OtDev + (Otcng) * dt C*******************Stocks********************** do 213 i=1,7 Gsta(i)=Gsta(i)+(Grecrg(i)-Gdis(i))*dt 213 continue C***************River Stations************** do 313 i=1,7 Rsta(i)=Rsta(i)+(inflow(i)-outflow(i))*dt 313 continue C----------------------------------------- strt = strt+1 go to 333 end if C**************LOOP TO TOP************** 777 continue C////The following part writes out the indicators** do 779 i=1,6 write(*,*) order(i) C write(*,*) inflow(i) 779 continue C// if station 7, then write total if (which.eq.8) then do 707 i=1,365 write(20,FMT='(F12.4)') tcon(i) 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.f)') 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 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)= TotLd(i)/(flow(i)*3.8) 901 continue end if 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,RO(365),inflow(365),tparLd(365) real parLd(365), parCon(365),parG(365),parR(365) do 912 i=1,365 parG(i)=Gdis(i)*3.8*GFac parR(i)=RO(i)*3.8*RFac parLd(i)=parG(i)+parR(i) parCon(i)=parLd(i)/(inflow(i)*3.8) tparLd(i)=tparLd(i)+parLd(i) 912 continue return end C*************************************************** C****************THE END*************************** C***************************************************