module OW_ctrillwei

use sep
use OWvelo_areal
use OW_type
use OW_fftw_areal
!  use OW_grabvel_areal
use OW_parms_areal

implicit none

integer,private					:: iwa,nwb
real,private					:: wb

contains

!---------------------------------------------------------
subroutine owctrlwei_init(iwa_in,wb_in,nwb_in)

integer						:: iwa_in,nwb_in
real						:: wb_in

iwa=iwa_in
wb=wb_in
nwb=nwb_in

end subroutine
!---------------------------------------------------------
subroutine owctrlwei(taper,reflect,source,vsal,vel,GG)

integer						:: iz,ikx,iref,ixw,nout
real						:: sqa,f,g
real						:: minv, maxv,vel(:,:),vsal,reflect(:,:),taper(:)
real, allocatable				:: vref(:),sqa0(:)
complex, allocatable				:: ws(:),gs0(:),gs(:,:)
complex, allocatable				:: tfld(:),refl(:)
complex						:: source(:,:,:),GG(:)

if (allocated(ws)) deallocate(ws);allocate(ws(img%xm%n));	ws=0
if (allocated(gs)) deallocate(gs);allocate(gs(img%xm%n,nref));	gs=0
if (allocated(gs0)) deallocate(gs0);allocate(gs0(img%xm%n));	gs0=0
if (allocated(sqa0)) deallocate(sqa0);allocate(sqa0(img%xm%n));	sqa0=0
if (allocated(vref)) deallocate(vref);allocate(vref(nref));	vref=0
if (allocated(tfld)) deallocate(tfld);allocate(tfld(img%xm%n));	tfld=0
if (allocated(refl)) deallocate(refl);allocate(refl(img%xm%n));	refl=0

!--------------------------------------------------
! Setting FFTW plans
!
call owfftw(tfld)

        depth:do iz=img%z%n,nwb+1,-1
!--------------------------------------------------
! Compute vref
!
            call owvref(nref,vsal,vel(iz,:),vref,minv,maxv)
!!!!!!!!!
!ATENCAO MUDAR DEPOIS DE TESTAR
!!!!!!!!
!minv=10000
!6096
!maxv=minv
!!!!!!
!!!!!
            tfld=cmplx(reflect(iz,:),0); call sfftw_execute(planF); refl=tfld/sqrt(1.*img%xm%n)
            tfld=ws; call sfftw_execute(planF); ws=tfld/sqrt(1.*img%xm%n)
            ws=ws+refl*source(1,iwa,1)
            
           vxz:if (minv .eq. maxv) then
                gs0=ws;

!      Operations in the wavenumber domain
!      phase shift

                wavenumber1:do ikx = 1,kxm%n

                   sqa = rec%w%all(iwa)**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 sfftw_execute(planI); ws=tfld/sqrt(1.*img%xm%n)

            else
                velref: do iref=1,nref                 
                    gs(:,iref)=ws
                    wavenumber2: do ikx = 1,kxm%n
!      phase shift
                        sqa = rec%w%all(iwa)**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 sfftw_execute(planI); gs(:,iref)=tfld/sqrt(1.*img%xm%n)

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

                do ikx=1,kxm%n
                    do iref=1,nref-1

                        if (vref(iref) .ne. vref(iref+1)) then
                            if (vel(iz,ikx) >= vref(iref) .and. vel(iz,ikx) < vref(iref+1)) then
                                f=(vel(iz,ikx) - vref(iref))/vel(iz,ikx)
                                g=1-f
                                ws(ikx)=g*gs(ikx,iref)+f*gs(ikx,iref+1)
                            end if
                            if (vel(iz,ikx) >= vsal) then ! ref(nref)) then
                                ws(ikx)=gs(ikx,nref)
                            end if
                        else
                            ws(ikx)=gs(ikx,iref)
                        end if                            
                    end do
                end do
           end if vxz
        end do depth
        ws=ws*taper
        GG=ws
!
deallocate(ws,gs0,gs,tfld)
end subroutine


end module
