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_shaper

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,nz),taper(kxm%n)
complex, allocatable				:: inset(:,:,:),grs(:,:,:)
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.
!--------------------------------------------------
!
! 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
!--------------------------------------------------
! 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

if (exist_file("green_s")) allocate(grs(kxm%n,sou%w%n,nz))

!$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.
           if (encode) then
              tflds(:,inode)=inset(:,iz,iw)
           else
              tflds(:,inode)=inset(:,iz,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
!
           ws(:,inode)=ws(:,inode)*taper
           if (exist_file("green_s")) grs(:,iw,iz)=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("green_s").or.exist_file("green_r")) then
   if (shaper) then
    do iz=1,nz
      call owshaper (grs(:,:,iz),1)
    end do
   end if
end if
if (exist_file("green_s").and.sign==-1) then
 call srite("green_s",grs,8*kxm%n*nz*sou%w%n)   
 deallocate(grs)
else if (exist_file("green_r").and.sign==1) then
 call srite("green_r",grs,8*kxm%n*nz*sou%w%n)   
 deallocate(grs)
end if
deallocate(inset)
end subroutine

end module
