C*******************************************************************************
C  POST.MAPLP.F77 -- A PROGRAM TO PROCESS OUTPUT FROM THE MODULAR GROUND
C  WATER MODEL BY MCDONALD AND HARBAUGH (1984).
C               BY   WES DANSKIN
C                    LAST REVISION 1/25/89
C
C  THE MODULAR MODEL WRITES AN UNFORMATTED FILE CONSISTING OF A HEADER RECORD
C  AND THEN CELL-BY-CELL VALUES FROM THE SOLUTION OUTPUT.  IN THE CASE OF HEADS
C  AND DRAWDOWNS THE OUTPUT IS A HEADER FOR EACH LAYER OF BUF(NCOL,NROW).  IN
C  THE CASE OF FLUXES A HEADER IS FOLLOWED BY BUFF(NCOL,NROW,NLAY).
C
C  VALUES IN THE UNFORMATTED RECORD BUFF(I,J,K) ARE INCREMENTED I,J,K.
C  PRINTING OF MAPS IS CONTROLLED BY MAPFLG (1=PRINT)
C*****************************************************************************
C  MAIN PROGRAM
C*****************************************************************************
      DIMENSION BUFF(40,180,2),BUF(40,180),IBND(40,180,2),TEXT(4),
     - SUM(2),MAPFLG(2),XLSD(40,180),HEAD1(40,180),HEAD2(40,180),
     - ACCUM(40,180,2),ASUM(2),T1984(40,180,2),ETMOD(40,180,2),
     - HEAD63U(40,180),HEAD63L(40,180),HEADSSU(40,180),HEADSSL(40,180),
     - HEAD84U(40,180),HEAD84L(40,180)
      CHARACTER*25 NAME
      CHARACTER*80 HEADNG
      COMMON /POST1/ BUFF,BUF,IBND,XLSD
      COMMON /POST2/ HEAD1,HEAD2,ACCUM,T1984,ETMOD,HEAD63U,HEAD63L
      COMMON /POST3/ HEADSSU,HEADSSL,HEAD84U,HEAD84L
C  DEFINE DIMENSIONS
C      NC=40
C      NR=180
C      NL=2
C  READ CONTROL PARAMETERS
      ISUM=1
      IOUT=6
      READ(30,902)HEADNG
      WRITE(IOUT,902)HEADNG
      READ(30,902)HEADNG
      WRITE(IOUT,902)HEADNG
      READ(30,903)NL,NR,NC,NYEARS
      WRITE(IOUT,903)NL,NR,NC,NYEARS
C  ZERO MATRICES AND VARIABLES
      CALL ZERO2D (BUF,NC,NR)
      CALL ZERO2D (XLSD,NC,NR)
      CALL ZERO2D (HEAD1,NC,NR)
      CALL ZERO2D (HEAD2,NC,NR)
      CALL ZERO3D (T1984,NC,NR,NL)
      CALL ZERO3D (ETMOD,NC,NR,NL)
      CALL ZERO3D (ACCUM,NC,NR,NL)
      CALL ZERO3D (BUFF,NC,NR,NL)
      DO 10 K=1,NL
      MAPFLG(K)=0
      DO 10 J=1,NR
      DO 10 I=1,NC
      IBND(I,J,K)=0
   10 CONTINUE
      SUM(1)=0.
      SUM(2)=0.
      ASUM(1)=0.
      ASUM(2)=0.
      TEXT(1)='    '
      TEXT(2)='    '
      TEXT(3)='    '
      TEXT(4)='    '
      KSTP=0
      KPER=0
      PERTIM=0.
      TOTIM=0.
      FLAG=0.
      IREAD=0
      IOUT=0
C--READ IN UPPER IBOUND ARRAY USED TO DEFINE MAP BOUNDARIES
      DO 100 J=1,NR
      READ(51,908) (IBND(I,J,1),I=1,NC)
  100 CONTINUE
C  READ IN TRANSPIRATION MATRIX FROM RAWSON
C  ONLY LAYER 1 IS NEEDED, BUT MODEL ET OUTPUT IS 3D MATRIX
      DO 110 J=1,NR
      READ(59,905) (T1984(I,J,1),I=1,NC)
  110 CONTINUE
C  CONVERT FROM IN/YR TO CFS/YR, AND CHANGE TO NEGATIVE
      DO 115 I=1,NC
      DO 115 J=1,NR
      DO 115 K=1,NL
      T1984(I,J,K)=T1984(I,J,K)/12.*2000.*2000./365./86400.*(-1.0)
  115 CONTINUE
C  READ AND CREATE LINE PRINTER MAPS FOR ALL MODEL VALUES
C  FOR SELECTED TIMESTEPS
      DO 800 KYEAR=1,NYEARS
