C*******************************************************************************
C  MAKE.PUMPAGE.F77  -- A PROGRAM TO CREATE PUMPAGE ARRAYS FOR EACH LAYER OF
C  THE VALLEYWIDE GROUND-WATER MODEL.  PUMPAGES ARE READ IN POSITIVE FOR
C  DISCHARGE AND OUTPUT POSITIVE FOR DISCHARGE.  VALUES ARE CONVERTED TO
C  NEGATIVE VALUES BY THE PROGRAM PRE.CONSTFLUX.F77
C
C               BY   WES DANSKIN
C                    LAST REVISION 11/14/89
C
C  BOTH STEADY STATE AND TRANSIENT CALIBRATION AS WELL AS MANAGEMENT SCENARIOS
C  WERE RUN WITH THIS PROGRAM.  TO FIND DIFFERENCES SEARCH FOR SCENARIO.
C*******************************************************************************
C  DIMENSIONS REQUIRED
C     30 = NUMBER OF YEARS
C    200 = NUMBER OF WELLS
C    180 = NUMBER OF MODEL ROWS
C     40 = NUMBER OF MODEL COLUMNS
C      2 = NUMBER OF MODEL LAYERS
C      9 = NUMBER OF WELL FIELDS
C*******************************************************************************
      DIMENSION AWELL(200),IL(200),IR(200),IC(200),IBOUND(180,40,2),
     - PERF(200,2),FLUX(180,40,2),TOTPUM(2),PERCNT(200,2),AVEPMP(200,2),
     - QMAX(200),AVETOT(200),NYRAVE(200),QALL(200,30),QMIN(200),
     - WFNAME(10),WFUP(10,30),WFLO(10,30),WFTOT(10,30),WFID(200),
     - WYNUM(30),WFAVE(10),QUP(200,30),QLO(200,30)
      CHARACTER*80 BLANK1,BLANK2,HEADN
      CHARACTER*20 TITLE
      CHARACTER*5 AWELL
      CHARACTER*3 AYES,AFLAG
      CHARACTER*2 WFNAME,WFID
      COMMON /BIG/ IBOUND,PERF,FLUX,AVEPMP,PERCNT,QMAX,AVETOT,NYRAVE,
     - QALL,QMIN,WFUP,WFLO,WFTOT,WFID,WFNAME,WYNUM,WFAVE,QUP,QLO
C  READ CONTROL PARAMETERS
      READ(29,911)HEADNG
      WRITE(6,911)HEADNG
      READ(29,911)HEADNG
      WRITE(6,911)HEADNG
      READ(29,912)NLAY,NROW,NCOL,NYEARS
      WRITE(6,912)NLAY,NROW,NCOL,NYEARS
C  DEFINE CONSTANTS AND SIZE OF GRID LOOPS
C      NYEARS=26
C      NROW=180
C      NCOL=40
C      NLAY=2
       NWF=10
C DEFINE YEARS FOR MANAGEMENT SCENARIOS
C      NYEARS=1
      NWELL=191
      IONE=1
      ITWO=2
      AYES='YES'
      WFNAME(1)=' ?'
      WFNAME(2)=' L'
      WFNAME(3)=' B'
      WFNAME(4)='BP'
      WFNAME(5)='TA'
      WFNAME(6)='TS'
      WFNAME(7)='IO'
      WFNAME(8)='SS'
      WFNAME(9)='BG'
      WFNAME(10)='LP'
C  ZERO VARIABLES AND MATRIXES
      IFLAG=0
      AFLAG='   '
      TITLE='                    '
      TOTPUM(1)=0.
      TOTPUM(2)=0.
      TPUMP=0.
      DO 10 I=1,NROW
      DO 10 J=1,NCOL
      DO 10 K=1,NLAY
      IBOUND(I,J,K)=0
      FLUX(I,J,K)=0.
   10 CONTINUE
      DO 20 I=1,NWELL
      IL(I)=0
      IR(I)=0
      IC(I)=0
      PERF(I,1)=0.
      PERF(I,2)=0.
      PERCNT(I,1)=0.
      PERCNT(I,2)=0.
      AVEPMP(I,1)=0.
      AVEPMP(I,2)=0.
      AVETOT(I)=0.
      NYRAVE(I)=0
      QMAX(I)=0.
      QMIN(I)=99999.
      DO 20 K=1,NYEARS
      QUP(I,K)=0.
      QLO(I,K)=0.
      QALL(I,K)=0.
      WYNUM(K)=0.
      DO 20 J=1,NWF
      WFAVE(J)=0.
      WFUP(J,K)=0.
      WFLO(J,K)=0.
      WFTOT(J,K)=0.
   20 CONTINUE
      DO 25 I=1,NYEARS
      WYNUM(I)=1962.25+I
   25 CONTINUE
