grnd Function

function grnd()
  real(8)            :: grnd    
  integer,parameter :: M = 397, MATA  = -1727483681 ! constant vector a
  integer,parameter :: LMASK =  2147483647          ! least significant r bits
  integer,parameter :: UMASK = -LMASK - 1           ! most significant w-r bits
  integer,parameter :: TMASKB= -1658038656, TMASKC= -272236544    
  integer,save       :: mag01(0:1)=[0,MATA] ! mag01(x) = x * MATA for x=0,1
  integer            :: kk,y
  ! These statement functions have been replaced with separate functions
  ! tshftu(y) = ISHFT(y,-11)
  ! tshfts(y) = ISHFT(y,7)
  ! tshftt(y) = ISHFT(y,15)
  ! tshftl(y) = ISHFT(y,-18)
  if (mti>=N) then               ! generate N words at one time
     if (mti==N+1) then          ! if sgrnd() has not been called,
        call sgrnd( defaultsd ) ! a default initial seed is used
     endif
     do kk=0,N-M-1
        y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
        mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1)))
     enddo
     do kk=N-M,N-2
        y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
        mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1)))
     enddo
     y=ior(iand(mt(N-1),UMASK),iand(mt(0),LMASK))
     mt(N-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1)))
     mti = 0
  endif
  y=mt(mti)
  mti = mti + 1 
  y=ieor(y,TSHFTU(y))
  y=ieor(y,iand(TSHFTS(y),TMASKB))
  y=ieor(y,iand(TSHFTT(y),TMASKC))
  y=ieor(y,TSHFTL(y))    
  if (y<0) then
     grnd=(dble(y)+2d0**32)/(2d0**32-1d0)
  else
     grnd=dble(y)/(2d0**32-1d0)
  endif
  return
contains
  integer function TSHFTU(y)
    integer,intent(in) :: y
    TSHFTU=ishft(y,-11)
    return
  end function TSHFTU
  integer function TSHFTS(y)
    integer,intent(in) :: y
    TSHFTS=ishft(y,7)
    return
  end function TSHFTS
  integer function TSHFTT(y)
    integer,intent(in) :: y
    TSHFTT=ishft(y,15)
    return
  end function TSHFTT
  integer function TSHFTL(y)
    integer,intent(in) :: y
    TSHFTL=ishft(y,-18)
    return
  end function TSHFTL
end function grnd