**************************************************************************** * BFI415.FOR - A COMPUTER PROGRAM FOR COMPUTING AN INDEX TO BASE FLOW * * COMPILED FOR THE PC USING MICROSOFT FORTRAN 4.10 * ALSO COMPILABLE USING THE UNIX F77 COMPILER * * April 2007, Version 4.15 * -------------------------------------------------------------------------- * Tony L. Wahl, Mail Code D-8560 | Kenneth L. Wahl, Mail Stop 406 * U.S. Bureau of Reclamation | U.S. Geological Survey * P.O. Box 25007 | Water Resources Division * Lakewood, Colorado 80225-0007 | Lakewood, Colorado 80225 * Phone: 303-445-2155 | Phone: 303-236-5950, ext. 220 * E-Mail: twahl@do.usbr.gov | E-mail: klwahl@usgs.gov * twahl@nyx.net | * * The latest version of the BFI program may be obtained on the * World Wide Web at http://www.usbr.gov/wrrl/twahl/bfi.html **************************************************************************** * * REFERENCES * ---------- * Additional information about the base-flow separation algorithm used * by this program can be found in the following references: * * Institute of Hydrology, 1980. Low Flow Studies: Wallingford, Oxon, * United Kingdom, Report No. 1, 21 p. * * Institute of Hydrology, 1980. Low Flow Studies: Wallingford, Oxon, * United Kingdom, Report No. 3, pp. 12-19. * * Swan, D. R. and R. Condie, 1983. "Computation of the Base Flow Index," * Water Resources Branch Inland Water Directorate, Environment Canada, * Ottawa, Ontario, 21 p. * * Wahl, Kenneth L. and Tony L. Wahl, 1988. "Effects of Regional Ground- * Water Level Declines on Streamflow in the Oklahoma Panhandle," * Proceedings of the American Water Resources Association Symposium * on Water-Use Data for Water Resources Management, Tucson, Arizona, * 11 p. (Please note: this paper misstates the criteria used for * determination of base-flow turning points.) * * Wahl, Kenneth L. and Tony L. Wahl, 1995. "Determining the Flow of Comal * Springs at New Braunfels, Texas," Texas Water '95, San Antonio, * Texas, August 16-17, 1995, 10 p. * * A comparison of this and many other automated base-flow separation * methods is contained in: * * Mau, David P., 1993, "Estimating Ground Water Recharge and Base Flow * from Streamflow Hydrographs for a Small Appalachian Mountain Basin," * Master of Science thesis, University of Colorado at Denver. * **************************************************************************** * * REVISION HISTORY * ---------------- * Version 1 - Original program used for 1988 Oklahoma Panhandle paper. * * Version 2 - Adds the ability to interpolate turning points for first * and last days when data are continuous at one gage from one year to * the next. Also, all turning points are output to a file, and base-flow * output is sent to both the screen and a file. Each year's base-flow * output shows the first and last turning point used for that year, and * indicates if those points were interpolated. If possible, the program * will use data from incomplete years to interpolate endpoints for * complete years. * * Version 3 (Aug. 1991) - Modified to facilitate easy use of turning * point periods other than 5 days. The GETQ subroutine was also modified to * account for the lack of leap-days in century years not divisible by 400. * (There is a leap-day in the year 2000, but not in 2100, 2200, or 2300). * * Version 3.1 (Nov. 1992) - Adds the capability to output daily base flow * and total flow values to a file. The interpolation of base-flow values * was also placed in a separate FUNCTION. * * Version 3.11 (May 1993) - Corrects bug that prevented printout of daily * base-flow value when last day of year was a base-flow turning point. * * Version 3.12 (June 1993) - Corrects bug that prevents printout of the daily * flow corresponding to the day of the last non-interpolated turning point in * each year. This bug was created by the fix in version 3.11. * * Version 3.13 (July 1993) - Corrects bug that caused program to abort * prematurely if N was evenly divisible into 366 and the last year in the * data file was a leap year. Also improved handling of incomplete or * scrambled data files. Fixed problem in which missing 3 cards caused math * errors or end-of-file errors. Program will now interpolate across years in * more situations than before and will output as much turning point and * base-flow data as possible, even in years where data are incomplete. * Simplified some of the interpolation logic. Added additional comments * to help users troubleshoot problems with missing data. Modified GETQ * and READCARD routines to monitor the month and week fields on 3-cards * and call out errors if cards are out of order or missing. This version * of the program was used for the Comal Springs paper. A special version * called BFICOMAL.FOR was also used by others for some analysis of Comal * Springs. It printed daily base flows to the nearest 1 cfs, and contained * a second modification that was subsequently included in version 3.14. * * Version 3.14 (June 1995) - Incorporates a modification originally * contained in a special version of 3.13 titled BFICOMAL.FOR. Daily base * flows to be printed to the daily flow output file are checked against the * corresponding daily total flow. If the computed base flow (interpolated * between turning points) is greater than the daily flow, then the base flow * for that day is printed out as being equal to the daily flow. Interpolated * base flows occasionally exceed the total flow when analyzing spring creeks * or other streams heavily dominated by base flow, or when analyzing records * with large numbers of zero-flow days. Random errors in daily discharge * measurements during extended base-flow periods can also produce this effect. * (Added IF statements before first two WRITE (DAYFILE,77) statements, and * a MIN function statement before third WRITE (DAYFILE,77) statement.) * Note that this does not affect the computation of the total volume of base * flow or the Base-Flow Index. Those are still computed using the * interpolated base flow, even though base flow on a given day might * occasionally exceed total flow. Added a printout of the program version * number on-screen and in the output files. * * Version 4.00 (July 1995) - Added new TEST_FOR_TP method, a variation of * the standard Institute of Hydrology method. Also modified standard * Institute of Hydrology method when one of adjacent minimums is zero. * It was previously impossible to have a minimum be a turning point if either * adjacent minimum had zero flow. Program now does the test in only one * direction if one adjacent minimum is zero. If both adjacent minima are zero * the point still cannot be a turning point. * Modified test so that all zero flow points would automatically be * turning points. This corrected a problem for years without any flow; * program previously never found a first turning point. Also added tests * on the calculation of the base-flow index and statistics to eliminate * divide-by-zero errors when dealing with years without any flow. * Added ability to process data files organized by water year or calendar * year. Program automatically detects data-file type and proceeds * accordingly. Modified date output in turning point and daily flows * files to give calendar year, month, and day (previously day was shown as * days from the start of the year), and eliminated gage ID from each line of * output files. Made modifications to make program calculate * the small increments of base flow on the first half day and last half * day of the year. The program computes the base flow as a definite * integral from day 1 through 365 in a normal, complete water year. This * only gives base flow for a sum of 364 days, with the first and last * half days being left out, if one assumes the reported average daily * discharge to be the instantaneous discharge at midday. The program * previously accounted for this by only counting the first and last days * as a half day when summing the total runoff, then calculating the BFI * from that result. The total annual runoff was then recalculated giving * all days of the year full weight and the base flow was back-calculated * (extrapolated) using the BFI and correctly calculated total runoff. * The program now uses the interpolated turning points going across years * to directly calculate the assumed instantaneous base flow at time t=0.5 days * and at time t=365.5 days (or 366.5 days in leap years). These times * correspond to the start of the first day and end of the last day, using the * previous assumption that the average daily flows are the instantaneous flows * at midday, and that t=1.000 corresponds to midday on the first day of the * year. The program can then calculate the base flow for the entire * year; all days receive equal weight in determining the base flow and * base-flow index. The program also now prints out a symbol ('!') next to * the annual base-flow volume if it could not be directly calculated due * to lack of data needed for the required interpolations. This gives the * user a more direct indication that the base-flow index and total annual * base-flow volume shown were obtained from an analysis that could not * include the entire water year. This could be determined previously * only by closely examining the last two columns of the printout, which * indicated the first and last turning points included in each year's * analysis. * * Version 4.01 (January 1996) - Added a UCASE function to convert user input * to uppercase and thus simplify the testing of user input. Added subroutine * to present a "parameter menu" to the user and allow them to change program * options and parameters using that menu. Added routines to deal with * output files that already exist, or an input file that does not exist. * Minor modifications to output text and prompts were made in February 1996. * COMMON blocks were removed March 2, 1996, and default summary output file * extension was changed to ".bfi". * * Version 4.02 (March 3, 1996) - Fixed error in turning point test for the * standard method. When Q2 was zero and Q0 was non-zero, it did not test * FRAC*Q1 against Q0, (it tested it against Q1 instead), so it always found * Q1 to be a turning point if FRAC was less than 1. Corrected problem with * calendar year printout in the output file. Calendar years were printing * as 1 greater than they should have been. My fix for this is kind of klunky * (doesn't fix the problem at the source, but corrects for it after the fact.) * * Version 4.1 (April 1998) - Added the GETWWWQ subroutine that can read * data from files in the rdb format used on USGS web pages. The first * version of GETWWWQ was written by Paul Weghorst of the Bureau of * Reclamation, for a Windows console version of BFI that he built using * Microsoft FORTRAN PowerStation. I made the following improvements: * * Removed non-standard FORTRAN syntax (Microsoft extensions) so that * the program can be used on Unix workstations and non-Microsoft compilers * * Reworked data parsing routines to handle data sets that have * flag character fields following the discharge values * * Modified leap-year handling to deal with century year exceptions * * Added date verifications as data are read. Web-formatted data files * don't actually indicate missing data, but rather just leave out lines * where data are missing and rely on user to verify that all data are * there. * * Modified to handle files that contain data from several gages * concatenated together into one file * * Modified the main program so that it automatically detects the * data file format and calls either GETQ or GETWWWQ as needed * * Updated to support data format provided by the NWIS-W software, * which varies slightly from the IL-SWR web retrieval software in * use at the time that GETWWWQ was first written. Differences are * in exact location of station name, gage id number, and date range * in header lines. BFI41 will handle both IL-SWR and NWIS-W formats. * * Added capability to auto-detect the format of dates and read file * accordingly. Both yyyy.mm.dd and mm/dd/yyyy options are now * available from USGS web pages. * * I also separated out the GETMINS routine so that it could be used by * both GETQ and GETWWWQ, and created an ADVANCEDATE routine that is * called from two different places in GETWWWQ. * * Finally, I fixed the bug in GETQ that was causing it to return the * wrong year (off by one) for data organized by calendar-year, and * cleaned up the klunky fix for this bug (in the main program) that * was used in version 4.02. * * Version 4.11 (January 2001) - Fixed problem in GETWWWQ subroutine that * caused erroneous results when reading a web-formatted input file * that was missing one or more entire years of data from the middle of * the file. Also modified main program to list missing years (all data * missing) in the output table. Previously they were listed as * incomplete if part of the data were present, but were not listed at * all if all data were missing. * * Version 4.12 (February 2001) - Added ability to read data files in the * newest NWIS-W format. Primary modifications in this version were in the * GETWWWQ routine. Added several variables to the GETWWWQ function call * from the main program. * * Version 4.13 (August 2003) - Fixed bug that prevented setting the * recession constant, 'K'. The code was still referring to it by an earlier * name 'R'. We changed 'R' to 'K' long ago in papers about the method. * * Version 4.14 (May 2006) - Fixed problem with reading web-format data files * that use the YYYYMMDD date format. * * Version 4.15 (April 2007) - Changed so that program recycles back to the * menu after completing calculations. This allows automating with * redirected input to process with varying options. This change does not * apply to the Windows version. **************************************************************************** * * DATA INPUT REQUIREMENTS AND NOTES * --------------------------------- * This program reads data from the standard WATSTORE 80-column card * format. It has also been tested with data obtained from Hydro-Data * CD-ROM disks that have been output using the card format. WATSTORE * files contain no decimal points, but the Hydro-Data files do include * a decimal point for discharges of 10 cfs or greater. In order to * properly read both data formats, an F7.0 format has been used. Also, * the Hydro-Data records use 999999 to fill positions of missing data and * to fill the unused data positions in the fourth card for each month. * The Hydro-Data files include H, Z, and N cards in addition to standard * 2 and 3 cards. The H and Z cards are ignored. If the N card (station name) * is present it is echoed to the output file. * * A base-flow index value can only be calculated for a complete water year. * (a single 2-card and 48 3-cards for each year). For single water years * some data from the start and end of the year will be unusable due to the * way the turning points are defined. For consecutive water years at the * same gage site, the program will interpolate through the break between * years, making use of as much of the record as possible. For those who * wish to use the program for portions of a year or for years with missing * data, the turning point and daily flow output files are available. * * To debug errors in data files, use the turning point and daily flow * output files to examine the data being read by the program. Missing * data will show as -99 in the output files. **************************************************************************** * * ARRAY DIMENSION NOTES * --------------------- * The program has been written to allow for base-flow separation using * flow periods other than 5 days. Changing the period to less than 5 days * requires larger arrays to hold the data for the 5-day minimum flows. * There are six arrays in the program and subroutines that are affected. * They are as follows: * * MAIN PROGRAM: SUBROUTINE GETQ: SUBROUTINE GETMINS: * DAY1MIN DAYMIN DAYMIN * DAY2MIN QMINS QMINS * QMINS1 * QMINS2 * * Each of these arrays must be dimensioned to (0:1+366/N), where N is the * length of the n-day period considered for each possible turning point. * Thus, for 5-day periods, these arrays would be dimensioned for (0:74). * In this version of the program, each of these arrays has been dimensioned * (0:367) to allow for values of N as small as 1. If you have difficulties * running the program due to memory occupied by these arrays, you should * modify the dimensioning of these arrays to match the requirements of the * lowest value of N that you wish to use. **************************************************************************** * * MAIN PROGRAM * ------------ * VARIABLES * --------- * A - character variable for answers to Y/N questions * ADDVOL - additional volume accumulated between successive turning points * BF2PRINT - daily base flow to print to daily flow output file * BFINDEX - base-flow index = (base flow / total flow) * BFCV - coefficient of variation of base-flow volume at specific gage * BFMEAN - mean of base-flow volume at specific gage * BFSDEV - standard deviation of base-flow volume at specific gage * BFSQR - sum of (base-flow volume)**2 at particular gage * BFSUM - sum of (base-flow volume) at particular gage * BFICV - coefficient of variation of base-flow index at specific gage * BFEXTRAP - string '!', indicates total base flow for the year was * extrapolated from base-flow index calculated with less than * full year's data (first or last turning point could not be * interpolated across end of year) * BFIMEAN - mean of base-flow index at specific gage * BFISDEV - standard deviation of base-flow index at specific gage * BFISQR - sum of (BFINDEX)**2 at particular gage * BFISUM - sum of BFINDEX's at particular gage * COMP1LETE - .TRUE. if year 1 data is complete * COMP2LETE - .TRUE. if year 2 data is complete * DAILYFILE - name of file for storage of daily flow and base-flow values * DATAFMT - string identifying format of input file * DATE1YR() - array containing year corresponding to each daily flow in year 1 * DATE1MTH() - array containing month for each daily flow in year 1 * DATE1DAY() - array containing day of the month for each daily flow in year 1 * DATE2YR() - array containing year corresponding to each daily flow in year 2 * DATE2MTH() - array containing month for each daily flow in year 2 * DATE2DAY() - array containing day of the month for each daily flow in year 2 * DAYFILE - unit number of output file for daily flows and base flows * DAYS - days between points used in interpolating turning points * DAY1MIN() - array of days of occurrence of each N-day minimum flow * for year 1 * DAY2MIN() - array of days of occurrence of each N-day minimum flow * for year 2 * DAYOFTP - day of the year for the current base-flow turning point * DAYS1YR - number of days in year 1 * DAYS2YR - number of days in year 2 * EOF - .TRUE. when end of input data file is reached * EOYQ - end-of-year base-flow (on day 365.5/366.5) * EOYDAY - end-of-year day (either 365.5 or 366.5, or -99) * FIRSTTP - day of first base-flow turning point of the year * FLAG - character variable to hold interpolation marker * FRAC - parameter (0.9) used in base-flow separation algorithm * FRSTPAS - .TRUE. if this is the first pass through main loop * GAGE1 - gage id number for year 1 * GAGE2 - gage id number for year 2 * I,J,II - counters * INFILE - unit number of input file * ISYM - symbol for interpolated turning points * LASTTP - day of the last base-flow turning point for the year * MDYSLASH - flag to indicate if date format IN GETWWWQ is mm/dd/yyyy * METHOD - integer indicating turning point test method (default=1). * (1=Standard Institute of Hydrology method; * 2=Modified method using 1-day recession constant) * M1SYM - Symbol used in FORMAT 56 statement to id chosen method * M2SYM - Symbol used in FORMAT 56 statement to id chosen method * N - N-day minimums parameter for base-flow separation (N=5 * in the standard implementation of this algorithm) * NCARD - character image of name card * NDAY1YR - Number of N-day periods per year (73 for N=5) in year 1 * NDAY2YR - Number of N-day periods per year (73 for N=5) in year 2 * NWISFMT - Format of NWIS-W data (explained in GETWWWQ) * NYEARS - number of years for which BFINDEX was calculated at each gage * OUTFILE - unit number of base-flow output file * PREVDAY - day of year on which previous turning point occurred * PREVQ - discharge at previous turning point * Q1() - array of mean daily discharges for year 1 * Q2() - array of mean daily discharges for year 2 * QMINS1() - array of N-day minimum flows for year 1 * QMINS2() - array of N-day minimum flows for year 2 * QQ1 - Q at first point for turning point interpolation * QQ3 - Q at second point for turning point interpolation * QTP - Q at a given turning point * QCV - coefficient of variation of runoff volume at specific gage * QMEAN - mean of runoff volume at specific gage * QSDEV - standard deviation of runoff volume at specific gage * QSQR - sum of (runoff volume)**2 at particular gage * QSUM - sum of (runoff volume) at particular gage * K1DAY - 1-day recession constant * PAUSE - variable to hold user response to pause prompt * RVARNAME - 1-character string holding name of the recession slope * parameter corresponding to chosen method ('f' or 'K') * RVAR - value of the recession slope parameter to be printed in * output file (either FRAC or K1DAY) * SEFLAG - integer flag to tell whether the complete year of data * was used to calculate base flow, including SOY and EOY points * SOYQ - start-of-year Q (instantaneous Q at t=0.5 days) * SOYDAY - start-of-year day (0.5 or -99) * THISLINE - string used to read first line of data file to determine format * TPFILE - unit number of turning point output file * TPQ1 - interpolated base flow on day 1 of water year (October 1) * TURNPT - .TRUE. if point is a turning point * TYPEF - holds interpolation marker for year's first turning point * TYPEL - holds interpolation marker for year's last turning point * VOLBF - volume of base flow * VOLQ - volume of total runoff * WATYR1 - water year for year 1 data * WATYR2 - water year for year 2 data *WWWGAGEFOUND- flag to indicate that GETWWWQ has already been called * once and a gage header was found. Second time it's called * it won't look for a gage header. * WWWDAY - last day read by GETWWWQ * WWWMONTH - last month read by GETWWWQ * WWWYEAR - last year read by GETWWWQ * XSYM - symbol for years where annual base flow is extrapolated * YEAR1BASIS - string, 'Calendar' or 'Water ' year --- year 1 * YEAR2BASIS - string, 'Calendar' or 'Water ' year --- year 2 * YYYYDOT - flag to indicate if date format in GETWWWQ is yyyy.mm.dd * YYYYDASH - flag to indicate if date format is yyyy-mm-dd * YYYYMMDD - flag to indicate if date format is yyyymmdd * MDYSLASH - flag to indicate if date format is mm/dd/yyyy * ZERO1 - holds no-flow symbol for year 1 * ZERO2 - holds no-flow symbol for year 2 * ZSYM - symbol for years with no-flow days * INPUTFILE - file name of input file (filenames limited to 255 characters) * OUTPUTFILE - file name of base-flow output file * TPOUTFILE - file name of turning point output file * * SUBROUTINES and FUNCTIONS (# INDICATES THOSE CALLED DIRECTLY BY MAIN PROGRAM) * ----------------------------------------------------------------------------- * # BASEFLOW - COMPUTES BASE-FLOW VOLUME FROM ONE TURNING POINT TO * THE NEXT * * # GETQ - READS ONE YEAR'S DAILY AVERAGE DISCHARGE DATA FROM * INPUT FILE SPECIFIED BY THE USER, AND CALLS GETMINS * * # GETWWWQ - READS ONE YEAR'S DATA FROM A FILE IN THE rdb FORMAT * USED ON USGS WORLD-WIDE-WEB SITES * * GETMINS - DETERMINES MINIMUM FLOWS IN EACH N-DAY PERIOD FOR * THE YEAR * * ADVANCEDATE - INCREMENTS CALENDAR DAY BY ONE FOR KEEPING TRACK * OF NEXT EXPECTED PIECE OF DATA * * READCARD - READS ONE CARD (LINE) FROM THE SPECIFIED INPUT DATA * FILE. THIS ROUTINE CHECKS FOR MISSING OR ERRONEOUS * DATA. * * # STATS - COMPUTES MEAN, STANDARD DEVIATION, AND COEFFICIENT * OF VARIATION FOR BASE-FLOW INDEX AT PARTICULAR GAGE * * # TEST_FOR_TP - PERFORMS LOGICAL TESTS TO DETERMINE IF A GIVEN DAY * IS A BASE-FLOW TURNING POINT, GIVEN 3 CONSECUTIVE * N-DAY MINIMUM FLOWS. * * # GETOPTIONS - SUBROUTINE TO GET BFI PARAMETERS FROM THE USER. * * # QINTERP - FUNCTION TO PERFORM STRAIGHT-LINE INTERPOLATION BETWEEN * TWO BASE-FLOW TURNING POINTS PLOTTED ON SEMI-LOG SCALES. * * UCASE - FUNCTION TO CONVERT FIRST CHARACTER OF A STRING TO * UPPERCASE. * * # EXISTS - LOGICAL FUNCTION TO INDICATE IF A FILE EXISTS * * # FILEEXISTS - SUBROUTINE TO DEAL WITH PROBLEMS PRESENTED WHEN ONE OF THE * OUTPUT FILES ALREADY EXISTS * **************************************************************************** PROGRAM BFI CHARACTER*1 ISYM, FLAG, TYPEF, TYPEL, ZERO1, ZERO2, ZSYM CHARACTER*1 BFEXTRAP,XSYM,M1SYM,M2SYM, PAUSE CHARACTER*1 RVARNAME CHARACTER*8 GAGE1, GAGE2, YEAR1BASIS, YEAR2BASIS CHARACTER*255 INPUTFILE, OUTPUTFILE, TPOUTFILE, DAILYFILE CHARACTER*78 NCARD CHARACTER*40 DATAFMT CHARACTER*80 THISLINE LOGICAL COMP1LETE, COMP2LETE, EOF, FRSTPAS, TURNPT,EXISTS LOGICAL WWWGAGEFOUND LOGICAL YYYYDOT, YYYYDASH, YYYYMMDD, MDYSLASH INTEGER DAYFILE, DAYS, DAYOFTP, DAYS1YR, DAYS2YR INTEGER DAY1MIN(0:367), DAY2MIN(0:367) INTEGER FIRSTTP, I, II, J, LASTTP, N, NDAY1YR, NDAY2YR INTEGER NYEARS, PREVDAY, WATYR1, WATYR2 INTEGER INFILE, OUTFILE, TPFILE, METHOD,SEFLAG INTEGER DATE1YR(366), DATE1MTH(366), DATE1DAY(366) INTEGER DATE2YR(366), DATE2MTH(366), DATE2DAY(366) INTEGER WWWDAY, WWWMONTH, WWWYEAR INTEGER NWISFMT REAL ADDVOL, BFINDEX, BFICV, BFIMEAN, BFISDEV,EOYQ,FRAC REAL SOYQ,SOYDAY,PREVQ, QMINS1(0:367), QMINS2(0:367) REAL BFISUM, BFISQR, QQ1, QQ3, QTP, TPQ1, VOLBF, VOLQ REAL Q1(366), Q2(366), EOYDAY REAL BFCV, BFMEAN, BFSDEV, BFSUM, BFSQR REAL QCV, QMEAN, QSDEV, QSUM, QSQR REAL BF2PRINT REAL K1DAY, RVAR *---- DECLARE FUNCTION DATA TYPES REAL QINTERP *---- GET INPUT FILE NAME AND CALL GETOPTIONS *---- These unit numbers are for the Data General and MS-DOS machines INFILE = 1 OUTFILE = 2 TPFILE=3 DAYFILE=4 *---- THESE UNIT NUMBERS ARE FOR THE PRIME C INFILE = 5 C OUTFILE = 6 C TPFILE = 7 C DAYFILE = 8 WRITE (*,3) 3 FORMAT (///////////////////////// 6X, +'________________________________' //6X, +' BFI ' //6X, +'A Computer Program for Computing' / 6X, +' an Index to Base Flow ' //6X, +' Version 4.15 - April 2007 ' / 6X, +'________________________________' /////) 4 WRITE (*,*) 'Enter name of input file:' READ (*,5) INPUTFILE 5 FORMAT (A255) IF (.NOT.EXISTS(INPUTFILE)) THEN WRITE (*,*) '*** ERROR with file: ',INPUTFILE(1:68) WRITE (*,*) 'The file does not exist.' WRITE (*,*) GOTO 4 ENDIF 10 CALL GETOPTIONS (INPUTFILE,OUTPUTFILE,TPOUTFILE,DAILYFILE, + ZSYM,ISYM,XSYM,METHOD,N,FRAC,K1DAY,TPFILE,DAYFILE) OPEN (INFILE,FILE=INPUTFILE,STATUS='OLD') IF (EXISTS(OUTPUTFILE)) CALL FILEEXISTS(OUTPUTFILE) IF (EXISTS(OUTPUTFILE)) THEN OPEN (OUTFILE,FILE=OUTPUTFILE,STATUS='OLD') ELSE OPEN (OUTFILE,FILE=OUTPUTFILE,STATUS='NEW') ENDIF IF (TPFILE.NE.0) THEN IF (EXISTS(TPOUTFILE)) CALL FILEEXISTS(TPOUTFILE) IF (EXISTS(TPOUTFILE)) THEN OPEN (TPFILE,FILE=TPOUTFILE,STATUS='OLD') ELSE OPEN (TPFILE,FILE=TPOUTFILE,STATUS='NEW') ENDIF ELSE TPOUTFILE=' ' ENDIF IF (DAYFILE.NE.0) THEN IF (EXISTS(DAILYFILE)) CALL FILEEXISTS(DAILYFILE) IF (EXISTS(DAILYFILE)) THEN OPEN (DAYFILE,FILE=DAILYFILE,STATUS='OLD') ELSE OPEN (DAYFILE,FILE=DAILYFILE,STATUS='NEW') ENDIF ELSE DAILYFILE=' ' ENDIF *---- DETERMINE WHETHER DATA ARE IN WATSTORE OR WEB/rdb FORMAT *---- AN H, Z, N, 2, OR 3 CARD INDICATES WATSTORE FORMAT *---- A '#' INDICATES AN RDB HEADER, OR A '<' INDICATES AN *---- HTML TAG AT THE TOP OF AN RDB HEADER 37 READ (INFILE,38,END=320) THISLINE 38 FORMAT (A80) IF (THISLINE(1:1) .EQ. 'H' .OR. + THISLINE(1:1) .EQ. 'Z' .OR. + THISLINE(1:1) .EQ. 'N' .OR. + THISLINE(1:1) .EQ. '2' .OR. + THISLINE(1:1) .EQ. '3' ) THEN DATAFMT = 'WATSTORE' BACKSPACE (INFILE) ELSEIF (THISLINE(1:1) .EQ. '#' .OR. + THISLINE(1:1) .EQ. '<' ) THEN DATAFMT = 'Web/rdb (NWIS-W)' BACKSPACE (INFILE) ELSE GOTO 37 *-------- KEEP LOOPING UNTIL WE FIGURE OUT WHAT KIND OF FILE IT IS ENDIF WRITE (OUTFILE,40) INPUTFILE, DATAFMT, OUTPUTFILE, TPOUTFILE, + DAILYFILE 40 FORMAT (33X,'Input file = ',T47,A30/, + 32X,'File format = ',T47,A40/, + 22X,'Base-flow output file = ',T47,A30/, + 18X,'Turning point output file = ',T47,A30/, + 1X,'Daily base flow and total flow output file = ',T47,A30/) 53 IF (METHOD.EQ.1) THEN M1SYM='*' M2SYM=' ' RVARNAME='f' RVAR=FRAC ELSE M2SYM='*' M1SYM=' ' RVARNAME='K' RVAR=K1DAY ENDIF WRITE (OUTFILE,56) M1SYM,M2SYM,METHOD,N,RVARNAME,RVAR IF (TPFILE.NE.0) WRITE (TPFILE,56) M1SYM,M2SYM,METHOD,N, + RVARNAME,RVAR IF (DAYFILE.NE.0) WRITE (DAYFILE,56) M1SYM,M2SYM,METHOD,N, + RVARNAME,RVAR 56 FORMAT (1X,'Program Version = BFI 4.15'//, + 1X,'AVAILABLE SEPARATION METHODS:'/, + 1X,A1,' 1 = STANDARD Institute of Hydrology method'/, + 1X,' (N-day avg. recession test; uses "N" and "f")'/, + 1X,A1,' 2 = MODIFIED method'/, + 1X,' (1-day recession constant adjusted for number of ', + 'days'/, + 1X,' between points; uses "N" and "K")'//, + 1X,'BASE-FLOW SEPARATION PARAMETERS'/ + 3X,'METHOD = ',I2/ + 3X,'N = ',I2/ + 3X,A1,' = ',F10.6) NYEARS = 0 BFISUM = 0. BFISQR = 0. BFSUM = 0. BFSQR = 0. QSUM = 0. QSQR = 0. GAGE1 = ' ' FRSTPAS = .TRUE. EOF = .FALSE. *---- READ ONE YEAR'S DATA *---- SEND THE NO FLOW SYMBOL TO GETQ USING PARAMETER ZERO2 60 ZERO2 = ZSYM IF (DATAFMT .EQ. 'WATSTORE') THEN CALL GETQ (Q2,QMINS2,DAY2MIN,WATYR2,GAGE2,NCARD,ZERO2,EOF, + COMP2LETE,DAYS2YR,DATE2YR,DATE2MTH,DATE2DAY,YEAR2BASIS, + INFILE,N) ELSE *-------- DATAFMT IS 'Web/rdb ' CALL GETWWWQ (Q2,QMINS2,DAY2MIN,WATYR2,GAGE2,NCARD,ZERO2,EOF, + COMP2LETE,DAYS2YR,DATE2YR,DATE2MTH,DATE2DAY,YEAR2BASIS, + INFILE,N,WWWGAGEFOUND,WWWDAY,WWWMONTH,WWWYEAR,YYYYDOT, + YYYYDASH,YYYYMMDD,MDYSLASH,NWISFMT) ENDIF *---- CALCULATE THE NUMBER OF N-DAY PERIODS IN THE YEAR OF DATA JUST READ. *---- IF THERE IS A FRACTION OF AN N-DAY PERIOD LEFT AT THE END OF A YEAR *---- IT IS LUMPED IN WITH THE LAST N-DAY PERIOD. NDAY2YR = DAYS2YR / N IF (FRSTPAS) GOTO 160 *---- DETERMINE TURNING POINTS AND COMPUTE VOLUME OF BASE FLOW BETWEEN *---- SUCCESSIVE TURNING POINTS FOR DATA IN ARRAYS FOR YEAR 1 VOLBF = 0. FIRSTTP = DAYS1YR LASTTP = 0 *---- TACK YEAR 2's FIRST QMIN ONTO END OF YEAR 1 IF ((WATYR2.EQ.WATYR1+1).AND.(GAGE2.EQ.GAGE1)) THEN QMINS1(NDAY1YR+1) = QMINS2(1) DAY1MIN(NDAY1YR+1) = DAY2MIN(1) ELSE QMINS1(NDAY1YR+1) = -99 ENDIF *---- IF AN INTERPOLATION WAS PERFORMED THE LAST TIME THROUGH, RECALL THE *---- RESULT AS THE FIRST TURNING POINT, AND RETRIEVE THE END-OF-YEAR Q *---- AS THIS YEAR'S START-OF-YEAR Q TYPEF = ' ' SOYQ = -99. SOYDAY = -99. SEFLAG = 0 IF (TPQ1.GE.0.) THEN FIRSTTP=1 PREVDAY=1 PREVQ=TPQ1 FLAG=ISYM TYPEF=ISYM SOYQ=EOYQ SOYDAY=0.5 ENDIF *---- LOCATE SUCCESSIVE TURNING POINTS THROUGHOUT THE YEAR DO 80 I=1,DAYS1YR/N TURNPT = .FALSE. CALL TEST_FOR_TP(FRAC,QMINS1(I-1),QMINS1(I),QMINS1(I+1), + DAY1MIN(I-1),DAY1MIN(I),DAY1MIN(I+1),K1DAY,METHOD,TURNPT) IF (TURNPT) THEN FIRSTTP = MIN(FIRSTTP,DAY1MIN(I)) LASTTP = MAX(LASTTP,DAY1MIN(I)) DAYOFTP = DAY1MIN(I) QTP = QMINS1(I) *-------- IF THE FIRST TURNING POINT WAS INTERPOLATED, AND A TURNING POINT *-------- WAS NOT SUBSEQUENTLY FOUND ON DAY 1 OF THE YEAR, THEN OUTPUT THE *-------- INTERPOLATED TURNING POINT WITH THE INTERPOLATION FLAG. OTHERWISE *-------- RESET TYPEF TO ' ' AND DON'T PRINT THE TURNING POINT. IT WILL GET *-------- PRINTED LATER. IF (TPQ1.GE.0) THEN IF (DAY1MIN(I).EQ.1) THEN TYPEF = ' ' ELSE IF (TPFILE.NE.0) WRITE (TPFILE,70) + DATE1YR(PREVDAY),DATE1MTH(PREVDAY), + DATE1DAY(PREVDAY),PREVQ,FLAG ENDIF TPQ1 = -99 FLAG = ' ' ENDIF *-------- WRITE THE TURNING POINT TO THE TURNING POINT OUTPUT FILE IF (TPFILE.NE.0) WRITE (TPFILE,70)DATE1YR(DAYOFTP), + DATE1MTH(DAYOFTP),DATE1DAY(DAYOFTP),QTP,FLAG 70 FORMAT (I4,I6,I6,F10.2,A2) *-------- WRITE THE INDIVIDUAL DAILY BASE FLOWS AND TOTAL FLOWS *-------- TO THE DAILY FLOW OUTPUT FILE. IF ((DAYFILE.NE.0).AND.(DAYOFTP.NE.FIRSTTP)) THEN DO 75, II=PREVDAY,DAYOFTP-1 BF2PRINT = QINTERP(PREVDAY,PREVQ,DAYOFTP,QTP,FLOAT(II)) IF (BF2PRINT.GT.Q1(II)) BF2PRINT=Q1(II) WRITE (DAYFILE,77)DATE1YR(II),DATE1MTH(II), + DATE1DAY(II),BF2PRINT,Q1(II) 75 CONTINUE 77 FORMAT (I4,I6,I6,F10.2,F12.2) ENDIF *-------- IF YEAR 1 IS COMPLETE THEN CALCULATE BASE FLOW IF (COMP1LETE) THEN *---------- IF THERE WAS AN INTERPOLATION AT THE START OF THIS YEAR, *---------- CALCULATE THE BASE FLOW FOR THE FIRST HALF-DAY *---------- OF THE YEAR (FROM SOYDAY TO DAY 1, WHICH IS CURRENTLY PREVDAY) IF (SOYQ.NE.-99.) THEN CALL BASEFLOW (SOYDAY,FLOAT(PREVDAY),SOYQ,PREVQ,ADDVOL) VOLBF = VOLBF + ADDVOL SEFLAG = SEFLAG + 1 SOYQ = -99. ENDIF IF ((DAYOFTP.NE.FIRSTTP).AND.(DAYOFTP.NE.PREVDAY)) THEN CALL BASEFLOW (FLOAT(PREVDAY),FLOAT(DAYOFTP),PREVQ,QTP, + ADDVOL) VOLBF = VOLBF + ADDVOL ENDIF ENDIF PREVQ = QTP PREVDAY = DAYOFTP ENDIF 80 CONTINUE *---- IF LAST TURNING POINT IN YEAR 1 WAS VALID AND GAGES AND WATER YEARS *---- MATCH, INTERPOLATE TURNING POINTS FOR LAST DAY OF YEAR 1 AND FIRST *---- DAY OF YEAR 2 TPQ1 = -99 TYPEL = ' ' EOYQ = -99. EOYDAY = -99. IF ((GAGE1.EQ.GAGE2).AND.(WATYR1+1.EQ.WATYR2).AND.(PREVQ.GE. + 0)) THEN QMINS2(0)=QMINS1(NDAY1YR) DAY2MIN(0)=DAY1MIN(NDAY1YR) *------ SEARCH THROUGH YEAR 2 DATA TO FIND FIRST TURNING POINT DO 90 J=1,(DAYS2YR/N)-1 TURNPT = .FALSE. CALL TEST_FOR_TP(FRAC,QMINS2(J-1),QMINS2(J),QMINS2(J+1), + DAY2MIN(J-1),DAY2MIN(J),DAY2MIN(J+1),K1DAY,METHOD,TURNPT) IF (TURNPT) THEN IF (QMINS2(J).GE.0.) THEN TYPEL=ISYM *------------ INTERPOLATE DAYS=DAY2MIN(J)-(PREVDAY-DAYS1YR) QQ1=PREVQ QQ3=QMINS2(J) QTP=QINTERP(0,PREVQ,DAYS,QMINS2(J),FLOAT(DAYS1YR-PREVDAY)) TPQ1=QINTERP(0,PREVQ,DAYS,QMINS2(J), + FLOAT(DAYS1YR+1-PREVDAY)) EOYQ=QINTERP(0,PREVQ,DAYS,QMINS2(J), + FLOAT(DAYS1YR)+.5+FLOAT(PREVDAY)) EOYDAY=DAYS1YR+.5 *------------ CALCULATE LAST BASE-FLOW INCREMENT FOR YEAR 1 DAYOFTP= DAYS1YR LASTTP=DAYS1YR CALL BASEFLOW (FLOAT(PREVDAY),EOYDAY,PREVQ,EOYQ,ADDVOL) SEFLAG=SEFLAG+1 VOLBF = VOLBF + ADDVOL IF (DAYOFTP.NE.PREVDAY) THEN IF (TPFILE.NE.0) WRITE (TPFILE,70) + DATE1YR(DAYOFTP),DATE1MTH(DAYOFTP), + DATE1DAY(DAYOFTP),QTP,TYPEL IF (DAYFILE.NE.0) THEN DO 85, II=PREVDAY,DAYOFTP-1 BF2PRINT = QINTERP(PREVDAY,PREVQ,DAYOFTP,QTP, + FLOAT(II)) IF (BF2PRINT.GT.Q1(II)) BF2PRINT=Q1(II) WRITE (DAYFILE,77)DATE1YR(II),DATE1MTH(II), + DATE1DAY(II),BF2PRINT,Q1(II) 85 CONTINUE ENDIF ELSE TYPEL=' ' ENDIF ENDIF GOTO 100 ENDIF 90 CONTINUE 100 ENDIF *---- PRINT DAILY FLOW VALUES FOR LAST TURNING POINT IN YEAR 1 BF2PRINT=MIN(QTP,Q1(LASTTP)) IF (DAYFILE.NE.0) WRITE (DAYFILE,77)DATE1YR(LASTTP), + DATE1MTH(LASTTP),DATE1DAY(LASTTP),QTP,Q1(LASTTP) *---- DETERMINE VOLUME OF TOTAL FLOW BETWEEN FIRST TURNING POINT AND *---- LAST TURNING POINT (FIRSTTP,LASTTP) THIS VOLUME IS IN CFS-DAYS IF (COMP1LETE) THEN VOLQ=0. DO 110 I=FIRSTTP,LASTTP IF (((I.EQ.FIRSTTP).AND.(SOYDAY.EQ.-99)).OR. + ((I.EQ.LASTTP ).AND.(EOYDAY.EQ.-99))) THEN VOLQ=VOLQ+(Q1(I)/2.) ELSE VOLQ=VOLQ+Q1(I) ENDIF 110 CONTINUE *---- CALCULATE BASE-FLOW INDEX IF (VOLQ.EQ.0) THEN BFINDEX = 0 ELSE BFINDEX = VOLBF/VOLQ ENDIF BFISUM = BFISUM + BFINDEX BFISQR = BFISQR + (BFINDEX**2) NYEARS = NYEARS + 1 *---- COMPUTE VOLUME OF RUNOFF AND BASE FLOW FOR THE FULL YEAR VOLQ = 0. DO 120 I=1,DAYS1YR VOLQ=VOLQ+Q1(I) 120 CONTINUE IF ((SEFLAG.LT.2).AND.(XSYM.EQ.'!').AND.((FIRSTTP.GT.1).OR. + (LASTTP.LT.DAYS1YR))) THEN BFEXTRAP=XSYM ELSE BFEXTRAP=' ' ENDIF VOLQ = VOLQ * (86400./43560.) VOLBF = VOLQ * BFINDEX QSUM = QSUM + VOLQ QSQR = QSQR + (VOLQ**2) BFSUM = BFSUM + VOLBF BFSQR = BFSQR + (VOLBF**2) *-------- OUTPUT BASE-FLOW RESULTS TO SCREEN AND OUTPUT FILE WRITE (*,130)WATYR1,BFINDEX,ZERO1,VOLBF,BFEXTRAP,VOLQ, + FIRSTTP,TYPEF,LASTTP,TYPEL WRITE (OUTFILE,130)WATYR1,BFINDEX,ZERO1,VOLBF,BFEXTRAP, + VOLQ,FIRSTTP,TYPEF,LASTTP,TYPEL 130 FORMAT (1X,I4,F14.3,A2,F14.0,A2,F14.0,I10,A2,I10,A2) *------ INCOMPLETE YEAR ELSE WRITE (*,140) WATYR1 WRITE (OUTFILE,140) WATYR1 140 FORMAT (1X,I4,5X,'Incomplete year. Base flow cannot be ', + 'determined.') ENDIF *---- IF GAGE HAS CHANGED SINCE LAST TIME, OR END OF FILE WAS REACHED, *---- CALL STATS AND RESET ACCUMULATORS IF ((GAGE2.NE.GAGE1).OR.(YEAR2BASIS.NE.YEAR1BASIS).OR.(EOF)) THEN IF (NYEARS.GE.5) THEN CALL STATS(NYEARS,BFISUM,BFISQR, BFIMEAN, BFISDEV, BFICV) CALL STATS(NYEARS,BFSUM,BFSQR,BFMEAN,BFSDEV,BFCV) CALL STATS(NYEARS,QSUM,QSQR,QMEAN,QSDEV,QCV) WRITE (*,150)NYEARS,YEAR1BASIS,GAGE1,BFMEAN,BFSDEV,BFCV, + QMEAN,QSDEV,QCV,BFIMEAN,BFISDEV,BFICV WRITE (OUTFILE,150)NYEARS,YEAR1BASIS,GAGE1,BFMEAN,BFSDEV, + BFCV,QMEAN,QSDEV,QCV,BFIMEAN,BFISDEV,BFICV 150 FORMAT (/1X,'Statistics for ',I3,A9, ' years at gage ',A8, + /1X,'--------------------------------------------------' + /1X,' ', + 'STANDARD COEFFICIENT'/ + 1X,' MEAN DEVIATION', + ' OF VARIATION'/ + 1X,'-------------------------------------------------', + '--------------------'/ + 1X,'BASE FLOW (AC-FT) ',F10.1,6X,F10.1,6X,F10.3/ + 1X,'TOTAL RUNOFF (AC-FT) ',F10.1,6X,F10.1,6X,F10.3/ + 1X,'BASE-FLOW INDEX ',F10.3,6X,F10.3,6X,F10.3) ENDIF IF (ZSYM.NE.' ') THEN WRITE (*,155) ZSYM,ISYM,XSYM WRITE (OUTFILE,155) ZSYM,ISYM,XSYM ENDIF 155 FORMAT (/1X,A1,' - Denotes years with one or more zero flow ', + 'days',/,1X,A1,' - Indicates interpolated turning point', + /,1X,A1,' - Indicates whole year could not be analyzed, but', + /,1X,' total annual base flow is extrapolated using the', + /,1X,' calculated Base-Flow Index and total annual runoff') NYEARS = 0 BFISUM = 0. BFISQR = 0. BFSUM = 0. BFSQR = 0. QSUM = 0. QSQR = 0. ENDIF 160 FRSTPAS = .FALSE. *---- MOVE YEAR 2 DATA INTO ARRAYS FOR YEAR 1 DO 170 J=1,DAYS2YR Q1(J)=Q2(J) DATE1YR(J)=DATE2YR(J) DATE1MTH(J)=DATE2MTH(J) DATE1DAY(J)=DATE2DAY(J) 170 CONTINUE *---- TACK END OF YEAR ON AS START OF NEXT YEAR IF ((WATYR2.EQ.WATYR1+1).AND.(GAGE2.EQ.GAGE1)) THEN QMINS2(0)=QMINS1(NDAY1YR) DAY2MIN(0)=DAY1MIN(NDAY1YR) ELSE QMINS2(0) = -99 TPQ1 = -99 IF ((GAGE2.EQ.GAGE1).AND.(WATYR2.GT.WATYR1+1)) THEN DO 175 I=WATYR1+1,WATYR2-1 WRITE (*,173) I WRITE (OUTFILE,173) I 173 FORMAT (1X,I4,5X,'Missing year.') 175 CONTINUE ENDIF ENDIF DO 180 I=0,(DAYS2YR/N) QMINS1(I)=QMINS2(I) DAY1MIN(I)=DAY2MIN(I) 180 CONTINUE WATYR1=WATYR2 ZERO1=ZERO2 COMP1LETE=COMP2LETE DAYS1YR=DAYS2YR NDAY1YR=NDAY2YR *---- IF THIS IS A NEW GAGE, PRINT NEW HEADER LINES IF ((GAGE2.NE.GAGE1).OR.(YEAR2BASIS.NE.YEAR1BASIS)) THEN IF (TPFILE.NE.0) THEN WRITE (TPFILE,190)GAGE2 190 FORMAT (/'==============================================='/ + /'Gage ',A8) IF (ISYM.NE.' ') WRITE (TPFILE,193) ISYM 193 FORMAT ('(',A1,' indicates interpolated turning point)') WRITE (TPFILE,195) 195 FORMAT (/'<-- Calendar --> Base Flow',/, + 'Year Month Day (cfs) ',/, + '---------------------------') ENDIF IF (DAYFILE.NE.0) THEN WRITE (DAYFILE,197)GAGE2 197 FORMAT (/'=================================================' + //,'Gage ',A8,//, + '<-- Calendar --> Base Flow Total Flow',/, + 'Year Month Day (cfs) (cfs) ',/, + '---------------------------------------') ENDIF WRITE (*,200)GAGE2,NCARD,YEAR2BASIS WRITE (OUTFILE,200)GAGE2,NCARD,YEAR2BASIS 200 FORMAT (/////// + 1X,'====================================================', + '========================='/ + 1X,'Base-Flow Index for gage ',A8/,1X,A78/, + 1X,A8, ' Base-Flow Base Flow Total Runoff', + ' | Day of Turning Point |',/, + 1X,'Year Index (acre-ft) (acre-ft) ', + ' | [First] [Last] |',/, + 1X,'----------------------------------------------------', + '-------------------------') ENDIF GAGE1=GAGE2 YEAR1BASIS=YEAR2BASIS *---- IF EOF HAS BEEN REACHED, QUIT, ELSE GO BACK TO 60 IF (.NOT.EOF) GOTO 60 300 CLOSE (INFILE) CLOSE (OUTFILE) IF (TPFILE.NE.0) CLOSE (TPFILE) IF (DAYFILE.NE.0) CLOSE (DAYFILE) WRITE (*,*) WRITE (*,*) 'Press "Enter" to continue, or "Q" to quit.' READ (*,310) PAUSE 310 FORMAT (A1) IF ((.NOT.(PAUSE.EQ.'Q')).AND.(.NOT.(PAUSE.EQ.'q'))) GOTO 10 STOP 320 WRITE (INFILE,*) '*** ERROR: Could not identify input file type.' WRITE (INFILE,*) ' for file: ', INPUTFILE WRITE (*,*) '*** ERROR: Could not identify input file type.' WRITE (*,*) ' for file: ', INPUTFILE GOTO 300 END SUBROUTINE STATS (N, SUM, SUMSQR, XBAR, SDEV, CV) ***************************************************************** * COMPUTE MEAN, STANDARD DEVIATION, AND COEFFICIENT OF VARIATION * * VARIABLES * CV - coefficient of variation of x * N - number of data items * SDEV - standard deviation of x * SUM - sum of x * SUMSQR - sum of x**2 * XBAR - mean of x ***************************************************************** INTEGER N REAL CV, SDEV, SUM, SUMSQR, XBAR *--------------------------------------------------------------------------- XBAR = SUM/FLOAT(N) SDEV = SQRT((SUMSQR-((SUM**2)/FLOAT(N)))/FLOAT(N-1)) IF (XBAR.EQ.0) THEN CV = 0 ELSE CV = SDEV/XBAR ENDIF RETURN END SUBROUTINE GETQ (Q,QMINS,DAYMIN,WATERYR,GAGEID,NCARD,ZERO,EOF, + COMPLETE,DAYSNYR,DATEYR,DATEMTH,DATEDAY,YEARBASIS,INFILE,N) ***************************************************************************** * READ ONE YEAR'S AVERAGE DAILY FLOWS FROM THE INPUT FILE IN WATSTORE FORMAT, * AND, IF SUCCESSFUL, CALL GETMINS TO DETERMINE N-DAY MINIMUMS * * VARIABLES * ACTMONTH - actual month (January=1, December=12) * BOGUSYR - bogus parameter to get back from READCARD when I * want to discard the water year value (only really need * and want it at the start...it can be off by one after that) * CARDTYPE - leftmost character of data card * COMPLETE - .TRUE. if no missing data * DATEYR() - array containing year corresponding to each daily flow * DATEMTH() - array containing month for each daily flow * DATEDAY() - array containing day of the month for each daily flow * DAYS() - array holding number of data items to be expected on the * fourth card of each month * DAYSCY() - same as DAYS(), but on a calendar-year basis * DAYMIN() - array of days of year on which N-day minimums occur * DAYSNYR - number of days in current year * DOY - day of the year * EOF - .TRUE. if end of file was reached, before any data was read * EOFFLAG - end of file flag returned from subroutine READCARD * FEBRUARY - month number for February (5 for water-year basis; otherwise 2) * FLOW() - array holding discharges returned by READCARD * GAGEID - gage identification number from the 2 card * GAGE - gage id number from the 3 card (should match GAGEID) * I,K,L - counters * INFILE - unit number of input file * MON - DO loop counter for months (October=1,September=12) * MONTHREAD - month value read from 3 card by READCARD * MTHOFWY - month of the water year (October=1; September=12) * N - parameter for base-flow separation * NCARD - string image of name card * Q() - array of mean daily discharges * QMINS() - array of N-day minimum flows * READDAYS - number of days data on current input card, sent to READCARD * WATERYR - water year for current year's data * WEEK - "week" value read from the 3-card by READCARD * YEARBASIS - string, 'Calendar' or 'Water ' year * YRREAD - actual year read from the 3-card by READCARD * ZSYM - symbol for years with no-flow days * ZERO - no-flow days marker * * SUBROUTINES USED: * ----------------- * READCARD, GETMINS **************************************************************************** CHARACTER*1 ZERO, CARDTYPE, ZSYM CHARACTER*8 GAGE,GAGEID,YEARBASIS CHARACTER*78 NCARD LOGICAL COMPLETE, EOF, EOFFLAG INTEGER DAYMIN(0:367), DAYS(12), DAYSCY(12), DAYSNYR, DOY INTEGER I, K, L, MON, MONTHREAD, MTHOFWY INTEGER READDAYS, WATERYR, WEEK, YRREAD, FEBRUARY INTEGER BOGUSYR INTEGER INFILE, N, DATEYR(366), DATEMTH(366), DATEDAY(366) REAL FLOW(8), Q(366), QMINS(0:367) *---- DEFINE DAYS IN 4TH LINE OF DATA FOR EACH MONTH DATA DAYS /7,6,7,7,4,7,6,7,6,7,7,6/ *---- NOTE: WATER YEAR STARTS WITH OCTOBER = MONTH 1 *---- DEFINE DAYS IN 4TH LINE OF DATA FOR EACH MONTH DATA DAYSCY /7,4,7,6,7,6,7,7,6,7,6,7/ *---- NOTE: CALENDAR YEAR STARTS WITH JANUARY = MONTH 1 *--------------------------------------------------------------------------- ZSYM=ZERO DO 400 I=1,366 Q(I) = -99 400 CONTINUE NCARD = ' ' ZERO = ' ' DOY = 0 READDAYS = 8 DAYSNYR = 365 *---- READ THE 2 CARD FOR THIS YEAR'S DATA 410 CALL READCARD (CARDTYPE,NCARD,COMPLETE,EOFFLAG,READDAYS,GAGEID, + WATERYR,FLOW,YRREAD,MONTHREAD,WEEK,INFILE) IF (EOFFLAG) GOTO 445 IF (CARDTYPE.NE.'2') GOTO 410 *---- READ FIRST 3 CARD TO DETERMINE CALENDAR OR WATER YEAR BASIS CALL READCARD (CARDTYPE,NCARD,COMPLETE,EOFFLAG,READDAYS,GAGEID, + WATERYR,FLOW,YRREAD,MONTHREAD,WEEK,INFILE) IF (EOFFLAG) GOTO 445 BACKSPACE (INFILE) IF (MONTHREAD.EQ.1) THEN YEARBASIS='Calendar' ELSE YEARBASIS='Water ' ENDIF IF (YEARBASIS.EQ.'Water ') THEN FEBRUARY=5 ELSE FEBRUARY=2 DO 411 I=1,12 DAYS(I)=DAYSCY(I) 411 CONTINUE ENDIF *---- READ DAILY DISCHARGES. COMPLETE = .TRUE. DO 440 MON=1,12 DO 430 K=1,4 READDAYS=8 IF (K.EQ.4) THEN *---------- DETERMINE DAYS IN 4th LINE FOR EACH MONTH IF (MON.EQ.FEBRUARY) THEN IF (((MOD(WATERYR,4).EQ.0).AND.(MOD(WATERYR,100).NE.0)) + .OR.(MOD(WATERYR,400).EQ.0)) THEN DAYS(FEBRUARY) = 5 DAYSNYR = 366 ELSE DAYS(FEBRUARY) = 4 ENDIF ENDIF READDAYS = DAYS(MON) ENDIF 412 CALL READCARD (CARDTYPE,NCARD,COMPLETE,EOFFLAG,READDAYS,GAGE, + BOGUSYR,FLOW,YRREAD,MONTHREAD,WEEK,INFILE) IF (EOFFLAG) GOTO 445 *-------- CHECK FOR UNEXPECTED END OF DATA OR UNEXPECTED CHANGE IN THE *-------- GAGE ID NUMBER. THERE SHOULD BE 48 3-CARDS FOR EACH YEAR OF *-------- DATA. IF A CARD IS MISSING, OR THE GAGE ID CHANGES BEFORE *-------- ALL 48 CARDS HAVE BEEN READ, THE PROGRAM MARKS THE YEAR *-------- INCOMPLETE. IF ((CARDTYPE.NE.'3').OR.(GAGE.NE.GAGEID)) THEN COMPLETE=.FALSE. BACKSPACE (INFILE) GOTO 445 ENDIF IF (YEARBASIS.EQ.'Water ') THEN MTHOFWY=MOD(MONTHREAD+2,12)+1 ELSE MTHOFWY=MONTHREAD ENDIF IF ((MTHOFWY.NE.MON).OR.(WEEK.NE.K)) THEN COMPLETE=.FALSE. IF ((MTHOFWY.LE.12).AND.(WEEK.LE.4)) THEN IF (MTHOFWY.GT.MON) THEN BACKSPACE (INFILE) ELSEIF ((MTHOFWY.EQ.MON).AND.(WEEK.GT.K)) THEN BACKSPACE (INFILE) ELSE GOTO 412 ENDIF ENDIF DO 415 L=1,READDAYS DOY=DOY+1 415 CONTINUE GOTO 430 ENDIF *-------- PUT FLOWS INTO ARRAY Q() AND DATES IN DATEx() ARRAYS DO 420 L=1,READDAYS DOY=DOY+1 Q(DOY) = FLOW(L) DATEYR(DOY) = YRREAD DATEMTH(DOY) = MONTHREAD DATEDAY(DOY) = (WEEK-1)*8+L 420 CONTINUE 430 CONTINUE 440 CONTINUE 445 CALL GETMINS (QMINS,DAYMIN,N,DAYSNYR,Q,ZERO,ZSYM) IF (.NOT. EOFFLAG) RETURN *---- END OF INPUT FILE WAS REACHED 450 IF (DOY.EQ.0) THEN *------ ONLY TELL THE MAIN PROGRAM WE'VE REACHED END OF FILE IF WE HIT IT *------ RIGHT AT THE START. OTHERWISE WE WANT IT TO GO AHEAD AND TRY TO READ *------ NEXT YEAR, AND THEN DO THE LAST BFI CALCULATIONS. EOF = .TRUE. ELSE COMPLETE = .FALSE. ENDIF RETURN END SUBROUTINE GETWWWQ (Q,QMINS,DAYMIN,WATERYR,GAGEID,NCARD,ZERO,EOF, + COMPLETE,DAYSNYR,DATEYR,DATEMTH,DATEDAY,YEARBASIS,INFILE,N, + GAGEFOUND, THISDAY,THISMONTH,THISYEAR,YYYYDOT,YYYYDASH, + YYYYMMDD,MDYSLASH,NWISFMT) **************************************************************************** * READ ONE YEAR'S AVERAGE DAILY FLOWS FROM AN INPUT FILE IN WORLD-WIDE-WEB * rdb FORMAT, AND IF SUCCESSFUL, CALL GETMINS TO DETERMINE N-DAY MINIMUMS * * VARIABLES * --------- * ALLFOUND - flag to keep track of whether all data were properly read * (used at end of routine to set COMPLETE = .TRUE. if there * are no missing data) * COMPLETE - .TRUE. if no missing data * DATEYR() - array containing year corresponding to each daily flow * DATEMTH() - array containing month for each daily flow * DATEDAY() - array containing day of the month for each daily flow * DATESTRING - temporarily holds date when reading header * DAYMIN() - array of days of year on which N-day minimums occur * DAYSMTH() - days per month * DAYSNYR - number of days in current year * DOY - day of the year * EOF - .TRUE. if end of file was reached, before any data was read * FLOW() - array holding discharges returned by READCARD * GAGEFOUND - .TRUE. if gage id number was found this time or previously * GAGEID - gage identification number * I,J,K - counters * IBY,IBM,IBD - beginning year, month, day for data file * IDEL - counter to find delimiter when interpreting date format * IEY,IEM,IED - ending year, month, day for data file * INFILE - unit number of input file * INFUTURE - flag to indicate type of date mismatch * MDYSLASH - flag to indicate if date format is mm/dd/yyyy * N - parameter for base-flow separation * NCARD - string image of name card * NWISFMT - integer indicating which NWIS-W format is used * 1 = pre-2001 format (3 columns of data: date, value, flag) * 2 = new NWIS-W (5 columns: agency, station no., date, * value, flag) * Q() - array of mean daily discharges * QDATASTR - holds string represenation of discharge value * QMINS() - array of N-day minimum flows * STARTNEWGAGE - flag indicating we're reading first data for a new gage * THISDAY - expected day value for next line * THISMONTH - expected month value for next line * THISYEAR - expected year value for next line * WATERYR - water year for current year's data * YEARBASIS - string, 'Calendar' or 'Water ' year * YYYYDOT - flag to indicate if date format is yyyy.mm.dd * YYYYDASH - flag to indicate if date format is yyyy-mm-dd * YYYYMMDD - flag to indicate if date format is yyyymmdd * MDYSLASH - flag to indicate if date format is mm/dd/yyyy * ZSYM - symbol for years with no-flow days * ZERO - no-flow days marker variable (returned equal to ZSYM if * a no-flow day is detected)...the symbol gets passed in * through this * * SUBROUTINES USED: * -------------------- * GETMINS, ADVANCEDATE **************************************************************************** CHARACTER*120 THISLINE, PREVLINE CHARACTER*1 ZERO, ZSYM CHARACTER*8 GAGEID,YEARBASIS CHARACTER*78 NCARD CHARACTER*10 DATESTRING CHARACTER*12 QDATASTR LOGICAL COMPLETE, EOF, GAGEFOUND, INFUTURE LOGICAL STARTNEWGAGE, ALLFOUND LOGICAL YYYYDOT, YYYYDASH, YYYYMMDD, MDYSLASH INTEGER DAYMIN(0:367), DAYSNYR, DOY INTEGER DAYSMTH(12) INTEGER I, J, K, IDEL, ICH INTEGER WATERYR INTEGER INFILE, N, DATEYR(366), DATEMTH(366), DATEDAY(366) INTEGER IBY,IBM,IBD,IEY,IEM,IED INTEGER THISDAY, THISMONTH, THISYEAR INTEGER NWISFMT REAL Q(366), QMINS(0:367) DATA DAYSMTH /31,28,31,30,31,30,31,31,30,31,30,31/ *---- NOTE: JANUARY = MONTH 1 500 ZSYM=ZERO DO 505 I=1,366 Q(I) = -99 505 CONTINUE NCARD = ' ' ZERO = ' ' COMPLETE = .FALSE. ALLFOUND = .TRUE. DOY = 0 STARTNEWGAGE = .NOT. GAGEFOUND *---- GET HEADER INFORMATION IF WE HAVEN'T GOTTEN IT PREVIOUSLY IF (.NOT.GAGEFOUND) THEN NWISFMT = 0 510 READ (INFILE,515,END=580) THISLINE 515 FORMAT (A80) IF (THISLINE(1:1).EQ.'#') THEN *-------- NEXT IF...ENDIF BLOCK OF STATEMENTS IS FOR LATEST NWIS-W FORMAT *-------- (RELEASED SOMETIME IN 2001?) IF (PREVLINE(3:29).EQ.'Sites in this file include:') THEN NWISFMT = 2 *---------- NEXT LINE LISTS STATION NUMBER AND DESCRIPTION I = 1 516 I = I + 1 NCARD = THISLINE(I:I+77) IF (NCARD(1:1) .EQ. ' ') GOTO 516 I = 0 517 I = I + 1 GAGEID = NCARD(I:I+7) ICH = ICHAR(GAGEID(1:1)) IF (((ICH.LT.48).OR.(ICH.GT.57)).AND.(I.LT.78)) GOTO 517 IF (I.LT.78) GAGEFOUND = .TRUE. ENDIF IF (THISLINE(3:17).EQ.'Station name :') THEN NWISFMT = 1 *---------- FIND STATION NAME. IT HAS MOVED BACK AND FORTH AS WEB DATA *---------- FORMAT HAS EVOLVED. SOMETIMES AT POSITION 19, SOMETIMES AT 20. I = 18 520 I = I + 1 NCARD = THISLINE(I:I+77) IF (NCARD(1:1) .EQ. ' ') GOTO 520 ENDIF IF (THISLINE(3:17).EQ.'Station number:') THEN NWISFMT = 1 *---------- FIND STATION NUMBER. IT HAS MOVED BACK AND FORTH AS WEB DATA *---------- FORMAT HAS EVOLVED. SOMETIMES AT POSITION 19, SOMETIMES AT 20. I = 18 525 I = I + 1 GAGEID = THISLINE(I:I+7) IF (GAGEID(1:1) .EQ. ' ') GOTO 525 GAGEFOUND = .TRUE. ENDIF IF (PREVLINE(7:24).EQ.'Date Range In File') THEN *---------- THIS IS ONE OF THE EARLY NWIS-W FORMATS (PRE-2001), AND *---------- THE CURRENT LINE DESCRIBES THE DATE RANGE. *---------- USE IT TO DETERMINE HOW THE DATA ARE ORGANIZED *---------- (BY WATER YEAR OR BY CALENDAR YEAR). *---------- MUST ALSO DETERMINE WHICH DATE FORMAT IS USED YYYYDOT = .FALSE. YYYYDASH = .FALSE. YYYYMMDD = .FALSE. MDYSLASH = .FALSE. DO 530 IDEL=1,26 IF (THISLINE(IDEL:IDEL) .EQ. '.') THEN *-------------- DATA FORMAT IS YYYY.MM.DD YYYYDOT = .TRUE. DATESTRING = THISLINE(IDEL-4:IDEL+5) GOTO 535 ELSEIF (THISLINE(IDEL:IDEL) .EQ. '/') THEN *-------------- DATA FORMAT IS MM/DD/YYYY MDYSLASH = .TRUE. DATESTRING = THISLINE(IDEL-2:IDEL+7) GOTO 535 ENDIF 530 CONTINUE *---------- FIGURE OUT WHETHER WE HAVE WATER-YEAR OR CALENDAR-YEAR DATA 535 IF (YYYYDOT) THEN READ (DATESTRING,540) IBY,IBM,IBD,IEY,IEM,IED 540 FORMAT (I4,1x,I2,1x,I2,1x,I4,1x,I2,1x,I2) ELSEIF (MDYSLASH) THEN READ (DATESTRING,541) IBM,IBD,IBY,IEM,IED,IEY 541 FORMAT (I2,1x,I2,1x,I4,1x,I2,1x,I2,1x,I4) ENDIF IF (IBM.EQ.1) THEN YEARBASIS='Calendar' ELSE YEARBASIS='Water ' ENDIF ENDIF PREVLINE = THISLINE ENDIF *------ KEEP LOOPING UNTIL WE'VE READ THE FIRST NON-HEADER LINE IF (THISLINE(1:1).EQ.'#') GOTO 510 *------ IF THIS IS A NEWER NWIS-W DATA FILE, WE WILL HAVE GONE THROUGH *------ ALL THAT AND STILL WON'T KNOW THE DATE FORMAT. WE WILL HAVE *------ TO READ THE FIRST LINE OF DATA TO DETERMINE IT. IF (NWISFMT.EQ.2) THEN BACKSPACE(INFILE) 542 READ (INFILE,515,END=549) THISLINE IF (THISLINE(1:6).EQ.'agency') THEN *---------- WE'RE ALMOST THERE. READ TWO MORE LINES. THEN BACKSPACE *---------- ONE SO THAT WE ARE READY TO READ THE FIRST LINE OF DATA *---------- AGAIN AS NUMERIC DATA. READ (INFILE,515,END=549) THISLINE READ (INFILE,515,END=549) THISLINE BACKSPACE (INFILE) *---------- EXTRACT THE STRING FOLLOWING THE GAGE ID. *---------- IT WILL BE THE DATE. I = 0 543 I = I + 1 IF (THISLINE(I:I+7) .EQ. GAGEID) THEN I = I + 7 544 I = I + 1 ICH = ICHAR(THISLINE(I:I)) *------------ REPEAT UNLESS WE FIND THE START OF THE DATE, SIGNIFIED *------------ BY ANY OF THE FOLLOWING: 0123456789 IF ((ICH.LT.45).OR.(ICH.GT.57)) GOTO 544 J = I 545 J = J + 1 ICH = ICHAR(THISLINE(J:J)) *------------ REPEAT AS LONG AS WE FIND ANY OF THE FOLLOWING: .-/0123456789 IF ((ICH.GE.45).AND.(ICH.LE.57)) GOTO 545 J = J - 1 ELSE GOTO 543 ENDIF DATESTRING = THISLINE(I:J) *---------- DETERMINE DATE FORMAT *---------- POSSIBILITIES ARE: *---------- yyyy.mm.dd *---------- yyyy-mm-dd *---------- yyyymmdd *---------- mm/dd/yyyy YYYYDOT = .FALSE. YYYYDASH = .FALSE. YYYYMMDD = .FALSE. MDYSLASH = .FALSE. IF (DATESTRING(5:5).EQ.'.') THEN YYYYDOT = .TRUE. ELSEIF (DATESTRING(5:5).EQ.'-') THEN YYYYDASH = .TRUE. ELSEIF (DATESTRING(3:3).EQ.'/') THEN MDYSLASH = .TRUE. ELSE YYYYMMDD = .TRUE. ENDIF *---------- DETERMINE IF IT'S WATER-YEAR OR CALENDAR-YEAR DATA 546 FORMAT (I4,1x,I2,1x,I2) 547 FORMAT (I4,I2,I2) 548 FORMAT (I2,1x,I2,1x,I4) IF ((YYYYDOT).OR.(YYYYDASH)) THEN READ (DATESTRING,546) IBY,IBM,IBD ELSEIF (YYYYMMDD) THEN READ (DATESTRING,547) IBY,IBM,IBD ELSEIF (MDYSLASH) THEN READ (DATESTRING,548) IBM,IBD,IBY ENDIF IF (IBM.EQ.1) THEN YEARBASIS='Calendar' ELSE YEARBASIS='Water ' ENDIF ELSE GOTO 542 ENDIF ENDIF 549 IF (.NOT. GAGEFOUND) GOTO 580 *---- THERE IS A BLANK LINE AFTER THE 'Date Discharge' HEADER IN *---- THE SAMPLE FILE I HAVE, BUT IT MAY NOT ALWAYS BE THERE. THE *---- DATE-CHECKING CODE THAT FOLLOWS WILL CAUSE US TO ADVANCE PAST *---- IT, IF IT IS ENCOUNTERED, SO I WON'T MAKE ANY SPECIAL EFFORT *---- TO READ PAST IT HERE, SINCE THAT WOULD CAUSE PROBLEMS IF IT'S *---- NOT THERE. ENDIF *---- READ DAILY FLOWS 550 DOY = DOY + 1 555 READ (INFILE,515,END=580) THISLINE *---- A '#' IS THE START OF AN RDB HEADER LINE, AND A '<' *---- IS AN HTML TAG THAT MIGHT BE AT THE TOP OF AN RDB FILE IF (THISLINE(1:1).EQ.'#' .OR. THISLINE(1:1).EQ.'<') THEN *------ WE'VE FOUND THE START OF ANOTHER GAGE *------ (OR AT LEAST THE END OF THIS ONE) BACKSPACE(INFILE) GAGEFOUND = .FALSE. IF (DOY.EQ.1) THEN *---------- THIS WAS THE FIRST LINE OF DATA. RESTART THIS ROUTINE FROM *---------- THE TOP, AS IF WE ARE READING IN A BRAND NEW GAGE HEADER. GOTO 500 ELSE *---------- GOTO THE END OF THIS ROUTINE, WHERE WE LEAVE *---------- THIS YEAR MARKED INCOMPLETE GOTO 580 ENDIF ENDIF *---- IF THIS IS A NEW-FORMAT NWIS-W FILE (NWISFMT=2) THEN READ *---- PAST THE FIRST TWO ITEMS OF DATA ON THE LINE *---- (AGENCY CODE AND GAGE ID). BUT, VERIFY THAT GAGE ID *---- HAS NOT CHANGED. IF IT HAS, SKIP THIS LINE BY GOING BACK TO 555. *---- HOPEFULLY, NEXT LINE WILL BE FOR THE CORRECT GAGE AGAIN. IF NOT, *---- WE'LL CONTINUE READING UNTIL WE FIND THE GAGE AGAIN, OR UNTIL WE *---- FIND THE BEGINNING OF ANOTHER HEADER FOR A NEW GAGE. IF (NWISFMT.EQ.2) THEN *------ READ THROUGH THE AGENCY CODE TO THE FIRST SPACE OR TAB CHARACTER I = 0 556 I = I + 1 ICH = ICHAR(THISLINE(I:I)) IF ((ICH.NE.9).AND.(ICH.NE.32)) GOTO 556 *------ READ THROUGH THE SPACES OR TAB(S) TO THE START OF THE GAGE ID 557 I = I + 1 ICH = ICHAR(THISLINE(I:I)) IF ((ICH.EQ.9).OR.(ICH.EQ.32)) GOTO 557 *------ NEXT EIGHT CHARACTERS ARE THE GAGE ID. IF IT DOESN'T MATCH *------ VALUE WE EXPECT, SKIP THIS LINE BY GOING BACK TO 555 TO READ *------ THE NEXT LINE OF DATA. IF (THISLINE(I:I+7).NE.GAGEID) GOTO 555 I = I + 7 *------ READ THROUGH THE SPACES OR TAB(S) TO THE START OF THE DATE 558 I = I + 1 ICH = ICHAR(THISLINE(I:I)) IF ((ICH.EQ.9).OR.(ICH.EQ.32)) GOTO 558 *------ MAKE THISLINE EQUAL TO THE REMAINDER OF THE STRING. FROM *------ HERE ON IT CAN THEN BE PROCESSED IDENTICALLY TO THE OLD *------ NWIS-W FORMAT. THISLINE = THISLINE(I:120) ENDIF *---- MAKE SURE IT LOOKS LIKE A VALID LINE, WITH PROPER DATE DELIMITERS IF (YYYYDOT) THEN IF (.NOT. (THISLINE(5:5).EQ.'.'.AND.THISLINE(8:8).EQ.'.') ) + GOTO 555 ELSEIF (YYYYDASH) THEN IF (.NOT. (THISLINE(5:5).EQ.'-'.AND.THISLINE(8:8).EQ.'-') ) + GOTO 555 ELSEIF (MDYSLASH) THEN IF (.NOT. (THISLINE(3:3).EQ.'/'.AND.THISLINE(6:6).EQ.'/') ) + GOTO 555 ENDIF *---- MUST STRIP OFF ANY TRAILING TAB OR TRAILING TAB AND FLAG CHARACTER *---- (SUCH AS AN 'e' FOR AN ESTIMATED VALUE). I TRIED USING UNFORMATTED *---- READS, BUT THE DATES CAN'T BE READ UNFORMATTED BECAUSE THE DECIMALS *---- BETWEEN THE DATE FIELDS MAKE FORTRAN THINK THE DATE IS NOT AN INTEGER. *---- THE FORMATTED READ WORKS FINE FOR THE DATES, BUT WILL BOMB ON THE *---- DISCHARGE DATA IF THE TRAILING TAB AND FLAG ARE NOT REMOVED. *---- THIS IS A PROBLEM BECAUSE THE DISCHARGE FIELD IS NOT _REALLY_ *---- 8 CHARACTERS WIDE. TO BE PRECISE, THE LINE OF DATA CONTAINS A *---- 10-CHARACTER DELIMITED DATE, FOLLOWED BY A SINGLE TAB CHARACTER, *---- THEN THE DISCHARGE VALUE AS A FLOATING POINT NUMBER. THE DISCHARGE *---- VALUE MAY USE AS MANY AS 8 CHARACTERS, BUT PROBABLY IS LESS. IF THERE *---- IS NO FLAG ON THE MEASUREMENT, THIS IS IT, BUT IF THERE IS A FLAG, THEN *---- THE DISCHARGE VALUE IS FOLLOWED BY ANOTHER TAB AND A FLAG CHARACTER. *---- IF YOU JUST READ AN F8.0 FOR THE DISCHARGE IT WORKS FINE WHEN THE LINE *---- ENDS IMMEDIATELY AFTER THE DISCHARGE VALUE, BUT IF THERE IS A TRAILING *---- TAB AND CHARACTER, AND THE DISCHARGE VALUE WAS NOT A FULL 8 CHARACTERS *---- LONG THEN THE TAB AND CHARACTER ARE PART OF THE 8 CHARACTERS ANALYZED *---- BY THE F8.0 FORMATTED READ, AND IT CAN'T HANDLE THEM. *---- IF WE CAN FIND TWO TABS (ASCII 9'S), TRUNCATE PRIOR TO THE SECOND TAB I = 0 DO 560 J=1,LEN(THISLINE) IF (ICHAR(THISLINE(J:J)).EQ.9) THEN I=I+1 IF (I.EQ.1) THEN DATESTRING = THISLINE(1:J-1) K = J+1 ENDIF ENDIF IF (I.EQ.2) THEN QDATASTR = THISLINE(K:J-1) GOTO 561 ENDIF 560 CONTINUE 561 IF (I.NE.2) THEN *------ IT MUST NOT HAVE BEEN DELIMITED (MAYBE WAS DOWNLOADED AND *------ SAVED WITH TABS CONVERTED TO SPACES) *------ USE THE FIRST 10 CHARACTERS FOR THE DATE DATESTRING = THISLINE(1:10) *------ USE THE REMAINDER FOR THE DATA, BUT STRIP OFF *------ ANY TRAILING FLAG CHARACTER BECAUSE IT WILL MESS UP *------ THE FORMATTED READ. MUST DO THIS BECAUSE WE DON'T *------ KNOW FOR SURE HOW THE TABS GOT CONVERTED TO SPACES. *------ MOST EDITORS WOULD REPLACE EACH TAB WITH 8 SPACES, *------ BUT THERE ARE ALSO OTHER POSSIBILITIES. I = 0 J = 0 DO 562 K=11,LEN(THISLINE) ICH = ICHAR(THISLINE(K:K)) IF ((I.EQ.0).OR.(J.EQ.0)) THEN *---------- LOOK FOR A NON-SPACE IF (ICH.NE.32) THEN *------------ LOOK FOR A NUMERIC CHARACTER IF ((ICH.EQ.46).OR.((ICH.GE.48).AND.(ICH.LE.57))) THEN *---------------- IF WE HAVEN'T FOUND THE START, THEN THIS IS IT IF (I.EQ.0) I = K ELSE *---------------- THIS MIGHT BE END OF DATA SINCE IT'S NON-NUMERIC IF (J.EQ.0) J = K-1 ENDIF ENDIF ENDIF 562 CONTINUE *------ IF WE DIDN'T FIND ANYTHING, THEN JUST TRY TO USE EVERYTHING *------ THAT FOLLOWED THE DATE. IF (I.EQ.0) I = 11 IF (J.EQ.0) J = LEN(THISLINE) QDATASTR = THISLINE(I:J) ENDIF *---- FINALLY, READ THE DATA USING FORMATTED READS IF ((YYYYDOT).OR.(YYYYDASH)) THEN READ (DATESTRING,563) DATEYR(DOY), DATEMTH(DOY), DATEDAY(DOY) 563 FORMAT (I4,1x,I2,1x,I2) ELSEIF (YYYYMMDD) THEN READ (DATESTRING,564) DATEYR(DOY), DATEMTH(DOY), DATEDAY(DOY) 564 FORMAT (I4,I2,I2) ELSEIF (MDYSLASH) THEN READ (DATESTRING,565) DATEMTH(DOY), DATEDAY(DOY), DATEYR(DOY) 565 FORMAT (I2,1x,I2,1x,I4) ENDIF READ (QDATASTR,566) Q(DOY) 566 FORMAT (F24.0) IF ( (Q(DOY).LT.0) .OR. (Q(DOY).EQ.999999) ) THEN Q(DOY) = -99. ALLFOUND = .FALSE. ENDIF IF (DOY.EQ.1) THEN IF (YEARBASIS.EQ.'Water ') THEN *---------- EXPECT TO GET DATA STARTING ON OCTOBER 1 OF THE *---------- CURRENT WATER YEAR IF (DATEMTH(DOY).GT.9) THISYEAR = DATEYR(DOY) IF (DATEMTH(DOY).LT.10) THISYEAR = DATEYR(DOY) - 1 THISMONTH = 10 THISDAY = 1 ELSE *---------- EXPECT TO GET DATA STARTING ON JANUARY 1 OF THE *---------- CURRENT CALENDAR YEAR THISYEAR = DATEYR(DOY) THISMONTH = 1 THISDAY = 1 ENDIF IF (YEARBASIS.EQ.'Water ') THEN WATERYR = THISYEAR + 1 ELSE WATERYR = THISYEAR ENDIF IF (((MOD(WATERYR,4).EQ.0).AND.(MOD(WATERYR,100).NE.0)) + .OR.(MOD(WATERYR,400).EQ.0)) THEN DAYSNYR = 366 DAYSMTH(2) = 29 ELSE DAYSNYR = 365 DAYSMTH(2) = 28 ENDIF ENDIF *---- MAKE SURE THE DATA WE READ IS FOR THE CORRECT DAY. *---- IF NOT, MARK ANY MISSED DAYS WITH THE MISSING DATA FLAG *---- (Q = -99), OR TRY READING THE NEXT LINE UNTIL WE FIND THE *---- CORRECT DATE. IF CORRECT DATE IS NEVER FOUND, THIS *---- WILL EXHAUST THE INPUT FILE, OR REACH A HEADER THAT *---- MARKS BEGINNING OF ANOTHER SECTION OF DATA (POSSIBLY *---- FOR A DIFFERENT GAGE). 570 IF (DATEYR(DOY) .NE. THISYEAR .OR. + DATEMTH(DOY) .NE. THISMONTH .OR. + DATEDAY(DOY) .NE. THISDAY ) THEN *------ IF DATE READ FROM FILE IS IN THE FUTURE (I.E., SOME LINES *------ WERE MISSING), THEN ADVANCE THE DATE COUNTERS UNTIL THEY *------ MATCH THE DATE IN THE FILE, AND MARK INTERVENING DAYS AS *------ MISSING. INFUTURE = .FALSE. IF (DATEYR(DOY).GT.THISYEAR) INFUTURE = .TRUE. IF (DATEYR(DOY).EQ.THISYEAR .AND. + DATEMTH(DOY).GT.THISMONTH ) INFUTURE = .TRUE. IF (DATEYR(DOY).EQ.THISYEAR .AND. + DATEMTH(DOY).EQ.THISMONTH .AND. + DATEDAY(DOY).GT.THISDAY ) INFUTURE = .TRUE. IF (INFUTURE) THEN *---------- MOVE DATA VALUES FORWARD A DAY AND MARK PREVIOUS FLOW AS MISSING DOY = DOY + 1 Q(DOY) = Q(DOY-1) Q(DOY-1) = -99 ALLFOUND = .FALSE. DATEYR(DOY) = DATEYR(DOY-1) DATEMTH(DOY) = DATEMTH(DOY-1) DATEDAY(DOY) = DATEDAY(DOY-1) DATEYR(DOY-1) = THISYEAR DATEMTH(DOY-1) = THISMONTH DATEDAY(DOY-1) = THISDAY *---------- ADVANCE DATE COUNTERS CALL ADVANCEDATE (THISDAY, THISMONTH, THISYEAR, DAYSMTH) IF (DOY .LE. DAYSNYR) THEN GOTO 570 ELSE BACKSPACE(INFILE) ENDIF ELSE *------ DATE READ FROM FILE IS IN THE PAST (MAYBE THERE WERE *------ DUPLICATE LINES), LOOP BACK AND READ ANOTHER LINE IN *------ HOPES THAT THE FILE WILL CATCH UP. GOTO 555 ENDIF ENDIF *---- ADVANCE DATE COUNTERS TO NEXT EXPECTED DATE CALL ADVANCEDATE (THISDAY, THISMONTH, THISYEAR, DAYSMTH) *------ LOOP UNTIL WE REACH END OF YEAR IF (DOY .LT. DAYSNYR) THEN GOTO 550 ELSEIF (ALLFOUND) THEN COMPLETE = .TRUE. ENDIF GOTO 590 *---- FOLLOWING LINE IS THE ERROR TRAP IN CASE END OF FILE IS REACHED *---- BACK OFF THE END OF FILE, SO WE CAN HIT IT AGAIN WHEN WE TRY *---- TO READ NEXT YEAR'S DATA. 580 BACKSPACE (INFILE) 590 CALL GETMINS (QMINS,DAYMIN,N,DAYSNYR,Q,ZERO,ZSYM) IF (COMPLETE) RETURN *---- END OF INPUT FILE WAS REACHED, OR THE DATE WERE NOT COMPLETE *---- BUT, ONLY TELL THE MAIN PROGRAM WE'VE REACHED END OF FILE IF WE HIT IT *---- RIGHT AT THE START. OTHERWISE WE WANT IT TO GO AHEAD AND TRY TO READ *---- NEXT YEAR, AND THEN DO THE LAST BFI CALCULATIONS. IF (DOY.LE.1) THEN EOF = .TRUE. ENDIF RETURN END SUBROUTINE GETMINS (QMINS,DAYMIN,N,DAYSNYR,Q,ZERO,ZSYM) **************************************************************************** * SUBROUTINE TO DETERMINE MINIMUMS FOR EACH N-DAY SEQUENCE AND PUT RESULTS * IN ARRAYS QMINS() AND DAYMIN() * IF THERE IS A FRACTION OF AN N-DAY PERIOD LEFT AT THE END OF A YEAR * IT IS LUMPED IN WITH THE LAST N-DAY PERIOD. * * VARIABLES * --------- * DAYMIN() - array of days of year on which N-day minimums occur * DAYN - temporary...day of water year on which N-day minimum occurs * DAYSNYR - number of days in current year * I,J,JJ - counters * N - parameter for base-flow separation, span of interval in days * Q() - array of mean daily discharges * QMIN - temporary...minimum Q in N-day sequence * QMINS() - array of N-day minimum flows * ZERO - no-flow days marker (returned equal to ZSYM if * a no-flow day is detected) * ZSYM - symbol to use for years with no-flow days **************************************************************************** INTEGER I, J, JJ INTEGER N, DAYSNYR, DAYN INTEGER DAYMIN(0:367) CHARACTER*1 ZERO, ZSYM REAL Q(366), QMIN, QMINS(0:367) DO 590 J=1,DAYSNYR/N QMIN = 999999 IF (J.EQ.(DAYSNYR/N)) THEN JJ=DAYSNYR-(J-1)*N ELSE JJ=N ENDIF IF (J.LT.(DAYSNYR/N)) THEN DO 585 I=1,JJ IF (Q(N*(J-1)+I).LT.QMIN) THEN QMIN = Q(N*(J-1)+I) DAYN = N*(J-1)+I ENDIF 585 CONTINUE ELSE DO 588 I=1,JJ IF (Q(N*(J-1)+I).LE.QMIN) THEN QMIN = Q(N*(J-1)+I) DAYN = N*(J-1)+I ENDIF 588 CONTINUE ENDIF IF (QMIN.EQ.0.) ZERO = ZSYM QMINS(J) = QMIN DAYMIN(J) = DAYN 590 CONTINUE RETURN END SUBROUTINE ADVANCEDATE (THISDAY, THISMONTH, THISYEAR, DAYSMTH) ********************************************************************* * ADVANCES DATE COUNTERS FOR USE IN KEEPING TRACK OF CALENDAR DATES * AS DATA ARE READ IN FROM THE INPUT FILE. * * VARIABLES: * ---------- * THISDAY - day to advance * THISMONTH - month * THISYEAR - year * DAYSMTH - array containing number of days per month for this year ********************************************************************* INTEGER THISDAY, THISMONTH, THISYEAR, DAYSMTH(12) THISDAY = THISDAY + 1 IF (THISDAY .GT. DAYSMTH(THISMONTH) ) THEN THISDAY = 1 THISMONTH = THISMONTH + 1 IF (THISMONTH .GT. 12) THEN THISMONTH = 1 THISYEAR = THISYEAR + 1 ENDIF ENDIF RETURN END SUBROUTINE READCARD (CARDTYPE,NCARD,COMPLETE,EOF,DAYS,GAGE, + WATERYR,FLOW,IYEAR,MONTH,IWEEK,INFILE) **************************************************************************** * SUBROUTINE TO READ ONE CARD FROM THE INPUT FILE AND CHECK IF DATA IS * COMPLETE * * VARIABLES * CARDTYPE - character in leftmost column of data card * COMPLETE - .TRUE. if no data missing * DAYS - number of days data to be expected on card * DUM() - string versions of flow data - used to detect missing data * EOF - .TRUE. if end of input file is encountered * FLOW() - average daily flows read from input file * GAGE - gage id string * I - counter * INFILE - unit number for input file * IWEEK - line number read from input file * IYEAR - year read from input file * JUNK - string used in detecting missing data * L - counter * MARK() - .TRUE. if data is missing * MONTH - month read from input file * NCARD - string image of the name card * THISLINE - string to hold line of data until we figure out what it is * WATERYR - water year corresponding to IYEAR **************************************************************************** CHARACTER*1 CARDTYPE CHARACTER*7 DUM(8) CHARACTER*8 GAGE CHARACTER*14 JUNK CHARACTER*78 NCARD CHARACTER*80 THISLINE LOGICAL COMPLETE, EOF, MARK(8) INTEGER DAYS, I, IYEAR, IWEEK, INFILE, MONTH INTEGER WATERYR, L REAL FLOW(DAYS) *--------------------------------------------------------------------------- EOF=.FALSE. DO 600 I=1,DAYS MARK(I)=.TRUE. 600 CONTINUE *---- READ THE CARD AS A STRING READ (INFILE,605,END=660) THISLINE 605 FORMAT (A80) *---- PARSE INTO SEPARATE STRINGS READ (THISLINE,610) CARDTYPE,GAGE,JUNK,(DUM(I),I=1,8) 610 FORMAT (A1,1x,A8,A14,8A7) IF (CARDTYPE.EQ.'N') THEN READ (THISLINE,620) NCARD 620 FORMAT (A78) ENDIF *---- IF IT'S A 3 CARD, THEN CHECK TO SEE IF ANY DATA IS MISSING IF (CARDTYPE.EQ.'3') THEN DO 630 L=1,DAYS IF((LEN(DUM(L)).NE.7).OR.(DUM(L).EQ.' ').OR. + (DUM(L).EQ.' 999999')) THEN COMPLETE=.FALSE. MARK(L)=.FALSE. ENDIF 630 CONTINUE *------ PARSE DATA CARD FOR REAL READ (THISLINE,640) CARDTYPE,GAGE,IYEAR,MONTH,IWEEK, + (FLOW(I),I=1,DAYS) 640 FORMAT (A1,1x,A8,6x,I4,I2,I2,8F7.0) IF (MONTH.LT.10) THEN WATERYR = IYEAR ELSE WATERYR = IYEAR + 1 ENDIF *------ CHECK FOR NEGATIVE FLOWS (MISSING DATA) *------ MARK ALL MISSING DATA POINTS WITH FLOW(I) = -99 DO 650 I=1,DAYS IF ((FLOW(I).LT.0.).OR.(.NOT.(MARK(I)))) THEN FLOW(I) = -99 COMPLETE=.FALSE. ENDIF 650 CONTINUE ENDIF RETURN *---- READ STATEMENT DETECTED END OF FILE 660 EOF=.TRUE. BACKSPACE (INFILE) RETURN END SUBROUTINE BASEFLOW (DAY1,DAY2,Q1,Q2,BF) **************************************************************************** * SUBROUTINE TO INTEGRATE BASE FLOW FROM DAY1 TO DAY2. * * IF POSSIBLE, USE LOGARITHMIC BASE-FLOW RECESSION TO COMPUTE * VOLUME OF BASE FLOW. IF NOT POSSIBLE, (Q1=Q2 OR Q1=0 OR Q2=0) * ASSUME LINEAR VARIATION. UNITS OF BASE-FLOW VOLUME ARE CFS-DAYS. * NOTE THAT THE FORTRAN LOG FUNCTION IS THE NATURAL LOGARITHM. * * VARIABLES * BF - base-flow volume from DAY1 to DAY2 (cfs-days) * DAY1 - day of water year for first turning point * DAY2 - day of water year for second turning point * Q1 - average daily base flow at first turning point * Q2 - average daily base flow at second turning point * SLOPE - slope of base-flow recession line (log cycles per day) **************************************************************************** REAL DAY1,DAY2,Q1,Q2,BF,SLOPE *--------------------------------------------------------------------------- IF ((Q1.EQ.Q2).OR.(Q1.EQ.0.).OR.(Q2.EQ.0.)) THEN BF=(((Q2+Q1)/2.)*(DAY2-DAY1)) ELSE SLOPE = LOG(Q2/Q1)/(DAY2-DAY1) BF=(Q1/SLOPE)*(EXP(SLOPE*(DAY2-DAY1))-1.) ENDIF RETURN END SUBROUTINE TEST_FOR_TP (FRAC,Q0,Q1,Q2,D0,D1,D2,K1DAY,METHOD, + TURNPT) **************************************************************************** * TEST TO SEE IF Q1 IS A TURNING POINT. * * VARIABLES * F - temporary variable * FRAC - ratio used in base-flow separation (default = 0.9) * D0 - day corresponding to Q0 * D1 - day corresponding to Q1 * D2 - day corresponding to Q2 * METHOD - integer defining which turning point test to use * Q0 - preceeding N-day minimum * Q1 - possible turning point * Q2 - following N-day minimum * K1DAY - 1-day recession constant * TURNPT - .TRUE. if Q1 is a turning point **************************************************************************** LOGICAL TURNPT INTEGER D0,D1,D2,METHOD REAL F,FRAC,Q0,Q1,Q2,K1DAY *--------------------------------------------------------------------------- TURNPT=.FALSE. IF (Q1) 700,710,720 *---- Q1 LESS THAN ZERO 700 TURNPT=.TRUE. RETURN *---- Q1 EQUAL ZERO 710 TURNPT=.TRUE. RETURN *---- Q1 GREATER THAN ZERO 720 IF (METHOD.EQ.1) THEN F=FRAC*Q1 IF ((Q0*Q2).NE.0) THEN IF ((F.LE.Q0) .AND. (F.LE.Q2)) TURNPT=.TRUE. ELSEIF (Q0.EQ.0) THEN IF (F.LE.Q2) TURNPT=.TRUE. ELSE IF (F.LE.Q0) TURNPT=.TRUE. ENDIF ELSE IF ((Q0*Q2).NE.0) THEN IF ((Q1*(K1DAY**(D1-D0)).LE.Q0) .AND. + (Q1*(K1DAY**(D2-D1)).LE.Q2)) TURNPT=.TRUE. ELSEIF (Q0.EQ.0) THEN IF (Q1*(K1DAY**(D2-D1)).LE.Q2) TURNPT=.TRUE. ELSE IF (Q1*(K1DAY**(D1-D0)).LE.Q0) TURNPT=.TRUE. ENDIF ENDIF RETURN END SUBROUTINE GETOPTIONS (INPUTFILE,OUTPUTFILE,TPOUTFILE,DAILYFILE, + ZSYM,ISYM,XSYM,METHOD,N,FRAC,K1DAY,TPFILE,DAYFILE) **************************************************************************** * THIS SUBROUTINE PRESENTS DEFAULT BFI PARAMETERS AND ALLOWS THE USER TO * MODIFY THOSE PARAMETERS USING AN ON-SCREEN MENU. * * VARIABLES: * INPUTFILE - name of the input file (255 chars., max) * OUTPUTFILE - name of the output file * TPOUTFILE - name of the turning point output file * DAILYFILE - name of the daily flow output file * ZSYM - the zero-flow symbol (empty if symbols aren't used) * ISYM - the interpolation symbol (empty if symbols aren't used) * XSYM - the extrapolation symbol (empty if symbols aren't used) * METHOD - integer to identify turning point test method * 1 = standard Institute of Hydrology method, using N and FRAC * 2 = modified method using N and K1DAY * N - partition length in days * FRAC - fraction used in turning point test for standard method * K1DAY - one-day recession constant used in modified method * TPFILE - unit number of turning point output file * DAYFILE - unit number of daily flow output file * CHOICE - user response to menu * BLANK - empty string for testing input * I - counter for decomposing input file name *DEFOUTPUTFILE - default name for output file * DEFTPOUTFILE - default name for turning point output file * DEFDAILYFILE - default name for daily flow (hydrograph) file * NEWFILE - new output filename supplied by user * * SUBROUTINES AND FUNCTIONS USED: * ------------------------------- * UCASE, EXISTS **************************************************************************** CHARACTER*255 INPUTFILE,OUTPUTFILE,TPOUTFILE,DAILYFILE CHARACTER*1 ZSYM,ISYM,XSYM INTEGER METHOD,N REAL FRAC,K1DAY INTEGER TPFILE,DAYFILE CHARACTER*14 NOFILE CHARACTER*1 CHOICE CHARACTER*255 BLANK INTEGER I CHARACTER*255 DEFOUTPUTFILE,DEFTPOUTFILE,DEFDAILYFILE,NEWFILE *---- FUNCTIONS CALLED FROM THIS ROUTINE CHARACTER*1 UCASE LOGICAL EXISTS NOFILE = '<--- none --->' BLANK = ' ' BLANK = BLANK//BLANK//BLANK//BLANK//BLANK//' ' *---- SET DEFAULT PARAMETER VALUES *------ BUILD DEFAULT OUTPUTFILE NAME BY REMOVING ANY EXTENSION FROM THE *------ INPUTFILE NAME AND ADDING '.bfi' IN IT'S PLACE. *-------- IF INPUTFILE DOES NOT CONTAIN '.', THEN JUST ADD '.bfi' IF (N.EQ.0) THEN IF (INDEX(INPUTFILE,'.').EQ.0) THEN OUTPUTFILE = INPUTFILE//'.bfi' DEFTPOUTFILE = INPUTFILE//'.tp' DEFDAILYFILE = INPUTFILE//'.q' ELSE *------------- FIND THE LAST OCCURRENCE OF '.' IN INPUTFILE DO 800, I=LEN(INPUTFILE),1,-1 IF (INPUTFILE(I:I).EQ.'.') THEN OUTPUTFILE = INPUTFILE(1:I-1)//'.bfi' DEFTPOUTFILE = INPUTFILE(1:I-1)//'.tp' DEFDAILYFILE = INPUTFILE(1:I-1)//'.q' GOTO 801 ENDIF 800 CONTINUE 801 ENDIF ENDIF DEFOUTPUTFILE = OUTPUTFILE TPOUTFILE = NOFILE DAILYFILE = NOFILE *----------SYMBOLS IN OUTPUT ARE TURNED OFF BY DEFAULT ZSYM = ' ' ISYM = ' ' XSYM = ' ' IF (METHOD.EQ.0) METHOD = 1 IF (N.EQ.0) N = 5 IF (FRAC.EQ.0) FRAC = 0.9 IF (K1DAY.EQ.0) K1DAY = 0.97915 *---- DISPLAY MENU WITH CURRENT OPTIONS AND PARAMETERS 805 WRITE (*,810) 810 FORMAT (/////////////////////////1X, +'_________________________', +'______________________________', +'________________________'//1X, +32X,'BFI Program Options',29X/1X, +'_________________________', +'______________________________', +'________________________'/) WRITE (*,820) INPUTFILE WRITE (*,830) OUTPUTFILE WRITE (*,840) TPOUTFILE WRITE (*,850) DAILYFILE 820 FORMAT (3X,'I) Input file:',T35,A40) 830 FORMAT (3X,'O) Output file:',T35,A40) 840 FORMAT (3X,'T) Turning point file:',T35,A40) 850 FORMAT (3X,'H) Hydrograph file:',T35,A40) IF (ZSYM.EQ.'*') THEN WRITE (*,*) ' S) Symbols: ', + 'WILL BE used in output' ELSE WRITE (*,*) ' S) Symbols: ', + 'will NOT be used in output' ENDIF IF (METHOD.EQ.1) THEN WRITE (*,855) 855 FORMAT (1X, ' M) Method: Institute ', + 'of Hydrology (Standard Method)') ELSE WRITE (*,*) + ' M) Method: Modified Method' ENDIF WRITE (*,860) N 860 FORMAT (1X, ' N) N: ',T35,I3) IF (METHOD.EQ.1) THEN WRITE (*,870) FRAC 870 FORMAT (3X, 'F) f: ',T35,F8.6//) ELSE WRITE (*,880) K1DAY 880 FORMAT (3X, 'K) K: ',T35,F8.6//) ENDIF *-------- GET USER'S CHOICE WRITE (*,890) 890 FORMAT (7X,'Choose option to change, or "Q"', + ' to quit without running BFI.'/7X,'Press to', + ' use values shown above.'/) *-------- IF THE USER PRESSES ENTER, CHOICE WILL BE EQUAL TO ' ' *-------- BECAUSE THE READ USES AN A1 FORMAT. READ (*,895) CHOICE 895 FORMAT (A1) CHOICE = UCASE(CHOICE) IF (.NOT.((CHOICE.EQ.'I').OR. + (CHOICE.EQ.'O').OR. + (CHOICE.EQ.'T').OR. + (CHOICE.EQ.'H').OR. + (CHOICE.EQ.'S').OR. + (CHOICE.EQ.'M').OR. + (CHOICE.EQ.'N').OR. + (CHOICE.EQ.' ').OR. + (CHOICE.EQ.'Q') )) THEN IF (METHOD.EQ.1) THEN IF (.NOT.(CHOICE.EQ.'F')) GOTO 805 ELSE IF (.NOT.(CHOICE.EQ.'K')) GOTO 805 ENDIF ENDIF *-------- PROCESS THE USER'S RESPONSE WRITE (*,*) WRITE (*,*) IF (CHOICE.EQ.' ') THEN GOTO 999 ELSEIF (CHOICE.EQ.'I') THEN 901 WRITE (*,*) 'Enter name of input file' READ (*,905) INPUTFILE 905 FORMAT (A255) IF (.NOT.EXISTS(INPUTFILE)) THEN WRITE (*,*) '*** ERROR: ',INPUTFILE, ' does not exist.' WRITE (*,*) GOTO 901 ENDIF GOTO 805 ELSEIF (CHOICE.EQ.'O') THEN WRITE (*,*) 'Enter name of output file for program results,' WRITE (*,*) 'or press to accept default filename:' WRITE(*,907) DEFOUTPUTFILE READ (*,905) NEWFILE IF (NEWFILE.EQ.BLANK) THEN OUTPUTFILE=DEFOUTPUTFILE ELSE OUTPUTFILE=NEWFILE ENDIF GOTO 805 ELSEIF (CHOICE.EQ.'T') THEN WRITE(*,*)'Enter name of output file for turning points,' IF (TPOUTFILE.EQ.NOFILE) THEN WRITE(*,*)'or press to accept default filename:' WRITE(*,907) DEFTPOUTFILE ELSE WRITE(*,*) + 'or press to disable turning point output.' ENDIF READ (*,905) NEWFILE IF (NEWFILE.EQ.BLANK) THEN IF (TPOUTFILE.EQ.NOFILE) THEN TPOUTFILE=DEFTPOUTFILE ELSE TPOUTFILE=NOFILE ENDIF ELSE TPOUTFILE=NEWFILE ENDIF GOTO 805 ELSEIF (CHOICE.EQ.'H') THEN WRITE (*,*) 'Enter output file for daily hydrograph data,' IF (DAILYFILE.EQ.NOFILE) THEN WRITE(*,*) 'or press to accept default filename:' WRITE(*,907) DEFDAILYFILE ELSE WRITE(*,*) + 'or press to disable hydrograph output.' ENDIF 907 FORMAT (1X,'Default = ',A68) READ (*,905) NEWFILE IF (NEWFILE.EQ.BLANK) THEN IF (DAILYFILE.EQ.NOFILE) THEN DAILYFILE=DEFDAILYFILE ELSE DAILYFILE=NOFILE ENDIF ELSE DAILYFILE=NEWFILE ENDIF IF (DAILYFILE.EQ.BLANK) DAILYFILE=NOFILE GOTO 805 ELSEIF (CHOICE.EQ.'S') THEN IF (ZSYM.EQ.'*') THEN ZSYM = ' ' ISYM = ' ' XSYM = ' ' ELSE ZSYM = '*' ISYM = '^' XSYM = '!' ENDIF GOTO 805 ELSEIF (CHOICE.EQ.'M') THEN 909 WRITE (*,910) 910 FORMAT(/////////////////////////6X, +'Two methods of base-flow separation are available.'/6X, +'The standard method is the Institute of Hydrology Method.'/6X, +'==========================================================='//6X +,'The Institute of Hydrology turning point test is:'/6X, +'-----------------------------------------------------------'/6X, +'Given three adjacent N-day minimum flows, (Q0,Q1,Q2)'//11X, +'Q1 is a turning point IF:'/16X, +'Q1*f <= Q0 AND Q1*f <= Q2'//6X, +'The Modified Method uses the following turning point test:'/6X, +'-----------------------------------------------------------'/6X, +'Given three adjacent N-day minimum flows, (Q0,Q1,Q2)'/6X, +'with N01 and N12 being the number of days separating the'/6X, +'occurrence of Q0 and Q1, and Q1 and Q2, respectively:'//11X, +'Q1 is a turning point IF:'/16X, +'Q1*(K^N01) <= Q0 AND Q1*(K^N12) <= Q2'//6X, +'Choose method (1=Standard, 2=Modified):') READ (*,*) METHOD IF ((METHOD.NE.1).AND.(METHOD.NE.2)) GOTO 909 GOTO 805 ELSEIF (CHOICE.EQ.'N') THEN 920 WRITE (*,*) 'Enter a value for "N" (1-365 days):' READ (*,*) N IF ((N.LT.1).OR.(N.GT.365)) GOTO 920 GOTO 805 ELSEIF (CHOICE.EQ.'F') THEN WRITE (*,*) 'Enter a value for "f" (typically < 1.0):' READ (*,*) FRAC GOTO 805 ELSEIF (CHOICE.EQ.'K') THEN WRITE (*,*) 'Enter a value for "K", the 1-day recession', + ' constant:' READ (*,*) K1DAY GOTO 805 ELSEIF (CHOICE.EQ.'Q') THEN STOP ENDIF *----------CLEANUP AND FINISH SUBROUTINE 999 CONTINUE *----------IF ANY OF THE FILENAMES ARE '', set the associated unit *----------NUMBER TO ZERO IF (TPOUTFILE.EQ.NOFILE) TPFILE=0 IF (DAILYFILE.EQ.NOFILE) DAYFILE=0 RETURN END REAL FUNCTION QINTERP (D1,Q1,D3,Q3,D2) **************************************************************************** * PERFORM STRAIGHT-LINE INTERPOLATION BETWEEN TWO BASE-FLOW DATA * POINTS PLOTTED ON SEMI-LOG PAPER. THE POINTS (D1,Q1) AND (D3,Q3) * SPECIFY THE LINE ON SEMI-LOG PAPER. QINTERP IS COMPUTED AT DAY * D2. IF Q1=Q3 OR IF Q1=0 OR IF Q3=0, THEN LOGARITHMIC INTERPOLATION * CAN NOT BE USED. IN THESE CASES THE FLOW WILL BE ASSUMED TO VARY * LINEARLY FROM DAY D1 TO DAY D3. * * VARIABLES * D1 - day 1 * D2 - day 2 (real) * D3 - day 3 * Q1 - discharge on day 1 * Q3 - discharge on day 3 * QINTERP - discharge on day 2, determined by interpolation between * day 1 and day 3 **************************************************************************** INTEGER D1, D3 REAL D2, Q1, Q3 *--------------------------------------------------------------------------- IF (D2-FLOAT(D1).LT.0.01) THEN QINTERP=Q1 ELSEIF (FLOAT(D3)-D2.LT.0.01) THEN QINTERP=Q3 ELSEIF ((Q1.EQ.-99).OR.(Q3.EQ.-99)) THEN QINTERP = -99 ELSEIF ((Q1.EQ.Q3).OR.(Q1.EQ.0.).OR.(Q3.EQ.0.)) THEN QINTERP=Q1+((Q3-Q1)/FLOAT(D3-D1))*(D2-FLOAT(D1)) ELSE QINTERP=EXP(LOG(Q1)+(LOG(Q3/Q1)/FLOAT(D3-D1))* + (D2-FLOAT(D1))) ENDIF RETURN END CHARACTER*1 FUNCTION UCASE(STR) ************************************************************************** * FUNCTION TO CONVERT FIRST CHARACTER OF A STRING TO UPPERCASE, TO * SIMPLIFY IF-THEN STATEMENTS INVOLVING USER RESPONSES TO PROMPTS. ************************************************************************** CHARACTER*(*) STR *------------------------------------------------------------------------- UCASE = STR(1:1) IF (ICHAR(UCASE).GE.97) UCASE = CHAR(ICHAR(UCASE)-32) RETURN END LOGICAL FUNCTION EXISTS(FILESPEC) ************************************************************************** * FUNCTION TO INDICATE IF FILESPEC ALREADY EXISTS ON THE SYSTEM. ************************************************************************** CHARACTER*(*) FILESPEC *------------------------------------------------------------------------- INQUIRE (FILE=FILESPEC,EXIST=EXISTS) RETURN END SUBROUTINE FILEEXISTS(FILESPEC) ************************************************************************** * SUBROUTINE TO DEAL WITH THE SITUATION WHERE ONE OF THE OUTPUT FILES * ALREADY EXISTS. ************************************************************************** CHARACTER*255 FILESPEC CHARACTER*1 A LOGICAL EXISTS *---- FUNCTIONS CALLED: CHARACTER*1 UCASE 1000 WRITE (*,1005) FILESPEC 1005 FORMAT (1X,'*** Problem with file: ',A55/1X, + '*** This file already exists.') WRITE (*,*) 1010 WRITE (*,*) 'Do you wish to overwrite this file? (Y/N)' READ (*,1020) A 1020 FORMAT (A1) A=UCASE(A) IF ((A.NE.'Y').AND.(A.NE.'N')) GOTO 1010 IF (A.EQ.'N') THEN WRITE (*,*) 'Please enter a new filename.' READ (*,1030) FILESPEC 1030 FORMAT (A255) IF (EXISTS(FILESPEC)) GOTO 1000 ENDIF RETURN END