module defconst
    implicit none
    real(8),parameter::pi=3.14159265358979323846D0
end module defconst

program f90_RAND
    use defconst
    implicit none
    integer::nd=100000
    real(8),allocatable::datax(:),datay(:)
    integer::i
    real(8)::z1,z2
    character::linebuf*100
    integer :: clock             ! for random_seed()
    integer :: nrand             ! for random_seed()
    integer,allocatable::seed(:) ! for random_seed()

    call random_seed(size=nrand)
    allocate(seed(nrand))
    call system_clock(count=clock)
    seed = clock
    call random_seed(put = seed)

    allocate(datax(1:nd),datay(1:nd))

    call random_number(datax)
    i=0
    do
        call BMM(z1,z2)
        i=i+1;datay(i)=z1
        if(i==nd)exit
        i=i+1;datay(i)=z2
        if(i==nd)exit
    end do

    open(12,file='inp_RU.txt',status='replace')
    write (12,'(a)') trim('random_number')
    do i=1,nd
        write(linebuf,*) datax(i)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    close(12)

    open(12,file='inp_RN.txt',status='replace')
    write (12,'(a)') trim('Box-Muller')
    do i=1,nd
        write(linebuf,*) datay(i)
        call del_spaces(linebuf)
        write (12,'(a)') trim(linebuf)
    end do
    close(12)

stop

contains

    subroutine BMM(z1,z2)
        !Box-Mullwe's method
        real(8),intent(out)::z1,z2
        real(8)::x,y

        call random_number(x)
        call random_number(y)
        z1=sqrt(-2.0D0*log(x))*cos(2.0D0*pi*y)
        z2=sqrt(-2.0D0*log(x))*sin(2.0D0*pi*y)
    end subroutine

    subroutine del_spaces(s)
        character (*), intent (inout) :: s
        character (len=len(s)) tmp
        integer i, j
        j = 1
        do i = 1, len(s)
            if (s(i:i)==' ') cycle
            tmp(j:j) = s(i:i)
            j = j + 1
        end do
        s = tmp(1:j-1)
    end subroutine del_spaces

end program f90_RAND