module oneway_mod
use sep
use myfft_mod
use ctaper_mod
use sort_mod
use vel_mod
use findsr_mod
!use global_defi_mod
!use global_mod

implicit none

contains
!==============================================
subroutine oneway(data,x, xs, xr, z, t, v, vnum, source, nxpad, imin_z, imax_z, jz, jx, jxr, jxs, reflect,restart,tape,opt,greens)

real    , dimension(:,:,:)      , pointer :: data
complex , dimension(:,:,:,:) ,allocatable :: G_hess
complex , dimension(:,:,:)      , pointer :: fdata
complex , dimension(:,:)        , pointer :: img
real    , dimension(:,:)        , pointer :: v
real    , dimension(:)       ,allocatable :: vtmp
real    , dimension(:,:)        , pointer :: reflect
real    , dimension(:,:)     ,allocatable :: vel
complex , dimension(:,:,:)   ,allocatable :: GG1
complex , dimension(:,:,:)   ,allocatable :: G1
complex , dimension(:,:)     ,allocatable :: green, green1
complex , dimension(:)          , pointer :: source
complex , dimension(:)       ,allocatable :: wavelet,wave,fwavelet
real    , dimension(:)          , pointer :: x, z, xs, xr, t
real    , dimension(:)          , pointer :: wall,kxall
real    , dimension(:)       ,allocatable :: gpts, gpts1,gpts2,igpts
real    , dimension(:)       ,allocatable :: vref
integer , dimension(:)       ,allocatable :: ind_vref        
complex , dimension(:)       ,allocatable :: tmp_wfld        
complex , dimension(:,:,:,:) ,allocatable :: wfld_d        
real    , dimension(:)       ,allocatable :: sqa        
real    , dimension(:,:,:)   ,allocatable :: tmp_wfld2        
real    , dimension(:,:,:)   ,allocatable :: tmp_wfld3        

real                    :: dx, dz, ox, oz  
real                    :: dxr, dxs, oxr, oxs 
real                    :: dt, ot 
real                    :: dw, ow, dkx, okx  
real                    :: tape

integer                 :: ix, iz, iw, it, itt, ik, iv
integer                 :: nx, nz
integer                 :: nxr, nxs
integer                 :: nw, nkx, nt
integer                 :: nxv
integer                 :: vnum, nref
integer                 :: nxpad, imax_z, imin_z
integer                 :: jz, jx, jxr, jxs
integer                 :: restart, what, opt
integer                 :: nsour2
integer                 :: iww, niw
integer                 :: i, j, k, ii, jj, kk
integer                 :: ir, is
integer                 :: ngpts
integer                 :: TID, OMP_GET_THREAD_NUM

integer                 :: ixp, ixg, ixs, ixr, iixs, iixr
integer                 :: sign
real                    :: pi
INTEGER, EXTERNAL       :: omp_get_num_procs 
INTEGER                 :: inode, nnode, i1, i2, thread
real                    :: xxs, xxr 

logical                 :: greens
character(len=256)      :: sou1_="wavelet1", wfld_d_="wfld_d", ghes_="ghes"
!=============================================================
pi=3.1415937; sign=1;    
!=============================================================
!=thread stuff
nnode=omp_get_num_procs(); thread=omp_get_thread_num(); write(*,*)'nnonde,thread'; write(*,*)nnode,thread
call omp_set_num_threads(floor(nnode*.98))
!call omp_set_num_threads(floor(nnode*1.))
!call omp_set_num_threads(floor(1.0))
!=================  Initial computations  ====================
!=model grid....vel grid...image grid...
        nx=size(x)  ;   nt=size(t)  ;   nz=size(z)  ;      
        dx=x(2)-x(1);   dt=t(2)-t(1);   dz=z(2)-z(1);    
        ox=x(1)     ;   ot=t(1)     ;   oz=z(1)     ;
!=============================================================
!=source/receiver axis
        nxs=size(xs) ; nxr=size(xr) ;
        if (nxs > 1)then; dxs=xs(2)-xs(1); oxs=xs(1); else; dxs=1; oxs=0.; endif;
        if (nxr > 1)then; dxr=xr(2)-xr(1); oxr=xr(1); else; dxr=1; oxr=0.; endif;
!=============================================================
!=Frequency axes
!=temporal freq.
        nw = floor(nt/2.); dw = 1./(nt*dt); ow=0;  allocate(wall(nw));
        do iw = 1, nw; wall(iw) = ow  + (iw - 1) * dw; enddo; wall = 2*pi * wall;
!==
!=spatial freq.
        nkx = nx+2*nxpad; dkx = 1./(nkx*dx); okx = 1./(2*dx)-dkx; allocate(kxall(nkx));
        do ik=1,nkx; kxall(ik) = mod(okx+(ik-1)*dkx,1./dx)-okx; enddo; kxall = 2*pi * kxall;
