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,igap)

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

ws=0.
refl=0.
   if (exist_file("xxx")) then
      call to_history("n1",kxm%n,"xxx"); call to_history("d1",sog%xm%d,"xxx"); call to_history("o1",xmin,"xxx")  
      call to_history("n2",1,"xxx") ; call to_history("d2",1.,"xxx") ; call to_history("o2",0.,"xxx")
      call to_history("n3",sog%z%n,"xxx") ; call to_history("d3",sog%z%d,"xxx") ; call to_history("o3",sog%z%o,"xxx")
      call to_history("n4",1,"xxx"); call to_history("d4",1,"xxx"); call to_history("o4",0,"xxx")
      call to_history("n5",1,"xxx"); call to_history("d5",1,"xxx"); call to_history("o5",0,"xxx")  
      call to_history("esize",4,"xxx")
      call srite("xxx",vel,4*sog%z%n*kxm%n)
   end if
!--------------------------------------------------
!
! Read rotated CIGs
!
intmp=0.
if (sign.eq.1) then
   call sreed_window("rec_cig",4,(/sog%xm%n,sog%sh%n,sog%z%n, nhorz/),&
                                 (/sog%xm%n,sog%sh%n,sog%z%n,     1/),&
                                 (/       0,       0,      0,ihor-1/),&
                                 (/       1,       1,      1,     1/),&
                                 4,intmp) 
else
   temp=0.;intmp=0.
   call sreed_window("sou_cig",4,(/sog%xm%n,sog%sh%n,sog%z%n, nhorz/),&
                                 (/sog%xm%n,sog%sh%n,sog%z%n,     1/),&
                                 (/       0,       0,      0,ihor-1/),&
                                 (/       1,       1,      1,     1/),&
                                  4,intmp) 
   do ih=1,sog%sh%n
      do ix=1,sog%xm%n
         temp(ix,ih,:)=intmp(ix,sog%sh%n-ih+1,:)
      end do
   end do
   intmp=temp
end if
! call to_history("n1",sog%xm%n,"zzz");call to_history("d1",sog%xm%d,"zzz");call to_history("o1",sog%xm%o,"zzz")
! call to_history("n2",sog%sh%n,"zzz");call to_history("d2",sog%sh%d,"zzz");call to_history("o2",sog%sh%o,"zzz")
! call to_history("n3",sog%z%n,"zzz");call to_history("d3",sog%z%d,"zzz");call to_history("o3",sog%z%o,"zzz")  
! call to_history("esize",4,"zzz");
! call to_history("n4",2,"zzz");
! call srite("zzz",intmp,4*sog%xm%n*sog%sh%n*sog%z%n)
!
!--------------------------------------------------
! Setting FFTW plans
tflds=0;
call owfftws(tflds,nodes)
if (encode) then
   allocate(inset(nx,nz,sou%w%n));inset=0
else
   allocate(inset(nx,nz,2));inset=0
end if
if (encode) then
   call owcombine(intmp_=intmp,inset_=inset,nodes_=nodes,rand1_=rand1,igap_=igap,ireal_=ireal)
else
   call owcombine(intmp_=intmp,inset_=inset,igap_=igap,ireal_=ireal)
end if
 call to_history("n1",nx,"zzz");call to_history("d1",sog%xm%d,"zzz");call to_history("o1",xmin,"zzz")
 call to_history("n2",sog%z%n,"zzz");call to_history("d2",sog%z%d,"zzz");call to_history("o2",sog%z%o,"zzz")  
 call to_history("n3",sou%w%n,"zzz");call to_history("d3",sou%w%d/(2*pi),"zzz");call to_history("o3",sou%w%o/(2*pi),"zzz")
 call to_history("esize",8,"zzz");
 call to_history("n4",2,"zzz");
 call srite("zzz",inset,8*kxm%n*sou%w%n*sog%z%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=nzmax,max(2,nzmin),-1
!
           tflds(:,inode)=0.
           if (encode) then
              tflds(:,inode)=inset(:,iz-nzmin+1,iw)
           else
              tflds(:,inode)=inset(:,iz-nzmin+1,1)
           end if
           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
!
        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()
deallocate(inset)
end subroutine
end module
