module OW_xwell1

use OWvel
use OW_type
use OW_image
use OW_fftw
use wei_myfft_mod
use OW_grabvel
use OW_parms
use OW_taper
use OW_wvfldIO2

implicit none
contains

!---------------------------------------------------------
subroutine owxwell(ireal,ixs,source,vel,taper) 

integer						:: iw,iz,ikx,iref,ixw,isp,ixs,ireal
integer						:: nv,ihso,nwb,node,inode
integer,allocatable					:: indref(:)
real						:: sqa,f,g
real						:: minv, maxv,vel(:,:),taper(:)
real, allocatable					:: vref(:),sqa0(:)
complex, allocatable				:: ws(:),gs0(:),gs(:,:),tfld(:),GSS(:,:)
complex						:: source(:)

integer, external 					:: omp_get_thread_num,omp_get_num_threads,omp_get_num_procs
!integer, external 					:: omp_set_num_threads

node=omp_get_num_threads()
node=8

!if (exist_file("GS")) then
!   allocate(GSS(img%z%n,kxm%n))
!   call to_history("n1",img%z%n,"GS");call to_history("d1",img%z%d,"GS");call to_history("o1",img%z%o,"GS")
!   call to_history("n2",img%xm%n,"GS");call to_history("d2",img%xm%d,"GS");call to_history("o2",img%xm%o,"GS")
!   call to_history("n3",sou%w%n,"GS");call to_history("d3",sou%w%d/(2*pi),"GS");call to_history("o3",sou%w%o/(2*pi),"GS")  
!   call to_history("esize",8,"GS")
!   allocate(GRR(img%z%n,kxm%n))
!   call to_history("n1",img%z%n,"GR");call to_history("d1",img%z%d,"GR");call to_history("o1",img%z%o,"GR")
!   call to_history("n2",img%xm%n,"GR");call to_history("d2",img%xm%d,"GR");call to_history("o2",img%xm%o,"GR")
!   call to_history("n3",sou%w%n,"GR");call to_history("d3",sou%w%d/(2*pi),"GR");call to_history("o3",sou%w%o/(2*pi),"GR")  
!   call to_history("esize",8,"GR")
!end if
!--------------------------------------------------

isp=int(((ixs-1)*sou%x%d+sou%x%o-sou%h%o)/sou%h%d)+1.5

write(0,*) isp

!--------------------------------------------------
! Setting FFTW plans
allocate(tfld(kxm%n));tfld=0
call owfftw(tfld)
deallocate(tfld)

!!$OMP parallel do private (inode,&
!                             iw,&
!                             iz,&
!                            ikx,&
!                           iref,&
!                            ixw,&
!                             ws,&
!                             gs,&
!                            gs0,&
!                           vref,&
!                           sqa0,&
!                            sqa,&
!                           minv,&
!                           maxv,&
!                           tfld,&
!                         indref,&
!                              f,&
!                              g) schedule(dynamic,1) num_threads(8)
! !$OMP DO schedule(dynamic,1) 

omega:do iw = 1,sou%w%n-subf

        if(allocated(ws)) deallocate(ws); allocate(ws(kxm%n)); ws=0.

        inode=omp_get_thread_num()+1; 

        ws(isp)=source(iw)

        depth:do iz=nwb+1,img%z%n

            if(allocated(tfld)) deallocate(tfld); allocate(tfld(kxm%n));tfld=0

!!--------------------------------------------------
!! Compute vref
!!
            if (maxval(vel(iz,:)).ge.vsal) then
              if(allocated(vref)) deallocate(vref); allocate(vref(nref+1)); vref=0
              if(allocated(gs)) deallocate(gs); allocate(gs(kxm%n,nref+1)); gs=0
              call owvref(nref,vel(iz,:),vref(1:nref),minv,maxv)
              vref(nref+1)=maxval(vel(iz,:))
            else if(maxval(vel(iz,:)).lt.vsal) then
              if(allocated(vref)) deallocate(vref); allocate(vref(nref)); vref=0
              if(allocated(gs)) deallocate(gs); allocate(gs(kxm%n,nref)); gs=0
              call owvref(nref,vel(iz,:),vref,minv,maxv)
            end if


