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_rndseed1
!use OW_pickbp

implicit none

integer,private				:: nodes,ihor,ireal,igapx,igapy
!complex,private,allocatable                       :: tflds(:,:,:)
contains

!--------------------------------------------------------------------
! subroutine owfftws_init()
! integer, external		:: omp_get_thread_num,omp_get_num_threads
! allocate (tflds(kxm%n,kym%n,nodes))
! call owfftws(tflds,nodes)
! end subroutine
!!--------------------------------------------------------------------
! subroutine owfftws_close()
! deallocate (tflds)
! call owfftdestroy(nodes)
! end subroutine
!--------------------------------------------------------------------
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(taper,vsal,vel,sign,vref1)

real,optional				:: vref1(:,:,:),vsal
integer					:: iw,iz,inode,endz,initz,nsamz
integer					:: sign,ihx,ihy,ix,iy,iwb,izb
real					:: start,finish,ww
real					:: vel(kxm%n,kym%n,nz),taper(kxm%n,kym%n)
!complex, allocatable			:: inset(:,:,:,:)
!real,allocatable				:: intmp(:,:,:,:,:)
!real, allocatable				:: randz(:)
complex					:: ws(kxm%n,kym%n,nodes)
complex					:: refl(kxm%n,kym%n,nodes),tflds(kxm%n,kym%n,nodes)
!complex,allocatable			:: GG(:,:,:)
integer, external 				:: omp_get_thread_num,omp_get_num_threads
real                                             :: intmp(sog%xm%n,sog%shx%n,sog%ym%n,sog%shy%n,nz),randz(sog%xm%n*sog%ym%n*wblock%ns)
!real                                              :: intmp(nnx+1,sog%shx%n,nny+1,sog%shy%n,nz),randz((nnx+1)*(nny+1)*wblock%ns)
complex                                           :: GG(kxm%n,kym%n,wblock%ns),inset(kxm%n,kym%n,nz,wblock%ns)

ws=0.
refl=0.

 call owfftws(tflds,nodes)
!--------------------------------------------------
!
! Read rotated CIGs
!
if (sign.eq.1) then
   intmp=0.
   call sreed_window("rec_cig",6,(/sog%xm%n,sog%shx%n,sog%ym%n,sog%shy%n,   nz     ,  nhorz/),&
                                 (/sog%xm%n,sog%shx%n,sog%ym%n,sog%shy%n,   nz     ,      1/),&
                                 (/       0,       0,        0,         0,   0     , ihor-1/),&
                                 (/       1,       1,        1,         1,   1     ,      1/),&
                                 4,intmp) 
else
   intmp=0.;
   call sreed_window("rec_cig",6,(/sog%xm%n,sog%shx%n,sog%ym%n,sog%shy%n,   nz     ,  nhorz/),&
                                 (/sog%xm%n,sog%shx%n,sog%ym%n,sog%shy%n,   nz     ,      1/),&
                                 (/       0,       0,        0,         0,   0     , ihor-1/),&
                                 (/       1,       1,        1,         1,   1     ,      1/),&
                                 4,intmp) 
end if
do iwb=1, wblock%nb
   randz=0.;GG=0.;inset=0.
   call init_random_seed3(iwb,randz)
   call owcombine(iwb=iwb,intmp_=intmp,inset_=inset,nodes_=nodes,rand1_=randz,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(:,:,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( inode, sign,  ww, 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( inode, sign,  ww, ws(:,:,inode), vel(:,:,iz-1), tflds(:,:,inode), vref1(:,:,iz-1))
              else
                 call owpspi( inode, sign,  ww, ws(:,:,inode), vel(:,:,iz-1), tflds(:,:,inode))
              end if
           end if
!
           ws(:,:,inode)=ws(:,:,inode)*taper
        end do
        GG(:,:,iw)=ws(:,:,inode)
   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
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=nreal*ngapx*ngapy
irxy=(ireal-1)*ngapx*ngapy+(igapx-1)*ngapy+igapy-1
if (iwb==1) then
   call srite_window(tag,4,(/nnx,nny,sou%w%n,nrxy/),&
                           (/nnx,nny,    nnw,   1/),&
                           (/  0,  0,      0,irxy/),&
                           (/  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,4,(/nnx,nny,sou%w%n,nrxy/),&
                           (/nnx,nny,    nnw,   1/),&
                           (/  0,  0,  initw,irxy/),&
                           (/  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
