g2bytes.F90

Go to the documentation of this file.

 
 
subroutine g2_gbytec(in, iout, iskip, nbits)
  implicit none
 
  character*1, intent(in) :: in(*)
  integer, intent(inout) :: iout(*)
  integer, intent(in) :: iskip, nbits
 
  interface
     subroutine g2_gbytesc(in, iout, iskip, nbits, nskip, n)
       character*1, intent(in) :: in(*)
       integer, intent(out) :: iout(*)
       integer, intent(in) :: iskip, nbits, nskip, n
     end subroutine g2_gbytesc
  end interface
     
  call g2_gbytesc(in, iout, iskip, nbits, 0, 1)
end subroutine g2_gbytec
 
subroutine g2_gbytec1(in, siout, iskip, nbits)
  implicit none
 
  character*1, intent(in) :: in(*)
  integer, intent(inout) :: siout
  integer, intent(in) :: iskip, nbits
 
  integer (kind = 4) :: iout(1)
 
  interface
     subroutine g2_gbytesc(in, iout, iskip, nbits, nskip, n)
       character*1, intent(in) :: in(*)
       integer, intent(out) :: iout(*)
       integer, intent(in) :: iskip, nbits, nskip, n
     end subroutine g2_gbytesc
  end interface
  
  call g2_gbytesc(in, iout, iskip, nbits, 0, 1)
  siout = iout(1)
end subroutine g2_gbytec1
 
!iteration.
subroutine g2_gbytescr(in, rout, iskip, nbits, nskip, n)
  implicit none
  character*1, intent(in) :: in(*)
  real (kind = 4), intent(out) :: rout(*)
  integer, intent(in) :: iskip, nbits, nskip, n
  integer (kind = 4) :: iout(n)
 
  interface
     subroutine g2_gbytesc(in, iout, iskip, nbits, nskip, n)
       character*1, intent(in) :: in(*)
       integer, intent(out) :: iout(*)
       integer, intent(in) :: iskip, nbits, nskip, n
     end subroutine g2_gbytesc
  end interface
 
  ! Unpack into integer array.
  call g2_gbytesc(in, iout, iskip, nbits, nskip, n)
 
  ! Transfer to real array.
  rout(1:n) = transfer(iout, rout(1:n), n)
end subroutine g2_gbytescr
 
subroutine g2_gbytesc(in, iout, iskip, nbits, nskip, n)
  implicit none
 
  character*1, intent(in) :: in(*)
  integer, intent(out) :: iout(*)
  integer, intent(in) :: iskip, nbits, nskip, n
  integer :: tbit, bitcnt
  integer, parameter :: ones(8) = (/ 1, 3, 7, 15, 31, 63, 127, 255 /)
 
  integer :: nbit, i, index, ibit, itmp
  integer, external :: g2_mova2i
 
  !     nbit is the start position of the field in bits
  nbit = iskip
  do i = 1, n
     bitcnt = nbits
     index = nbit / 8 + 1
     ibit = mod(nbit, 8)
     nbit = nbit + nbits + nskip
 
     !        first byte
     tbit = min(bitcnt, 8 - ibit)
     itmp = iand(g2_mova2i(in(index)), ones(8 - ibit))
     if (tbit .ne. 8 - ibit) itmp = ishft(itmp, tbit - 8 + ibit)
     index = index + 1
     bitcnt = bitcnt - tbit
 
     !        now transfer whole bytes
     do while (bitcnt .ge. 8)
        itmp = ior(ishft(itmp,8), g2_mova2i(in(index)))
        bitcnt = bitcnt - 8
        index = index + 1
     enddo
 
     !        get data from last byte
     if (bitcnt .gt. 0) then
        itmp = ior(ishft(itmp, bitcnt), iand(ishft(g2_mova2i(in(index)), &
             - (8 - bitcnt)), ones(bitcnt)))
     endif
 
     iout(i) = itmp
  enddo
 
end subroutine g2_gbytesc
 
