subroutine read_pck(lun,filnam,npixel,img,x,y,ier)
! expands the pck-compressed file 'filnam' into the array 'img'
! Fortran 90 version of Jan Pieter Abrahams's compression algorithm.
! W.Kabsch & K.Rohm 6-2002
!-------------------------------------------------------------------------
!Machine-dependent subroutine because the file was generated by C-program
!-------------------------------------------------------------------------
!Which iolength has a 1-byte, 4-byte, and 8-byte buffer on this machine?
!What happens to the last record when requesting more data than available?
!-------------------------------------------------------------------------
! f90 IRIX compiler on SGI            (computer=lego)
! iolength= 1  for a 1-byte record    can use 1-byte records
! iolength= 4  for a 4-byte record
! iolength= 8  for a 8-byte record
!  1   2   3   4   0   0   0   0      items in the last record are o.k.

! pgf90 compiler on INTEL PENTIUM III (computer=kai)
! iolength=            1  for a 1-byte record  can use 1-byte records
! iolength=            4  for a 4-byte record
! iolength=            8  for a 8-byte record
!  0   0   0   0   0   0   0   0      last record is all 0; not o.k. !!

! f90 compiler on alpha OSF/1         (computer=minos)
! iolength=           1  for a 1-byte record   cannot use 1-byte records!!
! iolength=           1  for a 4-byte record
! iolength=           2  for a 8-byte record
!  1   2   3   4   0   0   0   0      items in the last record are o.k.
!  0   0   0   0   0   0   0   0      -convert big_endian; not o.k. !!
!-------------------------------------------------------------------------

! read_pck is used in two modes:
! inquire mode:       if npixel=0, read_pck will only report x,y
! decompression mode: if npixel>=x*y, img contains the decompressed image
! lun       device unit number
! filnam    file name of the compressed image
! npixel    dimension of array img
! img       contains the decompressed image if npixel>=x*y
!  x        number of fast pixels in the image (returned from file header)
!  y        number of slow pixels in the image (returned from file header)
! ier       error flag
!            0  no error
!           -1  end-of-file
!           -2  cannot open or read file
!           -3  0<npixel<x*y  dimension of img too small

implicit none
intent(in)           :: lun,filnam,npixel
intent(out)          :: img,x,y,ier

integer              :: lun,ier,x,y
integer,parameter    :: bufsiz=1024
integer*1            :: buf(bufsiz),ibyte
integer*2            :: img(*)
integer*4            :: npixel,ibold,ib,ip,i,j,k,l,next,iword
logical              :: unbuffered,last
character(len=50)    :: filnam
character(len=35)    :: header,marker
character(len=1)     :: cbyte

integer              :: valids,spillbits,usedbits,total,bytrec
integer*4            :: setbits(0:32),window,spill,pixel,nextint,bitnum,pixnum
integer  ,parameter  :: bitdecode(0:7)=(/0, 4, 5, 6, 7, 8, 16, 32/)
data                    setbits                                           &
 /z'00000000',z'00000001',z'00000003',z'00000007',z'0000000F',z'0000001F',&
  z'0000003F',z'0000007F',z'000000FF',z'000001FF',z'000003FF',z'000007FF',&
  z'00000FFF',z'00001FFF',z'00003FFF',z'00007FFF',z'0000FFFF',z'0001FFFF',&
  z'0003FFFF',z'0007FFFF',z'000FFFFF',z'001FFFFF',z'003FFFFF',z'007FFFFF',&
  z'00FFFFFF',z'01FFFFFF',z'03FFFFFF',z'07FFFFFF',z'0FFFFFFF',z'1FFFFFFF',&
  z'3FFFFFFF',z'7FFFFFFF',z'FFFFFFFF'/

marker='CCP4 packed image, X:              '
x=0;y=0;ier=0;l=0;

!bytrec=iolength of a single-byte record for this computer
!bytrec=0 if the computer cannot use single-byte records
inquire(iolength=i)ibyte
inquire(iolength=j)iword
bytrec=0;if (i.lt.j)bytrec=i

last=.false.;unbuffered=.false.;ibold=-1;next=0;
inquire(iolength=i)buf
open(unit=lun,file=filnam,status='old',form='unformatted',                &
     access='direct',recl=i,action='read')

