!One-Way Modeling /Migration....modified from Yaxiuns
module oneway_mig_demig
  use type_mod
  use global_defi_mod
  use reflect_mod
  use ricker_mod
  use vel_mod
  use ctaper_mod
  use sep
  implicit none

  integer   :: unw , iww !usefil freq
contains
!=======================================================================================
  integer function oneway_init() result(stat)
    call from_param("adj", adj, .false. )
    call vel_ref_input_history()                                !read vel history
    if (.not. adj)then
        call wavelet_input_history()                            !read wavelethistory
        call param_init()                                       !initialize modeling oneway parameters 
        if (.not. only_green) call data_output_history()        !output the data history
        call data_array_alloc()                                 !allocate memory 
        call propagation_grid()                                 !propagation grid
        call velocity_grid()                                    !velocity grid
        if (.not. only_green) call data_io_init()               !data to zero
        if (lgreen    ) call green_output_history()             !output the greens history
        if (lgreen    ) call greens_io_init()                   !greens to zero
        if (wfld_d_out) call wfld_output_history(wfld_d_)       !output the wavefield history
        if (wfld_d_out) call wfld_io_init(wfld_d_)              !down wavefield to zero
        call k2_map()
    else
        call wavelet_input_history()                            !read wavelet history
        call param_init()                                       !initialize migration parameters 
        call mig_output_history()                               !output migration history
        call mig_array_alloc()                                  !allocate memory 
        call propagation_grid()                                 !propagation grid
        call velocity_grid()                                    !velocity grid
        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
        call k2_map()
    endif
    stat = 0
  end function
!=======================================================================================
  integer function oneway_clos() result(stat)
        call array_dealloc()
        stat = 0
  end function
!!=======================================================================================
  integer function oneway_oper() result(stat)
    integer                                        :: index_aw
    integer                                        :: ixg, iyg, niw
    real,    dimension(:)          ,allocatable    :: vref
    real,    dimension(:,:)        ,allocatable    :: vel_slice,ref_slice
    integer, dimension(:,:)        ,allocatable    :: ind_vref
    real                                           :: w 
    complex , dimension(:,:)       ,allocatable    :: tmp_wfld,split 
    real    , dimension(:,:)       ,allocatable    :: kz_dz
    real    , dimension(:,:,:)     ,allocatable    :: tmp_mig
    complex , dimension(:,:,:)     ,allocatable    :: green, green_s, green_r
    complex , dimension(:,:,:,:,:) ,allocatable    :: wfld_d 
    complex , dimension(:,:,:    ) ,allocatable    :: wfld_u 
    complex , dimension(:,:,:    ) ,allocatable    :: kz_dzv 
    real    , dimension(:)         ,allocatable    :: kz2  
    integer                                        :: isz, isx, isy, iisz, iisx, iisy
    integer                                        :: igz, igx, igy, iigz, iigx, iigy
    real                                           :: z, x, y  
    real                                           :: sx, sy, gx, gy
    real                                           :: kz,wkbj
    integer                                        :: icount 
    complex                                        :: phase_shaping 
    icount=0;    
    open(unit=10, file=report_, status='unknown')


    nnode =omp_get_num_procs(); thread=omp_get_thread_num();  nodes=min(nnode-1,nodes) 

    call omp_set_num_threads(nodes); write(*,*) 'using', nodes, 'of' , nnode

    allocate(vel_slice(pgrid%ax%n,pgrid%ay%n)); allocate(ref_slice(pgrid%ax%n,pgrid%ay%n))
    if (.not. adj) then
        do iz=1,vgrid%az%n; call data_io_vel(vel_slice, iz); call data_io_ref(ref_slice, iz); vel(iz,:,:)=vel_slice; refl(iz,:,:)=ref_slice; enddo
    else
        do iz=1,vgrid%az%n; call data_io_vel(vel_slice, iz); vel(iz,:,:)=vel_slice; enddo
        allocate(mig(igrid%az%n,igrid%ax%n,igrid%ay%n),tmp_mig(pgrid%az%n,pgrid%ax%n,pgrid%ay%n)); mig=0.; tmp_mig=0.
    endif

