PROGRAM READ_CBC
C     ******************************************************************
C     Modified from Arlen’s ZONEBUDGET Program for MODFLOW-96 to read 
C     to read cell-by-cell budget files – Modified by D.P. McAda 3/04
C
C     The original program is documented in USGS Open-File Report 90-392,
C     written by Arlen W. Harbaugh
C     ******************************************************************
C        SPECIFICATIONS:
      PARAMETER (NODDIM=300000,NTRDIM=12,NZDIM=25,MXCOMP=25,MXZWCZ=10)
C-----   NODDIM must be greater than or equal to the product of the
C-----          number of layers, rows, and columns in the model grid.
C-----   NTRDIM must be greater than or equal to the number of budget
C-----          terms, other than flow between zones, that will appear
C-----          the budget.  In the original model, there is a maximum
C-----          of 8 terms -- constant-head, storage, wells, rivers,
C-----          drains, recharge, general-head boundaries, and
C-----          evapotranspiration.
C-----   NZDIM  must be greater than or equal to the maximum zone
C-----          number.
C-----   MXCOMP is the maximum number of composite zones, and must be
C-----          less than 81.
C-----   MXZWCZ is the maximum number of numeric zones within each
C-----          composite zone.
      COMMON /BUFCOM/BUFF
      COMMON /ZONCOM/IZONE
      COMMON /CHCOM/ICH
      DIMENSION BUFF(NODDIM),IZONE(NODDIM),ICH(NODDIM),
     1          VBVL(2,NTRDIM,NZDIM),VBZNFL(2,0:NZDIM,0:NZDIM)
      DIMENSION ICOMP(MXZWCZ,MXCOMP),NZWCZ(MXCOMP)
      DOUBLE PRECISION VBVL,VBZNFL
      DIMENSION ITIME(2,10)
      CHARACTER*80 TITLE
      CHARACTER*80 NAME
      CHARACTER*16 VBNM(NTRDIM),TEXT
      CHARACTER*1 METHOD,IANS
