C=============================================================================C
C International Comprehensive Ocean-Atmosphere Data Set (ICOADS)  26 Apr 2010 C
C Filename:level: gsbytes.f:01B                            Fortran 77 program C
C Function: Bit-string packing/unpacking                     Author: S.Lubker C
C=============================================================================C
C Software Revision Information (original version):
C-----------------------------------------------------------------------3456789
C Software documentation for the (invariant) user-interface routines {gbytes,
C gbyte,sbytes,sbyte}:
C
C An implementation* of the Fortran subroutines designed by Jenne and
C Joseph (1974).  {gbytes} unpacks n b bit bytes from p into u, starting by
C skipping q bits in p, then skipping s bits between bytes.  {gbyte} unpacks
C one such byte.  {sbytes} packs n b bit bytes from u into p, starting by
C skipping q bits in p, then skipping s bits between bytes. {sbyte} packs one
C such byte.  Following are the constraints on q, b, s, and n:
C      q >= 0     number of bits to be skipped preceeding first byte in p
C  0 < b <= word  byte size
C      s >= 0     number of bits to be skipped between bytes
C      n >= 0     number of bytes to be packed/unpacked
C
C *endian independent
C
C Reference: Jenne, R.L., and D.H. Joseph, 1974: Techniques for the processing,
C storage, and exchange of data.  NCAR Tech. Note IA-93, National Center for
C Atmospheric Research, Boulder, Colo., 46 pp.
C
C Machine dependencies: Word size is 32 bits.
c-----------------------------------------------------------------------3456789
      subroutine gbytes(p,u,q,b,s,n)
      implicit integer(a-z)
      dimension p(*), u(n)
      q2 = q
      do i = 1, n
        call gbyte(p, u(i), q2, b)
        q2 = q2 + b + s
      enddo
      end
c-----------------------------------------------------------------------3456789
      subroutine gbyte(p,u,q,b)
      implicit integer(a-z)
      character p(*)
      dimension bits(3)
      u = 0
      q2 = q
      b2 = b
   10 byte = q2/8 + 1
      bits(1) = mod(q2, 8)
      bits(2) = min(b2, 8 - bits(1))
      bits(3) = 8 - bits(1) - bits(2)
      int4 = ichar(p(byte))
      call mvbits(int4, bits(3), bits(2), u, b2 - bits(2))
      q2 = q2 + bits(2)
      b2 = b2 - bits(2)
      if (b2.gt.0) goto 10
      end
c-----------------------------------------------------------------------3456789
      subroutine sbytes(p,u,q,b,s,n)
      implicit integer(a-z)
      dimension p(*), u(n)
      q2 = q
      do i = 1, n
        call sbyte(p, u(i), q2, b)
        q2 = q2 + b + s
      enddo
      end
c-----------------------------------------------------------------------3456789
      subroutine sbyte(p,u,q,b)
      implicit integer(a-z)
      character p(*)
      dimension bits(3)
      q2 = q
      b2 = b
   10 byte = q2/8 + 1
      bits(1) = mod(q2, 8)
      bits(2) = min(b2, 8 - bits(1))
      bits(3) = 8 - bits(1) - bits(2)
      int4 = ichar(p(byte))
      call mvbits(u, b2 - bits(2), bits(2), int4, bits(3))
      p(byte) = char(int4)
      q2 = q2 + bits(2)
      b2 = b2 - bits(2)
      if (b2.gt.0) goto 10
      end