C  READ IF TIMESTEP IS SELECTED FOR PRINTING
      READ(5,901)IMAP,IYEAR
      ITEST=1962+KYEAR
      IF(IYEAR.NE.ITEST) WRITE(97,913)IYEAR,KYEAR
      IOUT=31
      IF(IMAP.EQ.1) WRITE(IOUT,904)IYEAR
  200 CONTINUE
CC--LAND SURFACE DATA
C      IF(KYEAR.GT.1) GOTO 205
C      IREAD=60
C      IOUT=31
C      FLAG=1.
C      NAME='LAND SURFACE DATUM       '
C      CALL ZERO2D (BUF,NC,NR)
C      READ(IREAD,*)BUF
C      CALL SET2D (XLSD,BUF,NC,NR)
C      MAPFLG(1)=1
C      MAPFLG(2)=0
C      IF(IMAP.EQ.1) CALL MAP2D (BUF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR)
C  205 CONTINUE
CC--HEADS, LAYER 1 (UPPER)
C      IREAD=49
C      IOUT=31
C      FLAG=1.
C      NAME='HEADS -- UPPER LAYER     '
C      READ(IREAD)KSTP,KPER,PERTIM,TOTIM,TEXT,NCOL,NROW,ILAY
C      CALL ZERO2D (BUF,NC,NR)
C      READ(IREAD)BUF
C      CALL SET2D (HEAD1,BUF,NC,NR)
C      IF(ITEST.EQ.1963)CALL SET2D (HEAD63U,HEAD1,NC,NR)
C      MAPFLG(1)=1
C      MAPFLG(2)=0
C      IF(IMAP.EQ.1) CALL MAP2D (BUF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR)
C  210 CONTINUE
CC--HEADS, LAYER 2 (LOWER)
C      IREAD=49
C      IOUT=31
C      FLAG=1.
C      READ(IREAD)KSTP,KPER,PERTIM,TOTIM,TEXT,NCOL,NROW,ILAY
C      NAME='HEADS -- LOWER LAYER     '
C      CALL ZERO2D (BUF,NC,NR)
C      READ(IREAD)BUF
C      CALL SET2D (HEAD2,BUF,NC,NR)
C      IF(ITEST.EQ.1963)CALL SET2D (HEAD63L,HEAD2,NC,NR)
C      MAPFLG(1)=1
C      MAPFLG(2)=0
C      IF(IMAP.EQ.1) CALL MAP2D (BUF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR)
C  215 CONTINUE
CC--DRAWDOWN, LAYER 1
C      IREAD=50
C      IOUT=31
C      FLAG=9999999.99
C      READ(IREAD)KSTP,KPER,PERTIM,TOTIM,TEXT,NCOL,NROW,ILAY
C      NAME='DRAWDOWN -- UPPER        '
C      CALL ZERO2D (BUF,NC,NR)
C      READ(IREAD)BUF
C      MAPFLG(1)=1
C      MAPFLG(2)=0
C      IF(IMAP.EQ.1) CALL MAP2D (BUF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR)
C      IOUT=34
C      CALL SURF2A (BUF,IBND,NC,NR,IOUT,NAME)
C  220 CONTINUE
CC--DRAWDOWN, LAYER 2
C      READ(IREAD)KSTP,KPER,PERTIM,TOTIM,TEXT,NCOL,NROW,ILAY
C      NAME='DRAWDOWN -- LOWER        '
C      CALL ZERO2D (BUF,NC,NR)
C      READ(IREAD)BUF
C      MAPFLG(1)=1
C      MAPFLG(2)=0
C      IF(IMAP.EQ.1) CALL MAP2D (BUF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR)
C  225 CONTINUE
CC--LSD MINUS HEAD IN LAYER 1
C      IOUT=31
C      FLAG=11111.
C      NAME='LSD MINUS HEAD IN LAY1   '
C      MAPFLG(1)=1
C      MAPFLG(2)=0
C      CALL ZERO2D (BUF,NC,NR)
C      CALL DIFF2D (XLSD,HEAD1,BUF,NC,NR)
C      IF(IMAP.EQ.1) CALL MAP2D (BUF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR)
C      CONTINUE
CC--HEAD IN LAYER 2 BELOW LAYER 1
C      IOUT=31
C      FLAG=11111.
C      NAME='HEAD IN LAY1 MINUS LAY 2 '
C      MAPFLG(1)=1
C      MAPFLG(2)=0
C      CALL ZERO2D (BUF,NC,NR)
C      CALL DIFF2D (HEAD1,HEAD2,BUF,NC,NR)
C      IF(IMAP.EQ.1) CALL MAP2D (BUF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR)
C      CONTINUE
CC--DRAWDOWN PREDICTIONS (1963-1966,-1969,-1972)
C      IF(IMAP.EQ.1) THEN
C        IOUT=31
C        FLAG=11111.
C        NAME='DD PRED, 1963-84 LAYER 1 '
C        MAPFLG(1)=1
C        MAPFLG(2)=0
C        CALL ZERO2D (BUF,NC,NR)
C        CALL DIFF2D (HEAD1,HEAD63U,BUF,NC,NR)
C        CALL MAP2D (BUF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR)
C        IOUT=34
C        CALL SURF2 (BUF,IBND,NC,NR,IOUT,NAME)
C        CONTINUE
C        IOUT=31
C        FLAG=11111.
C        NAME='DD PRED, 1963-84 LAYER 2 '
C        MAPFLG(1)=1
C        MAPFLG(2)=0
C        CALL ZERO2D (BUF,NC,NR)
C        CALL DIFF2D (HEAD2,HEAD63L,BUF,NC,NR)
C        CALL MAP2D (BUF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR)
C        IOUT=34
C        CALL SURF2 (BUF,IBND,NC,NR,IOUT,NAME)
C      ENDIF
C      CONTINUE
CC--HEAD DIFFERENCE, RECENT STEADY STATE VS 1984 TRANSIENT CALIBRATION
CC--USED IF RUNNING SCENARIOS AND LOOKING BACK TO CALIBRATION RUN
CC--READ IN CALIBRATED HEADS, THEN MODEL RUN
C      IF(ITEST.EQ.1963)THEN
CC READ 1984 CALIBRATION DATASET AND SET ARRAYS
C        DO 250 IT=1,22
C        IREAD=52
C        READ(IREAD)KSTP,KPER,PERTIM,TOTIM,TEXT,NCOL,NROW,ILAY
C        CALL ZERO2D (BUF,NC,NR)
C        READ(IREAD)BUF
C        CALL SET2D (HEAD84U,BUF,NC,NR)
C        READ(IREAD)KSTP,KPER,PERTIM,TOTIM,TEXT,NCOL,NROW,ILAY
C        CALL ZERO2D (BUF,NC,NR)
C        READ(IREAD)BUF
C        CALL SET2D (HEAD84L,BUF,NC,NR)
C  250   CONTINUE
CC MAKE MAPS
C        IOUT=31
C        FLAG=11111.
C        NAME='HD DIF,1984 T VS S, LAY 1'
C        MAPFLG(1)=1
C        MAPFLG(2)=0
C        CALL ZERO2D (BUF,NC,NR)
C        CALL DIFF2D (HEAD84U,HEAD1,BUF,NC,NR)
C        CALL MAP2D (BUF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR)
C        IOUT=34
C        CALL SURF2 (BUF,IBND,NC,NR,IOUT,NAME)
C        CONTINUE
C        IOUT=31
C        FLAG=11111.
C        NAME='HD DIF,1984 T VS S, LAY 2'
C        MAPFLG(1)=1
C        MAPFLG(2)=0
C        CALL ZERO2D (BUF,NC,NR)
C        CALL DIFF2D (HEAD84L,HEAD2,BUF,NC,NR)
C        CALL MAP2D (BUF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR)
C        IOUT=34
C        CALL SURF2 (BUF,IBND,NC,NR,IOUT,NAME)
C      ENDIF
C      CONTINUE
CC--HEAD DIFFERENCE, 1984 TRANSIENT VS RECENT STEADY STATE
CC--USED IF RUNNING CALIBRATION RUN AFTER SCENARIOS
CC--READ IN PREDICTED HEADS, THEN MODEL RUN
C      IF(ITEST.EQ.1984)THEN
CC READ PREDICTED DATASET AND SET ARRAYS
C       IREAD=53
C        READ(IREAD)KSTP,KPER,PERTIM,TOTIM,TEXT,NCOL,NROW,ILAY
C        CALL ZERO2D (BUF,NC,NR)
C        READ(IREAD)BUF
C        CALL SET2D (HEADSSU,BUF,NC,NR)
C        READ(IREAD)KSTP,KPER,PERTIM,TOTIM,TEXT,NCOL,NROW,ILAY
C        CALL ZERO2D (BUF,NC,NR)
C        READ(IREAD)BUF
C        CALL SET2D (HEADSSL,BUF,NC,NR)
CC MAKE MAPS
C        IOUT=31
C        FLAG=11111.
C        NAME='HD DIF,1984 T VS S, LAY 1'
C        MAPFLG(1)=1
C        MAPFLG(2)=0
C        CALL ZERO2D (BUF,NC,NR)
C        CALL DIFF2D (HEAD1,HEADSSU,BUF,NC,NR)
C        CALL MAP2D (BUF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR)
C        IOUT=34
C        CALL SURF2 (BUF,IBND,NC,NR,IOUT,NAME)
C        CONTINUE
C        IOUT=31
C        FLAG=11111.
C        NAME='HD DIF,1984 T VS S, LAY 2'
C        MAPFLG(1)=1
C        MAPFLG(2)=0
C        CALL ZERO2D (BUF,NC,NR)
C        CALL DIFF2D (HEAD2,HEADSSL,BUF,NC,NR)
C        CALL MAP2D (BUF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR)
C        IOUT=34
C        CALL SURF2 (BUF,IBND,NC,NR,IOUT,NAME)
C      ENDIF
C      CONTINUE
C--STORAGE
      IREAD=55
      IOUT=32
      FLAG=111.11
      NAME='STORAGE                  '
      READ(IREAD)KSTP,KPER,TEXT,NCOL,NROW,NLAY
      CALL ZERO3D (BUFF,NC,NR,NL)
      READ(IREAD)BUFF
      CALL SUM3D (BUFF,NC,NR,NL,SUM)
