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,igapx,igapy
contains

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

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

ihor=ihor_in
ireal=ireal_in
igapx=igapx_in
igapy=igapy_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(intmp,taper,vsal,randz,vel,sign,vref1)

real,optional				:: vref1(:,:,:),vsal
integer					:: iw,iz,inode,endz,initz,nsamz
integer					:: sign,ihx,ihy,ix,iy,iwb,izb,iwi,iwf
real					:: start,finish,ww
real					:: vel(kxm%n,kym%n,nz),taper(kxm%n,kym%n)
complex,allocatable                               :: ws(:,:,:)
complex,allocatable                               :: refl(:,:,:),tflds(:,:,:)
integer, external 				:: omp_get_thread_num,omp_get_num_threads
real                                              :: intmp(:,:,:,:,:),randz(:)
complex,allocatable                               :: GG(:,:,:),inset(:,:,:,:)
complex,allocatable                               :: GGG(:,:,:,:)


allocate(GG(kxm%n,kym%n,wblock%ns),inset(kxm%n,kym%n,nz,wblock%ns))
allocate(refl(kxm%n,kym%n,nodes),tflds(kxm%n,kym%n,nodes))
allocate(ws(kxm%n,kym%n,wblock%ns))

!allocate(GGG(kxm%n,kym%n,nz,sou%w%n))

ws=0.
refl=0.

call owfftws(tflds,nodes)

do iwb=1, wblock%nb
   GG=0.;inset=0.
   iwi=(iwb-1)*sog%xm%n*sog%ym%n*wblock%ns+1
   iwf=iwi+sog%xm%n*sog%ym%n*wblock%blen(iwb)-1

   call owcombine(iwb=iwb,intmp_=intmp,inset_=inset,nodes_=nodes,rand1_=randz(iwi:iwf),igapx_=igapx,igapy_=igapy)

!$OMP parallel do private (inode,iw,iz,ww)  schedule(dynamic,1) num_threads(nodes) ordered
   do iw = 1,wblock%blen(iwb)
        inode=omp_get_thread_num()+1; 
        ww = wblock%os(iwb) + (iw-1)*wblock%ds
        ws(:,:,iw)=0;

        do iz=nz,2,-1
!
           tflds(:,:,inode)=0.;refl(:,:,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(:,:,iw)
           call owfftexes(inode)
           ws(:,:,iw)=tflds(:,:,inode)/sqrt(1.*kxm%n*kym%n)
           ws(:,:,iw)=ws(:,:,iw)+refl(:,:,inode)
           if (extrap(1:3).eq."ffd") then
              if (sou%w%all(iw).ne.0) then
                 call owffd( inode, sign,  ww, ws(:,:,iw), vel(:,:,iz-1), tflds(:,:,inode))
              else
                 ws(:,:,iw)=0.
              end if
           else if (extrap(1:4).eq."pspi") then
              if (present(vref1)) then
                 call owpspi( inode, sign,  ww, ws(:,:,iw), vel(:,:,iz-1), tflds(:,:,inode), vref1(:,:,iz-1))
              else
                 call owpspi( inode, sign,  ww, ws(:,:,iw), vel(:,:,iz-1), tflds(:,:,inode))
              end if
           end if
!
           ws(:,:,iw)=ws(:,:,iw)*taper
!           GGG(:,:,iz,wblock%indx(iwb)+iw-1)=ws(:,:,iw)
        end do
        GG(:,:,iw)=ws(:,:,iw)
   end do

!$OMP END PARALLEL DO

   if (sign.eq.1) then
      call srite("rec_out",GG(padtracesx+tapertrx+1:kxm%n-(padtracesx+tapertrx),&
                              padtracesy+tapertry+1:kym%n-(padtracesy+tapertry),1:wblock%blen(iwb)),&
                              8*(kxm%n-2*(padtracesx+tapertrx))*(kym%n-2*(padtracesy+tapertry))*wblock%blen(iwb)) 
    else
      call srite("sou_out",GG(padtracesx+tapertrx+1:kxm%n-(padtracesx+tapertrx),&
                              padtracesy+tapertry+1:kym%n-(padtracesy+tapertry),1:wblock%blen(iwb)),&
                              8*(kxm%n-2*(padtracesx+tapertrx))*(kym%n-2*(padtracesy+tapertry))*wblock%blen(iwb)) 
   end if

end do
!if (igapx==1.and.sign==1) then
! call to_history("n1",nx,"GGG");call to_history("d1",sog%shx%d,"GGG");call to_history("o1",xmin,"GGG")  
! call to_history("n2",ny,"GGG");call to_history("d2",sog%shy%d,"GGG");call to_history("o2",ymin,"GGG")  
! call to_history("n3",nz,"GGG");call to_history("d3",sog%z%d,"GGG");call to_history("o3",zmin,"GGG")  
! call to_history("n4",sou%w%n,"GGG");call to_history("d4",sou%w%d,"GGG");call to_history("o4",sou%w%o,"GGG")
! call to_history("n5",1,"GGG");	call to_history("d5",1.,"GGG");	call to_history("o5",0.,"GGG") 
!   call to_history("esize",8,"GGG")
!   call to_history("label1","cmp_x","GGG")
!   call to_history("label2","cmp_y","GGG")
!   call to_history("label3","depth","GGG")
!   call to_history("label4","frequency","GGG")
!   call srite("GGG",GGG,8*kxm%n*kym%n*nz*sou%w%n)
!end if
deallocate(GG,inset,ws)
deallocate(refl,tflds)

call owfftdestroy(nodes)
end subroutine
!!========================================
!subroutine wavefield_io(iwb,GG,tag)
!complex              :: GG(:,:,:)
!complex              :: temp(kxm%n,kym%n,size(GG,3))
!integer              :: iwb,initw
!character*7          :: tag
!integer              :: nnx,nny,nnw,nrxy,irxy
!write(0,*) " block #",iwb,"PREPARING TO WRITE ...",tag
!nnx=nx-2*(padtracesx+tapertrx)
!nny=ny-2*(padtracesy+tapertry)
!nnw=wblock%blen(iwb)
!nrxy=ngapx*ngapy
!irxy=(igapx-1)*ngapy+igapy-1
!if (iwb==1) then
!   call srite_window(tag,5,(/nnx,nny,sou%w%n,nrxy,nreal/),&
!                           (/nnx,nny,    nnw,   1,ireal/),&
!                           (/  0,  0,      0,irxy,    1/),&
!                           (/  1,  1,      1,   1,    1/),&
!                           8,temp) 
!else
!   temp=0
!   initw=(iwb-1)*wblock%blen(iwb-1)
!   call sseek(tag,0,)
!!write(0,*) "Wio 00"
!   call sreed_window(tag,5,(/nnx,nny,sou%w%n,nrxy,nreal/),&
!                           (/nnx,nny,    nnw,   1,ireal/),&
!                           (/  0,  0,  initw,irxy,    1/),&
!                           (/  1,  1,      1,   1,    1/),&
!                           8,temp) 
!!write(0,*) "Wio 01"
!!   write(0,*) size(temp,1),size(temp,2),size(temp,3)
!!   write(0,*) size(GG,1),size(GG,2),size(GG,3)
!   temp=temp+GG
!!write(0,*) "Wio 02"
!   call srite_window(tag,4,(/nnx,nny,sou%w%n,nrxy/),&
!                           (/nnx,nny,    nnw,   1/),&
!                           (/  0,  0,  initw,irxy/),&
!                           (/  1,  1,      1,   1/),&
!                           8,temp) 
!!write(0,*) "Wio 03"
!end if
!end subroutine
!========================================
end module
