program RaSe_
  use sep

  implicit none

  integer           :: n1,n2,n3,n4,n5, nsp
  integer           :: i1,i2,i3,i4,i5, isp
  real              :: d1,d2,d3,d4
  real              :: o1,o2,o3,o4
  complex, allocatable :: dat(:,:,:,:)
  integer, allocatable :: rn1(:), rn2(:)

  integer, dimension(:), allocatable :: seed,old
  integer                            :: i,k
  character(len=20)                  :: date, time, zone
  integer, dimension(8)              :: values
  integer                            :: s1,s2

  real  :: r

  call sep_init()

  call from_history("n1",n1,1)
  call from_history("n2",n2,1)
  call from_history("n3",n3,1)
  call from_history("n4",n4,1)
  call from_history("n5",n5,1)

  call from_history("d1",d1,1.)
  call from_history("d2",d2,1.)
  call from_history("d3",d3,1.)
  call from_history("d4",d4,1.)

  call from_history("o1",o1,0.)
  call from_history("o2",o2,0.)
  call from_history("o3",o3,1.)
  call from_history("o4",o4,1.)

  call from_param("nsp",nsp,1)

  call sep_close()

  allocate(dat(n1,n2,n3,n4))
  allocate(rn1(nsp), rn2(nsp))

  !!init the random number generator
  call random_seed(size=k)
  allocate(seed(k),old(k))
  call date_and_time(date,time,zone,values)

  s1=product(values(5:6)); s1=min(s1,31328)
  s2=product(values(5:7)); s2=min(s2,30081)
  seed(1:2) = (/s1, s2/)
  call random_seed(put = seed(1:k))

  do isp=1,nsp
     call random_number(r); rn1(isp)=ceiling((n1-1)*r); rn1(isp)=max(rn1(isp),1)
     call random_number(r); rn2(isp)=ceiling((n2-1)*r); rn2(isp)=max(rn2(isp),1)
  end do

  do i5=1,n5
     write(0,*) i5
     call sep_read(dat)
     do isp=1,nsp
        dat(rn1(isp),rn2(isp),:,:)=(0.,0.)
     end do
     call sep_write(dat)
  end do

  call exit(0)
end program RaSe_