subroutine g2_gbytec8(in, iout, iskip, nbits)
  implicit none
 
  character*1, intent(in) :: in(*)
  integer (kind = 8), intent(inout) :: iout(*)
  integer, intent(in) :: iskip, nbits
 
  interface
     subroutine g2_gbytesc8(in, iout, iskip, nbits, nskip, n)
       character*1, intent(in) :: in(*)
       integer (kind = 8), intent(out) :: iout(*)
       integer, intent(in) :: iskip, nbits, nskip, n
     end subroutine g2_gbytesc8
  end interface
     
  call g2_gbytesc8(in, iout, iskip, nbits, 0, 1)
end subroutine g2_gbytec8
 
subroutine g2_gbytec81(in, siout, iskip, nbits)
  implicit none
 
  character*1, intent(in) :: in(*)
  integer (kind = 8), intent(inout) :: siout
  integer, intent(in) :: iskip, nbits
  integer (kind = 8) :: iout(1)
 
  interface
     subroutine g2_gbytesc8(in, iout, iskip, nbits, nskip, n)
       character*1, intent(in) :: in(*)
       integer (kind = 8), intent(out) :: iout(*)
       integer, intent(in) :: iskip, nbits, nskip, n
     end subroutine g2_gbytesc8
  end interface
  
  call g2_gbytesc8(in, iout, iskip, nbits, 0, 1)
  siout = iout(1)
end subroutine g2_gbytec81
 
subroutine g2_gbytesc8(in, iout, iskip, nbits, nskip, n)
  implicit none
 
  character*1, intent(in) :: in(*)
  integer (kind = 8), intent(out) :: iout(*)
  integer, intent(in) :: iskip, nbits, nskip, n
  
  integer :: tbit, bitcnt
  integer, parameter :: ones(8) = (/ 1, 3, 7, 15, 31, 63, 127, 255 /)
 
  integer :: nbit, i, index, ibit, itmp
  integer (kind = 8) :: itmp8, itmp8_2, itmp8_3
  integer, external :: g2_mova2i
  integer (kind = 8), external :: g2_mova2i8
 
  ! Nbit is the start position of the field in bits.
  nbit = iskip
  do i = 1, n
     bitcnt = nbits
     index = nbit / 8 + 1
     ibit = mod(nbit, 8)
     nbit = nbit + nbits + nskip
 
     ! first byte
     tbit = min(bitcnt, 8 - ibit)
     itmp8 = iand(g2_mova2i8(in(index)), int(ones(8 - ibit), kind = 8))
     itmp = iand(g2_mova2i(in(index)), ones(8 - ibit))
     if (tbit .ne. 8 - ibit) itmp = ishft(itmp, tbit - 8 + ibit)
     if (tbit .ne. 8 - ibit) itmp8 = ishft(itmp8, tbit - 8 + ibit)
     index = index + 1
     bitcnt = bitcnt - tbit
 
     ! now transfer whole bytes
     do while (bitcnt .ge. 8)
        itmp = ior(ishft(itmp,8), g2_mova2i(in(index)))
        itmp8 = ior(ishft(itmp8,8), g2_mova2i8(in(index)))
        bitcnt = bitcnt - 8
        index = index + 1
     enddo
 
     ! get data from last byte
     if (bitcnt .gt. 0) then
        itmp = ior(ishft(itmp, bitcnt), iand(ishft(g2_mova2i(in(index)), - (8 - bitcnt)), ones(bitcnt)))
        itmp8_2 = ishft(g2_mova2i8(in(index)), int(-(8 - bitcnt), kind(8)))
        itmp8_3 = int(ones(bitcnt), kind(8))
        itmp8 = ior(ishft(itmp8, bitcnt), iand(itmp8_2, itmp8_3))
     endif
 
     iout(i) = itmp8
  enddo
end subroutine g2_gbytesc8
 
subroutine g2_sbytec(out, in, iskip, nbits)
  implicit none
 
  character*1, intent(inout) :: out(*)
  integer, intent(in) :: in(*)
  integer, intent(in) :: iskip, nbits
 
  interface
     subroutine g2_sbytesc(out, in, iskip, nbits, nskip, n)
       character*1, intent(out) :: out(*)
       integer, intent(in) :: in(n)
       integer, intent(in) :: iskip, nbits, nskip, n
     end subroutine g2_sbytesc
  end interface
     
  call g2_sbytesc(out, in, iskip, nbits, 0, 1)
