maths.F90

Go to the documentation of this file.

module maths_mod
  use datadefn_mod, only : default_precision
  implicit none
 
#ifndef TEST_MODE
  private
#endif
 
  integer, parameter :: k4b=selected_int_kind(9)
  integer(k4b), parameter :: ia=16807,im=2147483647
  integer(k4b), parameter :: iq=127773,ir=2836
 
  public random
contains
 
  real(kind=default_precision) function random(idum)    
    integer(k4b),intent(inout) :: idum
        
    real,save :: am
    integer(k4b), save :: ix=-1,iy=-1,k
 
    if (idum <=0 .or. iy < 0) then
      am = nearest(1.0,-1.0)/im
      iy=ior(ieor(888889999,abs(idum)),1)
      ix=ieor(777755555,abs(idum))
      idum=abs(idum)+1
    end if
    ix = ieor(ix,ishft(ix,13))
    ix = ieor(ix,ishft(ix,-17))
    ix = ieor(ix,ishft(ix,5))
    k=iy/iq
    iy=ia*(iy-k*iq)-ir*k
    if(iy < 0) iy = iy + im
    random = am*ior(iand(im,ieor(ix,iy)),1)
  end function random
end module maths_mod