external/io_grib2/g2lib/rdieee.F

   1       subroutine rdieee(rieee,a,num)
   2 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
   3 !                .      .    .                                       .
   4 ! SUBPROGRAM:    rdieee
   5 !   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-09
   6 !
   7 ! ABSTRACT: This subroutine reads a list of real values in
   8 !   32-bit IEEE floating point format.
   9 !
  10 ! PROGRAM HISTORY LOG:
  11 ! 2000-05-09  Gilbert
  12 !
  13 ! USAGE:    CALL rdieee(rieee,a,num)
  14 !   INPUT ARGUMENT LIST:
  15 !     rieee    - Input array of floating point values in 32-bit IEEE format.
  16 !     num      - Number of floating point values to convert.
  17 !
  18 !   OUTPUT ARGUMENT LIST:
  19 !     a        - Output array of real values.
  20 !
  21 ! REMARKS: None
  22 !
  23 ! ATTRIBUTES:
  24 !   LANGUAGE: Fortran 90
  25 !   MACHINE:  IBM SP
  26 !
  27 !$$$
  28
  29       real(4),intent(in) :: rieee(num)
  30       real,intent(out) :: a(num)
  31       integer,intent(in) :: num
  32
  33       integer(4) :: ieee
  34
  35       real,save :: two23
  36       real,save :: two126
  37       integer,save :: once=0
  38
  39       if ( once .EQ. 0 ) then
  40          once=1
  41          two23=scale(1.0,-23)
  42          two126=scale(1.0,-126)
  43       endif
  44
  45       do j=1,num
  46 !
  47 !  Transfer IEEE bit string to integer variable
  48 !
  49         ieee=transfer(rieee(j),ieee)
  50 !
  51 !  Extract sign bit, exponent, and mantissa
  52 !
  53         isign=ibits(ieee,31,1)
  54         iexp=ibits(ieee,23,8)
  55         imant=ibits(ieee,0,23)
  56         sign=1.0
  57         if (isign.eq.1) sign=-1.0
  58
  59         if ( (iexp.gt.0).and.(iexp.lt.255) ) then
  60           temp=2.0**(iexp-127)
  61           a(j)=sign*temp*(1.0+(two23*real(imant)))
  62
  63         elseif ( iexp.eq.0 ) then
  64           if ( imant.ne.0 ) then
  65             a(j)=sign*two126*two23*real(imant)
  66           else
  67             a(j)=sign*0.0
  68           endif
  69
  70         elseif ( iexp.eq.255 ) then
  71           a(j)=sign*huge(a(j))
  72
  73         endif
  74
  75       enddo
  76
  77       return
  78       end
  79