C  ACCUMULATED TOTALS FOR EACH LAYER
      ASUM(1)=ASUM(1)+SUM(1)
      ASUM(2)=ASUM(2)+SUM(2)
C  REZERO COUNTERS AND ACCUMULATORS
C      IF(IYEAR.EQ.1970) CALL ZERO3D (ACCUM,NC,NR,NL)
C      IF(IYEAR.EQ.1970) ASUM(1)=SUM(1)
C      IF(IYEAR.EQ.1970) ASUM(2)=SUM(2)
      CALL ACCUMM (BUFF,NC,NR,NL,ACCUM)
      MAPFLG(1)=1
      MAPFLG(2)=1
      IF(IMAP.EQ.1) CALL MAP3D (BUFF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
     -  IOUT,MAPFLG,IYEAR,ISUM)
      ISUM=1
      NAME='STORAGE -- ACCUMULATED   '
      IF(IMAP.EQ.1) CALL MAP3D (ACCUM,IBND,NC,NR,NL,NAME,TEXT,FLAG,ASUM,
     -  IOUT,MAPFLG,IYEAR,ISUM)
      ISUM=1
      CONTINUE
C--CONSTANT HEAD FLUX
      IREAD=55
      IOUT=32
      FLAG=111.11
      NAME='CONSTANT HEAD FLUX       '
      READ(IREAD)KSTP,KPER,TEXT,NCOL,NROW,NLAY
      CALL ZERO3D (BUFF,NC,NR,NL)
      READ(IREAD)BUFF
      CALL SUM3D (BUFF,NC,NR,NL,SUM)
      MAPFLG(1)=1
      MAPFLG(2)=1
      IF(IMAP.EQ.1) CALL MAP3D (BUFF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
     -  IOUT,MAPFLG,IYEAR,ISUM)
      CONTINUE