do
   if (l.eq.21)exit
   next=next+1
   if (unbuffered)go to 30
   ib=(next-1)/bufsiz
   ip=next-bufsiz*ib
   if (ib.ne.ibold)then
      if (last)go to 220
      ibold=ib
      read(lun,rec=1+ib,err=20)buf
   endif
10 ibyte=buf(ip)
   go to 40
20 close(lun)
   last=.true.
   if (bytrec.lt.1)go to 10
   open(unit=lun,file=filnam,status='old',form='unformatted',             &
	access='direct',recl=bytrec,action='read',err=230)
   unbuffered=.true.
30 read(lun,rec=next,err=220)ibyte
40 j=1+mod(next,35)
   cbyte=transfer(ibyte,cbyte)
   header(j:j)=cbyte    !note: equivalence  (cbyte,ibyte)
   if (next>34)then
      l=0
      do j=1,21
	 k=1+mod(next-35+j,35)
	 if (marker(j:j).ne.header(k:k))exit
	 l=j
      enddo
      if (l.eq.21)then
	 do j=22,35
	    k=1+mod(next-35+j,35)
	    marker(j:j)=header(k:k)
	 enddo
	 read(marker,'(21x,i5,4x,i5)')x,y
      endif
   endif
enddo
if (npixel<1)  go to 200
if (npixel<x*y)go to 240

next=next+1;
valids=0; spillbits = 0; total = x * y; window = 0; pixel = 0;

do
   if (pixel.ge.total)exit
   if (valids<6)then
      if (spillbits>0)then
	 window=ior(window,ishft(spill,valids))
	 valids=valids+spillbits
	 spillbits = 0
      else
	 next=next+1
	 if (unbuffered)go to 70
	 ib=(next-1)/bufsiz
	 ip=next-bufsiz*ib
	 if (ib.ne.ibold)then
	    if (last)go to 220
	    ibold=ib
	    read(lun,rec=1+ib,err=60)buf
	 endif
50       ibyte=buf(ip)
	 go to 80
60       close(lun)
	 last=.true.
	 if (bytrec.lt.1)go to 50
	 open(unit=lun,file=filnam,status='old',form='unformatted',       &
	      access='direct',recl=bytrec,action='read',err=230)
	 unbuffered=.true.
70       read(lun,rec=next,err=220)ibyte
80       spill=ibyte
	 if (spill.lt.0)spill=spill+256
	 spillbits = 8
      endif
   else
      pixnum = ishft(1,iand(window,setbits(3)))
      window = ishft(window,-3)
      bitnum = bitdecode(iand(window,setbits(3)))
      window = ishft(window,-3)
      valids = valids-6
      do
	 if ((pixnum.le.0).or.(pixel.ge.total))exit
	 if (valids < bitnum)then
	    if (spillbits > 0)then
	       window = ior(window,ishft(spill,valids))
	       if ((32 - valids) > spillbits)then
		  valids = valids+spillbits
		  spillbits = 0
	       else
		  usedbits = 32 - valids
		  spill = ishft(spill,-usedbits)
		  spillbits = spillbits-usedbits
		  valids = 32
	       endif
	    else
	       next=next+1
	       if (unbuffered)go to 110
	       ib=(next-1)/bufsiz
	       ip=next-bufsiz*ib
	       if (ib.ne.ibold)then
		  if (last)go to 220
		  ibold=ib
		  read(lun,rec=1+ib,err=100)buf
	       endif
90             ibyte=buf(ip)
	       go to 120
100            close(lun)
	       last=.true.
	       if (bytrec.lt.1)go to 90
	       open(unit=lun,file=filnam,status='old',form='unformatted', &
		    access='direct',recl=bytrec,action='read',err=230)
	       unbuffered=.true.
