C  Subroutine READBSQ -- Read and unpack band-sequential binary array
C
C  Purpose:
C     This subroutine extracts a specified window (rectangular subset)
C     from a band-sequential (BSQ) binary data array and places it into
C     an INTEGER*2 array specified by the calling program.  In a
C     band-sequential array, all data for the first band come first,
C     followed by all data for the second band, and so on until the last
C     band is complete.  Within each band, the data are organized by
C     rows.
C
C     The 'bands' of the array may represent spectral bands, geographic
C     layers, successive time intervals, or the like.
C
C     This subroutine may be used to unpack a two-dimensional binary
C     array by specifying inband = 1.
C
C     For multi-byte data values, this subroutine assumes that the input
C     array is in "big endian" format; i.e., the most-significant byte
C     of each value precedes the least-significant byte.  If the host
C     computer being used to read the data uses "little endian" format,
C     the compile-time parameter BIGEND should be changed -- see the
C     PARAMETER statements below.
C
C     This subroutine opens and closes the input file each time it is
C     called.  If it is to be called repeatedly from a main program
C     using different window specifications (e.g., to read in a small
C     number of rows at a time), it will be more efficient to move the
C     OPEN and CLOSE statements into the calling program.
C
C  Usage:
C     CALL READBSQ (filename, lun, inband, inrow, incol, inbyte,
C    +   iband, irow, icol, nband, nrow, ncol, outarray)
C
C  Arguments:
C     CHARACTER*(*) filename    Name of file containing the BSQ data
C                               array.
C     INTEGER*4     lun         Logical unit number to be used for
C                               accessing the input data.
C     INTEGER*4     inband	Number of bands in input data array; for
C                               two-dimensional array, inband = 1.
C     INTEGER*4     inrow	Number of rows in input data array.
C     INTEGER*4     incol	Number of columns per row in input data
C                               array.
C     INTEGER*4     inbyte  	Number of bytes per input data value.
C                               Allowed values are 1 and 2.
C     INTEGER*4     iband	Initial band for which data are desired.
C     INTEGER*4     irow	Initial row for which data are desired:
C                               rows are numbered from top to bottom
C                               (north to south), beginning with row 1.
C     INTEGER*4     icol	Initial column for which data are
C                               desired:  columns are numbered from left
C                               to right (west to east), beginning with
C                               column 1.
C     INTEGER*4     nband	Number of consecutive bands to be
C                               extracted from input data array
C     INTEGER*4     nrow	Number of consecutive rows to be
C                               extracted from each band of data array.
C     INTEGER*4     ncol	Number of consecutive columns to be
C                               extracted from each row of data array.
C     INTEGER*2     outarray(ncol,nrow,nband)
C                               Array to receive the data from input
C                               file.  If nband = 1, the array may be
C                               declared as outarray(ncol,nrow).
C
C  Compile-time parameters
C
C     BIGEND          Byte offset in memory from start of an INTEGER*2
C                     value to the most significant byte of the value.
C                     Set to 1 if host computer stores the least
C                     significant byte of integer data values at a
C                     lower memory address than the most significant
C                     byte ("little-endian"); set to 0 otherwise.
C     MAXINBYTE       Maximum number of bytes per row in input data
C                     array, MAXINBYTE >= incol.
C
C  Addtional variables and arrays used:
C
C     band            Current input band.
C     bytebuf         Buffer for moving bytes into integers.
C     col             Current input column.
C     inbuf           Buffer to receive one record of input data.
C                     declared as type CHARACTER.
C     intbuf	      Buffer for moving bytes into integers, declared
C                     INTEGER*2, equivalenced to  bytebuf  array.
C     row             Current input row.
C
C  Program history:
C     01-96   Initial version, R. A. White, PSU/ESSC,
C
C  **********  End of prolog for subroutine READBSQ  **********
C
C
      SUBROUTINE READBSQ (filename, lun, inband, inrow, incol, inbyte,
     + iband, irow, icol, nband, nrow, ncol, outarray)
C
C  Use BIGEND = 0 for big-endian hosts, such as Sun SPARC
C
      PARAMETER (BIGEND = 0)
C
C  Use BIGEND = 1 for little-endian hosts, such as VAX and Intel
C  80x86-compatible PCs
C
C     PARAMETER (BIGEND = 1)
C
      PARAMETER (MAXINBYTE = 16384)
C
C  In the following format statement, the repeat count must match the
C  value of MAXINBYTE
C
 1000 FORMAT (16384A1)
C
      CHARACTER*(*) filename
      INTEGER*4 lun, inband, inrow, incol, inbyte, iband, irow, icol
      INTEGER*4 nband, nrow, ncol
      INTEGER*2 outarray(ncol, nrow, nband)
C
      INTEGER*4 band, row, col
      INTEGER*2 intbuf
      CHARACTER*1 bytebuf(2)
      CHARACTER*1 inbuf(2*MAXINBYTE)
      EQUIVALENCE (intbuf, bytebuf)
C
      intbuf = 0
      OPEN (UNIT=lun, FILE=filename, ACCESS='DIRECT',FORM='FORMATTED',
     +      RECL=incol*inbyte)

      DO 100 band=iband,iband+nband-1
         DO 80 row=irow,irow+nrow-1
            READ (lun, REC=row+(band-1)*inrow, FMT=1000)
     +       (inbuf(i), i=1,incol*inbyte)
            DO 60 col=icol,icol+ncol-1
C
C  Unpack 1-byte data to 2-byte integers
C
               IF (inbyte .GT. 1) GOTO 20
                  bytebuf(2-BIGEND) = inbuf(col)
               GOTO 40
C
C  Extract two-byte data from input buffer           
C
 20            CONTINUE
                  bytebuf(BIGEND + 1) = inbuf(2*col-1)
                  bytebuf(2 - BIGEND) = inbuf(2*col)
 40            CONTINUE
               outarray(col-icol+1, row-irow+1, band-iband+1) = intbuf
 60         CONTINUE
 80      CONTINUE
 100  CONTINUE
      CLOSE (lun)

      RETURN
      END
