module OW_ctrillwei

use sep
use OWvelo_areal
use OW_type
use OW_fftw_areal
use OW_parms_areal
use OW_exec_time
use OW_pspi_areal
use OW_ffd_areal

implicit none

integer,private				:: nodes,ihor,ireal

contains

!--------------------------------------------------------------------
subroutine owctrlwei_init(ihor_in,ireal_in)

integer                                                     :: ihor_in,ireal_in
integer, external 					:: omp_get_thread_num,omp_get_num_threads

ihor=ihor_in
ireal=ireal_in
if (node.eq.9999) then
  !$OMP PARALLEL
  nodes = omp_get_num_threads()
  !$OMP END PARALLEL
else
  nodes=node
end if
write(0,*) "****************"
write(0,*) "* NUM.THREADS = ",nodes
write(0,*) "****************"

end subroutine

!---------------------------------------------------------
subroutine owctrlwei(taper,vsal,vel,GG,sign,vref1)

integer						:: igap
real,optional					:: vref1(:,:,:),vsal
integer						:: iw,iz,inode
integer						:: sign,ih,ix
real						:: start,finish
real						:: vel(kxm%n,nz),taper(kxm%n)
real						:: intmp(sog%xm%n,2,sog%z%n)
complex						:: ws(kxm%n,nodes),tflds(kxm%n,nodes),refl(kxm%n,nodes)
complex						:: GG(:,:)
complex, allocatable                                        :: wavefld(:,:,:)
integer, external 					:: omp_get_thread_num,omp_get_num_threads

ws=0.
refl=0.
!--------------------------------------------------
!
! Read rotated CIGs
!
intmp=0.
call sreed_window("cig",3,(/sog%xm%n,sog%sh%n,sog%z%n/),&
                              (/sog%xm%n,sog%sh%n,sog%z%n/),&
                              (/       0,      0, 0/),&
                              (/       1,      1, 1/),&
                              4,intmp(:,1,:)) 
!--------------------------------------------------
! Setting FFTW plans
tflds=0;
call owfftws(tflds,nodes)
if (exist_file("wvfld")) allocate(wavefld(kxm%n,nz,sou%w%n))
!$OMP parallel do private (inode,iw,iz) ordered schedule(dynamic,1) num_threads(nodes)
   do iw = 1,sou%w%n-subf
        inode=omp_get_thread_num()+1; 
        ws(:,inode)=0;
        do iz=nz,2,-1
!
           tflds(:,inode)=0.
           tflds(padtraces+tapertr+1:sog%xm%n+padtraces+tapertr,inode)=cmplx(intmp(:,1,iz),0.)
           call owfftexes(tflds(:,inode),inode)
           refl(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)
           tflds(:,inode)=0.
           tflds(:,inode)=ws(:,inode)
           call owfftexes(tflds(:,inode),inode)
           ws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)
           ws(:,inode)=ws(:,inode)+refl(:,inode)
!
           if (extrap(1:3).eq."ffd") then
              if (sou%w%all(iw).ne.0) then
                 call owffd( iw, inode, sign,  ws(:,inode), vel(:,iz-1), tflds(:,inode))
              else
                 ws(:,inode)=0.
              end if
           else if (extrap(1:4).eq."pspi") then
              if (present(vref1)) then
                 call owpspi( iw, inode, sign,  ws(:,inode), vel(:,iz-1), tflds(:,inode), vref1(:,:,iz-1))
              else
                 call owpspi( iw, inode, sign,  ws(:,inode), vel(:,iz-1), tflds(:,inode))
              end if
           end if
!
           ws(:,inode)=ws(:,inode)*taper

           if (exist_file("wvfld")) wavefld(:,iz,iw)=ws(:,inode)
        end do
        if (mod(iw,50)==1) write(0,5)  "DONE FREQUENCY # ",iw," OUT OF ",sou%w%n
5       format(A17,i5,A8,i5)
        GG(:,iw)=ws(:,inode)
   end do

!$OMP END PARALLEL DO
call owfftdestroy()
if (exist_file("wvfld")) then
 call to_history("n1",nx,"wvfld");call to_history("d1",sog%xm%d,"wvfld");call to_history("o1",xmin,"wvfld")  
 call to_history("n2",sog%z%n,"wvfld");call to_history("d2",sog%z%d,"wvfld");call to_history("o2",sog%z%o,"wvfld")
 call to_history("n3",sou%w%n,"wvfld");call to_history("d3",sou%w%d/(2*pi),"wvfld");call to_history("o3",sou%w%o/(2*pi),"wvfld")
 call to_history("n4",1,"wvfld");call to_history("d4",1.,"wvfld");call to_history("o4",1.,"wvfld")
 call to_history("esize",8,"wvfld")
 call srite("wvfld",wavefld,8*kxm%n*nz*sou%w%n)   
 deallocate(wavefld)
end if
end subroutine
end module
