program LCread
        
c  Example program which reads from LC95 data files and writes
c  demultiplexed data to files for external plotting or processing.
c  This would be version 6.0 of LCdump and its relatives
c  S.Constable 10/Nov/2000 
c  Modified 25 July 2002 for compass data

c  output file name root
      character*64 filename

c  header common
      integer*4 dir_count, dir_block, dir_start, data_start,
     * num_channel, sam_rate, data_type, slow_start,
     * slow_stop, slow_entry
      character*1 version(8)
      common /header/ dir_count, dir_block, dir_start, data_start,
     * num_channel, sam_rate, data_type, slow_start, slow_stop,
     * slow_entry, version

c  counter for record numbers in binary files
      integer*4 irec(4)

c  dimension the data arrays, leaving room for one block of over-run
      integer*4 idata(4,1249)
      real*4 fdata(4,1249)
      integer iskip

c common to indicate the presence of a PC
      logical lpc
      common /ifpc/ lpc

	logical ifdisk

c  set the pc flag
      lpc= .false.
c  set the unit number for the SIO file
      in_unit = 12

c  set the number of data per read
      ndata = 1000

c  open the SIO binary file- this will tell us how many channels we have
      call OpenSIO(in_unit)
      	
c  find out which blocks we want
      call GetRecs(jrec, iendrec, ns, iskip)
c	write(*,*) jrec, iendrec, ns, iskip
        
      write(*,*) ' Do you want binary(1) or ascii(2) data written?'
      read (*,*) idtype

c get the output filename:
      do i = 1, 64
        filename(i:i) = ' '
      end do
      write(*,*) ' Enter data filename root: '
      read (*,'(a)') filename
      if (filename(1:1) .eq. ' ') then
        stop
      end if
      j = index(filename, ' ')-1
      filename = filename(1:j)//'_c'
      j = j + 3
       
c open output files (one for each channel)
      do i = 1, num_channel
        filename(j:j) = char(48+i)
        write(*,*) ' Opening ', filename
        if (idtype .eq. 2) then
          open(12+i,file=filename(1:j))
        else
          open(12+i,file=filename(1:j),form='unformatted',
     *     access='direct',recl=ndata*4)
        end if
      end do

c  reset record counters for binary files
      do  i = 1,num_channel
        irec(i) = 0
      end do
c  keep track of total number of data points read
      numdat = 0
	 
c  get some data to play with
1     call ReadSIO(in_unit, jrec, iskip, ndata, idata, icount, ifdisk)
      do j=1,num_channel
c  copy the integer data into a floating point array
        do i=1,ndata
           fdata(j,i)=idata(j,i) 
        end do
c  binary data write
        if (idtype .eq. 1) then
           irec(j) = irec(j)+1
           write(unit=12+j,rec=irec(j)) (fdata(j,i),i=1,ndata)
        else
c  ascii data write
           write(unit=12+j,fmt='(f10.0)') (fdata(j,i), i=1, ndata)
        end if
      end do
      numdat = numdat + ndata
c      write(*,*) jrec, iskip, icount, ifdisk, numdat
	if (numdat .lt. ns) goto 1

      close (in_unit)
      do i = 1, num_channel
         close (12+i)
      end do
      stop

      end
     

c-----------------------------------------------------------------------
      subroutine GetRecs(jrec, iendrec, ns, iskip)
c-----------------------------------------------------------------------
 
c  GetRecs works out which block numbers we want to read
c  S.Constable 10/Nov/2000 
          
c  header common
      integer*4 dir_count, dir_block, dir_start, data_start,
     * num_channel, sam_rate, data_type, slow_start,
     * slow_stop, slow_entry
      character*1 version(8)
      common /header/ dir_count, dir_block, dir_start, data_start,
     * num_channel, sam_rate, data_type, slow_start, slow_stop,
     * slow_entry, version