C  READ IN IBOUND ARRAY FOR EACH LAYER
      DO 30 I=1,NROW
       READ(51,905)(IBOUND(I,J,1),J=1,NCOL)
       READ(52,905)(IBOUND(I,J,2),J=1,NCOL)
   30 CONTINUE
C  READ FLAG TO USE PERCENTAGES (0), OR PERFORATIONS (1)
      READ(5,901)IFLAG
      WRITE(6,901)IFLAG
C  READ IN WELL NUMBER, ROW, COLUMN, PERCENTAGE, DONE FLAG
C  PERFORATED INTERVAL (FT), HYDRAULIC CONDUCTIVITY
      READ(5,900)BLANK1,BLANK2
      WRITE(6,900)BLANK1,BLANK2
      READ(5,900)BLANK1,BLANK2
      WRITE(6,900)BLANK1,BLANK2
      DO 50 I=1,NWELL
C  READ STATEMENTS FOR USING PERCENTAGE DATA DIRECTLY OR
C  READ STATEMENTS FOR USING PERFORATION DATA TO CALCULATE PERCENTS
      READ(5,908)AWELL(I),IR(I),IC(I),PERCNT(I,1),PERCNT(I,2),AFLAG,
     - PERF(I,1),PERF(I,2),XKUP,XKLO,WFID(I)
      WRITE(6,908)AWELL(I),IR(I),IC(I),PERCNT(I,1),PERCNT(I,2),AFLAG,
     - PERF(I,1),PERF(I,2),XKUP,XKLO,WFID(I)
C  CHECK TO MAKE SURE NODES ARE ACTIVE
      IF(IBOUND(IR(I),IC(I),1).EQ.0) WRITE(95,914)IR(I),IC(I),IONE
      IF(IBOUND(IR(I),IC(I),2).EQ.0) WRITE(95,914)IR(I),IC(I),ITWO
C  CALCULATE PERCENTAGE OF WITHDRAWAL FROM EACH LAYER
      IF(IFLAG.EQ.1)THEN
        IF(AFLAG.EQ.AYES)THEN
          TOTUP=PERF(I,1)*XKUP
          TOTLO=PERF(I,2)*XKLO
          TOTPER=TOTUP+TOTLO
          PERCNT(I,1)=TOTUP/TOTPER*100.
          PERCNT(I,2)=TOTLO/TOTPER*100.
        ELSE
          PERCNT(I,1)=PERCNT(I,1)
          PERCNT(I,2)=PERCNT(I,2)
        ENDIF
      ELSE
          PERCNT(I,1)=PERCNT(I,1)
          PERCNT(I,2)=PERCNT(I,2)
      ENDIF
      WRITE(6,908)AWELL(I),IR(I),IC(I),PERCNT(I,1),PERCNT(I,2)
   50 CONTINUE
C  SKIP HEADER INFORMATION FOR MANAGEMENT SCENARIOS
      READ(41,900)BLANK
      READ(41,900)BLANK
      READ(41,900)BLANK
C  CALCULATE PUMPAGE FOR EACH YEAR
      DO 300 K=1,NYEARS
C  ZERO ANNUAL VARIABLES
      TOTPUM(1)=0.
      TOTPUM(2)=0.
      TPUMP=0.
      CALL ZERO3D(FLUX,NROW,NCOL,NLAY)
C  READ HEADER INFORMATION
      READ(40,901)IYEAR
C  CHECK THAT DATASET FOR CORRECT YEAR IS READ
      ITEST=1962+K
      IF(IYEAR.NE.ITEST) WRITE(95,913)IYEAR,K
C  READ HEADING OF MONTHS
      READ(40,900)BLANK
      READ(40,900)BLANK
C  READ YEARLY TOTALS IN ACFT FOR EACH PUMPING WELL
C  CALCULATE CFS AND FLUX IN EACH MODEL LAYER
      DO 100 I=1,NWELL
C READ STATEMENT FOR TRANSIENT CALIBRATION
      READ(40,902) Q
      WRITE(6,902) Q