!$OMP parallel do private (inode,&
			iw,&
			w,&
			index_aw,&
			ixg,&
			iyg,&
			wfld_d,&
			wfld_u,&
			iz,&
			vel_slice,&
			ind_vref,&
			vref,&
			nref,&
			green,&
			green_s,&
			green_r,&
			greens,&
			tmp_wfld,&
			kz2,&
			iv,&
			ikx,&
			iky,&
			ix,&
			iy,&
			data,&
			sx,&
			isx,&
			iisx,&
			sy,&
			isy,&
			iisy,&
			gx,&
			igx,&
			iigx,&
			gy,&
			igy,&
			iigy,&
			tmp_mig,&
			kz_dz,&
			kz_dzv,&
			split,&
			kz,&
			wkbj) schedule(dynamic,1) 

  do iw=1, dgrid%aw%n                                                            ! freq. count.
       w = dgrid%aw%o + (iw-1)*dgrid%aw%d; 
       inode=omp_get_thread_num(); 
       if(abs(csou(iw)) <= mcsou*.005 )cycle

       index_aw = 1.5 + (w - ggrid%aw%o)/dgrid%aw%d;
       allocate(wfld_d(imax_z, pgrid%ax%n, pgrid%ay%n, dgrid%asx%n, dgrid%asy%n)); wfld_d=0.;
       allocate(wfld_u(imax_z, pgrid%ax%n, pgrid%ay%n                          )); wfld_u=0.
       allocate(data(dgrid%agx%n, dgrid%agy%n, dgrid%asx%n, dgrid%asy%n)); data=0.
       allocate(tmp_mig(pgrid%az%n,pgrid%ax%n,pgrid%ay%n)); tmp_mig=0.
       if (allocated(split   )) deallocate(split   ); allocate(split(pgrid%ax%n,pgrid%ay%n))
       if (allocated(tmp_wfld)) deallocate(tmp_wfld); allocate(tmp_wfld(     pgrid%ax%n,pgrid%ay%n))
       if (allocated(kz_dz   )) deallocate(kz_dz   ); allocate(kz_dz   (     pgrid%ax%n,pgrid%ay%n ));kz_dz=0.
       !$OMP CRITICAL
                call data_io_data("reed", data, iw)
       !$OMP END CRITICAL
       if (pgrid%ay%n== 1 .or. pgrid%ax%n==1) then
       !!..2D: 45 deg. rotation i^(-1/2)
                phase_shaping = cmplx(cos(PI/4.), -sin(PI/4.)) / sqrt(w)           
       else
       !!..3D: 90 deg. rotation i^(-1)
                phase_shaping = cmplx(0., - 1. * w)       
       end if
       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)=conjg(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,imax_z;     ! depth loop
                         call vel_alloc(vel_slice,ind_vref,vref)
                         vel_slice=vel(iz,:,:)
                         call velocity(vref, vel_slice, ind_vref, nref, vnum);
                         if (allocated(green_s )) deallocate(green_s ); allocate(green_s (nref,pgrid%ax%n,pgrid%ay%n))
                         if (allocated(green_r )) deallocate(green_r ); allocate(green_r (nref,pgrid%ax%n,pgrid%ay%n))
                         if (allocated(kz2     )) deallocate(kz2     ); allocate(kz2     (nref                      ))
                         if (allocated(kz_dzv  )) deallocate(kz_dzv  ); allocate(kz_dzv  (pgrid%ax%n,pgrid%ay%n,nref));kz_dz=0.

                 if (iz/=imin_z)then
                         tmp_wfld=wfld_d(iz-1,:,:,isx,isy);call FFT2DX(tmp_wfld, 1); wfld_d(iz,:,:,isx,isy)=tmp_wfld;
                         tmp_wfld=wfld_u(iz-1,:,:        );call FFT2DX(tmp_wfld, 1); wfld_u(iz,:,:        )=tmp_wfld;

                         do iv=1, nref; green_s(iv,:,:)=wfld_d(iz,:,:,isx,isy);enddo
                         do iv=1, nref; green_r(iv,:,:)=wfld_u(iz,:,:        );enddo
                         do iv=1,nref                                                                           ! vel loop    
                                do ikx = 1, size(kx2)                                                           ! wavenumber loop
                                do iky = 1, size(ky2)                                                           !wavenumber loop
                                        kz2(iv) = w**2 / vref(iv)**2 - kx2(ikx) - ky2(iky);                     !phase shift
                                        if (kz2(iv) > 0) then                                                   !evanescence
                                                kz=sqrt(kz2(iv))
                                                green_s(iv,ikx,iky)=green_s(iv,ikx,iky)*exp(cmplx(0, +sign*pgrid%az%d*kz));
                                                green_r(iv,ikx,iky)=green_r(iv,ikx,iky)*exp(cmplx(0, -sign*pgrid%az%d*kz));
                                                wkbj= sqrt(kz_dz(ikx,iky)/kz);
                                        else
                                                green_s(iv,ikx,iky)=0.; green_r(iv,ikx,iky)=0.;
                                        endif
                                enddo                                                                           !end wavenumber loop
                                enddo
                                tmp_wfld=green_s(iv,:,:); call FFT2DX(tmp_wfld, -1); green_s(iv,:,:)=tmp_wfld
                                tmp_wfld=green_r(iv,:,:); call FFT2DX(tmp_wfld, -1); green_r(iv,:,:)=tmp_wfld

                                green_s(iv,:,:)=green_s(iv,:,:)*exp(cmplx(0,+sign*pgrid%az%d*w)*(1/vel_slice(:,:)-1/vref(iv)));
                                green_r(iv,:,:)=green_r(iv,:,:)*exp(cmplx(0,-sign*pgrid%az%d*w)*(1/vel_slice(:,:)-1/vref(iv)));

                                tmp_wfld=green_s(iv,:,:); call ctaper2d(tmp_wfld,tmp_wfld,pgrid,0,0,taper_pct,taper_type); green_s(iv,:,:)=tmp_wfld
                                tmp_wfld=green_r(iv,:,:); call ctaper2d(tmp_wfld,tmp_wfld,pgrid,0,0,taper_pct,taper_type); green_r(iv,:,:)=tmp_wfld
                        enddo
                 else
                         do iv=1, nref; green_s(iv,:,:)=wfld_d(iz,:,:,isx,isy);enddo
                         do iv=1, nref; green_r(iv,:,:)=wfld_u(iz,:,:        );enddo
                         do iv=1, nref                                                                           !vel. loop    
                                 do iky = 1, size(ky2)                                                           !wavenumber loop
                                 do ikx = 1, size(kx2)                                                           !wavenumber loop
                                         kz2(iv) = w**2 / vref(iv)**2 - kx2(ikx) - ky2(iky);                     !phase shift
                                         if (kz2(iv) >= 0) then                                                  !evanescence
                                                 kz=sqrt(kz2(iv)); kz_dzv(ikx,iky,iv)=kz;                        !wkbj
                                         else
                                         endif
                                 enddo
                                 enddo          
                          enddo
                  endif
                        do ix=1,pgrid%ax%n; do iy=1,pgrid%ay%n
                                wfld_d(iz,ix,iy,isx,isy)=green_s(ind_vref(ix,iy),ix,iy);
                                wfld_u(iz,ix,iy        )=green_r(ind_vref(ix,iy),ix,iy);
                                kz_dz(ix,iy)=kz_dzv(ix,iy,ind_vref(ix,iy))
                        enddo             ; enddo
                        tmp_mig(iz,:,:)=tmp_mig(iz,:,:)+real(wfld_u(iz,:,:)*conjg(wfld_d(iz,:,:,isx,isy))*phase_shaping);
                        if(wfld_d_out) wfld_out_d(:,:,:,iw,1) = wfld_d(:,pgrid%index_xmin:pgrid%index_xmax, &
                                                                         pgrid%index_ymin:pgrid%index_ymax, & 
                                                                         nint(dgrid%agx%n/2.),nint(dgrid%agy%n/2.));
                        if(wfld_u_out) wfld_out_u(:,:,:,iw,1) = wfld_u(:,pgrid%index_xmin:pgrid%index_xmax,pgrid%index_ymin:pgrid%index_ymax );
                enddo !end depth
       enddo
       enddo !end shot
       !$OMP CRITICAL
           write(*,*)'adding freq....', index_aw, sum(mig), sum(tmp_mig); mig=mig+tmp_mig; call mig_io_mig(mig, 1); 
       !$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
  write(*,*)pgrid%index_datpro_gymin,pgrid%index_datpro_gymax
  write(*,*)size(tmp_mig,1),size(tmp_mig,2),size(tmp_mig,3)
  write(*,*)size(mig,1),size(mig,2),size(mig,3)

  mig(:,:,:)=tmp_mig(:,:,:)
  write(*,*)'Max',maxval(tmp_mig),maxval(mig)
!  if (.not. nostack)  call sep_write(mig, mig_, esize=4)

  !mig(:,:,:)=tmp_mig(:,pgrid%index_datpro_gxmin:pgrid%index_datpro_gxmax, & 
  !                     pgrid%index_datpro_gymin:pgrid%index_datpro_gymax); call sep_write(mig, mig_, esize=4)
  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

!=========================================================================================
! compute data from all model points for one shot and one receiver
  subroutine data_compute(data, wfld_d)
    complex, dimension(:       ), pointer       :: wavelet
    complex, dimension(:,:,:,: ), pointer       :: data
    complex, dimension(:,:,:,:,:), pointer       :: wfld_d
    logical                                     :: add
    integer                                     :: isz, isx, isy
    integer                                     :: igz, igx, igy
    integer                                     :: iisz, iisx, iisy
    integer                                     :: iigz, iigx, iigy
    integer                                     :: iw
    real                                        :: z, x, y  
    real                                        :: sx , sy, gx, gy 
    complex                                     :: greens_mul 
    real, dimension(1)                          :: tmp

       write(0,*)'computing data frequency: ', iw
        do isx=1, dgrid%asx%n;        sx = dgrid%asx%o + (isx-1)*dgrid%asx%d;  iisx = 1.5 + (sx-pgrid%gxmin)/pgrid%ax%d
        do isy=1, dgrid%asy%n;        sy = dgrid%asy%o + (isy-1)*dgrid%asy%d;  iisy = 1.5 + (sy-pgrid%gymin)/pgrid%ay%d
                do igx=1, dgrid%agx%n; gx = dgrid%agx%o + (igx-1)*dgrid%agx%d;  iigx = 1.5 + (gx-pgrid%gxmin)/pgrid%ax%d
                do igy=1, dgrid%agy%n; gy = dgrid%agy%o + (igy-1)*dgrid%agy%d;  iigy = 1.5 + (gy-pgrid%gymin)/pgrid%ay%d
                        do iz = imin_z, pgrid%az%n;      
                        do ix=ggrid%index_modpro_minx, ggrid%index_modpro_maxx; 
                        do iy=ggrid%index_modpro_miny, ggrid%index_modpro_maxy; 
                                data(igx,igy,isx,isy)=data(igx,igy,isx,isy)+wfld_d(iz,ix,iy,iisx,iisy)*wfld_d(iz,ix,iy,iigx,iigy)*refl(iz,ix,iy);   ! convolve to generate data
                       end do; 
                       end do; 
                       end do; 
                end do; 
                end do; 
        end do; 
        end do
  end subroutine

!==========================================================================================
  subroutine data_array_alloc()
    integer :: memory
    allocate(wavelet(1))
    allocate(vel(vgrid%az%n,  vgrid%ax%n, vgrid%ay%n)); 
    if (.not. only_green .and. (.not. adj)) allocate(refl(rgrid%az%n, rgrid%ax%n, rgrid%ay%n)); 
    if (.not. only_green) allocate( data(dgrid%agx%n, dgrid%agy%n, dgrid%asx%n, dgrid%asy%n)); 
    if (.not. lfreq) allocate(tdata(dgrid%at%n )); 
    if (lgreen) allocate(greens(ggrid%ax_surf%n, ggrid%ay_surf%n, ggrid%ax_modl%n, ggrid%ay_modl%n, ggrid%az_modl%n)) ; 
    if (wfld_d_out) allocate(wfld_out_d(vgrid%az%n,vgrid%ax%n,vgrid%ay%n,dgrid%at%n,1)) ; 

    memory = 8*size(wavelet)    + & 
             4 * size(vel)      + 4 * size(refl)        + &
             8 * size(data)     + 8 * size(tdata)       + & 
             8 * (size(greens)) + 4 * size(wfld_out_d)
    memory = memory / 1024 / 1024
    write(0,*)"Memory allocated: ", memory, "Mb"
  end subroutine

!=======================================================================================
   subroutine mig_array_alloc()
    integer :: memory
    allocate(wavelet(1))
    allocate(vel(vgrid%az%n, vgrid%ax%n, vgrid%ay%n)); 
    allocate(mig(vgrid%az%n, vgrid%ax%n, vgrid%ay%n)); 
    if (wfld_d_out) allocate(wfld_out_d(vgrid%az%n,vgrid%ax%n,vgrid%ay%n,dgrid%at%n,1)) ; 
    if (wfld_u_out) allocate(wfld_out_d(vgrid%az%n,vgrid%ax%n,vgrid%ay%n,dgrid%at%n,1)) ; 

    memory = 8*size(wavelet) + & 
             4 * size(vel       )     + 4 * size(mig       ) +  & 
             4 * size(wfld_out_u)     + 4 * size(wfld_out_d)
    memory = memory / 1024 / 1024
    write(0,*)"Memory allocated: ", memory, "Mb"
  end subroutine
 
!=======================================================================================
  subroutine array_dealloc()
        if(associated(wavelet)) deallocate(wavelet)     ;       if(associated(vel    )) deallocate(vel)
        if(associated(refl   )) deallocate(refl)        ;       if(associated(data   )) deallocate(data)
        if(associated(wavelet)) deallocate(wavelet)     ;       if(associated(greens )) deallocate(greens)
        if(associated(wfld_out_d )) deallocate(wfld_out_d);     if(associated(wfld_out_u )) deallocate(wfld_out_u)
  end subroutine
!==========================================================================================
  subroutine param_init()
    real    :: xmin_modl_default, xmax_modl_default
    real    :: ymin_modl_default, ymax_modl_default
    real    :: surface_xmin_default, surface_xmax_default
    real    :: surface_ymin_default, surface_ymax_default
    sou_     = "sou"    ; rec_     = "rec"    ;  data_    = "dat"    ; tdata_   = "tdat"   ;       
    vel_     = "vel"    ; img_     = "img"    ;  ref_     = "ref"    ; mig_     = "mig";
    wfld_d_  = "wfld_d" ; wfld_u_  = "wfld_u" ;  hdr_sou_ = "hdr_sou"; hdr_rec_ = "hdr_rec";
    hessian_ = "hessian"; greens_  = "greens" ;  datumn_  = "datumn" ; report_ = "report"    ; 

    if (.not. adj)then
        call from_param("ngx", dgrid%agx%n, vgrid%ax%n)
        call from_param("dgx", dgrid%agx%d, vgrid%ax%d)
        call from_param("ogx", dgrid%agx%o, vgrid%ax%o)
        call from_param("ngy", dgrid%agy%n, vgrid%ay%n)
        call from_param("dgy", dgrid%agy%d, vgrid%ay%d)
        call from_param("ogy", dgrid%agy%o, vgrid%ay%o)

        call from_param("nsx", dgrid%asx%n, dgrid%agx%n)   !num. sources inline
        call from_param("dsx", dgrid%asx%d, dgrid%agx%d)   !samp. sources inline 
        call from_param("osx", dgrid%asx%o, dgrid%agx%o)   !orig. sources inline
        call from_param("nsy", dgrid%asy%n, dgrid%agy%n)   !              x-line
        call from_param("dsy", dgrid%asy%d, dgrid%agy%d)   !
        call from_param("osy", dgrid%asy%o, dgrid%agy%o)   !

    else
        call auxin(data_)
        call from_aux(data_, "n1", dgrid%agx%n)  
        call from_aux(data_, "d1", dgrid%agx%d) 
        call from_aux(data_, "o1", dgrid%agx%o) 

        call from_aux(data_, "n2", dgrid%agy%n)  
        call from_aux(data_, "d2", dgrid%agy%d) 
        call from_aux(data_, "o2", dgrid%agy%o) 

        !call from_aux(data_, "n3", dgrid%aw%n)  
        !call from_aux(data_, "d3", dgrid%aw%d)  
        !call from_aux(data_, "o3", dgrid%aw%o) 

        call from_aux(data_, "n4", dgrid%asx%n)  
        call from_aux(data_, "d4", dgrid%asx%d) 
        call from_aux(data_, "o4", dgrid%asx%o) 

        call from_aux(data_, "n5", dgrid%asy%n)  
        call from_aux(data_, "d5", dgrid%asy%d) 
        call from_aux(data_, "o5", dgrid%asy%o) 
    endif

    dgrid%nsh = dgrid%asx%n * dgrid%asy%n               !num  of sources..
    dgrid%ntr = dgrid%agx%n * dgrid%agy%n               !num. of traces/recievers

    write(0,*)"==============shot_x: ", dgrid%asx
    write(0,*)"==============shot_y: ", dgrid%asy
    write(0,*)"=============geoph_x: ", dgrid%agx
    write(0,*)"=============geoph_y: ", dgrid%agy
    write(0,*)"=============vgrid_x: ", vgrid%ax
    write(0,*)"=============vgrid_y: ", vgrid%ay
    write(0,*)"=============vgrid_z: ", vgrid%az

    call from_param("nxpad_beg", pgrid%nxpad_beg, 0)
    call from_param("nxpad_end", pgrid%nxpad_end, 0)
    call from_param("nypad_beg", pgrid%nypad_beg, 0)
    call from_param("nypad_end", pgrid%nypad_end, 0)

    call from_param("nx_modl", ggrid%ax_modl%n, vgrid%ax%n-2*pgrid%nxpad_beg)
    call from_param("dx_modl", ggrid%ax_modl%d, vgrid%ax%d)
    call from_param("ox_modl", ggrid%ax_modl%o, ggrid%surface_xmin)

    call from_param("ny_modl", ggrid%ay_modl%n, vgrid%ay%n)
    call from_param("dy_modl", ggrid%ay_modl%d, vgrid%ay%d)
    call from_param("oy_modl", ggrid%ay_modl%o, vgrid%ay%o)

    call from_param("nz_modl", ggrid%az_modl%n, vgrid%az%n)
    call from_param("dz_modl", ggrid%az_modl%d, vgrid%az%d)
    call from_param("oz_modl", ggrid%az_modl%o, vgrid%az%o)

    surface_xmin_default = dgrid%agx%o
    surface_xmax_default = dgrid%agx%o + (dgrid%agx%n-1)*dgrid%agx%d

    ggrid%index_modpro_minx=1.5+(ggrid%ax_modl%o-vgrid%ax%o)/vgrid%ax%d; ggrid%index_modpro_maxx=ggrid%index_modpro_minx+(ggrid%ax_modl%n-1)  
    ggrid%index_modpro_miny=1.5+(ggrid%ay_modl%o-vgrid%ay%o)/vgrid%ay%d; ggrid%index_modpro_maxy=ggrid%index_modpro_miny+(ggrid%ay_modl%n-1)   
    ggrid%index_modpro_minz=1.5+(ggrid%az_modl%o-vgrid%az%o)/vgrid%az%d; ggrid%index_modpro_maxz=min(ggrid%index_modpro_minz+ggrid%az_modl%n-1,vgrid%az%n)
    ggrid%az_modl%n=ggrid%index_modpro_maxz-ggrid%index_modpro_minz+1
    ggrid%ay_modl%n=ggrid%index_modpro_maxy-ggrid%index_modpro_miny+1
    ggrid%ax_modl%n=ggrid%index_modpro_maxx-ggrid%index_modpro_minx+1
    
    surface_ymin_default = vgrid%ay%o
    surface_ymax_default = vgrid%ay%o + (vgrid%ay%n-1)*vgrid%ay%d

    igrid%ax=vgrid%ax; igrid%ay=vgrid%ay;    igrid%az=vgrid%az;

!! assume receiver grids is same as velocity grid
    call from_param("surface_xmin", ggrid%surface_xmin, surface_xmin_default)
    call from_param("surface_xmax", ggrid%surface_xmax, surface_xmax_default)
    call from_param("surface_ymin", ggrid%surface_ymin, surface_ymin_default)
    call from_param("surface_ymax", ggrid%surface_ymax, surface_ymax_default)


    call from_param("report"      , report_      , 'report')
    call from_param("wfld_marine" , marine       , .false.)
    call from_param("nodes"       , nodes        ,      1 )
    call from_param("only_green"  , only_green   , .false.)
    call from_param("marine"      , marine       , .false.)
    call from_param("wfld_out_jt" , wfld_out_jt  , 10     )
    call from_param("sign"        , sign         , -1     )
    call from_param("radian"      , radian       , .false.)
    call from_param("spike"       , spike        , .false.)
    call from_param("vnum"        , vnum         , 2      )
    call from_param("nref"        , nref         , 2      )
    call from_param("lloyd"       , lloyd        , .false.)
    call from_param("taper"       , taper        , .true. )
    call from_param("taper_pct"   , taper_pct    , 0.1    )
    call from_param("taper_type"  , taper_type   , 2      )
    call from_param("wfld_d_out"  , wfld_d_out   , .false.)
    call from_param("wfld_u_out"  , wfld_u_out   , .false.)
    call from_param("green"       , lgreen       , .false.)
    call from_param("for"         , for          , .true. ) 
    call from_param("ricker"      , lricker      , .false.) 
    call from_param("freq"        , lfreq        , .true. )
    call from_param("nostack"     , nostack      , .false. )
    call from_param("max_z"       , max_z        , vgrid%az%o+vgrid%az%d*(vgrid%az%n-1))
    call from_param("min_z"       , min_z        , vgrid%az%o)
    imax_z=1.5+(max_z-vgrid%az%o)/vgrid%az%d 
    imin_z=1.5+(min_z-vgrid%az%o)/vgrid%az%d
    write(*,*)'max_z,min_z, imin_z,imax_z'
    write(*,*)max_z,min_z, imin_z,imax_z
    if (lgreen)then
            write(*,*)'Summary greens====================='
            write(*,*) 'imodx',ggrid%index_modpro_minx, ggrid%index_modpro_maxx
            write(*,*) 'imody',ggrid%index_modpro_miny, ggrid%index_modpro_maxy
            write(*,*) 'imodz',ggrid%index_modpro_minz, ggrid%index_modpro_maxz
            write(0,*) 'surf_x', ggrid%surface_xmin, ggrid%surface_xmax
            write(*,*) 'spike ',spike
            write(0,*) 'green_x',ggrid%ax_modl
            write(0,*) 'green_y',ggrid%ay_modl
            write(0,*) 'green_z',ggrid%az_modl
            write(*,*)'==================================='
    endif
  end subroutine
!========================================================================================== 
 subroutine vel_ref_input_history()
    integer             :: num_match
    character(len=1024) :: string
    vel_= "vel"; ref_ = "ref"
    call from_param("for" , for , .false. ) 

    !read the velocity from file   
    call auxin(vel_)
    call from_aux(vel_, "n1", vgrid%az%n)  
    call from_aux(vel_, "d1", vgrid%az%d) 
    call from_aux(vel_, "o1", vgrid%az%o) 

    call from_aux(vel_, "n2", vgrid%ax%n)
    call from_aux(vel_, "d2", vgrid%ax%d)
    call from_aux(vel_, "o2", vgrid%ax%o)

    call from_aux(vel_, "n3", vgrid%ay%n) 
    call from_aux(vel_, "d3", vgrid%ay%d) 
    call from_aux(vel_, "o3", vgrid%ay%o)

    !read the reflectivity from file if provided...if not generate reflectivity
    num_match =  exist_file(ref_)
    if ( num_match==1 .or. for ) then
        call auxin(ref_)
        call from_aux(ref_, "n1", rgrid%az%n)  
        call from_aux(ref_, "d1", rgrid%az%d) 
        call from_aux(ref_, "o1", rgrid%az%o) 

        call from_aux(ref_, "n2", rgrid%ax%n)
        call from_aux(ref_, "d2", rgrid%ax%d)
        call from_aux(ref_, "o2", rgrid%ax%o)

        call from_aux(ref_, "n3", rgrid%ay%n) 
        call from_aux(ref_, "d3", rgrid%ay%d) 
        call from_aux(ref_, "o3", rgrid%ay%o)

        if (rgrid%az%n /= vgrid%az%n .or. rgrid%ax%n /= vgrid%ax%n .or. rgrid%ay%n /= vgrid%ay%n) &
        call seperr("reflectivity and velocity grids must match...see vel_ref_input_history()")
    else
        if (.not. adj .and. .not. only_green) call reflec(vel_, ref_, vgrid) !generate reflectivity model
        rgrid%az=vgrid%az; rgrid%ax=vgrid%ax; rgrid%ay=vgrid%ay
    end if
    write(*,*)rgrid%az%o, rgrid%az%n, rgrid%az%d
    write(*,*)vgrid%az%o, vgrid%az%n, vgrid%az%d

  end subroutine
!========================================================================================== 
 subroutine wavelet_input_history()
    integer             :: iw, niw
    real                :: w 
    integer             :: num_match
    character(len=1024) :: string

    wavelet_="wavelet"; twavelet_="twavelet"

    call from_param("nt", dgrid%at%n, 1001 )
    call from_param("dt", dgrid%at%d, 0.004)
    call from_param("ot", dgrid%at%o, 0.   )
    call from_param("fund",fund   , 20.    )
    call from_param("flip", lflip , .true. ) 
    call from_param("tmax", tmax, dgrid%at%o + dgrid%at%d*(dgrid%at%n-1))
    call from_param("jw", jw , 1);   

!read the source wavelet from file    
    if (exist_file(wavelet_)) then
        call auxin(wavelet_)
        call from_aux(wavelet_, "n3", dgrid%aw%n); call from_aux(wavelet_, "d3", dgrid%aw%d); call from_aux(wavelet_, "o3", dgrid%aw%o)
    else
        call ricker(wavelet_, dgrid%at, fund, floor(dgrid%at%n*.1), lflip);
        call from_aux(wavelet_, "n3", dgrid%aw%n); call from_aux(wavelet_, "d3", dgrid%aw%d);  call from_aux(wavelet_, "o3", dgrid%aw%o)

    endif
    if (.not. radian) call data_radians_calc()

    allocate(csou(dgrid%aw%n),wavelet(1))
    do iw=1,dgrid%aw%n
            call data_io_wavelet(wavelet, iw)
            csou(iw)=wavelet(1)
    end do
    mcsou=maxval(abs(csou))
!=compute the useful freq.!!!!!!!!!!!!!!!!!!!!!!!!
!only take values >=0.5% of max value in wavelet
    iw=0.;
    !do niw=floor(dgrid%aw%n*0.1),dgrid%aw%n-1
    do niw=1,dgrid%aw%n-1
           !if (abs(csou(niw))<=mcsou*.001) then; iw=iw+1; if (iw>=floor(dgrid%aw%n*0.02)) exit; endif
           if (abs(csou(niw))<=mcsou*.005 .and. niw <= 0.1*dgrid%aw%n) iw=iw+1;
           if (abs(csou(niw))<=mcsou*.005 .and. niw >= 0.2*dgrid%aw%n) exit;
    enddo
    
    unw=niw-iw
    ggrid%aw%o=dgrid%aw%o+(iw-1)*dgrid%aw%d; ggrid%aw%d=dgrid%aw%d*jw; ggrid%aw%n=(niw-iw)/jw; 
    deallocate(wavelet)
    write(*,*)'---------------'
    write(*,*)'dgrid%aw & ggrid%aw '
    write(*,*)dgrid%aw
    write(*,*)ggrid%aw
    write(*,*)'useful freq/.......',iw,niw,unw
    write(*,*)'---------------'
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  end subroutine
!==========================================================================================
  subroutine data_output_history()
        if (lfreq)then
                write(*,*)'data===========>       freq!'
                call auxinout(data_)
                call auxputch("n1", "i", dgrid%agx%n, data_) ! 
                call auxputch("d1", "r", dgrid%agx%d, data_) !                       
                call auxputch("o1", "r", dgrid%agx%o, data_) !                     
                call auxputch("label1", "s", "geoph_x", data_)

                call auxputch("n2", "i", dgrid%agy%n, data_)
                call auxputch("d2", "r", dgrid%agy%d, data_)
                call auxputch("o2", "r", dgrid%agy%o, data_)
                call auxputch("label2", "s", "geoph_y", data_)

                call auxputch("n3", "i", dgrid%aw%n, data_) !
                call auxputch("d3", "r", dgrid%aw%d/(2*PI), data_) ! 
                call auxputch("o3", "r", dgrid%aw%o/(2*PI), data_) !
                call auxputch("label3", "s", "frequency", data_)

                call auxputch("n4", "i", dgrid%asx%n, data_)
                call auxputch("d4", "r", dgrid%asx%d, data_)
                call auxputch("o4", "r", dgrid%asx%o, data_)
                call auxputch("label4", "s", "shot_x", data_)

                call auxputch("n5", "i", dgrid%asy%n, data_)
                call auxputch("d5", "r", dgrid%asy%d, data_)
                call auxputch("o5", "r", dgrid%asy%o, data_)
                call auxputch("label5", "s", "shot_y", data_)

                call to_history('esize', 8, data_)

        else
                call auxinout(tdata_)
                call auxputch("n1", "i", dgrid%agx%n, tdata_)  
                call auxputch("d1", "r", dgrid%agx%d, tdata_)                        
                call auxputch("o1", "r", dgrid%agx%o, tdata_)                      
                call auxputch("label1", "s", "geoph_x", tdata_)

                call auxputch("n2", "i", dgrid%agy%n, tdata_)
                call auxputch("d2", "r", dgrid%agy%d, tdata_)
                call auxputch("o2", "r", dgrid%agy%o, tdata_)
                call auxputch("label2", "s", "geoph_y", tdata_)

                call auxputch("n3", "i", dgrid%at%n, tdata_) 
                call auxputch("d3", "r", dgrid%at%d, tdata_)  
                call auxputch("o3", "r", dgrid%at%o, tdata_) 
                call auxputch("label3", "s", "time", tdata_)

                call auxputch("n4", "i", dgrid%asx%n, tdata_)
                call auxputch("d4", "r", dgrid%asx%d, tdata_)
                call auxputch("o4", "r", dgrid%asx%o, tdata_)
                call auxputch("label4", "s", "shot_x", tdata_)

                call auxputch("n5", "i", dgrid%asy%n, tdata_)
                call auxputch("d5", "r", dgrid%asy%d, tdata_)
                call auxputch("o5", "r", dgrid%asy%o, tdata_)
                call auxputch("label5", "s", "shot_y", tdata_)

                call to_history('esize', 4, tdata_)
                write(*,*)'data=========>       time'
        endif
  end subroutine
!#####################################
  subroutine mig_output_history()
                call auxinout(mig_)
                call auxputch("n1", "i", igrid%az%n, mig_) ! 
                call auxputch("d1", "r", igrid%az%d, mig_) !                       
                call auxputch("o1", "r", igrid%az%o, mig_) !                     
                call auxputch("label1", "s", "z", mig_)

                call auxputch("n2", "i", igrid%ax%n, mig_) ! 
                call auxputch("d2", "r", igrid%ax%d, mig_) !                       
                call auxputch("o2", "r", igrid%ax%o, mig_) !                     
                call auxputch("label2", "s", "x", mig_)

                call auxputch("n3", "i", igrid%ay%n, mig_) ! 
                call auxputch("d3", "r", igrid%ay%d, mig_) !                       
                call auxputch("o3", "r", igrid%ay%o, mig_) !                     
                call auxputch("label3", "s", "y", mig_)

                if (nostack) then
                        call auxputch("n4", "i", ggrid%aw%n , mig_) !
                else
                        call auxputch("n4", "i", 1          , mig_) !
                endif
                call auxputch("d4", "r", ggrid%aw%d, mig_) !
                call auxputch("o4", "r", ggrid%aw%o, mig_) !
                call auxputch("label4", "s", "w", mig_)
                call to_history('esize', 4, mig_)
  end subroutine
!#####################################
 subroutine green_output_history()

        ggrid%ax_surf%o = ggrid%surface_xmin
        ggrid%ax_surf%d = vgrid%ax%d
        ggrid%ax_surf%n = 1.5 + (ggrid%surface_xmax - ggrid%surface_xmin) / ggrid%ax_surf%d
   
        ggrid%ay_surf%o = ggrid%surface_ymin
        ggrid%ay_surf%d = vgrid%ay%d
        ggrid%ay_surf%n = 1.5 + (ggrid%surface_ymax - ggrid%surface_ymin) / ggrid%ay_surf%d

        call auxinout(greens_)
        call auxputch("n1", "i", ggrid%ax_surf%n, greens_)
        call auxputch("d1", "r", ggrid%ax_surf%d, greens_)
        call auxputch("o1", "r", ggrid%ax_surf%o, greens_)

        call auxputch("n2", "i", ggrid%ay_surf%n, greens_)
        call auxputch("d2", "r", ggrid%ay_surf%d, greens_)
        call auxputch("o2", "r", ggrid%ay_surf%o, greens_)

        call auxputch("n3", "i", ggrid%ax_modl%n, greens_)
        call auxputch("d3", "r", ggrid%ax_modl%d, greens_)
        call auxputch("o3", "r", ggrid%ax_modl%o, greens_)

        call auxputch("n4", "i", ggrid%ay_modl%n, greens_)
        call auxputch("d4", "r", ggrid%ay_modl%d, greens_)
        call auxputch("o4", "r", ggrid%ay_modl%o, greens_)

        call auxputch("n5", "i", ggrid%az_modl%n, greens_)
        call auxputch("d5", "r", ggrid%az_modl%d, greens_)
        call auxputch("o5", "r", ggrid%az_modl%o, greens_)

        call auxputch("n6", "i",      ggrid%aw%n, greens_)
        call auxputch("d6", "r",      ggrid%aw%d/(2*PI), greens_)
        call auxputch("o6", "r",      ggrid%aw%o/(2*PI), greens_)

        call auxputch("esize" , "i",         8   , greens_)
        call auxputch("label1", "s", "surface_x" , greens_)
        call auxputch("label2", "s", "surface_y" , greens_)
        call auxputch("label3", "s", "model_x"   , greens_)
        call auxputch("label4", "s", "model_y"   , greens_)
        call auxputch("label5", "s", "model_z"   , greens_)
        call auxputch("label6", "s", "frequency" , greens_)
  end subroutine
!##################################################
  subroutine wfld_output_history(wfld_)
        character(len=256)        :: wfld_
        call auxinout(wfld_d_)
        call auxputch("n1","i",vgrid%az%n,wfld_); call auxputch("n2","i",vgrid%ax%n,wfld_); call auxputch("n3","i",dgrid%at%n/wfld_out_jt,wfld_); 
        call auxputch("d1","r",vgrid%az%d,wfld_); call auxputch("d2","r",vgrid%ax%d,wfld_); call auxputch("d3","r",dgrid%at%d*wfld_out_jt,wfld_);
        call auxputch("o1","r",vgrid%az%o,wfld_); call auxputch("o2","r",vgrid%ax%o,wfld_); call auxputch("o3","r",dgrid%at%o            ,wfld_);
        call auxputch("label1","s","nz",wfld_)
        call auxputch("label2","s","nx",wfld_)
        call auxputch("label3","s","t" ,wfld_)
        call auxputch("esize" ,"i", 4  ,wfld_)
  end subroutine

!===================================================
  subroutine data_io_init()
    complex, allocatable, dimension(:) :: ctmp
    integer                            :: i, nloop
    if(lfreq)then
        allocate(ctmp(dgrid%agx%n*dgrid%agy%n*dgrid%asx%n*dgrid%asy%n))
        ctmp = (0,0)
        nloop = dgrid%aw%n
        do i=1, nloop
                call sep_write(ctmp, data_, esize=8)
        end do
    else
        allocate(ctmp(dgrid%at%n))
        ctmp = 0
        nloop = dgrid%agx%n*dgrid%agy%n*dgrid%asx%n*dgrid%asy%n
        do i=1, nloop
                call sep_write(ctmp, data_, esize=4)
        end do
    endif
    deallocate(ctmp)
  end subroutine
!===================================================
  subroutine mig_io_init()
    complex, allocatable, dimension(:) :: ctmp
    integer                            :: i, nloop
        allocate(ctmp(igrid%ax%n*igrid%ay%n))
        ctmp = (0,0)
        if (nostack)then
               nloop = igrid%az%n*unw
        else
               nloop = igrid%az%n
        endif

        do i=1, nloop
                call sep_write(ctmp, mig_, esize=4)
        end do
        deallocate(ctmp)
  end subroutine
!===================================================
  subroutine greens_io_init()
    complex, allocatable, dimension(:) :: ctmp
    integer                            :: i, nloop

    allocate(ctmp(ggrid%ax_surf%n*ggrid%ay_surf%n*ggrid%ax_modl%n*ggrid%ay_modl%n))
    ctmp = (0,0)
    nloop = ggrid%az_modl%n*ggrid%aw%n
    write(*,*)"Initializing the greens functions"
    do i=1, nloop
            call sep_write(ctmp, greens_, esize=8)
    end do
    deallocate(ctmp)
  end subroutine
!===================================================
  subroutine wfld_io_init(wfld_)
    complex, allocatable, dimension(:) :: ctmp
    integer                            :: i, nloop
    character(len=256)                 :: wfld_

    allocate(ctmp(vgrid%ax%n*vgrid%ay%n*dgrid%at%n/wfld_out_jt))
    ctmp = (0,0)
    nloop = vgrid%az%n 
    write(*,*)"Initializing the wfld_d"
    do i=1, nloop
!            call sep_write(ctmp, wfld_, esize=4)
    end do

    deallocate(ctmp)
  end subroutine

!===================================================
subroutine data_io_data(flag, data, index_aw )
    complex, dimension(:,:,:,:  ), pointer :: data
    character(len=4)                    :: flag
    integer                             :: index_aw
    integer, dimension(5)               :: n_cube, j_wind, n_wind, f_wind
    integer                             :: n_dim_cube, err
    character :: y

    n_dim_cube = 5
    n_cube = (/dgrid%agx%n , dgrid%agy%n   , dgrid%aw%n    , dgrid%asx%n    ,  dgrid%asy%n/)
    j_wind = (/1           , 1             , 1             , 1              ,           1/)
    n_wind = (/dgrid%agx%n , dgrid%agy%n   , 1             , dgrid%asx%n    ,  dgrid%asy%n/)
    f_wind = (/0           , 0             , index_aw-1    , 0              ,  0          /)

!    write(*,*) n_cube
!    write(*,*) n_wind
!    write(*,*) f_wind
!    write(*,*) j_wind
    
    if (flag == "rite") then
     !  write(0,*)'writing data frequency: ', index_aw
       err = srite_window(data_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,data)
    else if (flag == "reed") then
       data = 0.
       write(0,*)'reading data frequency: ', index_aw
       err = sreed_window(data_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,data)
    else
       call seperr("you have to specify the flag for data_io_data()")
    end if
  end subroutine
!##################################################################
subroutine data_io_tdata(flag, tdata, index_agx, index_agy, index_asx, index_asy)
    complex, dimension(: ), pointer :: tdata
    character(len=4)                    :: flag
    integer                             :: index_agx, index_agy, index_asx, index_asy
    integer, dimension(5)               :: n_cube, j_wind, n_wind, f_wind
    integer                             :: n_dim_cube, err
    character :: y

    n_dim_cube = 5
    n_cube = (/dgrid%agx%n , dgrid%agy%n   , dgrid%at%n    , dgrid%asx%n    , dgrid%asy%n/)
    j_wind = (/1           , 1             , 1             , 1              ,           1/)
    n_wind = (/1           , 1             , dgrid%at%n    , 1              ,           1/)
    f_wind = (/index_agx-1 , index_agy-1   , 0             , index_asx-1    ,index_asy-1 /)

    if (flag == "rite") then
       err = srite_window(data_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,tdata)
    else if (flag == "reed") then
       data = 0.
       err = sreed_window(data_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,tdata)
    else
       call seperr("you have to specify the flag for data_io_data()")
    end if
  end subroutine

!===================================================
  subroutine data_io_greens(greens, index_ws)
    integer                                   :: index_ws
    integer, dimension(6)                     :: n_cube, j_wind, n_wind, f_wind
    integer                                   :: n_dim_cube, err
    complex, dimension(:,:,:,:,:), pointer    :: greens

    n_dim_cube = 6
    n_cube = (/ggrid%ax_surf%n, ggrid%ay_surf%n, ggrid%ax_modl%n, ggrid%ay_modl%n, ggrid%az_modl%n, ggrid%aw%n/)
    j_wind = (/1              , 1              , 1              , 1              , 1              , 1         /)
    n_wind = (/ggrid%ax_surf%n, ggrid%ay_surf%n, ggrid%ax_modl%n, ggrid%ay_modl%n, ggrid%az_modl%n, 1         /)
    f_wind = (/0              , 0              , 0              , 0              , 0              , index_ws-1/)
!    write(*,*) n_cube
!    write(*,*) n_wind
!    write(*,*) f_wind
!    write(*,*) j_wind
    write(0,*)'writing greens freq.: ', index_ws, 'of ', ggrid%aw%n, 'at node', omp_get_thread_num();  
    err = srite_window(greens_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,greens)
    if (err /= 0) call seperr("data_io_greens: error writing wavefield")

  end subroutine

!===================================================
  subroutine mig_io_mig(mig, index_ws)
    integer                                   :: index_ws
    integer, dimension(4)                     :: n_cube, j_wind, n_wind, f_wind
    integer                                   :: n_dim_cube, err
    real   , dimension(:,:,:), pointer    :: mig

    n_dim_cube = 4
    n_cube = (/igrid%az%n, igrid%ax%n, igrid%ay%n, unw         /)
    j_wind = (/1         , 1         , 1         , 1           /)
    n_wind = (/igrid%az%n, igrid%ax%n, igrid%ay%n, 1           /)
    f_wind = (/0         , 0         , 0         , index_ws-1  /)
!    write(*,*) n_cube
!    write(*,*) n_wind
!    write(*,*) f_wind
!    write(*,*) j_wind
    

    write(0,*)'writing migration freq.: ', index_ws, 'of ', unw, 'at node', omp_get_thread_num();  
    err = srite_window(mig_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,mig)
    if (err /= 0) call seperr("mig_io_mig: error writing mig")

  end subroutine
  
!===================================================
  subroutine data_io_wfld(wfld_, wfld, index_z)
    integer                                :: index_z
    integer, dimension(4)                  :: n_cube, j_wind, n_wind, f_wind
    integer                                :: n_dim_cube, err
    real , dimension(:,:,:,:), allocatable :: wfld
!   real , dimension(:,:), allocatable:: wfld
    character(len=256)                     :: wfld_

    n_dim_cube = 4
    n_cube = (/vgrid%az%n, vgrid%ax%n, vgrid%ay%n, dgrid%at%n/wfld_out_jt /)
    j_wind = (/1         , 1         , 1         , 1                      /)
    n_wind = (/vgrid%az%n, vgrid%ax%n, vgrid%ay%n, dgrid%at%n/wfld_out_jt /)
    f_wind = (/0         , 0        ,  0         , 0                      /)
!    n_wind = (/1         , vgrid%ax%n, dgrid%at%n/wfld_out_jt /)
!    f_wind = (/index_z-1 , 0        ,  0                      /)

    wfld = 0.
    err = srite_window(wfld_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,wfld)
    if (err /= 0) call seperr("data_io_greens: error reading Green's function")
  end subroutine

  !===================================================
  subroutine data_io_vel(vel, index_z)
    integer                              :: index_z
    integer, dimension(3)                :: n_cube, j_wind, n_wind, f_wind
    real   , dimension(:,:), allocatable :: vel
    integer                              :: n_dim_cube, err

    n_dim_cube = 3
    n_cube = (/vgrid%az%n, vgrid%ax%n        , vgrid%ay%n        /)
    j_wind = (/1         , 1                 , 1                 /)
    n_wind = (/1         , pgrid%ax%n        , pgrid%ay%n        /)
    f_wind = (/index_z-1 , pgrid%index_xmin-1, pgrid%index_ymin-1/)

    vel = 0.
    err = sreed_window(vel_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,vel)
    if (err /= 0) call seperr("data_io_vel: error reading vel slice")
  end subroutine
  !===================================================
!===================================================
  subroutine data_io_ref(ref, index_z)
    integer                              :: index_z
    integer, dimension(3)                :: n_cube, j_wind, n_wind, f_wind
    real   , dimension(:,:), allocatable :: ref
    integer                              :: n_dim_cube, err

    n_dim_cube = 3
    n_cube = (/rgrid%az%n, rgrid%ax%n        , rgrid%ay%n        /)
    j_wind = (/1         , 1                 , 1                 /)
    n_wind = (/1         , pgrid%ax%n        , pgrid%ay%n        /)
    f_wind = (/index_z-1 , pgrid%index_xmin-1, pgrid%index_ymin-1/)

    ref = 0.
    err = sreed_window(ref_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,ref)
    if (err /= 0) call seperr("data_io_ref: error reading ref slice")
  end subroutine

  !===================================================
  subroutine data_io_wavelet(wavelet, index_ws)
    complex, dimension(:), pointer :: wavelet
    integer                        :: index_ws
    integer, dimension(1)          :: n_cube, j_wind, n_wind, f_wind
    integer                        :: n_dim_cube, err

    n_dim_cube = 1
    n_cube = dgrid%aw%n
    j_wind = 1 
    n_wind = 1
    f_wind = index_ws-1

    wavelet = 0.
    err = sreed_window(wavelet_,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,wavelet)
    if (err .lt. 0) call seperr("io_wavelet: trouble reading source wavelet")
  end subroutine

!#####################################################

subroutine propagation_grid()
    real    :: pgrid_xmin, pgrid_xmax
    real    :: pgrid_ymin, pgrid_ymax
    real    :: pgrid_zmin, pgrid_zmax

    pgrid%ax%d = vgrid%ax%d; 
    pgrid%ay%d = vgrid%ay%d; 
    pgrid%az%d = vgrid%az%d
 
    pgrid%gxmin = ggrid%surface_xmin; pgrid%gxmax = ggrid%surface_xmax
    pgrid%gymin = ggrid%surface_ymin; pgrid%gymax = ggrid%surface_ymax


    !################################################
    pgrid_zmin=vgrid%az%o !min_z
    pgrid_zmax=max_z
    !################################################
    pgrid_xmin = pgrid%gxmin; pgrid_xmax = pgrid%gxmax
    pgrid_ymin = pgrid%gymin; pgrid_ymax = pgrid%gymax

    !!.. add the padding on both x and y direction
    pgrid_xmin = pgrid_xmin - pgrid%nxpad_beg*pgrid%ax%d
    pgrid_xmax = pgrid_xmax + pgrid%nxpad_end*pgrid%ax%d

    pgrid_ymin = pgrid_ymin - pgrid%nypad_beg*pgrid%ay%d
    pgrid_ymax = pgrid_ymax + pgrid%nypad_end*pgrid%ay%d

    write(*,*)'pgridx    min max:',  pgrid_xmin, pgrid_xmax
    !!.. make sure the propagation grid is within the velocity model 
!    pgrid_xmin = max(pgrid_xmin, vgrid%ax%o)
   ! pgrid_xmax = min(pgrid_xmax, vgrid%ax%o + (vgrid%ax%n-1)*vgrid%ax%d)
    pgrid_xmin = vgrid%ax%o
    pgrid_xmax = vgrid%ax%o + (vgrid%ax%n-1)*vgrid%ax%d

    pgrid_ymin = max(pgrid_ymin, vgrid%ay%o)
    pgrid_ymax = min(pgrid_ymax, vgrid%ay%o + (vgrid%ay%n-1)*vgrid%ay%d)
    !################################################
    !!.. propagation grid 
    pgrid%ax%o = pgrid_xmin;     pgrid%ay%o = pgrid_ymin
    pgrid%ax%n = 1.5 + (pgrid_xmax - pgrid_xmin)/pgrid%ax%d
    pgrid%ay%n = 1.5 + (pgrid_ymax - pgrid_ymin)/pgrid%ay%d
    pgrid%az%n = min(1.5 + (pgrid_zmax - pgrid_zmin)/pgrid%az%d, 1.*vgrid%az%n)

    !!.. the first and last index of surface x and y on the propagation grid (for output use)
    ggrid%index_surpro_xmin = 1.5 + (ggrid%surface_xmin - pgrid%ax%o) / pgrid%ax%d
    ggrid%index_surpro_xmax = 1.5 + (ggrid%surface_xmax - pgrid%ax%o) / pgrid%ax%d

    ggrid%index_surpro_ymin = 1.5 + (ggrid%surface_ymin - pgrid%ay%o) / pgrid%ay%d
    ggrid%index_surpro_ymax = 1.5 + (ggrid%surface_ymax - pgrid%ay%o) / pgrid%ay%d

    !!....first and last index of geoph. and source
    pgrid%index_datpro_gxmin = 1.5 + (pgrid%gxmin  - pgrid_xmin ) / pgrid%ax%d  
    pgrid%index_datpro_gymin = 1.5 + (pgrid%gymin  - pgrid_ymin ) / pgrid%ay%d  

    pgrid%index_datpro_gxmax = 1.5 + (pgrid%gxmax  - pgrid_xmin ) / pgrid%ax%d  
    pgrid%index_datpro_gymax = 1.5 + (pgrid%gymax  - pgrid_ymin ) / pgrid%ay%d 

    !!....first and last index of propagation grid in the velocity model
    pgrid%index_xmin = 1.5 + (pgrid_xmin  - vgrid%ax%o ) / pgrid%ax%d 
    pgrid%index_ymin = 1.5 + (pgrid_ymin  - vgrid%ay%o ) / pgrid%ay%d  

    pgrid%index_xmax = 1.5 + (pgrid_xmax  - vgrid%ax%o ) / pgrid%ax%d  
    pgrid%index_ymax = 1.5 + (pgrid_ymax  - vgrid%ay%o ) / pgrid%ay%d  

    !!...first index of shots in propagation grid
    pgrid%index_datpro_sx = 1.5 + (dgrid%asx%o  - pgrid_xmin ) / pgrid%ax%d  
    pgrid%index_datpro_sy = 1.5 + (dgrid%asy%o  - pgrid_ymin ) / pgrid%ay%d  
    write(*,*)'pgridx pad min max:',  pgrid_xmin, pgrid_xmax
    write(*,*)'index pgrid gx:', pgrid%gxmin, pgrid%gxmax, pgrid%index_datpro_gxmin, pgrid%index_datpro_gxmax  
  
  end subroutine
!######################################################################
  subroutine velocity_grid()
    real  :: vgrid_xmin, vgrid_xmax, vgrid_ymin, vgrid_ymax
    real  :: pgrid_xmin, pgrid_xmax, pgrid_ymin, pgrid_ymax

    vgrid_xmin = vgrid%ax%o
    vgrid_ymin = vgrid%ay%o
    vgrid_xmax = vgrid%ax%o + (vgrid%ax%n-1)*vgrid%ax%d
    vgrid_ymax = vgrid%ay%o + (vgrid%ay%n-1)*vgrid%ay%d

    pgrid_xmin = pgrid%ax%o
    pgrid_ymin = pgrid%ay%o
    pgrid_xmax = pgrid%ax%o + (pgrid%ax%n-1)*pgrid%ax%d
    pgrid_ymax = pgrid%ay%o + (pgrid%ay%n-1)*pgrid%ay%d

    if(pgrid_xmin < vgrid_xmin) then
       write(0,*)"pgrid_xmin, vgrid_xmin", pgrid_xmin, vgrid_xmin
       call seperr("propagation grid out of velocity grid: pgrid_xmin < vgrid_xmin")
    end if
    if(pgrid_xmax > vgrid_xmax) then
       write(0,*)"pgrid_xmax, vgrid_xmax", pgrid_xmax, vgrid_xmax
       call seperr("propagation grid out of velocity grid: pgrid_xmax > vgrid_xmax")
    end if
    if(pgrid_ymin < vgrid_ymin) then
       write(0,*)"pgrid_ymin, vgrid_ymin", pgrid_ymin, vgrid_ymin
       call seperr("propagation grid out of velocity grid: pgrid_ymin < vgrid_ymin")
    end if
    if(pgrid_ymax > vgrid_ymax) then
       write(0,*)"pgrid_ymax, vgrid_ymax", pgrid_ymax, vgrid_ymax
       call seperr("propagation grid out of velocity grid: pgrid_ymax > vgrid_ymax")
    end if

    !!.. computing the first and the last x,y index of the input velocity on the whole velocity grid
    vgrid%index_provel_xmin = 1.5 + (pgrid%ax%o - vgrid%ax%o)/vgrid%ax%d
    vgrid%index_provel_ymin = 1.5 + (pgrid%ay%o - vgrid%ay%o)/vgrid%ay%d
    vgrid%index_provel_xmax = vgrid%index_provel_xmin + vgrid%ax%n - 1
    vgrid%index_provel_ymax = vgrid%index_provel_ymin + vgrid%ay%n - 1

  end subroutine

 !#################################################
  subroutine k2_map()

    if(associated(kx2)) deallocate(kx2)
    if(associated(ky2)) deallocate(ky2)

    allocate(kx2(pgrid%ax%n))
    allocate(ky2(pgrid%ay%n))

    call k_map(kx2, pgrid%ax%n, pgrid%ax%d)
    call k_map(ky2, pgrid%ay%n, pgrid%ay%d)

  end subroutine

!###############################################
  subroutine k_map(k,n,d)
    real,  dimension(:) :: k
    integer, intent(in) :: n
    real,    intent(in) :: d
    integer             :: i,m
    k=0.
    if(n/=1) then
       m = n/2
       do i=1,m+1; k(i)=i-  1; end do
       do i=m+2,n; k(i)=i-n-1; end do
    else
       k=0. !! make sure k=0 if n=1
    end if

    k = k / (n*d)      !! 1/length
    k = 2 * PI * k     !! radians
    k = k * k          !! k squared
    
  end subroutine
!###############################################
  subroutine vel_alloc(vel_slice,ind_vref,vref)
     real,    dimension(:),        allocatable :: vref
     real,    dimension(:,:),      allocatable :: vel_slice
     integer, dimension(:,:),      allocatable :: ind_vref
     if (allocated(vel_slice)) deallocate(vel_slice); allocate(vel_slice(pgrid%ax%n,pgrid%ay%n))
     if (allocated(ind_vref )) deallocate(ind_vref ); allocate( ind_vref(pgrid%ax%n,pgrid%ay%n))
     if (allocated(vref)) deallocate(vref); allocate(vref(nref))
     vel_slice=0
     ind_vref=0
     vref=0
  end subroutine
!###############################################
  subroutine data_radians_calc()
  !!.. convert frequencies to radians
  if (.not.radian) then
     dgrid%aw%d = dgrid%aw%d * 2.0 * PI
     dgrid%aw%o = dgrid%aw%o * 2.0 * PI
  end if
  end subroutine

!###############################################
 subroutine wfld_out(wfld,wfld_)
        complex    , dimension(:,:,:,:,:), pointer     :: wfld
        complex    , dimension(:      ), allocatable :: tmp_wfld
        real       , dimension(:,:,:,:), allocatable :: tmp_wfld2
        real       , dimension(:,:,:,:), allocatable :: tmp_wfld3
        character (len=256) :: wfld_

        if (allocated(tmp_wfld )) deallocate (tmp_wfld ); allocate(tmp_wfld(dgrid%at%n))
        if (allocated(tmp_wfld2)) deallocate (tmp_wfld2); allocate(tmp_wfld2(vgrid%az%n,vgrid%ax%n,vgrid%ay%n,dgrid%at%n)) 
        if (allocated(tmp_wfld3)) deallocate (tmp_wfld3); allocate(tmp_wfld3(vgrid%az%n,vgrid%ax%n,vgrid%ay%n,dgrid%at%n/wfld_out_jt))

        nnode =omp_get_num_procs(); thread=omp_get_thread_num(); 
        nodes=floor(nnode*.98)
        call omp_set_num_threads(nodes); write(*,*) 'using', nodes, 'of' , nnode
        !$OMP parallel do private(it,itt,ix,iz,inode,tmp_wfld) schedule(dynamic,1) ordered 
        do iz=1,vgrid%az%n
        if (allocated(tmp_wfld )) deallocate (tmp_wfld ); allocate(tmp_wfld(dgrid%at%n))
                inode=omp_get_thread_num();  
                write(*,*)'iz.:', iz ,' of ',vgrid%az%n,' @ node ', inode         
                do iy=1,vgrid%ay%n 
                do ix=1,vgrid%ax%n 
                        tmp_wfld=wfld_out_d(iz,ix,iy,:,1); call FFT1DT(tmp_wfld, -1); tmp_wfld2(iz,ix,iy,:)=tmp_wfld
                enddo
                enddo
        enddo
        !$OMP end parallel do
        itt=1; 
        do it=1,dgrid%at%n-1,wfld_out_jt; 
                tmp_wfld3(:,:,:,itt)=tmp_wfld2(:,:,:,it); 
                write(*,*)'saved',maxval(abs(tmp_wfld2(:,:,:,it))),maxval(abs(tmp_wfld3(:,:,:,itt))) 
                itt=itt+1; 
        end do
                
!       call data_io_wfld (wfld_d_, tmp_wfld3, iz)

        call sep_write(tmp_wfld3, wfld_, esize=4)

        write(*,*)'written wavefields for one shot successfully'
        deallocate(tmp_wfld,tmp_wfld2,tmp_wfld3)
 end subroutine
!###############################################


  end module
