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