C READ STATEMENT FOR MANAGEMENT SCENARIOS
C      READ(41,917) Q
C      WRITE(6,902) Q
C  CONVERT AC-FT/YR TO CFS/YR
C  SKIP FOR MANAGEMENT SCENARIO
      Q=Q*0.00138127
C  SAVE Q IN BIG ARRAYS
      QUP(I,K)=PERCNT(I,1)/100.*Q
      QLO(I,K)=PERCNT(I,2)/100.*Q
      QALL(I,K)=Q
C  DETERMINE MAXIMUM PUMPAGE OF RECORD
      IF(ITEST.GE.1980)QMAX(I)=AMAX1(QMAX(I),Q)
      IF(ITEST.EQ.1984)QMIN(I)=AMIN1(QMIN(I),Q)
      IF(ITEST.EQ.1985)QMIN(I)=AMIN1(QMIN(I),Q)
C  ASSIGN PUMPAGE TO MODEL ROW/COL
      FLUX(IR(I),IC(I),1)=FLUX(IR(I),IC(I),1)+PERCNT(I,1)/100.*Q
      FLUX(IR(I),IC(I),2)=FLUX(IR(I),IC(I),2)+PERCNT(I,2)/100.*Q
C  COUNT NUMBER OF YEARS THAT WELL HAD PUMPAGE, AND CALCULATE TOTAL
C  PUMPAGE FOR THAT WELL FOR YEARS WITH PUMPAGE VALUE (NOTE SOME
C  YEARS A FLOWING WELL WOULD BE CAPABLE OF FLOWING, BUT WATER LEVELS
C  WOULD BE TOO LOW).
C      IF(Q.GT.0.) THEN
      IF((ITEST.GE.1970).AND.(ITEST.LE.1984)) THEN
        NYRAVE(I)=NYRAVE(I)+1
        AVETOT(I)=AVETOT(I)+Q
        AVEPMP(I,1)=AVEPMP(I,1)+PERCNT(I,1)/100.*Q
        AVEPMP(I,2)=AVEPMP(I,2)+PERCNT(I,2)/100.*Q
      ENDIF
C      ENDIF
  100 CONTINUE
C  CALCULATE TOTAL PUMPAGE FOR EACH LAYER
      CALL SUM3D(FLUX,NROW,NCOL,NLAY,TOTPUM)
C  WRITE OUT HEADER CARD AND TOTAL PUMPAGE
      WRITE(30,904) IYEAR,TOTPUM(1)
      WRITE(31,904) IYEAR,TOTPUM(2)
      TPUMP=TOTPUM(1)+TOTPUM(2)
      WRITE(90,906) IYEAR,TOTPUM(1),TOTPUM(2),TPUMP
      IF(K.EQ.1)WRITE(32,904) IYEAR,TOTPUM(1)
      IF(K.EQ.1)WRITE(33,904) IYEAR,TOTPUM(2)
      IF(K.EQ.22)WRITE(34,904) IYEAR,TOTPUM(1)
      IF(K.EQ.22)WRITE(35,904) IYEAR,TOTPUM(2)
      IF(K.EQ.23)WRITE(36,904) IYEAR,TOTPUM(1)
      IF(K.EQ.23)WRITE(37,904) IYEAR,TOTPUM(2)
C  WRITE OUT ARRAYS FOR EACH LAYER
      DO 200 I=1,NROW
      WRITE(30,903) (FLUX(I,J,1),J=1,NCOL)
      WRITE(31,903) (FLUX(I,J,2),J=1,NCOL)
      IF(K.EQ.1)WRITE(32,903) (FLUX(I,J,1),J=1,NCOL)
      IF(K.EQ.1)WRITE(33,903) (FLUX(I,J,2),J=1,NCOL)
      IF(K.EQ.22)WRITE(34,903) (FLUX(I,J,1),J=1,NCOL)
      IF(K.EQ.22)WRITE(35,903) (FLUX(I,J,2),J=1,NCOL)
      IF(K.EQ.23)WRITE(36,903) (FLUX(I,J,1),J=1,NCOL)
      IF(K.EQ.23)WRITE(37,903) (FLUX(I,J,2),J=1,NCOL)
  200 CONTINUE
C  READ MONTHLY TOTALS
      READ(40,900)BLANK
  300 CONTINUE
      WRITE(90,909)