110            read(lun,rec=next,err=220)ibyte
120            spill=ibyte
	       if (spill.lt.0)spill=spill+256
	       spillbits = 8
	    endif
	 else
	    pixnum = pixnum-1
	    if (bitnum.eq.0)then
	       nextint = 0
	    else
	       nextint = iand(window,setbits(bitnum))
	       valids  = valids-bitnum
	       window  = ishft(window,-bitnum)
	       if (iand(nextint,ishft(1,bitnum - 1)).ne.0) &
		  nextint= ior(nextint,not(setbits(bitnum)))
	    endif
	    i=0
	    if (pixel>0)i=img(pixel)
	    if (pixel>x)i=(i+img(pixel-x+2)+img(pixel-x+1)+img(pixel-x)+2)/4
	    pixel = pixel+1
	    img(pixel)=nextint+i ! =iand(nextint+i,z'0000FFFF')
	 endif
      end do
   endif
end do
200 ier=0
210 close(lun)
    return
220 ier=-1
    go to 210
230 ier=-2
    return
240 ier=-3
    go to 210
end subroutine read_pck
subroutine write_pck(lun,filnam,append,img,x,y,ier)
! Append packed integer*2 img(x*y) to 'filnam'. pack_wordimage
! Fortran 90 version of Jan Pieter Abrahams's compression algorithm.
! ("write_pck" was originally named "pack_wordimage").
! W.Kabsch & K.Rohm 6-2002

! lun       device unit number
! filnam    file name of the compressed image
! append    if .true.  append the compressed image to file 'filnam'
!           if .false. replace file 'filnam' by  the compressed image
! img       contains the decompressed image if npixel>=x*y
!  x        number of fast pixels in the image (returned from file header)
!  y        number of slow pixels in the image (returned from file header)
! ier       error flag
!            0  no error
!           -1  cannot write file
!           -2  cannot open  file

implicit   none
!external          :: diff_words,pack_chunk,bits
intent(in)        :: lun,filnam,append,img,x,y
intent(out)       :: ier
integer           :: lun,x,y,ier
integer*2         :: img(*)
character(len=*)  :: filnam
logical           :: append

integer           :: chunksiz,packsiz,nbits,next_nbits,tot_nbits, &
		     idiff,last,i,bits,bitmark,target_offset,oldline
integer,parameter :: DIFFBUFSIZ=16384,PACKBUFSIZ=16384,IOBUFSIZ=4096
integer*4         :: diffbuf(0:DIFFBUFSIZ-1),done,oldbytes,ibold,next,ib,ip
integer*1         :: packbuf(0:PACKBUFSIZ-1),iobuf(IOBUFSIZ),iheader(37)
character(len=80) :: zeile
character(len=37) :: header
logical           :: ex,offen
integer*1,allocatable  :: oldfile(:)

bitmark=0;target_offset=0

inquire(file=filnam,exist=ex,opened=offen)
if (offen)go to 110  !file should not be open if it exists

oldbytes=0
if (ex.and.append)then
   open(unit=lun,file=filnam,status='old',form='formatted',    &
	access='sequential',err=120)
   rewind lun
   oldline=0
   do
      zeile=' '
      read(lun,'(A80)',err=10,end=10)zeile
      oldline=oldline+1
   enddo
10 if (oldline>0)then
      oldbytes=oldline*80
      allocate (oldfile(oldbytes))
      rewind lun
      next=0
      do i=1,oldline
	 read(lun,'(A80)')zeile
	 do ip=1,80
	    next=next+1
	    oldfile(next)=ichar(zeile(ip:ip))
	 enddo
      enddo
   endif
   close(unit=lun,status='delete')
endif

write(header,"(' CCP4 packed image, X: ',i4,', Y: ',i4,1X)")x,y
iheader(1)=10;iheader(37)=10 ! ASCII values for \n
do i=2,36
   iheader(i)=ichar(header(i:i))
enddo

inquire(iolength=i)iobuf
open(unit=lun,file=filnam,status='replace',form='unformatted',    &
     access='direct',recl=i,err=120)
next=0;ibold=0
if (oldbytes>0)then
   do i=1,oldbytes
      ib=next/IOBUFSIZ;next=next+1;ip=next-ib*IOBUFSIZ
      if (ib.ne.ibold)then
	 write(lun,rec=1+ibold,err=110)iobuf
	 ibold=ib
      endif
      iobuf(ip)=oldfile(i)
   enddo
   deallocate (oldfile)
endif

