!One-Way Modeling /Migration....modified from Yaxiun's
module oneway_mod
  use history_mod 
  use io_mod
  use array_mod
  use param_mod 
  use wavelet_mod
  use myfft_mod
  use map_mod
  use grid_mod
  use vel_mod
  use ctaper_mod
  use ssf_mod
  implicit none

contains
!=======================================================================================
  integer function oneway_init() result(stat)
        call      mig_param()                                   !initialize migration parameters 
        call  wavelet_input()                                   !read wavelet history
        call      data_grid()                                   !get data grid     
        call  velocity_grid()                                   !get velocity grid     
        call  migration_grid()                                  !get velocity grid     
        call mig_output_history()                               !output migration history
        call mig_array_alloc()                                  !allocate memory 
        call mig_io_init()                                      !mig to zero
        if (wfld_d_out) call wfld_output_history(wfld_d_)       !output down wavefield history
        if (wfld_u_out) call wfld_output_history(wfld_u_)       !output up wavefield history
        if (wfld_d_out) call wfld_io_init(wfld_d_)              !down wavefield to zero
        if (wfld_u_out) call wfld_io_init(wfld_u_)              !up   wavefield to zero
    stat = 0
  end function
!=======================================================================================
  integer function oneway_clos() result(stat)
        call array_dealloc()
        stat = 0
  end function
!!=======================================================================================
  integer function oneway_oper() result(stat)
    icount=0;    
    call omp_set_num_threads(nodes); write(*,*) 'using', nodes, 'of' , nnode
!$OMP parallel do private(inode, iw, w, index_aw, wfld_d, wfld_u, iz, vel_slice, ind_vref, vref, nref, iv, ikx, iky, ix, iy, data, sx, isx, iisx, sy, isy, iisy, gx, igx, iigx, gy, igy, iigy, tmp_mig, wfld_slice_u, wfld_slice_d) schedule(dynamic,1) 
  do iw=1, dgrid%aw%n                                                            ! freq. count.
        inode=omp_get_thread_num();
        w = ggrid%aw%o + (iw-1)*ggrid%aw%d; index_aw=1.5 + (w-dgrid%aw%o)/ggrid%aw%d
        write(* ,*)'iw.:', iw , ' of ',ggrid%aw%n,' @ node ', inode      
        write(10,*)'iw.:', iw , ' of ',ggrid%aw%n,' @ node ', inode     
        do isx=1, dgrid%asx%n;    ! shot loop
                sx = dgrid%asx%o + (isx -1) * dgrid%asx%d; index_sx=1.5 + (sx-vgrid%ax%o)/vgrid%ax%d
        do isy=1, dgrid%asy%n;
                sy = dgrid%asy%o + (isy -1) * dgrid%asy%d; index_sy=1.5 + (sy-vgrid%ay%o)/vgrid%ay%d
                call mig_pro_grid(index_sx,index_sy,pgrid)
                call k2_map(pgrid);
                stat  = ssf_init(kx2, ky2, vgrid%az%d)
                if (allocated(wfld_slice_d)) deallocate(wfld_slice_d);  allocate(wfld_slice_d(pgrid%ax%n,pgrid%ay%n))
                if (allocated(wfld_slice_u)) deallocate(wfld_slice_u);  allocate(wfld_slice_u(pgrid%ax%n,pgrid%ay%n))
                if (allocated(tmp_wfld1   )) deallocate(tmp_wfld1   ); allocate(tmp_wfld1(pgrid%ax%n,pgrid%ay%n))
                if (allocated(tmp_wfld2   )) deallocate(tmp_wfld2   ); allocate(tmp_wfld2(pgrid%ax%n,pgrid%ay%n))
                if (pgrid%ay%n== 1 .or. pgrid%ax%n==1) then
                                phase_shaping = cmplx(cos(PI/4.), -sin(PI/4.)) / sqrt(w)          !!..2D: 45 deg. rotation i^(-1/2)
                else
                                phase_shaping = cmplx(0., - 1. * w)                               !!..3D: 90 deg. rotation i^(-1)
                end if
!
!       !$OMP CRITICAL
!                call data_io_data("reed", data, iw)
!       !$OMP END CRITICAL
!       do isx=1,dgrid%asx%n;  sx = dgrid%asx%o + (isx-1)*dgrid%asx%d ; iisx = 1.5 + (sx - pgrid%ax%o)/pgrid%ax%d !hsot loop
!       do isy=1,dgrid%asy%n;  sy = dgrid%asy%o + (isy-1)*dgrid%asy%d ; iisy = 1.5 + (sy - pgrid%ay%o)/pgrid%ay%d 
!                wfld_d(imin_z,iisx,iisy,isx,isy)=(csou(iw));
!                wfld_u(imin_z,pgrid%index_datpro_gxmin:pgrid%index_datpro_gxmax, &
!                         pgrid%index_datpro_gymin:pgrid%index_datpro_gymax)=data(:,:,isx,isy)
!                do iz=imin_z,pgrid%az%n;   ! depth loop
!                        if (iz /=imin_z) then 
!                                call vel_alloc(vel_slice,ind_vref,vref); vel_slice = vel(iz,:,:); call velocity(vref, vel_slice, ind_vref, nref, vnum);
!                                tmp_wfld_d=wfld_slice_d; wfld_slice_d=0.
!                                tmp_wfld_u=wfld_slice_u; wfld_slice_u=0.
!                                do iv=1, nref                                                                           !vel. loop  
!                                        tmp_wfld__d=tmp_wfld_d; 
!                                        tmp_wfld__u=tmp_wfld_u;
!                                        stat = ssf_oper(adj, -1, w, vref(iv), vel_slice, tmp_wfld__d)
!                                        stat = ssf_oper(adj, +1, w, vref(iv), vel_slice, tmp_wfld__u)
!                                        where (ind_vref == iv)
!                                             wfld_slice_d = wfld_slice_d + tmp_wfld__d
!                                             wfld_slice_u = wfld_slice_u + tmp_wfld__u
!                                        end where
!                                enddo
!                       else
!                           wfld_slice_d=wfld_d(1,:,:,isx,isy); 
!                           wfld_slice_u=wfld_u(1,:,:        ); 
!                       endif    
!                       tmp_mig(iz,:,:)=tmp_mig(iz,:,:)+real(wfld_slice_u(:,:)*conjg(wfld_slice_d(:,:))*phase_shaping);
!                enddo !end depth
        enddo
        enddo !end shot
        STOP
!       !$OMP CRITICAL
!       if (nostack)then; 
!           mig=tmp_mig; call mig_io_mig(mig, index_aw);  
!       else; 
!           write(*,*)'adding freq:', index_aw, 'sum_mig:',sum(mig), 'sum tmp_mig:',sum(tmp_mig); mig=mig+tmp_mig; call mig_io_mig(mig, 1); 
!       endif
!       !$OMP END CRITICAL
!       deallocate(wfld_d, wfld_u, tmp_mig)
!       endif
       enddo
!$OMP end parallel do

!       if (.not. adj) then; if(wfld_d_out) call wfld_out(wfld_out_d,wfld_d_); 
!       else
!       mig(:,:,:)=tmp_mig(:,:,:)
!       if(wfld_d_out) call wfld_out(wfld_out_d,wfld_d_); if(wfld_u_out) call wfld_out(wfld_out_u,wfld_u_);
!       stat = 0
!       deallocate(mig,tmp_mig)    
!       endif

       stat=0;

  end function

!!!!!!!!!!!!!!!!  ###########
  end module
