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 ! Recent versions of the PGI compilers apparently still do not fully support
  36 ! the use of all intrinsics in parameter statements, though this is part of
  37 ! the F2003 standard.
  38 !      real, parameter :: two23=scale(1.0,23)
  39 !      real, parameter :: two126=scale(1.0,126)
  40       real :: two23
  41       real :: two126
  42
  43       two23=scale(1.0,23)
  44       two126=scale(1.0,126)
  45
  46       alog2=alog(2.0)
  47
  48       do j=1,num
  49         ieee=0
  50
  51         if (a(j).eq.0.) then
  52           ieee=0
  53           rieee(j)=transfer(ieee,rieee(j))
  54 !       write(6,fmt='(f20.10,5x,b32)') a,a
  55 !       write(6,fmt='(f20.10,5x,b32)') rieee,rieee
  56           cycle
  57         endif
  58
  59 !
  60 !  Set Sign bit (bit 31 - leftmost bit)
  61 !
  62         if (a(j).lt.0.0) then
  63           ieee=ibset(ieee,31)
  64           atemp=abs(a(j))
  65         else
  66           ieee=ibclr(ieee,31)
  67           atemp=a(j)
  68         endif
  69 !
  70 !  Determine exponent n with base 2
  71 !
  72         if ( atemp .ge. 1.0 ) then
  73            n = 0
  74            do while ( 2.0**(n+1) .le. atemp )
  75               n = n + 1
  76            enddo
  77         else
  78            n = -1
  79            do while ( 2.0**n .gt. atemp )
  80               n = n - 1
  81            enddo
  82         endif
  83 !        n=floor(alog(atemp)/alog2)
  84         !write(6,*) ' logstuff ',alog(atemp)/alog2
  85         !write(6,*) ' logstuffn ',n
  86         iexp=n+127
  87         if (n.gt.127) iexp=255     ! overflow
  88         if (n.lt.-127) iexp=0
  89         !      set exponent bits ( bits 30-23 )
  90         call mvbits(iexp,0,8,ieee,23)
  91 !
  92 !  Determine Mantissa
  93 !
  94         if (iexp.ne.255) then
  95           if (iexp.ne.0) then
  96             atemp=(atemp/(2.0**n))-1.0
  97           else
  98             atemp=atemp*two126
  99           endif
 100           imant=nint(atemp*two23)
 101         else
 102           imant=0
 103         endif
 104         !      set mantissa bits ( bits 22-0 )
 105         call mvbits(imant,0,23,ieee,0)
 106 !
 107 !  Transfer IEEE bit string to real variable
 108 !
 109         rieee(j)=transfer(ieee,rieee(j))
 110 !       write(6,fmt='(f20.10,5x,b32)') a,a
 111 !       write(6,fmt='(f20.10,5x,b32)') rieee,rieee
 112
 113       enddo
 114
 115       return
 116       end
 117