C C This code has been slightly modified by the Integrated Decision Support Group C Colorado State University to allow it to use actual days in the month as well as C allow the input files to support running data prior to the historical period as C well as after the historical period. The main computational methodology has not C been modified by IDS (www.ids.colostate.edu). C C 1. Added option to use the actual days in the month rather than the average C days per month (see the USEAVGDAILY flag). C 2. Years that are leap years will have 29 days in February rather than C using every fourth year from the start of the simulation. C 3. Extra parameters were added to separate out projected data from C historical data. C C More changes made (3/2/04): C 1. Added option to calculate stream depletion using assumption of C average pumping and the theory of superposition. C C C ORIGINAL USGS SDF CODE FROM HRS CONSULTANTS (MODIFICATIONS BY C FRANK HAGAR 8/28/84 TO RUN ON VAX AND ALSO ALLOW ASSIGNMENT OF INPUT AND C OUTPUT FILE NAMES) RECEIVED BY NCWCD IN 1984 C C 1ST DATA CARD% FORMAT(7I5) C COLUMNS DATA ITEM(S) C 1-5 N=NUMBER OF MONTHS FROM DATAS THRU OUTE (LE 576) C 6-10 DATAS=FIRST YEAR OF PUMPING (1956 FOR SACWSD) C 11-15 DATAE=MOST RECENT YEAR FOR WHICH COMPLETE MONTHLY PUMPING C RECORDS ARE AVAILABLE C 16-20 OUTE=LAST PROJECTED YEAR OF PUMPING =(OUTE-DATAS+1 MUST C NOT EXCEED 48) C 21-25 YR1P=1ST YEAR IN WHICH DEPLETION RESULTS ARE DESIRED. C PROGRAM WILL PRINT MONTHLY DEPLETIONS BEGINING C WITH THIS YEAR AND FINISHING WITH OUTE C 26-30 NOLOC=TOTAL OF NUMBER OF MODEL WELLS AND IRRIGATION RETURN C PLOTS. C 31-35 CODE=A SINGLE DIGIT DESCRIBING DATA AND FORMAT DESIRED FOR C OUTPUT AS FOLLOWS% C 0 = OUTPUT is TABLES of the monthly PUMPING/RECHARGE C and depletions/accretions for each site. Also C Includes summary tables for all sites. C 1 = OUTPUT is TABLES of the TOTAL river depletions/ C accretions and net impact upon the river. C 2 = OUTPUT is TABLES with NO HEADINGS of the TOTAL C RIVER DEPLETIONS AND ACCRETIONS. C C A 2ND DATA CARD% FORMAT(F5.1,6A4,1X,2A1,2A4,10X,A1,A2,A4)(THIS IS THE C FIRST DATA CARD OF A DATA CARD SET FOR A PARTICULAR WELL SITE OR C IRRIGATION RETURN PLOT) C COLUMNS DATA ITEM(S) C 1-5 SDF=STREAM DEPLETION FACTOR FOR SITE OR PLOT DESCRIBED IN C THIS DATA CARD SET C 6-29 LOC=ADDRESS OR IDENTIFICATION NUMBER OF PLOT OR SITE C 31-40 TYPEW=CODE THE WORD WITHDRAWAL OR THE WORD RETURN FOR TYPE C OF PUMPING, BEGIN IN COL 31 C 51-57 WATER=CODE THE WORD SHALLOW OR THE WORD DEEP FOR TYPE OF C WATER, BEGIN IN COL 51 C NEXT (DATAE-DATAS+1) DATA CARDS% FORMAT(12F8.3) C EACH CARD REPRESENTS A SINGLE YEAR OF PUMPING WITH 12 MONTHLY PUMPING C RATES. THE MONTHLY RATES ARE IN UNITS OF ACRE FEET/DAY MONTHLY AVER- C AGES WITH THE RATES ARRANGED CHRONOLOGICALLY BY MONTH (NOV,DEC,JAN,ETC) C the following data lines are not used in Colorado MOA reregulation modeling =20 C NEXT CARD% FORMAT(12F8.3) THIS CARD REPRESENTS PROJECTED PUMPING RATES C FOR THE OPERATING YEAR FOLLOWING DATAE. THE DATA ON THIS C CARD AND THE NEXT ARE USED IN STRAIGHT LINE INTERPOLATION TO OBTAIN C PROJECTED PUMPING RATES FOR INTERVINING YEARS C NEXT CARD% FORMAT(12F8.3) THIS CARD CONTAINS PROJECTED PUMPING RATES FOR C THE OPERATING YEAR OUTE. THIS IS THE LAST DATA CARD IN THE SET DES- C CRIBING ONE PARTICULAR WELL SITE OR RETURN PLOT. C B REPEAT A UNTIL THE NUMBER OF DATA SETS FOR WELL SITES AND RETURN PLOTS C EQUALS NOLOC. C C NO OTHER DATA CARDS ARE REQUIRED. C USE MSFLIB INTEGER IN, OUT, MON, AQUIFTYP INTEGER DATAS, DATAE, OUTE, CODE, OUTS INTEGER NMONTHS, NOLOC, TYPEW, AQUIFER, SHOWIN INTEGER PROJTYPE, PROJSYR, PROJEYR, YEARTYPE INTEGER PREPROJTYPE, PREPROJSYR, PREPROJEYR INTEGER PRESKIP, POSKIP, UNITS, SDF_MODE REAL SDF, UNITCONV, PI CHARACTER*80 ARGMET CHARACTER*100 TXTSYNTH,PRETXTSYNTH C Flag indicating whether the average days in a month should be used C when converting monthly pumping to daily (otherwise the actual C days in the month will be used). INTEGER USEAVGDAILY C Number of days in the current month. REAL NDAYS INTEGER IARGC, M, Y C When using non-calendar year types, this will be the number to add C to get to the equivalent calendar month. INTEGER offset CHARACTER CTYPE*10, CAQUIFER*10 CHARACTER*45 FILIN,FILOU CHARACTER*90 LOC EXTERNAL ERF C Patterson: added to allow the user to ignore pumping or recharge C data after a certain year. INTEGER USE_IGNORE_YEAR, IGNORE_YEAR C C DIM_NY = MAXIMUM NUMBER OF YEARS ALLOWED C DIM_NM = MAXIMUM NUMBER OF MONTHS C DIM_ND = MAXIMUM NUMBER OF DAYS C INTEGER DIM_NM, DIM4, DIM_NY, DIM_ND C Be sure to modify these in subroutine PRNT as well! PARAMETER (DIM_NM = 5000) PARAMETER (DIM_NY = 400) PARAMETER (DIM_ND = DIM_NM*31) PARAMETER (DIM4 = 4*DIM_NM) PARAMETER (PI = 3.141592654) DIMENSION A(4),AD(4),QSD(DIM_NM), 1 Q(DIM_NM),ITM(DIM_NM),QR(4,DIM_NM),QSDR(4,DIM_NM), : SDFQ(DIM_ND),SDFV(DIM_ND),DAYSINMONTH(12), : IDELT(DIM_NM,DIM_NM) C 1 Q(DIM_NM),ITM(4000),QR(4,DIM_NM),QSDR(4,DIM_NM),SDFQ(40089) DATA QSD/ DIM_NM*0.000/ DATA QR/DIM4*0.0/,QSDR/DIM4*0.0/,A/4*0.0/,AD/4*0.0/ DATA DAYSINMONTH(1)/31/,DAYSINMONTH(2)/28/,DAYSINMONTH(3)/31/, $ DAYSINMONTH(4)/30/,DAYSINMONTH(5)/31/,DAYSINMONTH(6)/30/, $ DAYSINMONTH(7)/31/,DAYSINMONTH(8)/31/,DAYSINMONTH(9)/30/, $ DAYSINMONTH(10)/31/DAYSINMONTH(11)/30/,DAYSINMONTH(12)/31/ C-----Handle command line argument, prompt the user if argument is C-----not provided IF (IARGC().EQ.0) THEN PRINT *, 'ENTER THE INPUT DATA FILE NAME (45 CHARACTERS): ' PRINT * READ(*,11) FILIN 11 FORMAT(A45) PRINT *, 'ENTER THE OUTPUT FILE NAME (45 CHARACTERS): ' PRINT * READ(*,11) FILOU ELSE CALL GETARG(1,argmet) FILIN = argmet CALL GETARG(2,argmet) FILOU = argmet ENDIF PRINT *, "Running..." IN = 10 OUT = 20 OPEN(IN,FILE=FILIN,STATUS='OLD') OPEN(OUT,FILE=FILOU,STATUS='UNKNOWN') !OPEN(11, FILE="sdf.log",STATUS='UNKNOWN') C C READ AND WRITE INPUT DATA C C READ(IN,*) N,DATAS,DATAE,OUTE,YR1P,NOLOC,CODE C C DATAS = First year of pumping for dataset C DATAE = Last year for which pumping is available C NPRESIMYRS = Number of pre-synthesized years C NSIMYRS = Number of synthesized years C OUTS = First year of output C OUTE = Last year of output C NOLOC = Number of wells C CODE = Code for the type format for the output to be used C PREPROJTYPEW = Type of project to be used C 0 = Average (enter PREPROJSYR, PREPROJEYR) start and end year C to use for averaging. C 1 = Given year (enter PREPROJSYR) as given year C 2 = Historical Cycle (enter PREPROJSYR, PREPROJEYR) start and C end year of historical cycle. C PREPROJSYR = projected start year C PREPROJEYR = projected end year C PROJTYPEW = Type of project to be used C 0 = Average (enter PROJSYR, PROJEYR) start and end year C to use for averaging. C 1 = Given year (enter PROJSYR) as given year C 2 = Historical Cycle (enter PROJSYR, PROJEYR) start and C end year of historical cycle. C PROJSYR = projected start year C PROJEYR = projected end year C YEARTYPE = Type of year C 0 = Calendar year C 1 = Irrigation (November to October) C 2 = USGS (October to September) C UNITS = Units of input C 0 = Daily Average for a Month Ac-Ft/Day C 1 = Total Monthly value Ac-Ft C USEAVGDAILY = Flag indicating whether to use the actual days in the C to convert from monthly to daily pumping, or use the average C months in a year (30.417) C 0 = Use the actual days in the month to make the conversion. C 1 = Use 30.417 to convert from monthly to average daily C pumping. C USE_IGNORE_YEAR = Flag for whether well pumping should be ignored after C a given year. C 0 = Don't use IGNORE_YEAR parameter C 1 = Use IGNORE_YEAR parameter. C IGNORE_YEAR = The last year that pumping will be considered in the C calculations. C C SDF_MODE = Flag indicating whether to assume all pumping occurs in the C middle of the month or constant daily. C 0 = Use middle of the month C 1 = Use cosntant pumping. C First line is version number. Modify this in the future to check C for version strings. READ(IN,*) READ(IN,*) DATAS,DATAE,NPRESIMYRS,NSIMYRS,OUTS,OUTE,NOLOC,CODE, : PREPROJTYPE,PREPROJSYR,PREPROJEYR,PROJTYPE,PROJSYR, : PROJEYR,YEARTYPE,UNITS,USEAVGDAILY, : USE_IGNORE_YEAR,IGNORE_YEAR, SDF_MODE C If we aren't going to be using the ignore year, then set it to C zero. This will help reduce the number of parameters that C need to be passed to subroutines. IF (USE_IGNORE_YEAR .EQ. 0) THEN IGNORE_YEAR = DATAE ENDIF C Build synthesized year output text string. IF (NPRESIMYRS > 0) THEN IF(PREPROJTYPE .EQ. 1) THEN WRITE(PRETXTSYNTH,700) DATAS+PREPROJSYR,DATAS+PREPROJEYR 700 FORMAT('AVERAGE OF HISTORICAL DATA FROM',I5,' TO ',I5) ELSEIF(PREPROJTYPE .EQ.2) THEN WRITE(PRETXTSYNTH,710) DATAS+PREPROJSYR 710 FORMAT('HISTORICAL DATA YEAR',I5) ELSEIF(PREPROJTYPE .EQ.3) THEN WRITE(PRETXTSYNTH,720) DATAS+PREPROJSYR, : DATAS+PREPROJSYR + PREPROJEYR 720 FORMAT('HISTORICAL CYCLE FROM YEAR',I5,' TO ',I5) ELSEIF(PREPROJTYPE .EQ.4) THEN WRITE(PRETXTSYNTH,730) 730 FORMAT('USER-DEFINED DATA') ENDIF ENDIF IF (NSIMYRS > 0) THEN IF(PROJTYPE .EQ. 1) THEN WRITE(TXTSYNTH,700) DATAS+PROJSYR,DATAS+PROJEYR ELSEIF(PROJTYPE .EQ.2) THEN WRITE(TXTSYNTH,710) DATAS+PROJSYR ELSEIF(PROJTYPE .EQ.3) THEN WRITE(TXTSYNTH,720) DATAS+PROJSYR, DATAS+PROJSYR+PROJEYR ELSEIF(PROJTYPE .EQ.4) THEN WRITE(TXTSYNTH,730) ENDIF ENDIF C C If UNITS = 0 then pumping in right units conversion is one. C If UNITS = 1 then pumping in monthly total divide by average C number of days in a month. C IF ( UNITS .EQ. 0) THEN UNITCONV = 1.0 ELSE UNITCONV = 30.417 ENDIF CIDS NOYRS= OUTE-DATAS+1 CIDS NOYRS= SIMEYR-OUTS+1 C C If the output is longer than the input + synthesized data then C the implication is that you have extra years with zero pumping. C This is reflected by having NEXTRAYRS greater than zero. C IF( OUTE .GT. (DATAE + NSIMYRS)) THEN NOYRS = (DATAE + NSIMYRS) - OUTS + 1 NEXTRAYRS = OUTE - (DATAE + NSIMYRS) ELSE NOYRS = OUTE - OUTS + 1 NEXTRAYRS = 0 ENDIF C C Check to make sure that the input span doesn not exceed DIM_NY years C IF(NOYRS + NEXTRAYRS .GT. DIM_NY) THEN C open(unit=1, file='ERR1') PRINT *, "ERROR: input contains more than ", DIM_NY, : " years of data." CALL EXIT(-1) C STOP "ERROR: input contains more than ", DIM_NY, C : " years of input." ENDIF CIDS CALL OUTPUT(OUT,SIMBYR,DATAS,NOYRS,OUTS,OUTE) CALL OUTPUT(OUT,DATAS,DATAE,NPRESIMYRS,NSIMYRS,OUTS,OUTE, : PREPROJSYR,PREPROJEYR,PROJSYR,PROJEYR, : PRETXTSYNTH,TXTSYNTH,USEAVGDAILY,IGNORE_YEAR) C WRITE(OUT,410) DATAE, DATAS, NOYRS, OUTS, OUTE C WRITE(OUT,*) DATAE, DATAS, NOYRS, OUTS, OUTE C CIDS Original way to calculate the cumulative number of days. This C computation was updated by IDS and is now done later in the code. C Part of the update was to allow the code to handle leap years. C C *** Determine values of ITM for arrays 14 - 636 *** C *** Calculates the begining day of each month for C *** all the months in the simulation C *** This loops first do three years of 365 and then C *** one year of 366 (loop 104). C C DO 101 I=1,14 C DO 101 I=1,DIM_NY/4 C DO 102 J=1,3 C DO 102 K=1,12 C KK=((I*4-3)*12+(J-1)*12+K+3) C ITM(KK)=ITM(KK-12)+365 C 102 CONTINUE C DO 104 K=1,12 C KK=((I*4-3)*12+36+K+3) C ITM(KK)=ITM(KK-12)+366 C 104 CONTINUE C 101 CONTINUE C INITIALIZE THE VARIABLES QR, AND QSDR FOR ITERATION C ADDED BY NENGJIN C DO 111 IH= 1,4 DO 112 IP=1,DIM_NM QR(IH,IP) =0.0 QSDR(IH,IP)=0.0 112 CONTINUE 111 CONTINUE C C *** For individual sites READ in PUMPING or RECHARGE *** C DO 152 IJK=1,NOLOC READ(IN,*) SDF, TYPEW, SHOWIN C C SHOWIN tells the program to use or ignore a well C IF SHOWIN is 0 ignore the well otherwise use it C IF( SHOWIN .NE. 0) THEN READ(IN,350) LOC WRITE(*, *) LOC C C SET THE TYPEW TO PUMPING AND AQUIFER TO SHALLOW C IF( TYPEW .EQ. 1 ) THEN CTYPE = 'RECHARGE' ELSE CTYPE = 'PUMPING' ENDIF C C FOR NOW WE ARE HARD WIRED TO ALWAYS USE A SHALLOW AQUIFER C AQUIFER = 1 IF( AQUIFER .EQ. 1 ) THEN CAQUIFER = 'SHALLOW' AQUIFTYP = 1 ELSE CAQUIFER = 'DEEP' AQUIFTYP = 2 ENDIF IF (CODE .EQ. 0) WRITE (OUT,699) IF (CODE .EQ. 0) WRITE(OUT,690) LOC, SDF, + CTYPE,CAQUIFER C C Compute the # of months (NMONTHS) C ELSE READ(IN,*) ENDIF PRESKIP = (OUTS - DATAS + NPRESIMYRS) POSKIP = (DATAE + NSIMYRS) - OUTE CIDS N = ((OUTE - SIMBYR)+1)*12 NMONTHS = NOYRS*12 + (NEXTRAYRS*12) C C IDS New algorithm for building the ITM array. C Build ITM array based on year type. offset = 0 if (YEARTYPE .eq. 1) then offset = 10 elseif (YEARTYPE .eq. 2) then offset = 9 endif ITM(1) = DAYSINMONTH(offset + 1) DO 104 K=2,NMONTHS m = mod(K-1+offset, 12)+1 ! m == month in calendar year mode y = (K-1)/12 ITM(K) = ITM(K-1) + DAYSINMONTH(m) C Make a check for leap year IF (MOD(DATAS+y, 4) .EQ. 0) THEN IF (m .EQ. 2) THEN ITM(K) = ITM(K) + 1 ENDIF ENDIF 104 CONTINUE NODAYS=ITM(NMONTHS) C C *** Read PUMPING or RECHARGE from 1st year to PRESENT C ALSO 1st year of PROJECTED Pumping or Recharge. C DO 65 J = 1,NMONTHS Q(J) = 0.0 QSD(J) = 0.0 65 CONTINUE C C Skip the number of historical years prior to the begining C of the simulation (simulation years must equal the output years). C DO 67 J=1,PRESKIP READ(IN,*) 67 CONTINUE C C Read the number of years to use for the simulation C (simulation years must equal the output years). C IF(SHOWIN .NE. 0) THEN NDAYS = DAYSINMONTH(offset+1) DO 70 J=1,NOYRS IF ((USE_IGNORE_YEAR .EQ. 1) .AND. : (OUTS + J - 1 .GT. IGNORE_YEAR)) THEN PRINT *, "Skipping historical year ", OUTS+J READ(IN,*) ELSE READ(IN,*) (Q(I+(J-1)*12),I=1,12) DO 71 MON = 1, 12 CIDS Q(MON+(J-1)*12) = Q(MON+(J-1)*12) / UNITCONV IF ( UNITS .EQ. 0) THEN IF (USEAVGDAILY .EQ. 1) THEN Q(MON+(J-1)*12) =Q(MON+(J-1)*12) * UNITCONV ELSEIF (USEAVGDAILY .EQ. 0) THEN Q(MON+(J-1)*12) =Q(MON+(J-1)*12) * NDAYS C Set the number of days to the next month's. NDAYS = ITM((J-1)*12+MON+1)-ITM((J-1)*12+MON) ENDIF ENDIF 71 CONTINUE ENDIF 70 CONTINUE ELSE DO 75 J=1,NOYRS READ(IN,*) 75 CONTINUE ENDIF C C Skip the number of historical years after the end C of the simulation (simulation years must equal the output years). C DO 68 J=1,POSKIP READ(IN,*) 68 CONTINUE C C If the output is longer than then input plus the synthesized data C then we need to add some zero pumping to the end. C IF( SHOWIN .NE. 0) THEN IF( NEXTRAYRS .GT. 0) THEN DO 72 J=NOYRS+1,NOYRS+NEXTRAYRS DO 73 I = 1, 12 Q(I+(J-1)*12) = 0.0 73 CONTINUE 72 CONTINUE ENDIF C C C **** CALCULATE UNIT SDF RESPONSE CURVE *** C ! (50,FILE='sdfq.out',STATUS='unknown') DO 130 J=1,NODAYS TIME=FLOAT(J) Z=0.5*SQRT(SDF/TIME) SDFV(J) = ((SDF/TIME/2.0) + 1.0) * ( 1.0 - ERF(Z)) : - (Z * (2.0*EXP(-(SDF/TIME/4.0)))/SQRT(PI)) IF (Z.LE.3.00) THEN SDFQ(J)=1.-ERF(Z) ELSE SDFQ(J)=0.0 ENDIF IF (SDF .GT. 0) then SDFV(J) = SDFV(J) * TIME / SDF ENDIF !write(50,1000) J,J/sdf,sdfq(J) C1000 format(I5,2x,f5.2,2x,f5.3) 130 CONTINUE close(50) C C CALCULATE STREAM DEPLETION USING INDICATED PUMPING SCHEDULE AND C THE THEORY OF SUPERPOSITION C IT1=0 Q1=0.0 I=1 Q2=Q(1) IDAY = 0 NDAYS = ITM(1) DO 131 I=1,NMONTHS if (I .GT. 1) then NDAYS = ITM(I) - ITM(I-1) endif IF (I .EQ. 1) THEN IF (SDF_MODE .EQ. 0) THEN IDELT(1,I) = (NDAYS - 15) ELSE IDELT(1,I) = NDAYS ENDIF ELSE IDELT(1,I) = IDAY + NDAYS ENDIF IDAY = IDELT(1,I) 131 CONTINUE DO 132 J=2,NMONTHS IDAY = 0 NDAYS = ITM(J) - ITM(J-1) DO 133 I = 1,NMONTHS IF (I .GT. 1) THEN NDAYS = ITM(J+I-1) - ITM(J+I-2) ENDIF IF (I .EQ. 1) THEN IF (SDF_MODE .EQ. 0) THEN IDELT(J,I) = (NDAYS - 15) ELSE IDELT(J,I) = NDAYS ENDIF ELSE IDELT(J,I) = IDAY + NDAYS ENDIF IDAY = IDELT(J,I) 133 CONTINUE 132 CONTINUE C IF ASSUMPTION IS THAT PUMPING OCCURS DURING THE MIDDLE OF THE MONTH C (SDF_MODE .EQ. 0) IF (SDF_MODE .EQ. 0) THEN DO 141 J=1,NMONTHS c$$$ write(11,*)"Loop 141: J = ",J,"IDELT = ",IDELT(1,J) IF (IDELT(1,J) .LT. 1) THEN WRITE(*, *) "IDELT OUT OF RANGE: ",IDELT(1,J) ENDIF QSD(J)=QSD(J)-(Q1-Q2)*SDFQ(IDELT(1,J)) 141 CONTINUE Q1=Q2 DO 150 I=2,NMONTHS Q2=Q(I) CIDS CODE WAS CHANGED SO THAT LOOP 140 REPRESENTS THE DEPLETIONS FOR CIDS FOR THE PUMPING DURING THE CURRENT MONTH DO 140 J=I,NMONTHS M = mod(J-1+offset, 12)+1 ! m == month in calendar year mode c$$$ write(11,*) "Loop 141: I = ", I, "J = ", J, c$$$ $ "IDELT = ",IDELT(I,J-I+1) IF (IDELT(I,J-I+1) .LT. 1) THEN WRITE(*,*) "IDELT OUT OF RANGE: ",IDELT(I,J-I+1) ENDIF QSD(J)=QSD(J)+(Q2*SDFQ(IDELT(I,J-I+1))) 140 CONTINUE CIDS LOOP 143 REPRESENTS THE SUPERPOSITION (IMAGINARY WELL THAT IS ASSUME CIDS TO RECHARGE AT THE SAME RATE AS PUMPING WELL OFFSET BY ONE MONTH) DO 143 J=I,NMONTHS c$$$ write(11,*) "Loop 143: I = ", I, "J = ", J, c$$$ $ "IDELT = ",IDELT(I-1,J-I+1) IF (IDELT(I-1,J-I+1) .LT. 1) THEN WRITE(*,*)"IDELT OUT OF RANGE: ",IDELT(I-1,J-I+1) ENDIF QSD(J)=QSD(J)-(Q1*SDFQ(IDELT(I-1,J-I+1))) 143 CONTINUE Q1=Q2 150 CONTINUE ELSE C C CALCULATE STREAM DEPLETION USING ASSUMPTION OF AVERAGE PUMPING AND C THE THEORY OF SUPERPOSITION C DO 151 I=1,NMONTHS CUM_TIME = 0 QSD(I)=0 DO 145 J=1,I IF (J .GT. 1) THEN D_IN_J = ITM(J) - ITM(J-1) ELSE D_IN_J = ITM(J) ENDIF S = SDFV(ITM(I)-CUM_TIME); IF (ITM(I)-CUM_TIME - D_IN_J .GE. 1) then s = s - SDFV(ITM(I)-CUM_TIME - D_IN_J) ENDIF IF (SDF .GT. 0) THEN S = S * SDF/D_IN_J ENDIF QSD(I) = QSD(I) + S * Q(J) CUM_TIME = CUM_TIME + D_IN_J 145 CONTINUE 151 CONTINUE C C We have calculated the cumulative total, so subtract the previous. C IF (SDF .GT. 0) THEN DO I=NMONTHS,2,-1 QSD(I) = QSD(I) - QSD(I-1) END DO END IF ENDIF C C *** Add INDIVIDUAL SITE to TOTAL of other sites. ***** C *** IF CODE = 0, WRITE to output INDIVIDUAL site. *** C IF(TYPEW .EQ. 1) THEN C C ***** RECHARGE Sites *** C DO 123 MON=1,NMONTHS Q(MON) = TYPEW*(Q(MON)) QSD(MON) = TYPEW*(QSD(MON)) QR(AQUIFTYP,MON)=Q(MON)+QR(AQUIFTYP,MON) QR(4,MON)=QR(4,MON)+Q(MON) QSDR(4,MON)=QSDR(4,MON)+QSD(MON) QSDR(AQUIFTYP,MON)=QSDR(AQUIFTYP,MON)+QSD(MON) 123 CONTINUE C IF (CODE .EQ. 0) THEN WRITE (OUT,600) CALL PRNT (OUT, NOYRS+NEXTRAYRS, CODE, DATAS, OUTS, Q, $ offset) WRITE (OUT,610) CALL PRNT (OUT, NOYRS+NEXTRAYRS, CODE, DATAS, OUTS, $ QSD, offset) ENDIF ELSE C C *** PUMPING SITES *** C DO 154 MON=1,NMONTHS Q(MON) = TYPEW*(Q(MON)) QSD(MON) = TYPEW*(QSD(MON)) QR(3,MON)=QR(3,MON)+Q(MON) QSDR(3,MON)=QSDR(3,MON)+QSD(MON) QR(4,MON)=QR(4,MON)+Q(MON) QSDR(4,MON)=QSDR(4,MON)+QSD(MON) 154 CONTINUE C IF (CODE .EQ. 0) THEN WRITE (OUT,625) CALL PRNT (OUT, NOYRS+NEXTRAYRS, CODE, DATAS, OUTS, Q, $ offset) WRITE (OUT,635) CALL PRNT (OUT, NOYRS+NEXTRAYRS, CODE, DATAS, OUTS, $ QSD, offset) ENDIF ENDIF ENDIF 152 CONTINUE C C *** WRITE TOTALS to the Output Files *** C DO 280 ITEM = 1,4 WRITE (OUT,699) DO 200 N1 = 1,NMONTHS Q(N1) =QR(ITEM,N1) QSD(N1) = QSDR(ITEM,N1) 200 CONTINUE C IF( ITEM .EQ. 1) THEN IF (CODE .NE. 2) THEN WRITE (OUT,660) WRITE (OUT,600) CALL PRNT (OUT, NOYRS+NEXTRAYRS, CODE, DATAS, OUTS, Q, $ offset) IF (NOYRS .GT. 15) WRITE (OUT,699) WRITE (OUT,660) WRITE (OUT,610) ENDIF CALL PRNT (OUT, NOYRS+NEXTRAYRS, CODE, DATAS, OUTS, QSD, $ offset) ELSEIF( ITEM .EQ. 3 ) THEN 230 IF (CODE .NE. 2) THEN WRITE (OUT,660) WRITE (OUT,625) CALL PRNT (OUT, NOYRS+NEXTRAYRS, CODE, DATAS, OUTS, Q, $ offset) IF (NOYRS .GT. 15 ) WRITE (OUT,699) WRITE (OUT,660) WRITE (OUT,635) ENDIF CALL PRNT (OUT, NOYRS+NEXTRAYRS, CODE, DATAS, OUTS, QSD, $ offset) ELSEIF(ITEM .EQ. 4) THEN 240 IF (CODE .NE. 2) THEN WRITE (OUT,650) CALL PRNT (OUT, NOYRS+NEXTRAYRS, CODE, DATAS, OUTS, QSD, $ offset) ENDIF ENDIF 280 CONTINUE PRINT *, "Finished." C C *** FORMAT STATEMENTS *** C 361 FORMAT(12F8.3) 350 FORMAT(A90) 352 FORMAT(20X,F8.1,20X, 6A4, 4X,2A1,2A4,15X,A1,A2,A4) c 410 FORMAT(54X,'STREAM DEPLETION' c 1/45X,'CALCULATED BY U.S.G.S. SDF METHOD'///// c 210X,'THE LAST YEAR OF HISTORICAL DATA IS',I5,'.' c 3//10X,'NUMBER OF YEARS OF SYNTHESIZED DATA ARE ',I5, c 4//10X,'SYNTHESIZED DATA WAS GENERATED USING ',A50, c 4//10X,'PUMPING BEGINS IN' c 3 ,I5,' AND CONTINUES',I3,' YEARS.'//10X,'PRINT OUT OF DEPLETIONS c 4 IS FOR YEARS ',I4,' TO ',I4,' INCLUSIVE.') 411 FORMAT(//54X,'STREAM DEPLETION' 1/45X,'CALCULATED BY U.S.G.S. SDF METHOD'///// 210X,'THE LAST YEAR OF RECORD IS',I5,'.'//10X,'PUMPING BEGINS IN' 3 ,I5,' AND CONTINUES',I3,' YEARS.'//10X,'PRINT OUT OF DEPLETIONS 4 IS FOR YEARS ',I4,' TO ',I4,' INCLUSIVE.'//10X, 'DATA FOR WEL 5LS ARE%'//15X,'STREAM DEPLETION FACTOR',13X,'LOCATION',19X,'TYPE', 620X,'WATER'/) C C 600 FORMAT(/'NET DIVERSION INTO RECHARGE SITES'/' ( acre feet )'/1X) 610 FORMAT(/'ACCRETION TO THE STREAM RESULTING FROM RECHARGE'/ 1 ' ( acre feet )'/1X) 625 FORMAT(/'DIVERSIONS FROM WELLS'/' ( acre feet ) '/1X) 635 FORMAT(/'DEPLETIONS OCCURRING FROM THE STREAM DUE TO WELLS'/ 1 ' ( acre feet )'/1X) 650 FORMAT(/'NET IMPACT UPON THE STREAM'/' Depletions + Accretions'/ 1 ' ( acre feet )'/1X) 660 FORMAT('SUMMARY OF ALL LOCATIONS') 690 FORMAT('LOCATION: ',A60/'SDF FACTOR: ',F9.1/ 1 'TYPE OF SITE: ',A10/ 2 'DEPTH:',11X,A10/1X) 699 FORMAT(' ') 412 FORMAT(I5) C *** CLOSE FILES *** C CALL FILES(2,1,1) C C CALL EXIT CLOSE(IN) CLOSE(OUT) STOP END C C FUNCTION ERF(X) AX=ABS(X) T=1.0/(1.0+.3275911*AX) D=EXP(-X**2.) P=1.0-D*T*((((1.061405429*T-1.453152027)*T+1.421413741)*T- 1 0.284496736)*T+0.254829592) IF(X)2,1,1 1 ERF=P 2 RETURN END C *************************************************************************= C C THIS SUBROUTINE WRITES OUT THE DATA IN TABULAR FORM C C *************************************************************************= C SUBROUTINE PRNT(OUT, YRS, CODE, IYR, JYR, Q, m_offset) C INTEGER DIM_NM, DIM4, DIM_NY PARAMETER (DIM_NM = 5000) PARAMETER (DIM_NY = 400) PARAMETER (DIM_ND = DIM_NM*31) C PARAMETER (DIM_NM = 2500) C PARAMETER (DIM_NY = 200) PARAMETER (DIM4 = 4*DIM_NM) INTEGER OUT, CODE, YRS, m_offset DIMENSION SUM(12), AVG(12), YRT(DIM_NY), Q(DIM_NM), MON(12) DATA MON/ 'JAN','FEB','MAR','APR','MAY', 1 'JUN','JUL','AUG','SEP','OCT','NOV','DEC'/ C IF (CODE .NE. 2) WRITE (20,80) (MON(mod(M+m_offset, 12)+1), $ M = 0,11) CIDS XNYR = YRS - (JYR - IYR) TOT = 0.0 DO 20 I = 1,12 SUM(I) = 0.0 DO 10 J = 1,YRS CIDS IK = IYR + J -1 CIDS IF (IK .GE. JYR) THEN SUM(I) = SUM(I) + Q(12*(J-1)+I) CIDS ENDIF 10 CONTINUE CIDS AVG(I) = SUM(I) / XNYR AVG(I) = SUM(I) / YRS TOT = TOT + SUM(I) 20 CONTINUE C CIDS ATOT = TOT / XNYR ATOT = TOT / YRS DO 30 J = 1,YRS YRT(J) = 0.0 DO 30 I = 1,12 YRT(J) = YRT(J) + Q(12*(J-1)+I) 30 CONTINUE C 35 DO 40 J = 1,YRS C C ONLY APPLICABLE WHEN THE SIMULATION YEARS ARE DIFFERENT THAN THE PRINT C YEARS FOR NOW THEY ARE THE SAME. C CIDS IY = IYR + J - 1 CIDS IF (IY .LT. JYR) GO TO 40 CIDS WRITE(OUT,50) IY,(Q(12*(J-1)+I), I=1,12),YRT(J) C Mod DAP added in the monthly conversion to test old-style model. Replace all 30.417 to restore WRITE(OUT,50) (JYR+J-1),(Q(12*(J-1)+I), I=1,12), 1 YRT(J) 40 CONTINUE WRITE(OUT,60) (AVG(I), I=1,12),ATOT C 50 FORMAT(I4, 12(F12.4, 1X), 1X, F12.4) 60 FORMAT(/' AVG', 13(F12.4, 1X)) 70 FORMAT(/' AVG', 12(F12.4, 1X)) 80 FORMAT ('YEAR', 12(3X,A3,3X), 4X,'TOTAL') C 999 RETURN END SUBROUTINE OUTPUT(OUT,DATAS,DATAE,NPRESIMYRS,NSIMYRS,OUTS,OUTE, : PREPROJSYR,PREPROJEYR,PROJSYR,PROJEYR,PRETXTSYNTH,TXTSYNTH, : USEAVGDAILY,IGNORE_YEAR) INTEGER OUT, DATAS, DATAE, NSIMYRS, OUTS, OUTE, PROJSYR, PROJEYR INTEGER NPRESIMYRS, PREPROJSYR, PREPROJEYR INTEGER N, USEAVGDAILY CHARACTER*100 PRETXTSYNTH,TXTSYNTH INTEGER PUMP_YEARS C Patterson: IGNORE_YEAR == DATAE if USE_IGNORE_YEAR is off, C so this should be safe. IF(IGNORE_YEAR+N .LT. OUTE) THEN PUMP_YEARS = IGNORE_YEAR+N-OUTS ELSE PUMP_YEARS = OUTE-OUTS END IF C If the output start years includes no historical data, then C reduce the number of synthesized years accordingly. N = NPRESIMYRS + NSIMYRS IF(DATAE .LT. OUTS) THEN N = NSIMYRS + NPRESIMYRS - (OUTS-DATAE) IF(N .LT. 0) N = 0 ENDIF WRITE(OUT,410) DATAS,DATAE,NPRESIMYRS,NSIMYRS,PRETXTSYNTH, : TXTSYNTH,OUTS,PUMP_YEARS,OUTS,OUTE IF (USEAVGDAILY .EQ. 1) THEN WRITE(OUT,511) ELSE WRITE(OUT,512) ENDIF 410 FORMAT(54X,'STREAM DEPLETION' 1/45X,'CALCULATED BY U.S.G.S. SDF METHOD'///// 210X,'THE FIRST YEAR OF HISTORICAL DATA IS',I5,'.' 2//10X,'THE LAST YEAR OF HISTORICAL DATA IS',I5,'.' 3//10X,'NUMBER OF YEARS OF PREHISTORIC DATA ARE ',I5, 3//10X,'NUMBER OF YEARS OF POSTHISTORIC DATA ARE ',I5, 4//10X,'PREHISTORIC DATA WAS GENERATED USING ',A50, 4//10X,'POSTHISTORIC DATA WAS GENERATED USING ',A50, 4//10X,'PUMPING BEGINS IN' 3 ,I5,' AND CONTINUES',I3,' MORE YEARS.'//10X, 4'PRINT OUT OF DEPLETIONS IS FOR YEARS ',I4,' TO ',I4, 5' INCLUSIVE.') 511 FORMAT(10X,'MONTHLY TO DAILY CONVERSIONS USED ', 1 'AVERAGE DAYS PER MONTH (30.417).') 512 FORMAT(10X,'MONTHLY TO DAILY CONVERSIONS USED ', 1 'ACTUAL DAYS IN THE MONTH.') RETURN END