end subroutine g2_sbytec
 
subroutine g2_sbytec1(out, in, iskip, nbits)
  implicit none
 
  character*1, intent(inout) :: out(*)
  integer, intent(in) :: in
  integer, intent(in) :: iskip, nbits
  integer :: ain(1)
 
  interface
     subroutine g2_sbytesc(out, in, iskip, nbits, nskip, n)
       character*1, intent(out) :: out(*)
       integer, intent(in) :: in(n)
       integer, intent(in) :: iskip, nbits, nskip, n
     end subroutine g2_sbytesc
  end interface
 
  ain(1) = in
  call g2_sbytesc(out, ain, iskip, nbits, 0, 1)
end subroutine g2_sbytec1
 
subroutine g2_sbytescr(out, rin, iskip, nbits, nskip, n)
  implicit none
  character*1, intent(out) :: out(*)
  real, intent(in) :: rin(n)
  integer, intent(in) :: iskip, nbits, nskip, n
  integer :: in(n)
 
  interface
     subroutine g2_sbytesc(out, in, iskip, nbits, nskip, n)
       character*1, intent(out) :: out(*)
       integer, intent(in) :: in(n)
       integer, intent(in) :: iskip, nbits, nskip, n
     end subroutine g2_sbytesc
  end interface
 
  ! Transfer real array to integer array.
  in(1:n) = transfer(rin, in(1:n), n)
 
  ! Pack into character array.
  call g2_sbytesc(out, in, iskip, nbits, nskip, n)
end subroutine g2_sbytescr
 
subroutine g2_sbytesc(out, in, iskip, nbits, nskip, n)
  implicit none
 
  character*1, intent(out) :: out(*)
  integer, intent(in) :: in(n)
  integer, intent(in) :: iskip, nbits, nskip, n
  
  integer :: bitcnt, tbit
  integer, parameter :: ones(8)=(/ 1,  3,  7, 15, 31, 63, 127, 255/)
  integer :: nbit, i, itmp, index, ibit, imask, itmp2, itmp3
  integer, external :: g2_mova2i
 
  ! number bits from zero to ...
  ! nbit is the last bit of the field to be filled
  nbit = iskip + nbits - 1
  !print *, 'nbit', nbit, 'nbits ', nbits, 'nskip', nskip, 'n', n
  do i = 1, n
     itmp = in(i)
     bitcnt = nbits
     index = nbit / 8 + 1
     ibit = mod(nbit, 8)
     nbit = nbit + nbits + nskip
     !print *, 'i', i, 'itmp', itmp, 'bitcnt', bitcnt, 'index', index, 'ibit', ibit, 'nbit', nbit
     
     ! make byte aligned
     if (ibit .ne. 7) then
        tbit = min(bitcnt, ibit + 1)
        imask = ishft(ones(tbit), 7 - ibit)
        itmp2 = iand(ishft(itmp, 7 - ibit),imask)
        itmp3 = iand(g2_mova2i(out(index)), 255 - imask)
        out(index) = char(ior(itmp2, itmp3))
        bitcnt = bitcnt - tbit
        itmp = ishft(itmp, -tbit)
        index = index - 1
     endif
 
     ! now byte aligned
 
     ! do by bytes
     do while (bitcnt .ge. 8)
        out(index) = char(iand(itmp, 255))
        !print '(z2.2, x, z2.2, x, z2.2)', out(index), itmp, iand(itmp, 255)
        itmp = ishft(itmp, -8)
        bitcnt = bitcnt - 8
        index = index - 1
     enddo
 
     ! Do left over bits.
     if (bitcnt .gt. 0) then
        itmp2 = iand(itmp, ones(bitcnt))
        !print '(z2.2, x, z2.2)', ones(bitcnt), itmp2
        itmp3 = iand(g2_mova2i(out(index)), 255 - ones(bitcnt))
        out(index) = char(ior(itmp2, itmp3))
     endif
  enddo
