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

implicit none

integer,private				:: nodes,ihor,ireal,igapx,igapy,iset

contains

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

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

ihor=ihor_in
ireal=ireal_in
igapx=igapx_in
igapy=igapy_in
iset=iset_in
!$OMP PARALLEL
nodes = omp_get_num_threads()
!$OMP END PARALLEL
write(0,*) "****************"
write(0,*) "* NUM.THREADS = ",nodes
write(0,*) "****************"


end subroutine

!---------------------------------------------------------
subroutine owctrlwei(taper,vsal,vel,sign,randz,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(ncigx,sog%shx%n,ncigy,sog%shy%n,nz),randz(:)
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",7,(/ncigx,sog%shx%n,   ncigy,sog%shy%n,  nzimg,  nhorz,   nset/),&
                                 (/ncigx,sog%shx%n,   ncigy,sog%shy%n,     nz,      1,      1/),&
                                 (/    0,        0,       0,        0,nzmin-1, ihor-1, iset-1/),&
                                 (/    1,        1,       1,        1,      1,      1,      1/),&
                                 4,intmp) 
else
   intmp=0.;
   call sreed_window("rec_cig",7,(/ncigx,sog%shx%n,   ncigy,sog%shy%n,  nzimg,  nhorz,   nset/),&
                                 (/ncigx,sog%shx%n,   ncigy,sog%shy%n,     nz,      1,      1/),&
                                 (/    0,        0,       0,        0,nzmin-1, ihor-1, iset-1/),&
                                 (/    1,        1,       1,        1,      1,      1,      1/),&
                                 4,intmp) 
end if
do iwb=1, wblock%nb
   do izb=zblock%nb,1,-1
write(0,*) iwb," from ",wblock%nb, " and ", izb, " from ", zblock%nb
      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)
!write(0,*)      finz,nz,zblock%blen(zblock%nb),zblock%nb,izb,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
      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 (ww.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
         call wavefield_io(iwb,GG,"rec_out","rite")
      else
         call wavefield_io(iwb,GG,"sou_out","rite")
      end if
   end do
end do
call owfftdestroy(nodes)

end subroutine
!========================================
subroutine wavefield_io(iwb,GG,tag,op)
complex              :: GG(:,:,:)
integer              :: iwb
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,7,(/nnx,nny,          sou%w%n,  nset,  nreal,1,1/),&
                          (/nnx,nny, wblock%blen(iwb),     1,      1,1,1/),&
                          (/  0,  0,(iwb-1)*wblock%ns,iset-1,ireal-1,0,0/),&
                          (/  1,  1,                1,     1,      1,1,1/),&
                          8,GG(ninx:nfnx,niny:nfny,1:wblock%blen(iwb))) 
else
  call srite_window(tag,7,(/nnx,nny,          sou%w%n,  nset,  nreal,1,1/),&
                          (/nnx,nny, wblock%blen(iwb),     1,      1,1,1/),&
                          (/  0,  0,(iwb-1)*wblock%ns,iset-1,ireal-1,0,0/),&
                          (/  1,  1,                1,     1,      1,1,1/),&
                          8,GG(ninx:nfnx,niny:nfny,1:wblock%blen(iwb))) 
end if
end subroutine
!========================================
subroutine gathers_io(GG,tag,op)
complex              :: GG(:,:,:)
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,7,(/nnx,nny, sou%w%n,  nset,  nreal,1,1/),&
                          (/nnx,nny, sou%w%n,     1,      1,1,1/),&
                          (/  0,  0,       0,iset-1,ireal-1,0,0/),&
                          (/  1,  1,       1,     1,      1,1,1/),&
                          8,GG) 
else
  call srite_window(tag,7,(/nnx,nny, sou%w%n,  nset,  nreal,1,1/),&
                          (/nnx,nny, sou%w%n,     1,      1,1,1/),&
                          (/  0,  0,       0,iset-1,ireal-1,0,0/),&
                          (/  1,  1,       1,     1,      1,1,1/),&
                          8,GG) 
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