!=============================================================
!=propagation grid...
        nxv=size(v(1,:));
        if (nxv<nkx)then 
            allocate (vel(nz,nx+2*nxpad)) 
            do iv=1,nxpad; vel(:,iv)=v(:,1); vel(:,nxpad+nxv+iv)=v(:,size(v(1,:))); enddo; 
            vel(:,nxpad+1:nxpad+nxv)=v;
        else
            allocate (vel(nz,nx)); vel=v;
        endif;
!=============================================================
!=pad and position wavelet velocities
    allocate(wavelet(nt)); wavelet(1:nt)=source;
    call FFT1DT(wavelet, 1); allocate(wave(nw)); wave=wavelet(1:nw);
    call ctaper1d(wave,wave,dw, nw, 0, 0.1, 2);
!    wave=wave*exp(cmplx(0.,wall*nw*dt))


write(*,*)'writing source wavelet.....'
call auxinout(sou1_); 
call auxputch("n1    ","i", nw,sou1_);call auxputch("d1   ","r",dw,sou1_); 
call auxputch("o1", "r", ow , sou1_); 
call auxputch("label1","s","t",sou1_);call auxputch("esize","i",8,sou1_); 
call srite(sou1_,wave , 8*nw)


!=compute the useful freq. 
    iww=0;
    do niw=floor(nw*0.1),nw-1
        if (abs(wave(niw))<=maxval(abs(wave))*.005)then; iww=iww+1; if (iww>=floor(nw*0.02)) exit; endif
    enddo
!=============================================================
!=get green surf. functions points
    allocate(gpts1(nxr),gpts2(nxr))
    gpts1=xr;  gpts2=xr; j=0;
    do is=1,nxs
        do jj=1,nxr
                if (xs(is)==gpts1(jj))then; j=0; exit; else;  j=1; endif
        enddo
        if (j/=0)then;
            if(allocated(gpts1))deallocate(gpts1); allocate(gpts1(size(gpts2)));   gpts1=gpts2;
            if(allocated(gpts2))deallocate(gpts2); allocate(gpts2(size(gpts1)+1)); 
            gpts2(1:size(gpts1))=gpts1; gpts2(size(gpts1)+1)=xs(is); 
        endif
   enddo
   ngpts=size(gpts2); allocate(igpts(ngpts)); call sort(gpts2,igpts,ngpts,1); allocate(gpts(ngpts)); gpts=gpts2 
!=============================================================
!=initialize data   
allocate (GG1(nt,nxr,nxs), wfld_d(nz,nx,nt,1)); GG1=0; wfld_d=0;
if (greens) then; write(*,*)'Plan to store greens'; if (allocated(G_hess)) deallocate(G_hess); allocate(G_hess(nz,nx,ceiling(niw/4.),ngpts)); endif
G_hess=cmplx(0.,0.);
!=============================================================
write(*,*)'Finished prelim stuff....starting the frequency loop.....'
write(*,*)'======= ====== ========= ======== ==== == = = = = == = ...'
!$OMP parallel do private(inode,iw,ixg,ixp,ix,iz,iv,ik,nref,sqa,tmp_wfld,G1,vtmp,green,green1,vref,ind_vref,iixs,iixr,ixs,ixr,xxs,xxr) schedule(dynamic,1) ordered  
do iw = 1,niw;      ! frequency loop
        inode=omp_get_thread_num();  
        write(*,*)'freq. num.: ', iw ,' of ',niw,' inode ', inode;           
        if(abs(wave(iw))<=maxval(abs(wave))*.001)cycle 
        if (allocated(G1)) deallocate(G1); allocate(G1(nz,nkx,ngpts));G1=0
        xxs = oxs + (2-1)*dxs; call findsr(ixs,xxs,gpts);
        do ixg=ixs,ixs !1,ngpts;    ! shot/geoph. pos. loop
                ixp=floor(gpts(ixg)-ox)/dx+1;  
