double precision function uni2ieee_dp (iarray,iwrd)

C Function converts Univac 72-bit floating point to IEEE 64-bit floating point

      integer*4     iarray(1)    !array of Univac numbers
      integer*4     iwrd         !word number of first half of Univac number
                                 ! where words are counted 0,1,2,... and a
                                 ! word is 36 bits (DP = 72 bits = 2 words)

      integer*4     jwrd         !32-bit word pointer in iarray (1,2,3,...)
      integer*4     jbit         !bit position in iarray (31,30,...0)

      integer*4     jsign        !sign bit for the number
      integer*4     jexp         !exponent for number
      integer*4     jman0        !more significant 30 mantissa bits
      integer*4     jman1        !less significant 30 mantissa bits

      integer*4     kexp         !temporary exponent for 2's complement
      integer*4     kman0        !temporary jman0 for 2's complement
      integer*4     kman1        !temporary jman1 for 2's complement

      integer*4     zero         !constant

      double precision  exp_fac  !factor corresponnding to Univac exponent
      double precision  mantissa !Univac mantissa in local binary
      double precision  sfac30   !constant

      integer*4     get_bits     !integer function


C Initializations

      jsign = 0
      jexp  = 0
      jman0 = 0
      jman1 = 0
      sfac30 = dfloat(2)**30
      zero  = 0


C Translate 36-bit pointer to 32-bit pointers

      jwrd = iwrd*36/32 + 1
      jbit = 31 - (iwrd*36 - 32*(jwrd-1))

c      write(*,'("init:  jwrd,jbit=           ",2i5)')
c     *     jwrd,jbit

C Extract components of Univac number

      jsign = get_bits(iarray,jwrd,jbit,1)
      jexp = get_bits(iarray,jwrd,jbit,11)
      jman0 = get_bits(iarray,jwrd,jbit,30)
      jman1 = get_bits(iarray,jwrd,jbit,30)

c      write(*,'("jsign: jwrd,jbit,jsign=     ",3i5)')
c     *     jwrd,jbit,jsign


c      write(*,'("jexp:  jwrd,jbit,jsign,jexp=",4i5)')
c     *     jwrd,jbit,jsign,jexp

c      write(*,'(4i10)') jsign,jexp,jman0,jman1


C Do two's complement for negative numbers

      if (jsign .eq. 1) then

        kexp = not(jexp)                     !One's complement
        call mvbits(zero,11,21,kexp,11)
        kman0 = not(jman0)
        call mvbits(zero,30,2,kman0,30)
        kman1 = not(jman1)
        call mvbits(zero,30,2,kman1,30)

        kman1 = kman1 + 1                    !Two's complement
        if (kman1 .eq. 2**30) then
          kman1 = 0
          kman0 = kman0 + 1
          if (kman0 .eq. 2**30) then
            kman0 = 0
            kexp = kexp + 1
          end if
        end if

        jexp = kexp
        jman0 = kman0
        jman1 = kman1

      end if


C Reconstruct double precision number in local binary

      mantissa = (dfloat(jman0) + dfloat(jman1)/sfac30)/sfac30
      exp_fac = dfloat(2)**(jexp-1024)
      uni2ieee_dp = exp_fac*mantissa


C Put on the correct sign

      if (jsign .eq. 1) uni2ieee_dp = - uni2ieee_dp

      return
      end