c  data start time
      integer*4 isy, ismo, isd, ish, ism, jrec
      real fss
      common /startime/ isy, ismo, isd, ish, ism, fss
       
c  common for number of data per block
      integer nblk, diskblks(2000)
      common /blk/ nblk, sampfreq, diskblks
      
      integer iskip

      write(*,*)
     *  ' Stop (0), Recover by Time (1) or by Block Number (2)?'
      read (*,*) iitype

      if (iitype.eq.0) then       
        close (in_unit)
        stop
      else if (iitype.eq.1) then
c  Peform the classic Y2K patch on the year number
        if (isy .le. 50) then
           isy = isy + 2000
        else
           isy = isy + 1900
        end if

        write(*,'(a16,i5,4i3,f7.3)') 
     * ' SIO data starts at: ',isy, ismo, isd, ish, ism, fss

        write(*,*) ' Enter required data start time:'//
     *        ' yr, mo, day, hour, min, sec'
        read (*,*) iy, imo, id, ih, im, is

        write(*,*) ' Enter total data time: days, hours, minutes'
        read (*,*) itd, ith, itm
         
c  Compute starting and ending blocks:

c  modified julian day of SIO start time:
        ismjd = mjd(isy, ismo, isd)
	
c  day second number of SIO data start
        it0 = int(fss) + 60*(ism + 60*ish)
	
c  modified julian day of data start time:
        imjd = mjd(iy, imo, id)
	
c  day second number of  required data start
        istart = is + 60*(im + 60*ih) 
	
c  total number of samples per channel from SIO start to data start
        ns = int(sampfreq*(istart-it0 + 86400.*(imjd-ismjd)))
	
c  total number of blocks/ch from data start (rounded down)
        nblocks = ns/nblk
	
c  so the block number in which the first datum lies is:
        jrec = data_start + nblocks*num_channel
	
c  and this is the number of data in that block to skip past:
        iskip = ns - nblocks*nblk
                 
c  total number of samples requied:                
        ns = 60*(itm + 60*(ith + 24*itd))*sampfreq
        nblks = int(ns/nblk)
        if (mod(nblks*nblk, ns) .ne. 0) nblks = nblks+1
        nblocks = nblks*num_channel
        iendrec = jrec + nblocks -1 
               
      else
         iskip=0  
         write(*,'(''Enter starting block number:  '',$)')
         read (*,*) jrec
         write(*,'(''Enter ending block number:  '',$)')
         read (*,*) iendrec
	   ns = nblk*(iendrec-jrec+1)/num_channel
      end if
           
      if (jrec .lt. data_start) then
         write(*,*) ' Sorry, start time or block preceeds data'
         stop
      end if
      
      write (*,*) ' Will try to read data blocks', jrec, iendrec   
      return
	end


c-----------------------------------------------------------------------
      subroutine OpenSIO(in_unit)
c-----------------------------------------------------------------------
 
c  opens an SIO binary LC95 data file and prepares life for reading it
c  S.Constable 9/Nov/2000 

      character*80 cfin
          
c  header common
      integer*4 dir_count, dir_block, dir_start, data_start,
     * num_channel, sam_rate, data_type, slow_start,
     * slow_stop, slow_entry
      character*1 version(8)
      common /header/ dir_count, dir_block, dir_start, data_start,
     * num_channel, sam_rate, data_type, slow_start, slow_stop,
     * slow_entry, version
       
c  common for number of data per block
      integer nblk, diskblks(2000)
      common /blk/ nblk, sampfreq, diskblks

c  flags for 16 or 24 bit
      logical l_24
      common /bits/ l_24

1     write(*,'(''Enter SIO binary image file name:  '',$)')
      read (*,'(a80)') cfin

c   open binary file and read header 
      open(unit=in_unit, file=cfin, form='unformatted', 
     * access='direct', recl=512, status='old', err = 888)
 
c  read the Disk Header
      call GetHeader(in_unit)      

      nch = num_channel