C--CELL FLOW RATES
      IREAD=55
      IOUT=32
      FLAG=111.11
      NAME='CELL FLOW -- RIGHT FACE  '
      READ(IREAD)KSTP,KPER,TEXT,NCOL,NROW,NLAY
      CALL ZERO3D (BUFF,NC,NR,NL)
      READ(IREAD)BUFF
      WRITE(93,999)ITEST
  999 FORMAT(I5)
      WRITE(91,998)BUFF(22,129,1),BUFF(23,129,1),
     - BUFF(22,129,2),BUFF(23,129,2)
  998 FORMAT(4F10.4)
      CALL SUM3D (BUFF,NC,NR,NL,SUM)
      MAPFLG(1)=1
      MAPFLG(2)=1
      IF(IMAP.EQ.1) CALL MAP3D (BUFF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
     -  IOUT,MAPFLG,IYEAR,ISUM)
      CONTINUE
      NAME='CELL FLOW -- FRONT FACE  '
      READ(IREAD)KSTP,KPER,TEXT,NCOL,NROW,NLAY
      CALL ZERO3D (BUFF,NC,NR,NL)
      READ(IREAD)BUFF
      WRITE(92,998)BUFF(23,128,1),BUFF(23,129,1),
     - BUFF(23,128,2),BUFF(23,129,2)
      CALL SUM3D (BUFF,NC,NR,NL,SUM)
      MAPFLG(1)=1
      MAPFLG(2)=1
      IF(IMAP.EQ.1) CALL MAP3D (BUFF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
     -  IOUT,MAPFLG,IYEAR,ISUM)
      CONTINUE
      NAME='CELL FLOW -- BOTTOM FACE '
      READ(IREAD)KSTP,KPER,TEXT,NCOL,NROW,NLAY
      CALL ZERO3D (BUFF,NC,NR,NL)
      READ(IREAD)BUFF
      WRITE(93,998)BUFF(23,129,1),BUFF(23,129,2)
      CALL SUM3D (BUFF,NC,NR,NL,SUM)
      MAPFLG(1)=1
      MAPFLG(2)=1
      IF(IMAP.EQ.1) CALL MAP3D (BUFF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
     -  IOUT,MAPFLG,IYEAR,ISUM)
      CONTINUE
