[数值问题] fortran随机数调用问题

module phase
  use imslf90
  implicit none
  real*8,parameter,private:: pi=3.141592653589793238462643d0
  integer,parameter,private:: N=233, NN=N*2, M=1, tran=2000, dot=2000, npara=101,D=0.001d0
  !integer,parameter,private:: mxparm=50, N=30, NN=N*2, M=1, tran=200, dot=200, npara=101 !Fig. 4
  real*8, parameter,private:: gamma=0.7d0, a=1.d0, b=144.d0/233.d0, beta=89.d0/144.d0, L=144.d0, Kappa=3.d0*pi
  !real*8, parameter,private:: gamma=0.7d0, a=1.d0, b=1.d0, beta=24.d0/30.d0, L=30.d0, Kappa=3.d0*pi !Fig. 4
  real*8, private:: changeF
  private derivs
  public plot 
  contains

  function ran()   !returns random number between 0 - 1  
    implicit none 
    integer , save :: flag = 0
    double precision :: ran 
    if(flag==0) then 
      call random_seed()
      flag = 1 
    endif 
    call random_number(ran)     ! built in fortran 90 random number function  
  end function ran
  
  function normal(mean,sigma) 
    implicit none 
    integer :: flag 
    double precision, parameter :: pi = 3.141592653589793239  
    double precision :: uuu1, uuu2, yyy1, yyy2, normal, mean, sigma 
    save flag 
    data flag /0/ 
    uuu1 = ran(); uuu2 = ran() 
    if (flag.eq.0) then 
      yyy1 = sqrt(-2.0d0*log(uuu1))*cos(2.0d0*pi*uuu2) 
      normal = mean + sigma*yyy1 
      flag = 1 
    else 
      yyy2 = sqrt(-2.0d0*log(uuu1))*sin(2.0d0*pi*uuu2) 
      normal = mean + sigma*yyy2 
      flag = 0 
    endif  
  end function normal
 
 subroutine plot()
  integer i,j,k
  real*8 y(NN),dydt(NN)
  real*8 hh,th,h,t
  real*8 dym(NN),dyt(NN),yt(NN),c(N)
  real*8 y0(NN)
  real*8 parameters(npara),averageVelocity(npara),parameters2(npara),averageVelocity2(npara),temp
  character*38 filenm
  ! set initial conditions
  t = 0.0d0
  !call random_seed()
  do i=1,N
    !call random_number(y(i))
    y0(i)  =(i-1)*b
	y0(i+N)=0.d0
  enddo


!绝热增加F
  do k=1,npara
  y=y0
  changeF= 0.01d0*(k-1)
  t=0.d0; h=0.001d0
  hh=h/2.d0
  temp=0.d0
  do i=1,dot
     call derivs(t,y,dydt)
     yt=y+h*dydt
	 t=t+h
     call derivs(t,yt,dyt) 
     y=y+hh*(dydt+dyt) 
     do j=1,N
        temp=temp+y(j+N)
     enddo
  enddo
  temp=temp/N/dot
  print*, k,changeF,temp
  averageVelocity(k)=temp
  parameters(k)=changeF
  y0=y
  enddo


!输出结果
  filenm="eeee.txt"
  open (5555, file = filenm)
  write (5555,*) " x= [ "
  do i = 1 , npara 
     write(5555,'(f13.5)') parameters(i)
  enddo
  write (5555,*) " ]; "
  write (5555,*) " y= [ "
  do i = 1 , npara 
     write(5555,'(f13.5)') averageVelocity(i) 
  enddo
  write (5555,*) " ]; "
  write (5555,*) 
  write (5555,*) "figure(8)"
  write (5555,*) "plot(x,y,'-b^');"
  print*,filenm
  close (5555)
 endsubroutine plot

subroutine derivs(t,y,dydt)
implicit none
real*8 t
real*8 y(NN),dydt(NN),c(N)
integer i,j
real*8 u(N),v(N),w(-N+1:NN),PForce,f1,f2
  do i=1,N
     u(i)=y(i); v(i)=y(i+N)
	 c(i)=normal(5.0d0,2.0d0)
	 !为了处理周期性边界条件Eq.(3)，引入如下中间变量
	 w(i)=y(i); w(i+N)=y(i)+N*b; w(i-N)=y(i)-N*b
  enddo
  do i=1,N
     PForce=0.d0;f1=0.d0;f2=0.d0
     do j=1,M
	    f1=f1-Kappa*(1.d0-dexp(b+w(i)-w(i+j))) * dexp(b+w(i)-w(i+j))   ! when x_i<x_j
		f2=f2+Kappa*(1.d0-dexp(b-w(i)+w(i-j))) * dexp(b-w(i)+w(i-j))
        PForce=f1+f2
	 enddo
     dydt(i)  =v(i)+D*c(i)
	 dydt(i+N)=-gamma*v(i) - dsin(2.d0*pi/a*u(i))/2.d0 - dsin(2.d0*pi*beta/a*u(i))/2.d0 -PForce+changeF+D*c(i)
  enddo
end subroutine  derivs

endmodule phase



program phase_plot
  use phase
  implicit none
  real*8 t1,t2
  call cpu_time(t1)
  call plot()
  call cpu_time(t2)
  print*, "Time used:", t2-t1, "seconds."
endprogram phase_plot