# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
program randmap
  use non_lin_pick_mod
  use sep
  use random
  use sep3d_struct_mod
  use cgstep_mod
  use erf_mod
  use npolydiv
  use solver_mod
  use io_nfilt
  use weight
  use quantile_mod
  implicit none
  type(sep3d) :: data,model,smooth
  type(nfilter) :: filt
  real , pointer :: mbuf(:),junk(:),dbuf(:,:),map(:,:),wt(:),lookup(:)&
    &,ttt(:),wt2(:),mid(:),std(:)
  real :: eps,f
  integer :: seed,seed_it(2),ierr,ireal
  integer :: i,niter,i4,i3,nmap,ilast,ih,i1,i2,nreal,ih2,iloc(1)
  real :: dmap,omap,xp,scale,down,pct,mymax
  real(kind=dp) :: p
  logical :: verb,found
  call from_param("down",down,.75)
  call from_param("pct",pct,.675)
  call from_param("verb",verb,.false.)
  call init_sep3d("in",data,"INPUT")
  call from_param("nreal",nreal,1)
  if (sep3d_ndims(data).eq.3) then
    call sep3d_reshape(data,(/1,2,3,3/))
  end if
  call init_sep3d(model,"OUTPUT","FLOAT","REGULAR",(/data%n(1), data%n&
    &(3), data%n(4)/),(/data%o(1), data%o(3), data%o(4)/),(/data%d&
    &(1), data%d(3), data%d(4)/))
  call init_sep3d(smooth,"OUTPUT","FLOAT","REGULAR",(/data%n(1), data%n&
    &(3), data%n(4)/),(/data%o(1), data%o(3), data%o(4)/),(/data%d&
    &(1), data%d(3), data%d(4)/))
  call sep3d_write_description("out",model)
  call sep3d_write_description("smooth",smooth)
  call to_history("n4",nreal,"out")
  call to_history("n4",nreal,"smooth")
!allocate(dbuf(product(data%n(:sep3d_ndims(data)))))
  allocate(mbuf(product(model%n(:sep3d_ndims(model)))))
  allocate(mbuf(product(model%n(:sep3d_ndims(model)))))
  allocate(mbuf(product(model%n(:sep3d_ndims(model)))))
!INITIALIZE RANDOM GENERATOR
  call from_param("seed",seed,-1)
  if (seed.ne.-1) then
    seed_it(1)=seed
    if (seed_it(1).eq.0) then
      seed_it(1)=2000
    end if
    seed_it(2)=seed_it(1)/31328
    seed_it(1)= seed_it(1)-seed_it(2)*31328
    write(0,*) "seed it",seed_it
    call random_seed(put=seed_it)
  end if
  if (verb) then
    write(0,*) "before filter read"
  end if
  if (.not. read_nfilt("filt",filt)) then
    call seperr("trouble writing filter")
  end if
  call npolydiv_init(size(mbuf),filt)
  allocate(dbuf(data%n(1),data%n(2)))
  allocate(wt2(data%n(1)*data%n(3)*data%n(4)))
  allocate(mid(data%n(1)*data%n(3)*data%n(4)))
  allocate(std(data%n(1)*data%n(3)*data%n(4)))
  call from_param("nmap",nmap,101)
  omap=0.
  dmap=1./real(nmap-1)
  allocate(map(nmap,data%n(1)*data%n(3)*data%n(4)))
  wt=1
  ih=0
  ih2=0
  do i4=1,data%n(4)
    do i3=1,data%n(3)
      call sreed("in",dbuf,size(dbuf)*4)
      do i1=1,data%n(1)
        iloc=maxloc(dbuf(:,i1))
        mymax=dbuf(iloc(1),i1)*down
        wt2(ih)=mymax/down
        i=1
        ih= ih+1
        found=.false.
        do while (.not. found .and. iloc(1)-i >=1 .and. iloc(1)+i <=&
          & data%n(2))
          if (dbuf(i1,iloc(1)-i) < mymax.or.dbuf(i1,iloc(1)+i) < mymax&
            &) then
            found=.true.
          end if
          i= i+1
        end do
        mid(ih)=(iloc(1)-1)*data%d(2)+data%o(2)
        std(ih)=(data%d(2)*i/.675)
      end do
    end do
  end do
  wt2=wt2**4
  call weight_init(wt)
  call from_param("niter",niter,100)
  call from_param("eps",eps,1.)
  allocate(ttt(size(mbuf)))
  allocate(junk(size(mbuf)))
  do ireal=1,nreal
    do i=1,size(mbuf)
      mbuf(i)=random_normal()
    end do
    if (verb) then
      write(0,*) "realization",ireal,nreal
    end if
    call solver_prec(i_op,cgstep,npolydiv_lop,size(mbuf),junk,mbuf,niter&
      &,eps,wt=wt2)
! call solver_prec(i_op,cgstep,i_op,size(mbuf),junk,mbuf,niter,eps,wt=wt2)
! call solver_prec(i_op,cgstep,i_op,size(mbuf),junk,mbuf,niter,eps)
!ttt=abs(junk)
! call srite("smooth",junk,size(junk)*4)
!call from_param("scale",scale,.5)
!scale=scale/maxval(abs(junk))
!write(0,*) "COMPARE",scale,.3/maxval(abs(junk))
!junk=junk*scale
    ih=0
    do i2=1,data%n(3)*data%n(4)
      do i1=1,data%n(1)
        ih= ih+1
        junk(ih)=mid(ih)+junk(ih)
      end do
    end do
    call srite("out",junk,size(junk)*4)
  end do
end program
