subroutine wrteps
c	kim sharp/9 feb 88
c	reformats epsilon array epsmap
c	to compact array (MKG format)
c	old format: epsmap(65,65,65,3) real*4
c	new format: neps(5,65,65,3) integer*2
c	where first index 1-65 now compressed
c	into 1-5 plus offset into 16 bit words
c	compact format also contains oldmid, the
c	center of the protein in real coordinates
c	compaction is effected by storing
c	real eps (which take values of 0. and 1.)
c	as bits in a 16 bit word
c	access is via pointers idimx and ioffset
c	thus x arrary indices of reps 0-15 -> word 1
c	16-31 -> word 2 etc
c	
	include 'qdiffpar4.h'
	dimension neps(5,ngrid,ngrid,3)	
	dimension keps(5,ngrid,ngrid)	
c	compact fine epsilon map
	dimension idimx(ngrid),ip2(ngrid)
c	array of pointers to words
	dimension ioffset(ngrid)		
c	array of pointers to bit offsets
	integer*2  neps,i,j ,ioffset,keps
c 	integer  neps,ioffset
	character*80 filnam
c---------------------------------------------------------------------

      write(6,*)' setting up pointers...'
	do 9000 ix = 1, ngrid
	  idimx(ix) = ix/16 + 1
	  ioffset(ix) = mod(ix,16)
	  ip2(ix)=2**(mod(ix,16))
9000	continue
	  ip2(15)=-2**15
	  ip2(31)=-2**15
	  ip2(47)=-2**15
	  ip2(63)=-2**15
      write(6,*)' clearing bits...'
      do 9001 idir = 1, 3
        do 9002 iz=1,ngrid
          do 9003 iy=1,ngrid
            neps(1,iy,iz,idir) = 0
            neps(2,iy,iz,idir) = 0
            neps(3,iy,iz,idir) = 0
            neps(4,iy,iz,idir) = 0
            neps(5,iy,iz,idir) = 0
9003	    continue
9002	  continue
9001	continue

      write(6,*)' generating compact fine epsilon array...'
        do 9006 iz=1,ngrid
            do 9008 ix=1,ngrid
	    i=idimx(ix)
          do 9007 iy=1,ngrid
	    j1=ip2(ix)*iepsmp(ix,iy,iz,1) 
	    j2=ip2(ix)*iepsmp(ix,iy,iz,2) 
	    j3=ip2(ix)*iepsmp(ix,iy,iz,3) 
	    neps(i,iy,iz,1)=neps(i,iy,iz,1)+j1
	    neps(i,iy,iz,2)=neps(i,iy,iz,2)+j2
	    neps(i,iy,iz,3)=neps(i,iy,iz,3)+j3
9007	    continue
9008		continue
9006	  continue

        do 9016 iz=1,ngrid
          do 9017 iy=1,ngrid
            do 9014 ix=1,5
              keps(ix,iy,iz) = 0
9014		continue
            do 9018 ix=1,ngrid
              if(idebmap(ix,iy,iz).ne.0) then
c set bit
                keps(idimx(ix),iy,iz) = ibset(
     &          keps(idimx(ix),iy,iz),ioffset(ix))
              end if
9018		continue
9017	    continue
9016	  continue
	kmap = 1

      write(6,*)' writing to compact epsilon file'
	open(unit=17,form='unformatted')
	filnam = ' '
	inquire(17,name = filnam)
	write(6,*)' '
	write(6,*)'dielectric map written to file'
	write(6,*)filnam
	write(6,*)' '
	imap = 0
	write (17) kmap, scale, oldmid
	write (17) neps
	write (17) keps
	close (17)
      return
	end