program main

  use sep
  use off2ang2d1_sd_mod
  use solver_tiny_mod
  use solver_prc_mod
  use cgstep_mod
  use igrad1
  use ident_mod

  implicit none

  real, dimension(:), pointer :: modl(:), data(:), m(:), d(:),mask(:)
  integer                     :: nz, nhx, nax, nmx, imx, iax, ihx, iz, stat, style, mode
  integer		      :: niter,NAop
  real                        :: z0, hx0, mx0 
  real                        :: dz, dhx, dax, dmx
  real                        :: eps, axmin, axmax, hxmin, hxmax
  real, parameter             :: PI = 3.14159265359
  logical, parameter          :: T = .true. , F = .false.
  logical	              :: adj	
  call initpar()

  call from_history("n1", nf); call from_history("d1", df); call from_history("o1", of)
  call from_history("n2", nh); call from_history("d2", dh); call from_history("o2", oh)  
  call from_history("n3", ns); call from_history("d3", ds); call from_history("o3", os)

  call from_aux("wave","n1", nw); call from_aux("wave","d1", dw); call from_aux("wave","o1", ow)

  call from_aux("velo","n1", nz); call from_aux("velo","d1", dz); call from_aux("velo","o1", oz)
  call from_aux("velo","n2", nm); call from_aux("velo","d2", dm); call from_aux("velo","o2", om)  

  call from_param("hmin", axmin,-1200.)
  call from_param("hmax", axmax, 1200.)
  call from_param("dhx", dax, 30.)
  nhx = (hmax - hmin) / dhx + 1

  call from_param("style", style, 0)
  call from_param("niter", niter, 50)
  call from_param("eps", eps, 1.)
  call from_param("mode",mode,1)

  call to_history("n1", nz) ; call to_history("d1", dz) ; call to_history("o1", z0)
  call to_history("n2", nax); call to_history("d2", dax); call to_history("o2", axmin) 
  call to_history("n3", nmx); call to_history("d3", dmx); call to_history("o3", mx0)

  allocate(shot(nf,nh), wave(nf), velo(nz,nm), migr(nz,nm))

  call sreed("wave", wave, 8*nf)
  call sreed("velo", velo, 4*nz*nm)
!     call sreed("in", shot, 8*nf*nh)

  pi = 4*atan(1.)
  kd = 1./(nx*dx)
  ko = 1./(2*dx)-kd 			!k does not increase uniformly!

! Downward propagation of the source function
  do is=1,ns
     shot = 0
     sx = (is - 1)*ds + os
     do iw=1,nw
        w = (iw - 1)*dw + ow
        do iz = nz,1,-1
           v0 = sum(v(iz,:))/nx
!           call fft_1d_data(fft_1d_id, -1, nx, C(iz,:))
!Operations in the wavenumber domain
           do ik = 1, kn
              k = mod(ko+(ik-1)*kd,1./dx)-ko
              k = 2*pi * k
              
!---phase shift term
              sqa = w**2 / v0**2 - k**2
!---avoid evanescent energy
              if ( sqa > 0 ) then
                U(iw,ik) = U(iw,ik) * exp( i * dz * sqrt(sqa) ) + C(iz,ik)
              else
                U(iw,ik) = 0.
              end if
           end do
        end do 
        call fft_1d_data(fft_1d_id, 1, k%n, U(iw,:))
     end do 
        


           call off2ang2d1_sd_init(z0, hx0, dz, dhx, nz, nhx, nax, dax, axmin, 1)
           stat = off2ang2d1_sd_oper(adj, F, modl, data)
	end if	
        call srite("out", modl, 4*nz*nax)
     end do
  else
     call from_history("n1", nz) ; call from_history("d1", dz) ; call from_history("o1", z0)
     call from_history("n2", nax); call from_history("d2", dax); call from_history("o2", axmin) 
     call from_history("n3", nmx); call from_history("d3", dmx); call from_history("o3", mx0)
     call from_param("hmin", hxmin,-600.)
     call from_param("hmax", hxmax, 600.)
     call from_param("dhx", dhx, 25.)
     nhx = (hxmax - hxmin) / dhx + 1
     call to_history("n1", nz ); call to_history("d1", dz ); call to_history("o1", z0 )
     call to_history("n2", nhx); call to_history("d2", dhx); call to_history("o2", hxmin)  
     call to_history("n3", nmx); call to_history("d3", dmx); call to_history("o3", mx0)

     allocate(modl(nz*nax), data(nz*nhx), m(nz*nax), d(nz*nhx))

     do imx=1,nmx
        data = 0;modl = 0.
        call sreed("in", modl, 4*nz*nax)
!--------------------------------------
!	m=modl;modl=0.
!        do iax=1,nax
!	   do iz=1,nz
!	      modl((iz-1)*nax+iax)=m((iax-1)*nz+iz)
!	   end do
!	end do
!-----------------------------------------
	if (style.eq.0) then
           call off2ang2d1_sd_init(z0, hxmin, dz, dhx, nz, nhx, nax, dax, axmin, 0)
           stat = off2ang2d1_sd_oper(adj, F, modl, data)
	else
           call off2ang2d1_sd_init(z0, hxmin, dz, dhx, nz, nhx, nax, dax, axmin, 1)
           stat = off2ang2d1_sd_oper(adj, F, modl, data)
	end if

!-------------------------------------
!	d=data;data=0.
!       do ihx=1,nhx
!	   do iz=1,nz
!	      data((ihx-1)*nz+iz)=d((iz-1)*nhx+ihx)
!	   end do
!	end do
!---------------------------------------
        call srite("out", data, 4*nz*nhx)
     end do
  end if 

end program
