module OW_ctrillwei1

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 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,finz,izz
integer					:: sign,ihx,ihy,ix,iy,iwb,izb
real					:: start,finish,ww,a
real					:: vel(kxm%n,kym%n,nz),taper(kxm%n,kym%n)
complex					:: ws(kxm%n,kym%n,nodes)
complex					:: refl(kxm%n,kym%n,nodes),tflds(kxm%n,kym%n,nodes)
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)
complex                                           :: GG(kxm%n,kym%n,wblock%ns),inset(kxm%n,kym%n,zblock%ns,wblock%ns)
!complex                                           :: GGG((kxm%n-2*(padtracesx+tapertrx)),(kym%n-2*(padtracesy+tapertry)),wblock%ns) 

 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.;
   call init_random_seed3(iwb,randz)
   do izb=zblock%nb,1,-1
      GG=0.
      inset=0.
      initz=1
      if (izb.eq.1) initz=2
      finz=nz-zblock%blen(zblock%nb)-(zblock%nb-izb-1)*zblock%blen(izb)
      call owcombine(iwb=iwb,izb=izb,finz=finz,intmp_=intmp,inset_=inset,nodes_=nodes,rand1_=randz,igapx_=igapx,igapy_=igapy)
      if (izb.ne.zblock%nb) then

         if (sign.eq.1) then
            call wavefield_io(iwb,GG,"rec_out","reed")
         else
            call wavefield_io(iwb,GG,"sou_out","reed")
         end if

!         if (sign.eq.1) then
!            call sreed_window("rec_out",3,(/kxm%n-2*(padtracesx+tapertrx),kym%n-2*(padtracesy+tapertry),                 sou%w%n/),&
!                                          (/kxm%n-2*(padtracesx+tapertrx),kym%n-2*(padtracesy+tapertry),        wblock%blen(iwb)/),&
!                                          (/                            0,                            0,       (iwb-1)*wblock%ns/),&
!                                          (/                            1,                            1,                       1/),&
!     8,GG(padtracesx+tapertrx+1:kxm%n-(padtracesx+tapertrx),padtracesy+tapertry+1:kym%n-(padtracesy+tapertry),1:wblock%blen(iwb))) 
!         else
!            call sreed_window("sou_out",3,(/kxm%n-2*(padtracesx+tapertrx),kym%n-2*(padtracesy+tapertry),                 sou%w%n/),&
!                                          (/kxm%n-2*(padtracesx+tapertrx),kym%n-2*(padtracesy+tapertry),        wblock%blen(iwb)/),&
!                                          (/                            0,                            0,       (iwb-1)*wblock%ns/),&
!                                          (/                            1,                            1,                       1/),&
!     8,GG(padtracesx+tapertrx+1:kxm%n-(padtracesx+tapertrx),padtracesy+tapertry+1:kym%n-(padtracesy+tapertry),1:wblock%blen(iwb))) 
!         end if
      end if
