external/io_grib2/g2lib/mkieee.F

   1       subroutine mkieee(a,rieee,num)
   2 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
   3 !                .      .    .                                       .
   4 ! SUBPROGRAM:    mkieee
   5 !   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-09
   6 !
   7 ! ABSTRACT: This subroutine stores 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 mkieee(a,rieee,num)
  14 !   INPUT ARGUMENT LIST:
  15 !     a        - Input array of floating point values.
  16 !     num      - Number of floating point values to convert.
  17 !
  18 !   OUTPUT ARGUMENT LIST:
  19 !     rieee    - Output array of floating point values in 32-bit IEEE format.
  20 !
  21 ! REMARKS: None
  22 !
  23 ! ATTRIBUTES:
  24 !   LANGUAGE: Fortran 90
  25 !   MACHINE:  IBM SP
  26 !
  27 !$$$
  28
  29       real,intent(in) :: a(num)
  30       real(4),intent(out) :: rieee(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       alog2=alog(2.0)
  46
  47       do j=1,num
  48         ieee=0
  49
  50         if (a(j).eq.0.) then
  51           ieee=0
  52           rieee(j)=transfer(ieee,rieee(j))
  53 !       write(6,fmt='(f20.10,5x,b32)') a,a
  54 !       write(6,fmt='(f20.10,5x,b32)') rieee,rieee
  55           cycle
  56         endif
  57
  58 !
  59 !  Set Sign bit (bit 31 - leftmost bit)
  60 !
  61         if (a(j).lt.0.0) then
  62           ieee=ibset(ieee,31)
  63           atemp=abs(a(j))
  64         else
  65           ieee=ibclr(ieee,31)
  66           atemp=a(j)
  67         endif
  68 !
  69 !  Determine exponent n with base 2
  70 !
  71         n=floor(alog(atemp)/alog2)
  72         iexp=n+127
  73         if (n.gt.127) iexp=255     ! overflow
  74         if (n.lt.-127) iexp=0
  75         !      set exponent bits ( bits 30-23 )
  76         call mvbits(iexp,0,8,ieee,23)
  77 !
  78 !  Determine Mantissa
  79 !
  80         if (iexp.ne.255) then
  81           if (iexp.ne.0) then
  82             atemp=(atemp/(2.0**n))-1.0
  83           else
  84             atemp=atemp*two126
  85           endif
  86           imant=nint(atemp*two23)
  87         else
  88           imant=0
  89         endif
  90         !      set mantissa bits ( bits 22-0 )
  91         call mvbits(imant,0,23,ieee,0)
  92 !
  93 !  Transfer IEEE bit string to real variable
  94 !
  95         rieee(j)=transfer(ieee,rieee(j))
  96 !       write(6,fmt='(f20.10,5x,b32)') a,a
  97 !       write(6,fmt='(f20.10,5x,b32)') rieee,rieee
  98
  99       enddo
 100
 101       return
 102       end
 103