c  correct fractional sample rates on 24-bit systems
      if (sam_rate .eq. 31) then
        sampfreq = 31.25
      else if (sam_rate .eq. 62) then
        sampfreq = 62.5
      else 
        sampfreq = sam_rate
      end if 
         
c  work out what sort of data type we have and set number of samples per block         
      if (data_type.gt.1) then
            l_24 = .true.
            nblk = 166
            write(*,*) 'data_type is 24-bit'
      else 
            l_24 = .false.
            nblk = 249
            write(*,*) 'data_type is 16-bit'
      end if
         
c read the direcory
      call GetDirectory(in_unit)

c  get the compass data
      call GetComp(in_unit)
     
      return

888   write(*,*) ' Error opening image file (does not exist?)'
      goto 1
     
       
      end 
 
c-----------------------------------------------------------------------
      subroutine ReadSIO(in_unit,jrec,iskip,ndata,idat,icount,ifdisk)
c-----------------------------------------------------------------------

c  reads data from an SIO binary LC95 data file
c  S.Constable 9/Nov/2000 

c  on input:
c    jrec = starting block number
c    iskip = number of points to skip in first block (jrec)
c    ndata = number of data required
c  on output;
c    idat(4,ndata) = raw count data on up to 4 channels
c    icount = actual number of data in idat
c    jrec = block number of last read
c    iskip = data number in block of last read

c  common for number of data per block, sam freq., and disk write indices
      integer nblk, diskblks(2000)
      common /blk/ nblk, sampfreq, diskblks
      
c  flags for 16 or 24 bit
      logical l_24
      common /bits/ l_24
        
c  header common
      integer*4 dir_count, dir_block, dir_start, data_start,
     * num_channel, sam_rate, data_type, slow_start,
     * slow_stop, slow_entry
      character*1 version(8)
      common /header/ dir_count, dir_block, dir_start, data_start,
     * num_channel, sam_rate, data_type, slow_start, slow_stop,
     * slow_entry, version

      integer*2 ihead(12)
c  common in which data are passed back
      integer*4 idata24(256)
      common /data24/ idata24
	
      logical ifdisk
      integer iskip

c  arguments     
      integer*4 idat(4,1249)

c  We pack 16 bit and 24 bit data into 32 bit integers; to get them 
c  back to the right numbers we have to byte-shift them by division.
c  Also compute the length of a disk write in blocks (2 min. for 16-bit
c  and 4 min. for 24 bit)
      if (l_24) then
        idivid = 256
	  ndisk = 4.*60.*sampfreq/nblk
      else 
        idivid = 256*256
	  ndisk = 2.*60.*sampfreq/nblk
      end if
         
	ifdisk = .false.
c work out how many blocks of data we need
      nblocks = (ndata + iskip)/nblk
      if (mod(nblocks*nblk, ndata+iskip) .ne. 0) nblocks = nblocks+1
      iend = jrec+nblocks-1
      do ii=1,dir_count
          if((diskblks(ii).ge.jrec).and.(diskblks(ii).le.iend)
     *         .or.
     *   (diskblks(ii)+ndisk.ge.jrec).and.(diskblks(ii)+ndisk.le.iend))
     *    then 
             ifdisk = .true.
	    end if
	 end do
c set record pointer to starting block
      irec = jrec
      icount = 0
      is = iskip
    
      do k = 1, nblocks
         do i = 1, num_channel
            call GetaBlock(in_unit, irec, ihead)
            irec = irec + 1
c  check to make sure channel number in header meets expectations
            if (ihead(10)+1 .ne. i) then
                write(*,*) ' readsio: channels out of synch, i= ', i
                write(*,'(i8,6i4,''.'',i3.3,4i4)') l, 
     *           (ihead(j),j=8,3,-1), ((ihead(1)*256)+ihead(2)), 
     *           (ihead(j),j=9,12)
            end if

            do j=1,nblk-is
              idat(i,j+icount)=idata24(j+is)/idivid
            end do
         end do 
         icount = icount + nblk - is
         is = 0
      end do
      jrec = irec - num_channel
	iskip = nblk - (icount - ndata)
      return
      end