CC--WELLS AND STREAM RECHARGE
C      IREAD=63
C      IOUT=33
C      FLAG=111.11
C      NAME='WELLS AND STREAMS        '
C      READ(IREAD)KSTP,KPER,TEXT,NCOL,NROW,NLAY
C      CALL ZERO3D (BUFF,NC,NR,NL)
C      READ(IREAD)BUFF
C      CALL SUM3D (BUFF,NC,NR,NL,SUM)
C      MAPFLG(1)=1
C      MAPFLG(2)=0
C      IF(IMAP.EQ.1) CALL MAP3D (BUFF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR,ISUM)
C      CONTINUE
CC--DRAINS
C      IREAD=66
C      IOUT=33
C      FLAG=111.11
C      NAME='DRAINS -- SPRINGS        '
C      READ(IREAD)KSTP,KPER,TEXT,NCOL,NROW,NLAY
C      CALL ZERO3D (BUFF,NC,NR,NL)
C      READ(IREAD)BUFF
C      CALL SUM3D (BUFF,NC,NR,NL,SUM)
C      MAPFLG(1)=1
C      MAPFLG(2)=0
C      IF(IMAP.EQ.1) CALL MAP3D (BUFF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR,ISUM)
C      CONTINUE
CC--PRECIP AND ET RECHARGE (NET)
C      IREAD=64
C      IOUT=33
C      FLAG=111.11
C      NAME='PRECIP AND ET ON FANS    '
C      READ(IREAD)KSTP,KPER,TEXT,NCOL,NROW,NLAY
C      CALL ZERO3D (BUFF,NC,NR,NL)
C      READ(IREAD)BUFF
C      CALL SUM3D (BUFF,NC,NR,NL,SUM)
C      MAPFLG(1)=1
C      MAPFLG(2)=0
C      IF(IMAP.EQ.1) CALL MAP3D (BUFF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR,ISUM)
C      CONTINUE
C--EVAPOTRANSPIRATION
      IREAD=62
      IOUT=33
      FLAG=111.11
      NAME='EVAPOTRANSPIRATION       '
      READ(IREAD)KSTP,KPER,TEXT,NCOL,NROW,NLAY
      CALL ZERO3D (BUFF,NC,NR,NL)
      READ(IREAD)BUFF
      WRITE(94,998)BUFF(23,129,1),BUFF(23,129,2)
      CALL SET3D (ETMOD,BUFF,NC,NR,NL)
      CALL SUM3D (BUFF,NC,NR,NL,SUM)
      MAPFLG(1)=1
      MAPFLG(2)=0
      IF(IMAP.EQ.1) CALL MAP3D (BUFF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
     -  IOUT,MAPFLG,IYEAR,ISUM)
      CONTINUE
CC--DIFFERENCE BETWEEN GROENEVELD ET AND MODEL ET
C      IOUT=33
C      FLAG=111.11
C      NAME='MODEL ET MINUS RAWSON T  '
C      MAPFLG(1)=1
C      MAPFLG(2)=0
C      CALL DIFF3D (ETMOD,T1984,BUFF,NC,NR,NL)
C      CALL SUM3D (BUFF,NC,NR,NL,SUM)
C      IF(IMAP.EQ.1) CALL MAP3D (BUFF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
C     -  IOUT,MAPFLG,IYEAR,ISUM)
C      CONTINUE
C--RIVER
      IREAD=65
      IOUT=33
      FLAG=111.11
      NAME='OWENS RIVER              '
      READ(IREAD)KSTP,KPER,TEXT,NCOL,NROW,NLAY
      CALL ZERO3D (BUFF,NC,NR,NL)
      READ(IREAD)BUFF
      CALL SUM3D (BUFF,NC,NR,NL,SUM)
      MAPFLG(1)=1
      MAPFLG(2)=0
      IF(IMAP.EQ.1) CALL MAP3D (BUFF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,
     -  IOUT,MAPFLG,IYEAR,ISUM)
      CONTINUE
C--GENERAL HEAD BOUNDARY
      CONTINUE
C  END OF MAJOR DO-LOOP
  800 CONTINUE
C  END OF MAPPING
C  FORMAT STATEMENTS
  901 FORMAT(2I10)
  902 FORMAT(A80)
  903 FORMAT(4I10)
  904 FORMAT('******************** WATER YEAR = ',I10,
     -' ***********************************************')
  905 FORMAT(8F10.4)
  908 FORMAT(40I3)
  913 FORMAT('ERROR --INCORRECT YEAR IS READ IN ',2I10)
C
C  END OF MAIN PROGRAM
      STOP
      END
C*****************************************************************************
      SUBROUTINE MAP2D (BUF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,IOUT,
     -  MAPFLG,IYEAR)