end subroutine g2_sbytesc
 
subroutine g2_sbytec8(out, in, iskip, nbits)
  implicit none
 
  character*1, intent(inout) :: out(*)
  integer (kind = 8), intent(in) :: in(*)
  integer, intent(in) :: iskip, nbits
 
  interface
     subroutine g2_sbytesc8(out, in, iskip, nbits, nskip, n)
       character*1, intent(out) :: out(*)
       integer (kind = 8), intent(in) :: in(n)
       integer, intent(in) :: iskip, nbits, nskip, n
     end subroutine g2_sbytesc8
  end interface
  
  call g2_sbytesc8(out, in, iskip, nbits, 0, 1)
end subroutine g2_sbytec8
 
subroutine g2_sbytec81(out, sin, iskip, nbits)
  implicit none
 
  character*1, intent(inout) :: out(*)
  integer (kind = 8), intent(in) :: sin
  integer, intent(in) :: iskip, nbits
 
  integer (kind = 8) :: in(1)
  
  interface
     subroutine g2_sbytesc8(out, in, iskip, nbits, nskip, n)
       character*1, intent(out) :: out(*)
       integer (kind = 8), intent(in) :: in(n)
       integer, intent(in) :: iskip, nbits, nskip, n
     end subroutine g2_sbytesc8
  end interface
  
  in(1) = sin
  call g2_sbytesc8(out, in, iskip, nbits, 0, 1)
end subroutine g2_sbytec81
 
subroutine g2_sbytesc8(out, in, iskip, nbits, nskip, n)
  implicit none
 
  character*1, intent(out) :: out(*)
  integer (kind = 8), intent(in) :: in(n)
  integer, intent(in) :: iskip, nbits, nskip, n
  
  integer :: bitcnt, tbit
  integer, parameter :: ones(8)=(/ 1,  3,  7, 15, 31, 63, 127, 255/)
  integer :: nbit, i, index, ibit, imask, itmp1, itmp2, itmp3
  integer (kind = 8) :: itmp8, itmp8_2
  integer, external :: g2_mova2i
 
  ! number bits from zero to ...
  ! nbit is the last bit of the field to be filled
  nbit = iskip + nbits - 1
  !print *, 'nbit', nbit, 'nbits ', nbits, 'nskip', nskip, 'n', n
  do i = 1, n
     itmp8 = in(i)
     bitcnt = nbits
     index = nbit / 8 + 1
     ibit = mod(nbit, 8)
     nbit = nbit + nbits + nskip
     !print *, 'i', i, 'itmp8', itmp8, 'bitcnt', bitcnt, 'index', index, 'ibit', ibit, 'nbit', nbit 
 
     ! make byte aligned
     if (ibit .ne. 7) then
        tbit = min(bitcnt, ibit + 1)
        imask = ishft(ones(tbit), 7 - ibit)
        itmp1 = int(ishft(itmp8, int(7 - ibit, kind(8))), kind(4))
        itmp2 = iand(itmp1, imask)
        itmp3 = iand(g2_mova2i(out(index)), 255 - imask)
        out(index) = char(ior(itmp2, itmp3))
        bitcnt = bitcnt - tbit
        itmp8 = ishft(itmp8, -tbit)
        index = index - 1
     endif
 
     ! now byte aligned
 
     ! Process a byte at a time.
     do while (bitcnt .ge. 8)
        !print *, bitcnt, iand(itmp8, 255_8)
        out(index) = char(iand(itmp8, 255_8))
        itmp8 = ishft(itmp8, -8)
        bitcnt = bitcnt - 8
        index = index - 1
     enddo
 
     ! Take care of left over bits.
     if (bitcnt .gt. 0) then
        itmp8_2 = int(ones(bitcnt), kind(8))
        itmp2 = int(iand(itmp8, itmp8_2), kind(4))
        itmp3 = iand(g2_mova2i(out(index)), 255 - ones(bitcnt))
        out(index) = char(ior(itmp2, itmp3))
     endif
  enddo