do i=1,37
   ib=next/IOBUFSIZ;next=next+1;ip=next-ib*IOBUFSIZ
   if (ib.ne.ibold)then
      write(lun,rec=1+ibold,err=110)iobuf
      ibold=ib
   endif
   iobuf(ip)=iheader(i)
enddo

done=1
do
   if (done>x*y)exit
   call diff_words(img,x,y,diffbuf,DIFFBUFSIZ,done,last)
   idiff=0
   do
      if (idiff>last)exit
      packsiz=0
      chunksiz=1
      nbits=bits(diffbuf(idiff),1)
      do
	 if (packsiz.ne.0)exit
	 if (last.le.idiff+2*chunksiz)then
	    packsiz=chunksiz
	 else
	    next_nbits=bits(diffbuf(idiff+chunksiz),chunksiz)
	    tot_nbits=2*max(nbits,next_nbits)
	    if (tot_nbits.ge.nbits+next_nbits+6)then
	       packsiz=chunksiz
	    else
	       nbits=tot_nbits
	       if (chunksiz.eq.64)then
		  packsiz=128
	       else
		  chunksiz=2*chunksiz
	       endif
	    endif
	 endif
      enddo
      call pack_chunk(diffbuf(idiff),packsiz,nbits/packsiz,PACKBUFSIZ,  &
	 packbuf,bitmark,target_offset,lun,iobuf,IOBUFSIZ,ibold,next,ier)
      if (ier.ne.0)go to 110
      idiff=idiff+packsiz
   enddo
enddo
call pack_chunk(diffbuf,0,0,PACKBUFSIZ,packbuf,bitmark,target_offset,   &
		lun,iobuf,IOBUFSIZ,ibold,next,ier)
if (ier.ne.0)go to 110
close(unit=lun)
ier=0
return
110 ier=-1
    return
120 ier=-2
    return
end subroutine write_pck
integer function bits(chunk,n)
! Returns the number of bits neccesary to encode array chunk(0:n-1)
! The size in bits of one encoded element can be 0,4,5,6,7,8,16 or 32.
implicit none
integer    :: n,i
integer*4  :: chunk(0:*),maxsize

maxsize = abs(chunk(0))
do i=1,n
   maxsize = max(maxsize, abs(chunk(i-1)))
enddo
if         (maxsize.eq. 0)then; bits = 0
   else if (maxsize<    8)then; bits = 4 * n
   else if (maxsize<   16)then; bits = 5 * n
   else if (maxsize<   32)then; bits = 6 * n
   else if (maxsize<   64)then; bits = 7 * n
   else if (maxsize<  128)then; bits = 8 * n
   else if (maxsize<65536)then; bits =16 * n
   else                       ; bits =32 * n
endif
return
end function bits
subroutine pack_chunk(lng,nmbr,bitsize,PACKBUFSIZ,packbuf,bitmark,       &
		      target_offset,lun,iobuf,IOBUFSIZ,ibold,next,ier)
! Packs 'nmbr' LONGs starting at 'lng[0]' into a packed array of 'bitsize'
! sized elements. If the internal packbuf in which the array is packed is
! full, it is flushed to 'file', making room for more of the packed array.
! If (nmbr <= 0), the packbuf is flushed aswell.
implicit none
!external             :: pack_longs
integer              :: PACKBUFSIZ,IOBUFSIZ
integer              :: nmbr,bitsize,bitmark,target_offset,lun,k,ier
integer*1            :: packbuf(0:PACKBUFSIZ-1),iobuf(IOBUFSIZ)
integer*4            :: lng(0:*),descriptor(0:1),ib,ibold,ip,next,i,j
integer*4, parameter :: bitsize_encode(0:32)=(/                       &
			0, 0, 0, 0, 1, 2, 3, 4, 5, 0, 0,              &
			0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0,              &
			0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7/)

ier=-1
if (nmbr>0)then
   i=nmbr;j=0
   do
      if (i<2)exit
      i=i/2;j=j+1
   enddo
   descriptor(0) = j
   descriptor(1) = bitsize_encode(bitsize)
   if (target_offset > (PACKBUFSIZ - (130 * 4))) then
      do k=0,target_offset-1
	 ib=next/IOBUFSIZ;next=next+1;ip=next-ib*IOBUFSIZ
	 if (ib.ne.ibold)then
	    write(lun,rec=1+ibold,err=100)iobuf
	    ibold=ib
	 endif
	 iobuf(ip)=packbuf(k)
      enddo
      packbuf(0) = packbuf(target_offset)
      target_offset=0
   endif
   call pack_longs(descriptor,2,packbuf,target_offset,bitmark,3)
   call pack_longs(lng,nmbr,packbuf,target_offset,bitmark,bitsize)