C
C  CREATE LINE PRINTER MAP FOR LAND SURFACE, HEAD, AND DRAWDOWN.  MAP ONLY ONE
C  LAYER AT A TIME.
C
C*****************************************************************************
      DIMENSION BUF(40,180),IBND(40,180,2),NCPRN(40),TEXT(4),
     -  SUM(2),MAPFLG(2)
      CHARACTER*25 NAME
C      COMMON /POST1/ BUFF,BUF,IBND,XLSD
C      COMMON /POST2/ HEAD1,HEAD2,ACCUM
C--TEST IF MAP IS DESIRED
      IF(MAPFLG(1).NE.1) RETURN
C--CREATE PRINT LABELS
      DO 20 I=1,NC
        NCPRN(I)=I
   20 CONTINUE
C--SET BOUNDARY FLUX TO FLAG VALUE
      DO 40 J=1,NR
      DO 40 I=1,NC
      IF(IBND(I,J,1).NE.0) GOTO 40
      BUF(I,J)=FLAG
   40 CONTINUE
C--WRITE HEADING
      WRITE(IOUT,902)NAME,IYEAR,FLAG,TEXT
  902 FORMAT('0',2X,'*** MAP OF ',A25,5X,'WATER YEAR =',I5,/,
     - T6,F10.2,' INDICATES NODE IS OUT OF THE BASIN',10X,
     - '+ FLUX IN OR TOO LOW HEADS    - FLUX OUT OR TOO HIGH HEADS',/
     - ,T10,'MODEL SAYS..',4A4)
C--WRITE FIRST HALF OF MAP
      WRITE(IOUT,900) (NCPRN(I),I=1,20)
      DO 100 J=1,NR
      WRITE(IOUT,901) J,(BUF(I,J),I=1,20)
  100 CONTINUE
C--WRITE SECOND HALF OF MAP
      WRITE(IOUT,900) (NCPRN(I),I=21,40)
      DO 110 J=1,NR
      WRITE(IOUT,901) J,(BUF(I,J),I=21,40)
  110 CONTINUE
C--FORMAT STATEMENTS
  900 FORMAT(1X,'ROW',1X,20(2X,I2,2X))
  901 FORMAT(1X,I3,1X,20F6.1)
      RETURN
      END
C*****************************************************************************
      SUBROUTINE MAP3D (BUFF,IBND,NC,NR,NL,NAME,TEXT,FLAG,SUM,IOUT,
     -  MAPFLG,IYEAR,ISUM)
C
C  CREATE LINE PRINTER MAP FOR FLUXES.  MAP ALL LAYERS
C
C*****************************************************************************
      DIMENSION BUFF(40,180,2),IBND(40,180,2),NCPRN(40),TEXT(4),
     -  SUM(2),MAPFLG(2)
      CHARACTER*25 NAME
C      COMMON /POST1/ BUFF,BUF,IBND,XLSD
C      COMMON /POST2/ HEAD1,HEAD2,ACCUM
C--CREATE PRINT LABELS
      DO 20 I=1,NC
        NCPRN(I)=I
   20 CONTINUE
C--SET BOUNDARY FLUX TO FLAG VALUE
      DO 40 K=1,NL
      DO 40 J=1,NR
      DO 40 I=1,NC
      IF(IBND(I,J,1).NE.0) GOTO 40
      BUFF(I,J,K)=FLAG
   40 CONTINUE
C--CREATE A LINE PRINTER MAP FOR EACH LAYER THAT IS REQUESTED
      DO 300 K=1,NL
      IF(MAPFLG(K).NE.1) GOTO 300
C--WRITE HEADING
      WRITE(IOUT,903)NAME,IYEAR,K,FLAG,SUM(K),TEXT
  903 FORMAT('0',2X,'*** MAP OF ',A25,5X,'WATER YEAR =',I5,
     - ' FOR LAYER ',I2,/,
     - T6,F10.2,' INDICATES NODE IS OUT OF THE BASIN',10X,
     - '+ FLUX IN OR TOO LOW HEADS    - FLUX OUT OR TOO HIGH HEADS',/
     - ,T10,'TOTAL FLUX = ',
     - F10.2,5X,'MODEL SAYS..',4A4)
C--WRITE FIRSTHALF OF MAP
      WRITE(IOUT,900) (NCPRN(I),I=1,20)
      DO 100 J=1,NR
      IF(ISUM.EQ.1)WRITE(IOUT,901) J,(BUFF(I,J,K),I=1,20)
      IF(ISUM.EQ.2)WRITE(IOUT,902) J,(BUFF(I,J,K),I=1,20)
  100 CONTINUE
C--WRITE SECOND HALF OF MAP
      WRITE(IOUT,900) (NCPRN(I),I=21,40)
      DO 110 J=1,NR
      IF(ISUM.EQ.1)WRITE(IOUT,901) J,(BUFF(I,J,K),I=21,40)
      IF(ISUM.EQ.2)WRITE(IOUT,902) J,(BUFF(I,J,K),I=21,40)
  110 CONTINUE
  300 CONTINUE
