 module global 
   double precision, parameter :: pi = 4.d0 * atan(1.d0) 
   double precision a,s     ! some parameters 
   double precision time    ! the time
   double precision dt      ! the time step
   double precision dth     ! the time step**0.5
   character (len=10) seed_char
   integer seed
 end module global
!========================================================  
 program Milshtein 
!
! This code implement the Milshtein method to 
! integrate a one-dimensional stochastic flow 
! of the form 
!
!   dx
!   -- = f(x,t) + g(x) xi(t)
!   dt
!
! where xi(t) is a Gaussian stochastic process 
! with delta correlation and unit variance
!
! This code implements the linear equation with 
! multiplicative noise:
!
!   f(x,t) = a*x
!   g(x,t) = s*x
!
!======================================================== 
!
 use global
 integer N_steps          ! the total number of steps
 integer N_samples        ! the total number of samples
 double precision x       ! the state variable
 double precision f,g,dg  ! the model functions
 double precision x0      ! the initial conditions
 double precision T_tot   ! the total integration time
 character (len=20) arg

!======================================================== 
 if ( iargc() == 1 ) then
     call getarg(1,arg)
     read (arg,*) seed
 else 
     print*, ' Usage: milshtein.gfx [ seed ]'
     STOP
 end if
 write (seed_char,'(i10.10)') seed
!========================================================
! Model parameters

 a = -0.1d0
 s = 1.0d0

!======================================================== 
 do i=1, seed 
    x = rand()
 end do
 
 dt = 0.01
 dth = dsqrt(dt)
 N_samples = 100          ! 
 T_tot = 20. 
 N_steps = int(T_tot/dt)
 i_samp = int ( dfloat(N_steps) / dfloat(N_samples) )
 
 x0 = 0.01
 
!======================================================== 
! Integration

 x = x0
 do i=1, N_steps
    time = i*dt
    call mils_tstep(x)
    if ( mod(i,i_samp) == 0 ) call printout(x)
 end do

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 contains
 
 subroutine printout(x)
 use global 
 double precision x
 
 open (1,file='x_seed_'//seed_char//'.dat', &
         access='append',status='unknown')
 write (1,'(2f16.8)') time, x
 close (1)
 
 end subroutine 

!---------------------------------------------------------

 subroutine mils_tstep(x)
 use global 
 double precision x,r1,r2,u,f,g,dg
 
 r1 = rand()
 r2 = rand()
 u = sqrt(-2.0 * log(r1)) * cos (2.0 * pi * r2)
 
 call model(f,g,dg,x)
 x = x + g*dth*u + ( f + g*dg*u*u/2.0d0 )*dt
 
 end subroutine

!---------------------------------------------------------

 subroutine model(f,g,dg,x)
 use global 
 double precision x,f,g,dg
 
 f  = a*x
 g  = s*x 
 dg = s
 
 end subroutine

 end program  
 
 
 
 
 
 
 