C  CALCULATE AND OUTPUT ANNUAL PUMPAGE AVERAGES FOR EACH WELL
      DO 400 I=1,NWELL
      IF(NYRAVE(I).NE.0) THEN
        AVEPMP(I,1)=AVEPMP(I,1)/NYRAVE(I)
        AVEPMP(I,2)=AVEPMP(I,2)/NYRAVE(I)
        AVETOT(I)=AVETOT(I)/NYRAVE(I)
      ENDIF
      WRITE(90,910)AWELL(I),AVEPMP(I,1),AVEPMP(I,2),AVETOT(I),QMAX(I),
     - QMIN(I)
C  NO PUMPAGE RECORDED WRITE ERROR MESSAGE
      IF(NYRAVE(I).EQ.0)WRITE(90,916)
  400 CONTINUE
C  OUTPUT ALL ANNUAL PUMPAGE AND AVERAGES FOR EACH WELL, IN CFS
      WRITE(90,909)
      DO 500 I=1,NWELL
      WRITE(90,910)AWELL(I),AVEPMP(I,1),AVEPMP(I,2),AVETOT(I),QMAX(I),
     - QMIN(I)
      WRITE(90,915)(QALL(I,J),J=1,NYEARS)
  500 CONTINUE
C  CALCULATE WELL FIELD VALUES
C  ANNUAL PUMPAGE FOR WELL FIELDS, BY WATER YEAR
      DO 550 K=1,NYEARS
      DO 550 I=1,NWELL
C  CONVERT Q BACK TO AC-FT
      QU=QUP(I,K)/0.00138127
      QL=QLO(I,K)/0.00138127
      QT=QALL(I,K)/0.00138127
      DO 550 M=1,NWF
      IF(WFID(I).EQ.WFNAME(M))WFUP(M,K)=WFUP(M,K)+QU
      IF(WFID(I).EQ.WFNAME(M))WFLO(M,K)=WFLO(M,K)+QL
      IF(WFID(I).EQ.WFNAME(M))WFTOT(M,K)=WFTOT(M,K)+QT
  550 CONTINUE
C  CALCULATE AVERAGE FOR EACH WELL FIELD
      DO 590 M=1,NWF
      DO 580 K=1,NYEARS
      WFAVE(M)=WFAVE(M)+WFTOT(M,K)
  580 CONTINUE
      WFAVE(M)=WFAVE(M)/NYEARS
  590 CONTINUE
C  WRITE OUT WELL FIELD PUMPAGE
      WRITE(90,920)(WFNAME(M),M=1,NWF)
      WRITE(90,922)(WFAVE(M),M=1,NWF)
      DO 600 K=1,NYEARS
      WRITE(90,921) WYNUM(K),(WFTOT(M,K),M=1,NWF)
  600 CONTINUE
C  WRITE OUT WELL FIELD PUMPAGE IN ANOTHER FORMAT
      DO 700 M=1,NWF
      WRITE(90,918)WFNAME(M),WFAVE(M)
      WRITE(90,919)(WYNUM(K),WFUP(M,K),WFLO(M,K),WFTOT(M,K),K=1,NYEARS)
  700 CONTINUE
C  OUTPUT ANNUAL PUMPAGE FOR EACH WELL, IN ACFT
      DO 800 I=1,NWELL
C  CONVERT QALL TO ACFT
      DO 810 J=1,NYEARS
      QALL(I,J)=QALL(I,J)/0.00138127
  810 CONTINUE
      WRITE(90,910)AWELL(I)
      WRITE(90,923)(QALL(I,J),J=1,NYEARS)
  800 CONTINUE
C  CALCULATE AVERAGE PUMPAGE PER WELL FOR SPECIFIC TIME PERIOD
      IOUT=90
      TITLE='AVE 1984-88'
      CALL AVEQ(TITLE,AWELL,QALL,NWELL,NYEARS,IOUT)