else
   do k=0,target_offset
      ib=next/IOBUFSIZ;next=next+1;ip=next-ib*IOBUFSIZ
      if (ib.ne.ibold)then
	 write(lun,rec=1+ibold,err=100)iobuf
	 ibold=ib
      endif
      iobuf(ip)=packbuf(k)
   enddo
   write(lun,rec=1+ib,err=100)iobuf
endif
ier=0
100 return
end subroutine pack_chunk
subroutine pack_longs(lng,n,target,target_offset,bit,size)
implicit none
integer             :: n,bit,size,valids,i,temp,temp_bit
integer*1           :: target(0:*)
integer*4           :: lng(0:*),target_offset,mask,window,setbits(0:32)
data                   setbits                                            &
 /z'00000000',z'00000001',z'00000003',z'00000007',z'0000000F',z'0000001F',&
  z'0000003F',z'0000007F',z'000000FF',z'000001FF',z'000003FF',z'000007FF',&
  z'00000FFF',z'00001FFF',z'00003FFF',z'00007FFF',z'0000FFFF',z'0001FFFF',&
  z'0003FFFF',z'0007FFFF',z'000FFFFF',z'001FFFFF',z'003FFFFF',z'007FFFFF',&
  z'00FFFFFF',z'01FFFFFF',z'03FFFFFF',z'07FFFFFF',z'0FFFFFFF',z'1FFFFFFF',&
  z'3FFFFFFF',z'7FFFFFFF',z'FFFFFFFF'/

temp_bit =bit;
if (size > 0)then
   mask = setbits(size)
   do i=0,n-1
      window =iand(lng(i),mask)
      valids = size
      if (temp_bit.eq.0)then
	 target(target_offset)=window
      else
	 temp=target(target_offset)
	 target(target_offset)=ior(temp,ishft(window,temp_bit))
      endif
      window =ishft(window,temp_bit-8)
      valids = valids-(8-temp_bit)
      if (valids<0)then
	 temp_bit=temp_bit+size
      else
	 do
	    if (valids.le.0)exit
	    target_offset=target_offset+1
	    target(target_offset)=window
	    window =ishft(window,-8)
	    valids = valids-8
	 enddo
	 temp_bit= 8+valids
      endif
      if (valids.eq.0)then
	 temp_bit=0
	 target_offset=target_offset+1
      endif
   enddo
   bit=mod((bit+(size*n)),8)
endif
return
end subroutine pack_longs
subroutine diff_words(word,x,y,diffs,DIFFBUFSIZ,done,last)
! Calculates the difference of WORD-sized pixels of an image with the
! truncated mean value of four of its neighbours. 'x' is the number of fast
! coordinates of the image 'img', 'y' is the number of slow coordinates,
! 'diffs' will contain the differences, 'done' defines the index of the pixel
! where calculating the differences should start. A pointer to the last
! difference is returned. Maximally DIFFBUFSIZ differences are returned in
! 'diffs'.
implicit      none
integer*2  :: word(*)
integer    :: x,y,DIFFBUFSIZ
integer*4  :: diffs(0:DIFFBUFSIZ-1),done,last,j,tot

last=0; tot=x*y
if (done.eq.1)then
   diffs(last)=word(done);done=done+1;last=last+1
endif
do
   if ((done>x+1).or.(last.ge.DIFFBUFSIZ))exit
   j=word(done);diffs(last)=j-word(done-1);
   done=done+1;last=last+1
enddo
do
   if ((done>tot).or.(last.ge.DIFFBUFSIZ))exit
   j=word(done-1)+2
   diffs(last)=word(done)-(j+word(done-x-1)+word(done-x)+word(done-x+1))/4
   done=done+1;last=last+1
enddo
last=last-1
return
end subroutine diff_words