C--FORMAT STATEMENTS
  900 FORMAT(1X,'ROW',1X,20(2X,I2,2X))
  901 FORMAT(1X,I3,1X,20F6.2)
  902 FORMAT(1X,I3,1X,20F6.1)
      RETURN
      END
C*****************************************************************************
      SUBROUTINE ZERO2D(BUF,NC,NR)
C*****************************************************************************
      DIMENSION BUF(40,180)
C      COMMON /POST1/ BUFF,BUF,IBND,XLSD
C      COMMON /POST2/ HEAD1,HEAD2,ACCUM
C--ZERO 2-DIMENSIONAL BUFFER
      DO 100 J=1,NR
      DO 100 I=1,NC
        BUF(I,J)=0.
  100 CONTINUE
      RETURN
      END
C*****************************************************************************
      SUBROUTINE ZERO3D(BUFF,NC,NR,NL)
C*****************************************************************************
      DIMENSION BUFF(40,180,2)
C      COMMON /POST1/ BUFF,BUF,IBND,XLSD
C      COMMON /POST2/ HEAD1,HEAD2,ACCUM
C--ZERO 3-DIMENISONAL BUFFER
      DO 100 K=1,NL
      DO 100 J=1,NR
      DO 100 I=1,NC
        BUFF(I,J,K)=0.
  100 CONTINUE
      RETURN
      END
C*****************************************************************************
      SUBROUTINE SUM2D(BUF,NC,NR,SUM)
C*****************************************************************************
      DIMENSION BUF(40,180),SUM(2)
C      COMMON /POST1/ BUFF,BUF,IBND,XLSD
C      COMMON /POST2/ HEAD1,HEAD2,ACCUM
C--SUM A 2-DIMENSIONAL ARRAY FOR ONE LAYER
      SUM(1)=0.
      SUM(2)=0.
      DO 100 J=1,NR
      DO 100 I=1,NC
        SUM(1)=SUM(1)+BUF(I,J)
  100 CONTINUE
      RETURN
      END
C*****************************************************************************
      SUBROUTINE SUM3D(BUFF,NC,NR,NL,SUM)
C*****************************************************************************
      DIMENSION BUFF(40,180,2),SUM(2)
C      COMMON /POST1/ BUFF,BUF,IBND,XLSD
C      COMMON /POST2/ HEAD1,HEAD2,ACCUM
C--SUM A 3-DIMENSIONAL ARRAY FOR 2 LAYERS
      SUM(1)=0.
      SUM(2)=0.
      DO 100 K=1,NL
      DO 100 J=1,NR
      DO 100 I=1,NC
        SUM(K)=SUM(K)+BUFF(I,J,K)
  100 CONTINUE
      RETURN
      END
C*****************************************************************************
      SUBROUTINE DIFF2D (ARRAY1,ARRAY2,ARRAY3,NC,NR)
C
C  SUBTRACT TWO 2-DIMENSIONAL ARRAYS AND RETURN THE RESULTING ARRAY
C
C*****************************************************************************
      DIMENSION ARRAY1(40,180),ARRAY2(40,180),ARRAY3(40,180)
C      COMMON /POST1/ BUFF,BUF,IBND,XLSD
C      COMMON /POST2/ HEAD1,HEAD2,ACCUM
C--SUBTRACT ARRAYS
      DO 100 J=1,NR
      DO 100 I=1,NC
        ARRAY3(I,J)=ARRAY1(I,J)-ARRAY2(I,J)
  100 CONTINUE
      RETURN
      END
C*****************************************************************************
      SUBROUTINE DIFF3D (ARRAY1,ARRAY2,ARRAY3,NC,NR,NL)
C
C  SUBTRACT TWO 3-DIMENSIONAL ARRAYS AND RETURN THE RESULTING ARRAY
C
C*****************************************************************************
      DIMENSION ARRAY1(40,180,2),ARRAY2(40,180,2),ARRAY3(40,18,2)
C      COMMON /POST1/ BUFF,BUF,IBND,XLSD
C      COMMON /POST2/ HEAD1,HEAD2,ACCUM
C--SUBTRACT ARRAYS
      DO 100 J=1,NR
      DO 100 I=1,NC
      DO 100 K=1,NL
        ARRAY3(I,J,K)=ARRAY1(I,J,K)-ARRAY2(I,J,K)
  100 CONTINUE
      RETURN
      END
C*****************************************************************************
      SUBROUTINE SET2D (ARRAY1,ARRAY2,NC,NR)
C
C  SETS TWO 2-DIMENSIONAL ARRAYS EQUAL TO EACH OTHER
C
C*****************************************************************************
      DIMENSION ARRAY1(40,180),ARRAY2(40,180)