C     ------------------------------------------------------------------
C     set string for use if RCS ident command
      NAME =
     &’$Id: zonebdgt.f,v 1.0 1996/12/20 13:40:08 rsregan Exp rsregan $’
      NAME = 
     &’@(#)ZONEBDGT - Program for computing subregional water budgets’
      NAME = ‘@(#)           for results of MODFLOW simulations’
      NAME = ‘@(#)ZONEBDGT - USGS Open File Report 90-392, Harbaugh’
      NAME = ‘@(#)ZONEBDGT - Contact: h2osoft@usgs.gov’
      NAME = ‘@(#)ZONEBDGT - Version: 1.0x 1996/12/20’
C     ------------------------------------------------------------------
C
C-----DEFINE INPUT AND OUTPUT UNITS AND INITIALIZE OTHER VARIABLES
      INZN1=10
      INZN2=11
      INBUD=12
      IOUT=13
      K1=0
      K2=0
      NZONES=0
      MSUM=0
C
C-----TELL THE USER WHAT THIS PROGRAM IS
      WRITE(*,*)
      WRITE(*,*) ‘ Modified from ZONEBUDGET version 1.00 to read’
      WRITE(*,*) ‘ cell-by-cell flow data from the USGS Modular Ground-W
     1ater Flow Model.’
C

C-----OPEN A LISTING FILE
      WRITE(*,*)
3     WRITE(*,*) ‘ Enter the name of a LISTING FILE for writing budget f
     1ile contents:’
      READ(*,’(A)’) NAME
      OPEN(UNIT=IOUT,FILE=NAME,ERR=3)
      WRITE(IOUT,4)
4     FORMAT(1H ,’modified from ZONEBUDGET version 1.00’/
     1’ Program to read cell-by-cell flow data and format it for input t
     3o other programs.’)
C
C-----OPEN THE CELL-BY-CELL BUDGET FILE
      WRITE(*,*)
10    WRITE(*,*) ‘ Enter the name of the file containing CELL-BY-CELL BU
     1DGET TERMS:’
      READ(*,’(A)’) NAME
      OPEN(UNIT=INBUD,FILE=NAME,STATUS=’OLD’,FORM=’UNFORMATTED’,ERR=10)
      WRITE(IOUT,*)
      WRITE(IOUT,*) ‘ The cell-by-cell budget file is:’
      WRITE(IOUT,*) NAME
C
C-----READ GRID SIZE FROM BUDGET FILE AND REWIND
      READ(INBUD,END=2000) KSTP,KPER,TEXT,NCOL,NROW,NLAY
      REWIND(UNIT=INBUD)
      WRITE(*,*)
      WRITE(*,14) NLAY,NROW,NCOL
      WRITE(IOUT,*)
      WRITE(IOUT,14) NLAY,NROW,NCOL
14    FORMAT(1X,I3,’ layers’,I10,’ rows’,I10,’ columns’)
C
C-----CHECK TO SEE IF NODIM IS LARGE ENOUGH
      NODES=NCOL*NROW*NLAY
      IF(NODES.GT.NODDIM) THEN
         WRITE(*,*) ‘ PROGRAM ARRAYS ARE DIMENSIONED TOO SMALL’
         WRITE(*,*) ‘ PARAMETER NODDIM IS CURRENTLY’,NODDIM
         WRITE(*,*) ‘ CHANGE NODDIM TO BE’,NODES,’ OR GREATER’
         STOP
      END IF
C
C-----READ A TITLE TO BE PRINTED IN THE LISTING
      WRITE(*,*)
      WRITE(*,*) ‘ Enter a TITLE to be printed in the listing:’
      READ(*,’(A)’) TITLE
      WRITE(IOUT,*)
      WRITE(IOUT,’(1X,A)’) TITLE
C
c  get the pumping period and timestep you want
C
      write (*,*) ‘ Enter the pumping period you would like’
      read (*,*) KPERYES
      write (*,*) ‘ Enter the timestep you would like’
      read (*,*) KSTPYES
C
C-----READ BUDGET DATA AND ACCUMULATE AS LONG AS TIME REMAINS CONSTANT.
C-----WHEN TIME CHANGES, PRINT THE BUDGET, REINITIALIZE, AND START OVER
100   READ(INBUD,END=1000) KSTP,KPER,TEXT,NC,NR,NL
c   Write timestep#, pumping period#, flow term, col,row,lay
c
      write (iout,*)  KSTP,KPER,TEXT,NC,NR,NL
c
c   Read the budget terms
c
      read (INBUD) (BUFF(I),I=1,NODES)
c
c   If it’s the right pumping period, timestep, & flow across faces we want it
c
      open (unit=21,file=’flowritf’)
      open (unit=22,file=’flowfntf’)
      open (unit=23,file=’flowlowf’)
      open (unit=24,file=’flowlftf’)
      open (unit=25,file=’flowbckf’)
      if (kstp .eq. kstpyes .and. kper .eq. kperyes) then
c   see if it is one of the flow faces 
        if (text(1:6) .eq. ‘FLOW R’) then
               INFLOW=21
               write (*,*) ‘right’
            else if (text(1:6) .eq. ‘FLOW F’) then
               INFLOW=22
               write (*,*) ‘front’
            else if (text(1:6) .eq. ‘FLOW L’) then
               INFLOW=23
               write (*,*) ‘lower’
            else
               go to 100
         end if
c  if this is what we want, let’s write it
c
         call writit (buff,inflow,ncol,nrow,nlay)
c
        if (text(1:6) .eq. ‘FLOW R’) then
               INFLOW=24
               write (*,*) ‘left’
               call writleft (buff,inflow,ncol,nrow,nlay)
            else if (text(1:6) .eq. ‘FLOW F’) then
               INFLOW=25
               write (*,*) ‘back’
               call writbck (buff,inflow,ncol,nrow,nlay)
            else
               go to 100
         end if
c
      end if
      go to 100
C
C-----EMPTY BUDGET FILE
2000  WRITE(*,*) ‘CELL-BY-CELL FLOW TERM FILE WAS EMPTY’
      STOP
1000  WRITE(*,*) ‘END OF CELL-BY-CELL FLOW TERM FILE’
      STOP
      end
C
      subroutine writit (buff,inflow,nc,nr,nl)
      dimension buff(nc,nr,nl)
         do 105 ir=1,nr
            do 105 ic=1,nc
               iseq=(ir-1)*nc+ic
               write (inflow,200) ir,ic,iseq,(buff(ic,ir,il),il=1,nl)
 200           format( 2(i4,’,’),i10,’,’,5(e20.6,’,’),e20.6)
 105     continue

      end
c
      subroutine writleft (buff,inflow,nc,nr,nl)
      dimension buff(nc,nr,nl),zero(100)
      do 101 i=1,100
 101     zero(i)=0
         do 105 ir=1,nr
            iseq=(ir-1)*nc+1
            inc=1
               write (inflow,200) ir,inc,iseq,(zero(il),il=1,nl)
            do 105 ic=1,(nc-1)
               inc=ic+1
               iseq=(ir-1)*nc+inc
               do 106 il=1,nl
 106              buff(ic,ir,il)=0-buff(ic,ir,il)
               write (inflow,200) ir,inc,iseq,(buff(ic,ir,il),il=1,nl)
 200           format( 2(i4,’,’),i10,’,’,5(e20.6,’,’),e20.6)
 105     continue

      end
c
      subroutine writbck (buff,inflow,nc,nr,nl)
      dimension buff(nc,nr,nl),zero(100)
      do 101 i=1,100
 101     zero(i)=0
         do 105 ir=0,(nr-1)
            ibr=ir+1
            do 105 ic=1,nc
               iseq=(ibr-1)*nc+ic
              if (ibr .eq. 1) then
               write (inflow,200) ibr,ic,iseq,(zero(il),il=1,nl)
              else
c  Remember the flow on back face will be negitive of what is in buff
               do 106 il=1,nl
 106              buff(ic,ir,il)=0-buff(ic,ir,il)
               write (inflow,200) ibr,ic,iseq,(buff(ic,ir,il),il=1,nl)
              end if
 200           format( 2(i4,’,’),i10,’,’,5(e20.6,’,’),e20.6)
 105     continue

      end