!!---------------------------------------------------
!!
!            tfld=ws; call sfftw_execute(planF); ws=tfld/sqrt(1.*kxm%n);

            tfld=ws; call FFT1DX(tfld,1); ws=tfld;

           vxz:if (minv .eq. maxv .and. maxv.eq.vref(size(vref))) then

               if(allocated(gs0)) deallocate(gs0); allocate(gs0(kxm%n));gs0=ws
!
!      Operations in the wavenumber domain
!      phase shift
                wavenumber1:do ikx = 1,kxm%n
                   sqa = sou%w%all(iw)**2 / minv**2 - kxm%all(ikx)**2
                   evanescent1:if (sqa > 0) then
                       gs0(ikx)=gs0(ikx)*exp(cmplx(0,v%z%d*sqrt(sqa)))
                   else
                       gs0(ikx)=gs0(ikx)*exp(-sqrt(-sqa)*v%z%d)
                   end if evanescent1
                end do wavenumber1

            tfld=gs0; call FFT1DX(tfld,-1); ws=tfld;

!             tfld=gs0; call sfftw_execute(planI); ws=tfld/sqrt(1.*kxm%n)
!
            else
               if(allocated(sqa0)) deallocate(sqa0); allocate(sqa0(kxm%n));sqa0=0.
                velref: do iref=1,size(vref)                 
                    gs(:,iref)=ws;
                    wavenumber2: do ikx = 1,kxm%n
!      phase shift
                        sqa = sou%w%all(iw)**2 / vref(iref)**2 - kxm%all(ikx)**2;
                        evanescent2: if (sqa > 0) then
                            gs(ikx,iref)=gs(ikx,iref)*exp(cmplx(0,v%z%d*sqrt(sqa)));
                        else
                            gs(ikx,iref)=gs(ikx,iref)*exp(-sqrt(-sqa)*v%z%d)
                        end if evanescent2
                    end do wavenumber2

! Operations in x            

            tfld=gs(:,iref); call FFT1DX(tfld,-1); gs(:,iref)=tfld;

!             tfld=gs(:,iref); call sfftw_execute(planI); gs(:,iref)=tfld/sqrt(1.*kxm%n)

!      split step
                      gs(:,iref) = gs(:,iref) * exp(cmplx(0,v%z%d*sou%w%all(iw))*(1./vel(iz,:)-1./vref(iref)));
                end do velref

!                do iref=1,nref-1
!                interpwvfld: where (vel(iz,:) >= vref(iref) .and. vel(iz,:) < vref(iref+1))
!                     ws=gs(:,iref)
!                     end where interpwvfld
!                end do
!                where (vel(iz,:) >= vref(nref))
!                      ws=gs(:,nref)
!                end where

                do ikx=1,kxm%n
                    if(allocated(indref)) deallocate(indref); allocate(indref(1));indref=0
                    indref=maxloc(vref, mask = vref.lt.vel(iz,ikx)) 
                    if (indref(1).ne.0.and.indref(1).ne.size(vref)-1) then
                       f=(vel(iz,ikx) - vref(indref(1)))/vel(iz,ikx)
                       g=1-f
                       ws(ikx)=g*gs(ikx,indref(1))+f*gs(ikx,indref(1)+1)
                    else
                       ws(ikx)=gs(ikx,indref(1)+1)
                    end if
                end do
           end if vxz
        ws=ws*taper

        end do depth

!
!        !$OMP CRITICAL
           call owwvfld_io(ws,iw,ixs)!,ireal)
!        !$OMP END CRITICAL
          if (mod(iw,20)==1) write(0,5)  "DONE FREQUENCY # ",iw," OUT OF ",sou%w%n-subf," / SHOT # ",ixs
5                            format(A17,i5,2(A8,i5))


end do omega
!!$OMP END PARALLEL DO
end subroutine
end module