end subroutine g2_sbytesc8
 
subroutine rdieeec(cieee, a, num)
  implicit none
  
  character(len = 1), intent(in) :: cieee(*)
  real, intent(out) :: a(num)
  integer, intent(in) :: num
  real (kind = 4) :: rieee(num)
 
  interface
     subroutine rdieee(rieee, a, num)
       real(4), intent(in) :: rieee(num)
       real, intent(out) :: a(num)
       integer, intent(in) :: num
     end subroutine rdieee
  end interface
 
  rieee(1:num) = transfer(cieee(1:num * 4), rieee, num)
  call rdieee(rieee, a, num)
end subroutine rdieeec
 
subroutine rdieee(rieee, a, num)
  implicit none
  
  real(4), intent(in) :: rieee(num)
  real, intent(out) :: a(num)
  integer, intent(in) :: num
 
  integer(4) :: ieee
  real, parameter :: two23 = scale(1.0, -23)
  real, parameter :: two126 = scale(1.0, -126)
  integer :: iexp, imant, isign, j
  real :: sign, temp
 
  do j = 1, num
     ! Transfer IEEE bit string to integer variable.
     ieee = transfer(rieee(j), ieee)
 
     ! Extract sign bit, exponent, and mantissa.
     isign = ibits(ieee, 31, 1)
     iexp = ibits(ieee, 23, 8)
     imant = ibits(ieee, 0, 23)
     sign = 1.0
     if (isign .eq. 1) sign = -1.0
 
     if (iexp .gt. 0 .and. iexp .lt. 255) then
        temp = 2.0**(iexp - 127)
        a(j) = sign * temp * (1.0 + (two23 * real(imant)))
     elseif (iexp .eq. 0) then
        if (imant .ne. 0) then
           a(j) = sign * two126 * two23 * real(imant)
        else
           a(j) = sign * 0.0
        endif
     elseif (iexp .eq. 255) then
        a(j) = sign * huge(a(j))
     endif
  enddo
end subroutine rdieee
 
subroutine mkieee(a, rieee, num)
  implicit none
  
  real(4), intent(in) :: a(num)
  real(4), intent(out) :: rieee(num)
  integer, intent(in) :: num
 
  integer(4) :: ieee
  real, parameter :: two23 = scale(1.0,23)
  real, parameter :: two126 = scale(1.0,126)
  real :: alog2, atemp
  integer :: iexp, imant, j, n
 
  alog2 = alog(2.0)
 
  do j = 1, num
     ieee = 0
     if (a(j) .eq. 0.) then
        ieee = 0
        rieee(j) = transfer(ieee, rieee(j))
        cycle
     endif
 
     ! Set Sign bit (bit 31 - leftmost bit).
     if (a(j) .lt. 0.0) then
        ieee = ibset(ieee, 31)
        atemp = abs(a(j))
     else
        ieee = ibclr(ieee, 31)
        atemp = a(j)
     endif
 
     ! Determine exponent n with base 2.
     if (atemp .ge. 1.0) then
        n = 0
        do while (2.0**(n+1) .le. atemp)
           n = n + 1
        enddo
     else
        n = -1
        do while (2.0**n .gt. atemp )
           n = n - 1
        enddo
     endif
     iexp = n + 127
     if (n .gt. 127) iexp = 255 ! overflow
     if (n .lt. -127) iexp = 0
     call mvbits(iexp, 0, 8, ieee, 23)
 
     ! Determine Mantissa.
     if (iexp .ne. 255) then
        if (iexp .ne. 0) then
           atemp = (atemp / (2.0**n)) - 1.0
        else
           atemp = atemp * two126
        endif
        imant = nint(atemp * two23)
     else
        imant = 0
     endif
     ! set mantissa bits (bits 22-0).
     call mvbits(imant, 0, 23, ieee, 0)
 
     ! Transfer IEEE bit string to real variable.
     rieee(j) = transfer(ieee, rieee(j))
  enddo
end subroutine mkieee