C      COMMON /POST1/ BUFF,BUF,IBND,XLSD
C      COMMON /POST2/ HEAD1,HEAD2,ACCUM
C--SET ARRAYS EQUAL TO EACH OTHER
      DO 100 J=1,NR
      DO 100 I=1,NC
        ARRAY1(I,J)=ARRAY2(I,J)
  100 CONTINUE
      RETURN
      END
C*****************************************************************************
      SUBROUTINE SET3D (ARRAY1,ARRAY2,NC,NR,NL)
C
C  SETS TWO 3-DIMENSIONAL ARRAYS EQUAL TO EACH OTHER
C
C*****************************************************************************
      DIMENSION ARRAY1(40,180,2),ARRAY2(40,180,2)
C      COMMON /POST1/ BUFF,BUF,IBND,XLSD
C      COMMON /POST2/ HEAD1,HEAD2,ACCUM
C--SET ARRAYS EQUAL TO EACH OTHER
      DO 100 J=1,NR
      DO 100 I=1,NC
      DO 100 K=1,NL
        ARRAY1(I,J,K)=ARRAY2(I,J,K)
  100 CONTINUE
      RETURN
      END
C*****************************************************************************
      SUBROUTINE ACCUMM(BUFF,NC,NR,NL,ACCUM)
C*****************************************************************************
      DIMENSION BUFF(40,180,2),ACCUM(40,180,2)
C      COMMON /POST1/ BUFF,BUF,IBND,XLSD
C      COMMON /POST2/ HEAD1,HEAD2,ACCUM
C--ADD RECENT VALUES OF CELL-BY-CELL FLUXES TO ACCUMULATED VALUE
      DO 100 K=1,NL
      DO 100 J=1,NR
      DO 100 I=1,NC
        ACCUM(I,J,K)=ACCUM(I,J,K)+BUFF(I,J,K)
  100 CONTINUE
      RETURN
      END
C*****************************************************************************
      SUBROUTINE SURF2 (ARRAY,IBND,NC,NR,IOUT,NAME)
C
C  WRITE OUT SURFACE 2 FILE
C
C*****************************************************************************
      DIMENSION ARRAY(40,180),IBND(40,180,2)
      CHARACTER*25 NAME
C--WRITE OUT ARRAYS
      WRITE(IOUT,902)NAME
      DO 100 J=1,NR
      DO 100 I=1,NC
      IF(IBND(I,J,1).EQ.0)THEN
        WRITE(IOUT,901)
      ELSE
        WRITE(IOUT,900)ARRAY(I,J)
      ENDIF
  100 CONTINUE
  900 FORMAT(F10.1)
  901 FORMAT('ERROR--INACTIVE NODE')
  902 FORMAT(A25)
      RETURN
      END
C*****************************************************************************
      SUBROUTINE SURF2A(VALUE,IBND,NCOL,NROW,IOUT,NAME)
C*****************************************************************************
      DIMENSION VALUE(40,180),IBND(40,180,2)
      CHARACTER*25 NAME
C  DEFINE SYBMOL CODES
      N1=-1
      I1=1
      I2=2
      I3=3
      I4=4
      I5=5
      I6=6
C  WRITE OUT SURFACE2 IDXY FILE
      WRITE(IOUT,901)NAME
      DO 100 J=1,NROW
      DO 100 I=1,NCOL
      IF(IBND(I,J,1).EQ.0)THEN
        WRITE(IOUT,903)
        GOTO 90
      ENDIF
      IF(VALUE(I,J).LT.0.1)VALUE(I,J)=0.
      VALUE(I,J)=VALUE(I,J)*(-1.0)
      WRITE(IOUT,900)VALUE(I,J),I1,I1
C      IF(VALUE(I,J).LT.15.)
C     -  WRITE(IOUT,904)
C      IF(VALUE(I,J).GE.12.AND.VALUE(I,J).LT.18)
C     -  WRITE(IOUT,900)VALUE(I,J),I1,I3
C      IF(VALUE(I,J).GE.18.AND.VALUE(I,J).LT.24)
C     -  WRITE(IOUT,900)VALUE(I,J),I1,I4
C      IF(VALUE(I,J).GE.24.AND.VALUE(I,J).LT.99)
C     -  WRITE(IOUT,900)VALUE(I,J),I1,I5
C      IF(VALUE(I,J).EQ.99)
C     -  WRITE(IOUT,903)
C      IF(VALUE(I,J).GT.99.)
C     -  WRITE(IOUT,904)
   90 CONTINUE
  100 CONTINUE
  900 FORMAT(F10.2,2I10)
  901 FORMAT('*** SURFACE 2 IDXY LIST OF FLOATING POINT VALUES ***'
     - ' FOR ',A25)
  902 FORMAT('ERROR--ZERO VALUE')
  903 FORMAT('ERROR--INACTIVE NODE')
  904 FORMAT('ERROR--VALUE OUTSIDE RANGE')
      RETURN
      END