c----------------------------------------------------------------------
      subroutine GetaBlock(in_unit, nrec, ihead)
c----------------------------------------------------------------------

c Gets a block of LC-95 data from various devices and peels off the
c 12 header bytes (ihead in argument list).
c The main chunk of data is passed back through /data/, which will
c automatically perform the 2's complement conversion by declaring
c it character*512 in Getablock and integer*2 times 256 in the calling
c routine.
c (This was an easy way to handle the mix of unsigned integer
c information (the header bytes) and 2's compliment (the actual data))
c  S.Constable 9/Nov/2000 

c arguments
      integer*4 nrec
      integer*2 ihead(12)
c binary data buffer
      character*1 buffer(512), buffer24(1024), nul
      character*512 buff
      equivalence (buffer(1), buff(1:1))
      common /data24/ buffer24
      
      logical lpc
      common /ifpc/ lpc
c  flags for 16 or 24 bit
      logical l_24
      common /bits/ l_24
c  common for number of data per block
      integer nblk, diskblks(2000)
      common /blk/ nblk, sampfreq, diskblks
      
c reading from disk image file only; add 1 for Fortran files
      read(in_unit,rec=nrec+1,err=19) buff

c peel off the header data while we still have the block in 
c character format
      do i = 1,12
        ihead(i) = ichar(buffer(i))
	end do

c  we solve all our problems, (PC/MAC/SUN byte ordering, variable ADC width,
c  and conversion to 2's compiment) by packing the data bytes into a int*4
c  word, here represented as a character*1 array

      nul = char(0)
      if(l_24) then
         if (lpc) then
            do i = 0, nblk-1
               buffer24(4*i+1) = nul
               buffer24(4*i+2) = buffer(3*i+17)
               buffer24(4*i+3) = buffer(3*i+16)
               buffer24(4*i+4) = buffer(3*i+15)
            end do
         else
            do i = 0, nblk-1
               buffer24(4*i+1) = buffer(3*i+15)
               buffer24(4*i+2) = buffer(3*i+16)
               buffer24(4*i+3) = buffer(3*i+17)
               buffer24(4*i+4) = nul
             end do
          end if
       else 
          if (lpc) then
            do i = 0, nblk-1         
               buffer24(4*i+1)=nul
               buffer24(4*i+2)=nul
               buffer24(4*i+3)=buffer(2*i+16)
               buffer24(4*i+4)=buffer(2*i+15)
             end do
          else
             do i = 0, nblk-1
                buffer24(4*i+1) = buffer(2*i+15)
                buffer24(4*i+2) = buffer(2*i+16)
                buffer24(4*i+3) = nul
                buffer24(4*i+4) = nul
             end do

          end if
      end if

      return
19    write(*,*) ' Getablock: read error from image file'
      stop
      end
c-----------------------------------------------------------------------
      integer function mjd(i,j,k)
c-----------------------------------------------------------------------

c  on input:
c    i = year
c    j = month
c    k = day
c  on output:
c    mjd = modified Julian day

c  From a JD routine provided me by Duncan Agnew, who got it from somewhere
c  else, and altered by me to return Modified Julian Day
c  S.Constable 9/Nov/2000 (MJD 51857)

      mjd = -678927 + k + 1461*(i+(j-14)/12)/4 + 367*(j-2-12 *
     $           ((j-14)/12))/12 + (24002-12*i-j)/1200 
      return
      end
      
c----------------------------------------------------------------------
      subroutine GetDirectory (in_unit)
c----------------------------------------------------------------------

c This version for LCheapo-95 data structure and hard-wired for disk 
c images
c Copyright Steven Constable, Scripps Institution of Oceanography.
      
c read buffer
      character*1 buffer(512)
c local variables
      integer*4 ibuff(512)
c  header common
      integer*4 dir_count, dir_block, dir_start, data_start,
     * num_channel, sam_rate, data_type, slow_start,
     * slow_stop, slow_entry
      character*1 version(8)
      common /header/ dir_count, dir_block, dir_start, data_start,
     * num_channel, sam_rate, data_type, slow_start, slow_stop,
     * slow_entry, version

c data start time
      integer*4 isy, ismo, isd, ish, ism
      real fss
      common /startime/ isy, ismo, isd, ish, ism, fss

c  common for number of data per block
      integer nblk, nblock(2000)
      common /blk/ nblk, sampfreq, nblock

      write(*,*) ' Reading directory'
      write(*,'(a)') ' block no yr:mo:dy:hr:mn:seconds rate blksz 
     * reclen flg mux'

      ic = 0
      do 600 nb = dir_start, dir_block
        read(in_unit,rec=nb+1,err=19) (buffer(k), k=1,512)

        do 500 k = 1, 16
	    ic = ic+1
          if (ic. gt. dir_count) goto 700
          index = (k-1)*32
          do i = 1,32
            ibuff(i) = ichar(buffer(index+i))
          end do

          ims = ibuff(1)*256+ibuff(2)
          seconds = float(ibuff(3)) + float(ims)/1000.
          nblock(ic) = ((ibuff(9)*256 + ibuff(10))*256 + 
     *       ibuff(11))*256 + ibuff(12)
             
          lengthrec = ((ibuff(13)*256 + ibuff(14))*256 + 
     *       ibuff(15))*256 + ibuff(16)
          irate = ibuff(17)*256 + ibuff(18)
          nblocks = ibuff(19)*256 + ibuff(20)

          if (data_start .eq. nblock(ic)) then
              isy = ibuff(8)
              ismo = ibuff(7)
              isd = ibuff(6)
              ish = ibuff(5)
              ism = ibuff(4)
              fss = seconds
          end if
          write(*,'(i9,1x,5(i2,'':''),f6.3,i5,i6,i6,2(1x,i5))')
     *      nblock(ic),(ibuff(j),j=8,4,-1),
     *      seconds,irate, nblocks, lengthrec,ibuff(21),ibuff(22)
500     continue

600   continue
700   return
19    write(*,*) ' read error from image file'
      stop
      end

c----------------------------------------------------------------------
      subroutine GetHeader (in_unit)
c----------------------------------------------------------------------

c Copyright Steven Constable, Scripps Institution of Oceanography.
c The header is always in block 2 (= 3 in fortran binary file)
      
c read buffer
      character*1 buffer(512)

c  header common
      integer*4 dir_count, dir_block, dir_start, data_start,
     * num_channel, sam_rate, data_type, slow_start,
     * slow_stop, slow_entry
      character*1 version(8)
      common /header/ dir_count, dir_block, dir_start, data_start,
     * num_channel, sam_rate, data_type, slow_start, slow_stop,
     * slow_entry, version

      read (in_unit,rec=3,err=19) (buffer(k), k=1,512)

c here is where the 512-byte block of header information is decoded.

      i = 12
      dir_start = ichar(buffer(i+4)) + 256*(ichar(buffer(i+3)) +
     * 256*(ichar(buffer(i+2)) + 256*ichar(buffer(i+1))))
      write(*,*) 'current directory start:  ', dir_start

      i = i + 4
      write(*,*) 'current directory size:  ', ichar(buffer(i+4)) +
     * 256*(ichar(buffer(i+3)) +
     * 256*(ichar(buffer(i+2)) + 256*ichar(buffer(i+1))))

      i = i + 4
      dir_block = ichar(buffer(i+4)) + 256*(ichar(buffer(i+3)) +
     * 256*(ichar(buffer(i+2)) + 256*ichar(buffer(i+1))))
      write(*,*) 'current directory block:  ', dir_block

      i = i + 4
      dir_count = ichar(buffer(i+4)) + 256*(ichar(buffer(i+3)) +
     * 256*(ichar(buffer(i+2)) + 256*ichar(buffer(i+1))))
      write(*,*) 'current directory count:  ', dir_count

      i = i + 4
      slow_start = ichar(buffer(i+4)) +
     * 256*(ichar(buffer(i+3)) +
     * 256*(ichar(buffer(i+2)) + 256*ichar(buffer(i+1))))
      write(*,*) 'slow data start block:  ', slow_start

      i = i + 4
      write(*,*) 'slow data size:  ', ichar(buffer(i+4)) +
     * 256*(ichar(buffer(i+3)) +
     * 256*(ichar(buffer(i+2)) + 256*ichar(buffer(i+1))))

      i = i + 4
      slow_stop = ichar(buffer(i+4)) +
     * 256*(ichar(buffer(i+3)) +
     * 256*(ichar(buffer(i+2)) + 256*ichar(buffer(i+1))))
      write(*,*) 'slow data current block:  ', slow_stop
      
      i = i + 4
      slow_entry = ichar(buffer(i+4)) +
     * 256*(ichar(buffer(i+3)) +
     * 256*(ichar(buffer(i+2)) + 256*ichar(buffer(i+1))))
      write(*,*) 'slow data current entry:  ', slow_entry 

      i = i + 4
      write(*,*) 'log data start block:  ', ichar(buffer(i+4)) +
     * 256*(ichar(buffer(i+3)) +
     * 256*(ichar(buffer(i+2)) + 256*ichar(buffer(i+1))))

      i = i + 4
      write(*,*) 'log data size:  ', ichar(buffer(i+4)) +
     * 256*(ichar(buffer(i+3)) +
     * 256*(ichar(buffer(i+2)) + 256*ichar(buffer(i+1))))

      i = i + 4
      write(*,*) 'log data current block:  ', ichar(buffer(i+4)) +
     * 256*(ichar(buffer(i+3)) +
     * 256*(ichar(buffer(i+2)) + 256*ichar(buffer(i+1))))

      i = i + 4
      write(*,*) 'log data current byte:  ', ichar(buffer(i+4)) +
     * 256*(ichar(buffer(i+3)) +
     * 256*(ichar(buffer(i+2)) + 256*ichar(buffer(i+1))))

      i = i + 4
      data_start = ichar(buffer(i+4)) + 256*(ichar(buffer(i+3)) +
     * 256*(ichar(buffer(i+2)) + 256*ichar(buffer(i+1))))
      write(*,*) 'data start block:  ', data_start
      i = i + 4

      write(*,*) 'Disk number:  ', ichar(buffer(i+2)) + 
     * 256*ichar(buffer(i+1))
      i = i + 2

      do j = 1, 8
          version(j) = buffer(i+j)
      end do
      write(*,'(a22)') ' Program version: ', (version(j),j=1,4)
      i = i + 10

      write(*,'(95a)') ' Description:          ',
     * (buffer(j),j=i+1,i+80)

      i = i + 80
      sam_rate = ichar(buffer(i+2)) + 256*ichar(buffer(i+1))
      write(*,*) 'Sample rate:  ', sam_rate

      i = i + 2
      write(*,*) 'Start Channel:  ', ichar(buffer(i+2)) +
     * 256*ichar(buffer(i+1))

      i = i + 2
      num_channel =  ichar(buffer(i+2)) + 256*ichar(buffer(i+1))
      write(*,*) 'Number of Channels:  ', num_channel

      i = i + 2
      write(*,*) 'Slow data rate:  ', ichar(buffer(i+2)) +
     * 256*ichar(buffer(i+1))

      i = i + 2
      write(*,*) 'Slow data start channel:  ', ichar(buffer(i+2)) +
     * 256*ichar(buffer(i+1))

      i = i + 2
      write(*,*) 'Number of slow channels:  ', ichar(buffer(i+2)) +
     * 256*ichar(buffer(i+1))

      i = i + 2
      write(*,*) 'Data type (compressed/not):  ', ichar(buffer(i+2)) +
     * 256*ichar(buffer(i+1))
      data_type = ichar(buffer(i+2)) + 256*ichar(buffer(i+1))
      
      i = i + 2
      write(*,*) 'Disk size:  ', ichar(buffer(i+2)) + 
     * 256*ichar(buffer(i+1))

      i = i + 2
      write(*,*) 'RAM size:  ', ichar(buffer(i+2)) + 
     * 256*ichar(buffer(i+1))

      i = i + 2
      write(*,*) 'number of windows:  ', ichar(buffer(i+2)) +
     * 256*ichar(buffer(i+1))

      i = i + 2
c rest of block reserved for future use
      return
19    write(*,*) ' read error from image file'
      stop
      end

c----------------------------------------------------------------------
      subroutine GetComp(in_unit)
c----------------------------------------------------------------------
c Copyright Steven Constable, Scripps Institution of Oceanography.
c Reads compass data off the file
      
c read buffer
      character*1 buffer(512)


c  header common
      integer*4 dir_count, dir_block, dir_start, data_start,
     * num_channel, sam_rate, data_type, slow_start,
     * slow_stop, slow_entry
      character*1 version(8)
      common /header/ dir_count, dir_block, dir_start, data_start,
     * num_channel, sam_rate, data_type, slow_start, slow_stop,
     * slow_entry, version
     
      character*4 intbuf
      character*22 comstring

c check that the version is 8.08 or later      
      do j = 1, 4
        intbuf(j:j) = version(j)
      end do
     
      read(intbuf, '(f4.2)') rv
      if (rv .lt. 8.08) return

c check to see that the block numbers are reasonable
      if (slow_start .le. 2 .or. slow_start .ge. 32000) then
        write(*,*) ' No compass installed???'
        return
      end if

c write out the compass data      
      write(*,*) ' '
      write(*,*) ' tilt1, tilt2, compass, temp:'

c take care of whole blocks
      num = 0
      t1sum = 0.
      t2sum = 0.
      comsum = 0.
      tempsum = 0.
      do j = slow_start, slow_stop-1
      
        if (slow_start .ne. slow_stop) then 
          read (in_unit,rec=j+1,err=19) (buffer(k), k=1,512)
	  do i = 0,7
	    do k = 1,22
              comstring(k:k) = buffer(i*64+k)
	    end do
            write(*,'(x,22a)') comstring
	    num = num+1
	    read (comstring,*) t1, t2, com, temp
	    if (num .gt. 1) then
	      t1sum = t1sum + t1
              t2sum = t2sum + t2
              comsum = comsum + com
              tempsum = tempsum + temp
            end if
	  end do
	end if
      end do

c now do the partially filled block
      read (in_unit,rec=j+1,err=19) (buffer(k), k=1,512)      
      do i = 0,slow_entry-1
	do k = 1,22
          comstring(k:k) = buffer(i*64+k)
	end do
        write(*,'(x,22a)') comstring
        num = num+1
	read (comstring,*) t1, t2, com, temp
	if (i .ne. slow_entry-1 .and. num .gt. 1) then
	  t1sum = t1sum + t1
          t2sum = t2sum + t2
          comsum = comsum + com
          tempsum = tempsum + temp
        end if
      end do

      write(*,*) '   '
c write out averages(neglecting first and last measurements)
      write(*,*) ' averages:'
      write(*,*) t1sum/(num-2), t2sum/(num-2), 
     * comsum/(num-2), tempsum/(num-2)
      return
19    write(*,*) ' read error from image file'
      stop
      end