# 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&
    &(:)
  real ,pointer :: ttt(:),wt2(:),var(:,:,:),real(:),semb(:)
  real :: eps,f,mx
  integer :: seed,seed_it(2),ierr,ireal
  integer :: i,niter,i4,i3,nmap,ilast,ih,i1,i2,nreal,ih2
  real :: dmap,omap,pct,xp,scale
  real(kind=dp) :: p
  logical :: verb
  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("var",model)
  call sep3d_write_description("out",model)
  call sep3d_write_description("smooth",smooth)
  call to_history("n4",nreal,"out")
  call to_history("n4",nreal,"smooth")
  call to_history("n4",nreal,"var")
  allocate(semb(data%n(2)))
  do i=1,data%n(2)
    semb(i)=data%o(2)+data%d(2)*(i-1)
  end do
  write(0,*) "STYU 1",data%n(4)
!allocate(dbuf(product(data%n(:sep3d_ndims(data)))))
  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
  write(0,*) "STYU a"
  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(var(data%n(1),data%n(3),data%n(4)))
  allocate(dbuf(data%n(1),data%n(2),data%n(3),data%n(4)))
  allocate(wt2(data%n(1)*data%n(3)*data%n(4)))
  allocate(wt(data%n(1)*data%n(3)*data%n(4)))
  call from_param("nmap",nmap,101)
  omap=0.
  dmap=1./(nmap-1.)
  allocate(map(nmap,data%n(1)*data%n(3)*data%n(4)))
  wt=1
  ih=0
  ih2=0
  call sreed("in",dbuf,size(dbuf)*4)
  dbuf=dbuf**4
  do i4=1,data%n(4)
    do i3=1,data%n(3)
      do i1=2,data%n(1)
        ih2= ih2+1
        wt2(ih2)=maxval(dbuf(i1,:,i3,i4))
      end do
      do i1=1,data%n(1)
        dbuf(i1,:,i3,i4)=dbuf(i1,:,i3,i4)-min(.09,minval(dbuf(i1,:,i3&
          &,i4)))
      end do
      do i2=2,data%n(2)
        dbuf(:,i2,i3,i4)= dbuf(:,i2,i3,i4)+dbuf(:,i2-1,i3,i4)
      end do
      do i1=1,data%n(1)
        ih= ih+1
        wt(ih)=dbuf(i1,data%n(2),i3,i4)
        if (wt(ih)>.01) then
          dbuf(i1,:,i3,i4)=dbuf(i1,:,i3,i4)/wt(ih)
        end if
        map(1,ih)=data%o(2)
        ilast=1
        do i=2,nmap
          pct=omap+dmap*(i-1)
          do while (dbuf(i1,ilast,i3,i4) < pct)
            ilast= ilast+1
          end do
          if (ilast.eq.1) then
            map(i,ih)=data%o(2)
          else
            f=(pct-dbuf(i1,ilast-1,i3,i4))/(dbuf(i1,ilast,i3,i4)-dbuf&
              &(i1,ilast-1,i3,i4))
            map(i,ih)=data%o(2)+data%d(2)*(ilast-2+f)
          end if
        end do
      end do
    end do
  end do
  call srite("map.H",map,size(map)*4)
  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)))
  write(0,*) "what the 1"
!wt2=1
  do ireal=1,nreal
    write(0,*) "real ",ireal
    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)
    ttt=abs(junk)
    call srite("smooth",junk,size(junk)*4)
    call from_param("scale",scale,.5)
    scale=scale/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
        p=junk(ih)
        xp=erf(p)/2.+.5
        f=(xp-omap)/dmap+1.
        i=f
        f= f-i
        junk(ih)=(1.-f)*map(i,ih)+f*map(i+1,ih)
      end do
    end do
    var=0
    ih=0
    do i4=1,data%n(4)
      do i3=1,data%n(3)
        do i1=1,data%n(1)
          ih= ih+1
          do i2=1,data%n(2)
            var(i1,i3,i4)=(dbuf(i1,i2,i3,i4)*(semb(i2)-junk(ih)))**2
          end do
        end do
      end do
    end do
    var=var/data%n(2)
    call srite("var",var,size(var)*4)
    call srite("out",junk,size(junk)*4)
  end do
end program
