program greenf_Vz

  use sep, only: sep_init, sep_close, sep_read, sep_write, from_param
!  use ooio, only: axis, readaxis, writeaxis, chk_dim_corresp
  use ooio
!  use gr_post3dvz_mod, only: post3dvz_init, post3dvz_mixdom, clean_post3dvz, rft_fwd_axis1, rft_bck_axis1, post3dvz_fd,post3dvz_ffd
  use grf_vz_mod
!  use Shan_post3dvz_mod

  implicit none

  real   , allocatable, dimension (:,:)		 :: U ! input wavefield, 2-D, offset=0
  real   , allocatable, dimension (:,:)		 :: v ! interval velocity field
  complex, allocatable, dimension (:,:)	 	 :: C ! wavefield in Fourier domain
  real, allocatable, dimension (:,:,:)	 	 :: D ! Green_Fct.
  integer                               	 :: ix,iz,i,it
  integer                                	 :: op
  type(axis)                            	 :: x, z, t, w


  call sep_init()

  call from_param('operator',op,1)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Read/write axes
  t = readaxis(1, chk_even=.true.)
  if(t%o/=0) call seperr('For migration t must start at zero (o1=0 in input)')

  z = readaxis(1, from='vel')
  t = readaxis(1, from='in')
  x = readaxis(2, chk_even=.true.)

  call writeaxis(z, 1, axisname='z')
  call writeaxis(x, 2, axisname='x')
  call writeaxis(t, 3, axisname='t')

  write(0,*) z%n,x%n,t%n

  w = post3dvz_init(x, z, t)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Velocity stuff
  call chk_dim_corresp('vel' ,(/z, x/))
  call sep_close()
  allocate( v(z%n,x%n))
  call sep_read(v,'vel')
  if(any(v<=0.)) call seperr('zero or neg value in the velocity model!')
  v = v/2.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Input data
  allocate (U(t%n,x%n))
  call sep_read(U)
  allocate (C(w%n,x%n),D(z%n,x%n,t%n)); C=0! nt/2 frequencies, no zero freq
  C = rft_fwd_axis1(U)

  deallocate(U)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Modeling
!  time:do it = 1, t%n
!      depth:do iz = 1,z%n
!           if (op==1) then
!       call post3dvz_mixdom(C, v(iz,:),D) ! split-step Downward continuation step
               call post3dvz_mixdom(C, v, D) ! split-step Downward continuation step
!               C = rft_bck_axis1(C)
!          else if (op==2) then
!              call post3dvz_fd(C, v(iz,:),D) ! implicit finite-difference
!          else if (op==3) then
!              call post3dvz_ffd(C, v(iz,:),D) ! Foutier finite-difference
!           end if
           call sep_write(D) 
!           write(0,*) iz*100./z%n, '% migrated'
!      end do depth
!  end do time

  deallocate(C,v,D)
  call clean_post3dvz
end program greenf_Vz
