      PROGRAM MAPACK
C
C     PROGRAMMER: R.S. ZHOU
C	          BOSTON COLLEGE, 10/89
C
C     THIS ROUTINE PACKS ELECTRON DENSITY FILE INTO BYTE LAYERS
C     ALONG X,Y,Z RESPECTIVELY FOR 3D CONTOURING.
C
      CHARACTER FNAME*32, TITLE*64
      INTEGER*4 RANG(6)
      REAL*4    CELL(6), LINE(128)
      BYTE      BLOCK(128,128,128)

   1  WRITE(6,1000)
 1000 FORMAT(' ENTER MAP FILE: ',$)
      READ(5,1001,END=999) FNAME
 1001 FORMAT(A32)
      IF( FNAME.EQ.' ' ) GOTO 1
      OPEN( UNIT=1, FILE=FNAME, STATUS='OLD', FORM='UNFORMATTED' )
C
C  READ THE HEADER INFORMATION...
C
      READ(1,ERR=111,END=222) TITLE, CELL, MX,MY,MZ,RANG
      WRITE(6,2000) TITLE, CELL, MX,MY,MZ,RANG
 2000 FORMAT(' TITLE: ',A64,/' UNIT CELL: ',6F10.4,/
     +' GRID SIZE: ',3I5,' RANG:',6I5)
   3  WRITE(6,3000)
 3000 FORMAT(' ENTER OUTPUT FILENAME: ',$)
      READ(5,1001,END=999) FNAME
      IF( FNAME.EQ.' ' ) GOTO 3
      NX = RANG(2)-RANG(1)+1
      NY = RANG(4)-RANG(3)+1
      NZ = RANG(6)-RANG(5)+1
C
C  READ DENSITY FILE ONE PASS TO GET SCALE FACTOR AND OFFSET...
C
      SUMA = 0.0
      SUMB = 0.0
      DMAX =-1.0E37
      DMIN = 1.0E37
      DO I=1, NZ
         DO J=1, NY
             READ(1,ERR=111,END=222) (LINE(K),K=1,NX)
             DO K=1, NX
     		TEMP = LINE(K)
                IF( TEMP.GT.DMAX ) DMAX = TEMP
                IF( TEMP.LT.DMIN ) DMIN = TEMP
                SUMA = SUMA + TEMP
                SUMB = SUMB + TEMP*TEMP
             ENDDO
         ENDDO
      ENDDO
      K = NX*NY*NZ
      SIGMA = SQRT( (SUMB-SUMA*SUMA/K)/(K-1) )
      PRINT *, 'DMIN =', DMIN,';  DMAX =', DMAX,';  SIGMA=', SIGMA
      SCALE = 250.0/(DMAX-DMIN)
      OFFSET=-125.0-DMIN*SCALE
      PRINT *, 'SCALE=', SCALE,';   OFFSET=', OFFSET
C
C  WRITE THE OUTPUT FILE HEADER...
C
      I = NX*NY
      J = NX*NZ
      K = NY*NZ
      IREC = MAX( I, J, K ) / 4 + 1   !RECORDLENGTH IN 4-BYTES
      OPEN( UNIT=2, FILE=FNAME, STATUS='NEW', FORM='UNFORMATTED',
     +      ACCESS='DIRECT', RECL=IREC, RECORDTYPE='FIXED' )
      WRITE(2,REC=1) TITLE, MX,MY,MZ, RANG, CELL,
     +               SCALE, OFFSET, SIGMA, DMIN, DMAX
C
C  READ DENSITY LAYERS AND PACK INTO NEW FILE...
C
      REWIND(1)
      READ(1)
      DO I=1, NZ
         DO J=1, NY
             READ(1,ERR=111,END=222) (LINE(K),K=1,NX)
             DO K=1, NX
                BLOCK(K,J,I) = INT(SCALE*LINE(K)+OFFSET)
             ENDDO
         ENDDO
      ENDDO
      IREC = 1
      DO K=1, NX
         IREC = IREC + 1
         WRITE(2,REC=IREC) ((BLOCK(K,J,I),J=1,NY),I=1,NZ)
      ENDDO
      DO J=1, NY
         IREC = IREC + 1
         WRITE(2,REC=IREC) ((BLOCK(K,J,I),K=1,NX),I=1,NZ)
      ENDDO
      DO I=1, NZ
         IREC = IREC + 1
         WRITE(2,REC=IREC) ((BLOCK(K,J,I),K=1,NX),J=1,NY)
      ENDDO
      GOTO 999
  111 STOP 'OOPS: ERROR DURING READING INPUT FILE'
  222 STOP 'OOPS: EOF DURING READING INPUT FILE'
  999 CALL EXIT
      END