!$OMP parallel do private (inode,iw,iz,izz,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.
        ws(:,:,inode)=GG(:,:,iw)
        do iz=zblock%blen(izb),initz,-1
           izz=finz+iz-zblock%blen(izb)-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(:,:,izz), 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(:,:,izz), tflds(:,:,inode), vref1(:,:,izz))
              else
                 call owpspi( inode, sign,  ww, ws(:,:,inode), vel(:,:,izz), 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

         if (sign.eq.1) then
            call wavefield_io(iwb,GG,"rec_out","rite")
         else
            call wavefield_io(iwb,GG,"sou_out","rite")
         end if
!         call srite_window("rec_out",3,(/kxm%n-2*(padtracesx+tapertrx),kym%n-2*(padtracesy+tapertry),                 sou%w%n/),&
!                                       (/kxm%n-2*(padtracesx+tapertrx),kym%n-2*(padtracesy+tapertry),        wblock%blen(iwb)/),&
!                                       (/                            0,                            0,        (iwb-1)*wblock%ns/),&
!                                       (/                            1,                            1,                       1/),&
!     8,GG(padtracesx+tapertrx+1:kxm%n-(padtracesx+tapertrx),padtracesy+tapertry+1:kym%n-(padtracesy+tapertry),1:wblock%blen(iwb))) 
!      else
!         call srite_window("sou_out",3,(/kxm%n-2*(padtracesx+tapertrx),kym%n-2*(padtracesy+tapertry),                 sou%w%n/),&
!                                       (/kxm%n-2*(padtracesx+tapertrx),kym%n-2*(padtracesy+tapertry),        wblock%blen(iwb)/),&
!                                       (/                            0,                            0,        (iwb-1)*wblock%ns/),&
!                                       (/                            1,                            1,                       1/),&
!     8,GG(padtracesx+tapertrx+1:kxm%n-(padtracesx+tapertrx),padtracesy+tapertry+1:kym%n-(padtracesy+tapertry),1:wblock%blen(iwb))) 
!      end if


!      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
end do
call owfftdestroy(nodes)
end subroutine
!========================================
subroutine wavefield_io(iwb,GG,tag,op)
complex              :: GG(:,:,:)
complex              :: temp(kxm%n,kym%n,size(GG,3))
integer              :: iwb,initw,sign
character*7          :: tag
character*4          :: op
integer              :: nnx,nny,ninx,nfnx,niny,nfny,nrxy,irxy

nnx=nx-2*(padtracesx+tapertrx)
nny=ny-2*(padtracesy+tapertry)
ninx=padtracesx+tapertrx+1
nfnx=kxm%n-(padtracesx+tapertrx)
niny=padtracesy+tapertry+1
nfny=kym%n-(padtracesy+tapertry)
nrxy=ngapx*ngapy
irxy=(igapx-1)*ngapy+igapy-1
if(op.eq."reed") then
  call sreed_window(tag,5,(/nnx,nny,          sou%w%n,  nrxy,  nreal/),&
                          (/nnx,nny, wblock%blen(iwb),     1,      1/),&
                          (/  0,  0,(iwb-1)*wblock%ns,  irxy,ireal-1/),&
                          (/  1,  1,                1,     1,      1/),&
                          8,GG(ninx:nfnx,niny:nfny,1:wblock%blen(iwb))) 
else
  call srite_window(tag,5,(/nnx,nny,          sou%w%n,  nrxy,  nreal/),&
                          (/nnx,nny, wblock%blen(iwb),     1,      1/),&
                          (/  0,  0,(iwb-1)*wblock%ns,  irxy,ireal-1/),&
                          (/  1,  1,                1,     1,      1/),&
                          8,GG(ninx:nfnx,niny:nfny,1:wblock%blen(iwb))) 
end if
end subroutine
!========================================
end module

!if(exist_file("inset")) then
!   if (iwb==1.and.a==1) then
!          call to_history("n1",nx,"inset"); call to_history("d1",sog%shx%d,"inset"); call to_history("o1",xmin,"inset")
!          call to_history("n2",ny,"inset"); call to_history("d2",sog%shy%d,"inset"); call to_history("o2",ymin,"inset")
!          call to_history("n3",sog%z%n,"inset") ; call to_history("d3",sog%z%d,"inset") ; call to_history("o3",sog%z%o,"inset")
!          call to_history("n4",1,"inset"); call to_history("d4",1.,"inset"); call to_history("o4",0.,"inset")
!          call to_history("n5",1,"inset"); call to_history("d5",1.,"inset"); call to_history("o5",0.,"inset")  
!          call to_history("esize",4,"inset")
!   call srite_window("inset",3,(/kxm%n,kym%n,              nz/),&
!                                 (/kxm%n,kym%n,zblock%blen(izb)/),&
!                                 (/    0,    0,finz-zblock%blen(izb)/),&
!                                 (/    1,    1,                      1/),&
!                                 4,real(inset(:,:,1:zblock%blen(izb),1))) 
!   end if
!end if
