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_combine
use OW_pspi_areal
use OW_ffd_areal
!use OW_pickbp

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,rand1,vref1,igapx,igapy)

integer						:: igapx,igapy
real,optional					:: vref1(:,:,:),rand1(sog%xm%n,sog%ym%n,sou%w%n),vsal
integer						:: iw,iz,inode
integer						:: sign,ihx,ihy,ix,iy
real						:: start,finish
real						:: vel(kxm%n,kym%n,nz),taper(kxm%n,kym%n)
complex, allocatable				:: inset(:,:,:,:)
real						:: intmp(sog%xm%n,sog%shx%n,sog%ym%n,sog%shy%n,sog%z%n)
real, allocatable					:: temp(:,:,:,:,:)
complex						:: ws(kxm%n,kym%n,nodes),tflds(kxm%n,kym%n,nodes)
complex						:: refl(kxm%n,kym%n,nodes)
complex						:: GG(:,:,:)
integer, external 					:: omp_get_thread_num,omp_get_num_threads

ws=0.
refl=0.
!--------------------------------------------------
!
! Read rotated CIGs
!
intmp=0.
if (sign.eq.1) then
   call sreed_window("rec_cig",6,(/sog%xm%n,sog%shx%n,sog%ym%n,sog%shy%n,sog%z%n, nhorz/),&
                                 (/sog%xm%n,sog%shx%n,sog%ym%n,sog%shy%n,sog%z%n,     1/),&
                                 (/       0,       0,       0,        0,       0,ihor-1/),&
                                 (/       1,       1,       1,        1,       1,     1/),&
                                 4,intmp) 
else
   allocate(temp(sog%xm%n,sog%shx%n,sog%ym%n,sog%shy%n,sog%z%n))
   temp=0.;intmp=0.
   call sreed_window("sou_cig",6,(/sog%xm%n,sog%shx%n,sog%ym%n,sog%shy%n,sog%z%n, nhorz/),&
                                 (/sog%xm%n,sog%shx%n,sog%ym%n,sog%shy%n,sog%z%n,     1/),&
                                 (/       0,       0,       0,        0,       0,ihor-1/),&
                                 (/       1,       1,       1,        1,       1,     1/),&
                                 4,intmp) 
   do ihy=1,sog%shy%n
      do iy=1,sog%ym%n
         do ihx=1,sog%shx%n
            do ix=1,sog%xm%n
               temp(ix,ihx,iy,ihy,:)=intmp(ix,sog%shx%n-ihx+1,iy,sog%shy%n-ihy+1,:)
            end do
         end do
      end do
   end do
   intmp=temp
   deallocate(temp)
end if
!--------------------------------------------------
! Setting FFTW plans
tflds=0;
call owfftws(tflds,nodes)
allocate(inset(kxm%n,kym%n,nz,sou%w%n));inset=0
call owcombine(intmp_=intmp,inset_=inset,nodes_=nodes,rand1_=rand1,igapx_=igapx,igapy_=igapy,ireal_=ireal)
   if (exist_file("zzz")) then
       call to_history("n1",nx,"zzz");  call to_history("d1",sog%xm%d,"zzz");    call to_history("o1",xmin,"zzz")
       call to_history("n2",ny,"zzz"); call to_history("d2",sog%ym%d,"zzz");   call to_history("o2",ymin,"zzz")
       call to_history("n3",nz,"zzz");  call to_history("d3",sog%z%d,"zzz");    call to_history("o3",sog%z%o,"zzz")
       call to_history("n4",sou%w%n,"zzz"); call to_history("d4",sou%w%d/(2*pi),"zzz");   call to_history("o4",sou%w%o/(2*pi),"zzz")
       call to_history("n5",1,"zzz");   call to_history("d5",1.,"zzz");     call to_history("o5",0.,"zzz")  
       call to_history("esize",8,"zzz");
       call srite("zzz", inset,8*nx*ny*sou%w%n*nz)
   end if

!$OMP parallel do private (inode,iw,iz)  schedule(dynamic,1) num_threads(nodes) !ordered
   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(:,:,inode)=inset(:,:,iz,iw)
           call owfftexes(inode)
           refl(:,:,inode)=tflds(:,:,inode)/sqrt(1.*kxm%n*kym%n)
           tflds(:,:,inode)=0.
           tflds(:,:,inode)=ws(:,:,inode)
           call owfftexes(inode)
           ws(:,:,inode)=tflds(:,:,inode)/sqrt(1.*kxm%n*kym%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
        end do
        if (mod(iw,50)==1) write(0,5)  "DONE FREQUENCY # ",iw," OUT OF ",sou%w%n
5       format(A17,i5,A8,i5)
        
!        !$OMP ORDERED
        GG(:,:,iw)=ws(:,:,inode)
!        !$OMP END ORDERED
   end do

!$OMP END PARALLEL DO
call owfftdestroy(nodes)
deallocate(inset)
end subroutine
end module