!               write(*,*) 'ixp, ixg, gpts(ixg)',ixp, ixg, gpts(ixg)
                if (ixp>nkx) then; call seperr('source out of propagation grid'); endif;  
                G1(1,ixp+nxpad,ixg)=wave(iw); !%plot(1:nkx,real(G1(1,:,ixg)));axis 'tight';
                do iz=2,nz;     ! depth loop
                        if (allocated(vtmp)) deallocate(vtmp); allocate(vtmp(nkx))
                        if (allocated(ind_vref)) deallocate(ind_vref); allocate(ind_vref(nkx))
                        if (allocated(vref)) deallocate(vref); allocate(vref(vnum))
                        vtmp=vel(iz,:); 
                        call  velocity(vref,vtmp, ind_vref, vnum, nref); 
                        if (allocated(green )) deallocate(green ); allocate(green (nref,nkx))
                        if (allocated(green1)) deallocate(green1); allocate(green1(nref,nkx))
                        if (allocated(sqa   )) deallocate(sqa   ); allocate(sqa   (nref)    )
                        if (allocated(tmp_wfld)) deallocate(tmp_wfld); allocate(tmp_wfld(nkx))
                        tmp_wfld=G1(iz-1,:,ixg);
                        call FFT1DX(tmp_wfld, 1); G1(iz,:,ixg)=tmp_wfld;
                        do iv=1, nref; green(iv,:)=G1(iz,:,ixg);enddo
!% operations in the wavenumber domain
                        do iv=1,nref                    ! vel loop    
                           do ik = 1,nkx                ! wavenumber loop
                              sqa(iv) = wall(iw)**2 / vref(iv)**2 - kxall(ik)**2;  !phase shift 
                              if (sqa(iv) > 0) then                                !evanescence
                                 green(iv,ik)=green(iv,ik)*exp(cmplx(0,sign*dz*sqrt(sqa(iv))));
                              else
                                 green(iv,ik)=0.;
                             endif
                           enddo                           !end wavenumber loop
!% split-step         
                           tmp_wfld=green(iv,:); call FFT1DX(tmp_wfld, -1); green(iv,:)=tmp_wfld
                           green(iv,:)=green(iv,:)*exp(cmplx(0,sign*dz*wall(iw))*(1/vel(iz,:)-1/vref(iv)));!spl.step 
                           tmp_wfld=0; tmp_wfld=green(iv,:); call ctaper1d(tmp_wfld,tmp_wfld,dx, nkx, nxpad*0, tape, opt); green1(iv,:)=tmp_wfld
                        enddo
!                        
                        do ix=1,nkx;
                           G1(iz,ix,ixg)=green1(ind_vref(ix),ix);
                        enddo
                enddo        ! end depth loop
        enddo   ! end shot/receiver loop
!                    
!% convolve greens + reflectivity
        do iixs=1,nxs; xxs = oxs + (iixs-1)*dxs; call findsr(ixs,xxs,gpts);
                do iixr=1,nxr; xxr = oxr + (iixr-1)*dxr;  call findsr(ixr,xxr,gpts);
                     do ix=1,nx; do iz=1,nz;
                        GG1(iw,iixr,iixs)=GG1(iw,iixr,iixs)+G1(iz,ix+nxpad,ixs)*G1(iz,ix+nxpad,ixr)*reflect(iz,ix);   ! convolve to generate data
                     enddo;      enddo 
                enddo;
        enddo
        xxs = oxs + (2-1)*dxs; call findsr(ixs,xxs,gpts);
        wfld_d(:,:,iw,1)=G1(:,1+nxpad:nx+nxpad,ixs);
        !if (greens .and. mod(iw,4)==0) G_hess(:,:,floor(iw/4.),:)=G1(:,1+nxpad:nx+nxpad,:) 
if (allocated(G1))       deallocate(G1);
if (allocated(green ))   deallocate(green );
if (allocated(green1))   deallocate(green1);
if (allocated(sqa   ))   deallocate(sqa   );
if (allocated(tmp_wfld)) deallocate(tmp_wfld);
if (allocated(vtmp))     deallocate(vtmp); 
if (allocated(ind_vref)) deallocate(ind_vref); 
if (allocated(vref))     deallocate(vref);
end do                  ! end frequency loop
!$OMP end parallel do

if (greens) then
        niw=floor(niw/5.)
        call auxinout(ghes_)
        call auxputch("n1", "i", nz , ghes_); call auxputch("d1", "r", dz  , ghes_); call auxputch("o1", "r", oz , ghes_); 
        call auxputch("n2", "i", nx , ghes_); call auxputch("d2", "r", dx  , ghes_); call auxputch("o2", "r", ox , ghes_); 
        call auxputch("n3", "i", niw, ghes_); call auxputch("d3", "r", dw*5, ghes_); call auxputch("o3", "r", ow , ghes_); 
        call auxputch("n4", "i", nxr, ghes_); call auxputch("d4", "r", dxr , ghes_); call auxputch("o4", "r", oxr, ghes_); 
        call auxputch("label1", "s", "x"  , ghes_)
        call auxputch("label2", "s", "z"  , ghes_)
        call auxputch("label3", "s", "w"  , ghes_)
        call auxputch("label4", "s", "nxr", ghes_)
        call auxputch("esize" , "i",  8   , ghes_)