C  FORMAT STATEMENTS
  900 FORMAT(A80,A80)
  901 FORMAT(I10)
  902 FORMAT(80X,F10.0)
  903 FORMAT(8F10.3)
  904 FORMAT(I10,'WY  TOTAL PUMPAGE CFS = ',F10.2)
  905 FORMAT(40I3)
  906 FORMAT(I10,'WY  TOTAL PUMPAGE UPPER = ',F10.2,
     - ' TOTAL PUMPAGE LOWER = ',F10.2,
     - ' TOTAL PUMPAGE BOTH, CFS = ',F10.2)
  908 FORMAT(A5,5X,2I5,2F10.0,7X,A3,2F10.0,10X,2F10.0,10X,8X,A2)
  909 FORMAT('WELL NO       AVE UP    AVE LO   AVE TOT   MAX WY',
     - ' MIN 84-85')
  910 FORMAT(A5,5X,6F10.2)
  911 FORMAT(A80)
  912 FORMAT(4I10)
  913 FORMAT('ERROR --INCORRECT YEAR IS READ IN ',2I10)
  914 FORMAT('ERROR -- NODE IS NOT ACTIVE, I,J,K = ',3I5)
  915 FORMAT(16X,8F8.2,/,5(10F8.2,/))
  916 FORMAT('POSSIBLE ERROR = NUMBER OF YEARS WITH PUMPAGE = 0')
  917 FORMAT(60X,F10.2)
  918 FORMAT('WELLFIELD = ',A2,' AVERAGE 1963-88 PUMPAGE = ',F10.0,/,
     - 'WATER YEAR  AVE UP    AVE LO  AVE TOTAL')
  919 FORMAT(F10.2,3F10.0)
  920 FORMAT(10X,10('     ',A2,'   '))
  921 FORMAT(F10.2,10F10.0)
  922 FORMAT('  AVERAGE ',10F10.0)
  923 FORMAT(16X,8F8.0,/,5(10F8.0,/))
      STOP
      END
C*****************************************************************************
      SUBROUTINE SUM3D(BUF,NR,NC,NL,SUM)
C*****************************************************************************
      DIMENSION BUF(180,40,2),SUM(2)
C      COMMON /SUM3/BUF
C--SUM A 3-DIMENSIONAL ARRAY FOR EACH LAYER
      SUM(1)=0.
      SUM(2)=0.
      DO 100 I=1,NR
      DO 100 J=1,NC
        SUM(1)=SUM(1)+BUF(I,J,1)
        SUM(2)=SUM(2)+BUF(I,J,2)
  100 CONTINUE
      RETURN
      END
C*****************************************************************************
      SUBROUTINE ZERO3D(BUFF,NROW,NCOL,NLAY)
C*****************************************************************************
      DIMENSION BUFF(180,40,2)
C      COMMON /ZERO3/ BUFF
C--ZERO 3-DIMENISONAL BUFFER
      DO 100 K=1,NLAY
      DO 100 I=1,NROW
      DO 100 J=1,NCOL
        BUFF(I,J,K)=0.
  100 CONTINUE
      RETURN
      END
C*****************************************************************************
      SUBROUTINE AVEQ(TITLE,NAME,Q,NW,NY,IOUT)
C*****************************************************************************
      DIMENSION Q(200,30),AVE(200),NAME(200),QMX(200)
      CHARACTER*5 NAME
      CHARACTER*10 TITLE
C--ZERO ARRAY
      DO 10 I=1,NW
      AVE(I)=0.
      QMX(I)=0.
   10 CONTINUE
C--AVERAGE FOR A GIVEN PERIOD
      DO 100 I=1,NW
      DO 100 J=1,NY
        IYEAR=1962+J
        IF(IYEAR.GE.1985.AND.IYEAR.LE.1988)THEN
C        AVE(I)=AVE(I)+Q(I,J)
          IF(IYEAR.EQ.1985) AVE(I)=AVE(I)+Q(I,J)
          IF(IYEAR.EQ.1986) AVE(I)=AVE(I)+Q(I,J)
          IF(IYEAR.EQ.1988) AVE(I)=AVE(I)+Q(I,J)
          QMX(I)=AMAX1(Q(I,J),QMX(I))
        ENDIF
  100 CONTINUE
C--CALCULATE AVERAGE AND CONVERT TO AC-FT, IF NECESSARY
C--WRITE OUT VALUES
      WRITE(IOUT,902)TITLE
      DO 200 I=1,NW
C      AVE(I)=AVE(I)/3.*365.*86400./43560.
      AVE(I)=AVE(I)/3.
      WRITE(IOUT,900)NAME(I),AVE(I),QMX(I)
C      WRITE(IOUT,901)NAME(I),AVE(I),QMX(I)
  200 CONTINUE
C--FORMAT STATEMENTS
  900 FORMAT(A5,5X,2F10.0)
  901 FORMAT(A5,5X,2F10.2)
  902 FORMAT('AVERAGE VALUES FOR ',A10)
      RETURN
      END
