ungrib/src/ngl/g2/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(4),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         if ( atemp .ge. 1.0 ) then
  72            n = 0
  73            do while ( 2.0**(n+1) .le. atemp )
  74               n = n + 1
  75            enddo
  76         else
  77            n = -1
  78            do while ( 2.0**n .gt. atemp )
  79               n = n - 1
  80            enddo
  81         endif
  82 !        n=floor(alog(atemp)/alog2)
  83         !write(6,*) ' logstuff ',alog(atemp)/alog2
  84         !write(6,*) ' logstuffn ',n
  85         iexp=n+127
  86         if (n.gt.127) iexp=255     ! overflow
  87         if (n.lt.-127) iexp=0
  88         !      set exponent bits ( bits 30-23 )
  89         call mvbits(iexp,0,8,ieee,23)
  90 !
  91 !  Determine Mantissa
  92 !
  93         if (iexp.ne.255) then
  94           if (iexp.ne.0) then
  95             atemp=(atemp/(2.0**n))-1.0
  96           else
  97             atemp=atemp*two126
  98           endif
  99           imant=nint(atemp*two23)
 100         else
 101           imant=0
 102         endif
 103         !      set mantissa bits ( bits 22-0 )
 104         call mvbits(imant,0,23,ieee,0)
 105 !
 106 !  Transfer IEEE bit string to real variable
 107 !
 108         rieee(j)=transfer(ieee,rieee(j))
 109 !       write(6,fmt='(f20.10,5x,b32)') a,a
 110 !       write(6,fmt='(f20.10,5x,b32)') rieee,rieee
 111
 112       enddo
 113
 114       return
 115       end
 116