call srite(ghes_, G_hess , 8*nx*nz*nxr*niw)
if (allocated(G_hess)) deallocate(G_hess);
end if 

if (allocated(tmp_wfld2)) deallocate(tmp_wfld2); allocate(tmp_wfld2(nx,nz,nt))
if (allocated(tmp_wfld)) deallocate(tmp_wfld); allocate(tmp_wfld(nt));tmp_wfld=0

!call omp_set_dynamic(.false.)
!$OMP parallel do private(ix,iz,inode,tmp_wfld) schedule(dynamic,1) ordered 
do iz=1,nz 
        if (allocated(tmp_wfld)) deallocate(tmp_wfld); allocate(tmp_wfld(nt));tmp_wfld=0
        inode=omp_get_thread_num();  
        write(*,*)'iz.: ', iz ,' of ',nz,' inode ', inode;           
        do ix=1,nx
                tmp_wfld=wfld_d(iz,ix,:,1); call FFT1DT(tmp_wfld, -1); tmp_wfld2(iz,ix,:)=tmp_wfld 
        enddo
enddo
!$OMP end parallel do

if (allocated(tmp_wfld3)) deallocate(tmp_wfld3); allocate(tmp_wfld3(nz,nx,ceiling(nt/10.)))

itt=1
do it=1,nt,10; tmp_wfld3(:,:,itt)=tmp_wfld2(:,:,it); itt=itt+1; end do

call auxinout(wfld_d_)
call auxputch("n1", "i", nz , wfld_d_); call auxputch("n2", "i", nx , wfld_d_); call auxputch("n3", "i", ceiling(nt/10.), wfld_d_); 
call auxputch("d1", "r", dz , wfld_d_); call auxputch("d2", "r", dx , wfld_d_); call auxputch("d3", "r", dt*10 , wfld_d_);
call auxputch("o1", "r", oz , wfld_d_); call auxputch("o2", "r", ox , wfld_d_); call auxputch("o3", "r", ot , wfld_d_);
call auxputch("label1", "s", "nz" , wfld_d_)
call auxputch("label2", "s", "nx" , wfld_d_)
call auxputch("label3", "s", "t" , wfld_d_)
call auxputch("esize" , "i",  4  , wfld_d_)
call srite(wfld_d_,tmp_wfld3 , 4*ceiling(nt/10.)*nz*nx)
STOP
if (allocated(tmp_wfld)) deallocate(tmp_wfld); allocate(tmp_wfld(nt));tmp_wfld=0
if (associated(data)) deallocate(data); allocate(data(nt,nxr,nxs));data=0
do iixs=1,nxs
        do iixr=1,nxr
                tmp_wfld=GG1(:,iixr,iixs); 
                call FFT1DT(tmp_wfld, -1);
                data(:,iixr,iixs)=real(tmp_wfld)
        enddo
enddo


if (allocated(GG1)) deallocate(GG1);
if (allocated(vel)) deallocate(vel); 
if (allocated(wavelet)) deallocate(wavelet); 
if (allocated(wave)) deallocate(wave); 
if (allocated(gpts1)) deallocate(gpts1); 
if (allocated(gpts2)) deallocate(gpts2); 
if (greens) then
        niw=floor(niw/5.)
        call auxinout(ghes_)
        call auxputch("n1", "i", nz , ghes_); call auxputch("d1", "r", dz  , ghes_); call auxputch("o1", "r", oz , ghes_); 
        call auxputch("n2", "i", nx , ghes_); call auxputch("d2", "r", dx  , ghes_); call auxputch("o2", "r", ox , ghes_); 
        call auxputch("n3", "i", niw, ghes_); call auxputch("d3", "r", dw*5, ghes_); call auxputch("o3", "r", ow , ghes_); 
        call auxputch("n4", "i", nxr, ghes_); call auxputch("d4", "r", dxr , ghes_); call auxputch("o4", "r", oxr, ghes_); 
        call auxputch("label1", "s", "x"  , ghes_)
        call auxputch("label2", "s", "z"  , ghes_)
        call auxputch("label3", "s", "w"  , ghes_)
        call auxputch("label4", "s", "nxr", ghes_)
        call auxputch("esize" , "i",  8   , ghes_)
call srite(ghes_, G_hess , 8*nx*nz*nxr*niw)
if (allocated(G_hess)) deallocate(G_hess);
end if 

if (allocated(wfld_d)) deallocate(wfld_d);
if (allocated(tmp_wfld)) deallocate(tmp_wfld);
if (allocated(tmp_wfld2)) deallocate(tmp_wfld2);
if (allocated(tmp_wfld3)) deallocate(tmp_wfld3);

end subroutine
end module
