module OW_combine

use OW_type
use OW_parms_areal
use OW_exec_time

implicit none

contains

!------------------------------------------------------
subroutine owcombine(intmp_,inset_,nodes_,rand1_,igap_,ireal_)

real,allocatable					:: pick1(:,:)
integer						:: par1(2)
integer						:: igap_,ireal_
integer,optional					:: nodes_
real,optional					:: rand1_(:,:)
integer						:: iw,iz,inode,array(sog%xm%n,nreal)
integer						:: jxs,ixs,n,irest,kxs,lxs,nareal
real						:: start,finish
complex 						:: inset_(:,:,:)
real						:: intmp_(:,:,:)
real, allocatable					:: weight(:,:,:), temp(:,:)
integer, external 					:: omp_get_thread_num,omp_get_num_threads

if (encode.and.(rand.or.exist_file("InCode"))) then
   if (encz) then
      call owexec_time(start)
      !$OMP parallel do private (   iw,ixs,jxs,inode) schedule(dynamic,1) num_threads(nodes_)
      do iw = 1,sou%w%n-subf
         do ixs=igap_,sog%xm%n,ngap
            jxs=ixs+padtraces+tapertr
            do iz=nzmax,nzmin,-1
               inset_(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)=inset_(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)+ &
                                      cmplx(intmp_(ixs,:,iz),0)*exp(cmplx(0,eps*rand1_(iw,ixs)))
            end do
         end do
      end do
      !$OMP end parallel do
   else if (encrsxz) then
      call owexec_time(start)
      call piksrt(array)
      nareal=sog%xm%n/ngap
      irest=mod(sog%xm%n,ngap)
      !$OMP parallel do private (   iw,ixs,jxs,kxs,lxs,inode) schedule(dynamic,1) num_threads(nodes_)
      do iw = 1,sou%w%n-subf
         if (igap_.le.int((sog%xm%n-irest)/ngap)) then
            do ixs=1,ngap
               kxs=array((igap_-1)*nareal+ixs,ireal_)
               jxs=(igap_-1)*nareal+ixs
               lxs=kxs+padtraces+tapertr-sog%sh%n/2
               do iz=nzmax,nzmin,-1
                  inset_(iz-nzmin+1,lxs:lxs+sog%sh%n/2-1,iw)=inset_(iz-nzmin+1,lxs:lxs+sog%sh%n/2-1,iw)+ &
                                      cmplx(intmp_(kxs,:,iz),0)*exp(cmplx(0,eps*rand1_(iw,jxs)))
               end do
            end do
         else
            do ixs=1,irest
               kxs=array((igap_-1)*nareal+ixs,ireal_)
               jxs=(igap_-1)*nareal+ixs
               lxs=kxs+padtraces+tapertr-sog%sh%n/2
               do iz=nzmax,nzmin,-1
                  inset_(iz-nzmin+1,lxs:lxs+sog%sh%n/2-1,iw)=inset_(iz-nzmin+1,lxs:lxs+sog%sh%n/2-1,iw)+ &
                                      cmplx(intmp_(kxs,:,iz),0)*exp(cmplx(0,eps*rand1_(iw,jxs)))
               end do
            end do
         end if
      end do
      !$OMP end parallel do
   else if (encw) then
      !$OMP parallel do private (iw,iz,ixs,jxs,inode) schedule(dynamic,1) num_threads(nodes_)
      do iw = 1,sou%w%n-subf
         do ixs=igap_,sog%xm%n,ngap
            jxs=ixs+int((sog%xm%o+sog%sh%o-xmin)/sog%xm%d)+0.5
            do iz=nzmax,nzmin,-1
                  inset_(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)=inset_(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)+ &
                                         cmplx(intmp_(ixs,:,iz),0)*exp(cmplx(0,eps*rand1_(iw,ixs)))
            end do
         end do
      end do
      !$OMP end parallel do
   else if (ench) then
      !$OMP parallel do private (iw,iz,ixs,jxs,inode) schedule(dynamic,1) num_threads(nodes_)
      do iw = 1,sou%w%n-subf
         do ixs=igap_,sog%xm%n,ngap
            jxs=ixs+int((sog%xm%o+sog%sh%o-xmin)/sog%xm%d)+0.5
            do iz=nzmax,nzmin,-1
                  inset_(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)=inset_(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)+ &
                                         cmplx(intmp_(ixs,:,iz),0)*exp(cmplx(0,eps*rand1_(iw,1)))
            end do
         end do
      end do
      !$OMP end parallel do
      call owexec_time(finish)
      write(0,6)  "FINISHED PHASE-ENCODING ",int(finish-start),'(s)'
6     format(A24,i5,A4)
   end if
else if (encode.and.gold.or.encode.and.sgnpi) then
   call owexec_time(start)
   if (encz) then
      !$OMP parallel do private (iw,ixs,jxs,inode) schedule(dynamic,1) num_threads(nodes_)
      do iw = 1,sou%w%n-subf
         do ixs=igap_,sog%xm%n,ngap
            jxs=ixs+padtraces+tapertr
            do iz=nzmax,nzmin,-1
               inset_(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)=inset_(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)+ &
                                      cmplx(intmp_(ixs,:,iz),0)*exp(cmplx(0,eps*rand1_(iw,ixs)))
            end do
         end do
      end do
      !$OMP end parallel do
   else
      !$OMP parallel do private (iw,ixs,jxs,inode) schedule(dynamic,1) num_threads(nodes_)
      do iw = 1,sou%w%n-subf
         do ixs=igap_,sog%xm%n,ngap
            jxs=ixs+int((sog%xm%o+sog%sh%o-xmin)/sog%xm%d)+0.5
            do iz=nzmax,nzmin,-1
               inset_(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)=inset_(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)+ & 
                                      cmplx(intmp_(ixs,:,iz),0)*exp(cmplx(0,eps*rand1_(iw,ixs)))
            end do
         end do
      end do
      !$OMP end parallel do
   end if
   call owexec_time(finish)
   write(0,6)  "FINISHED PHASE-ENCODING ",int(finish-start),'(s)'
!
!else if (.not.encode) then
!   do ixs=igap,sog%xm%n,ngap
!      jxs=ixs+int((sog%xm%o+sog%sh%o-xmin)/sog%xm%d)+0.5
!      inset(:,jxs:jxs+sog%sh%n-1,1)=inset(:,jxs:jxs+sog%sh%n-1,1)+cmplx(intmp(nzmin:nzmax,:,ixs),0)
!   end do
end if
end subroutine 
!--------------------------------------------------------
SUBROUTINE PIKSRT(array)
  real                  :: r1(sog%xm%n,nreal),a,b
  real                  :: arr(sog%xm%n,nreal+1)
  integer               :: i,j,k,array(:,:)

  call Random_Number(r1)
  arr(:,1:nreal)=r1

 do k=1,nreal
  do i=1,sog%xm%n
     arr(i,nreal+1)=real(i)
  end do
  do j=2, sog%xm%n
    a=arr(j,k)
    b=arr(j,nreal+1)
    do i=j-1,1,-1
      if (arr(i,k)<=a) goto 10
      ARR(i+1,nreal+1)=arr(i,nreal+1)
      ARR(i+1,k)=arr(i,k)
    end do
    i=0
10  ARR(i+1,nreal+1)=b
    ARR(i+1,k)=a
  end do
  array(:,k)=int(arr(:,nreal+1))
end do
END subroutine
!---------------------------------------
end module
!            jxs=ixs+int((sog%xm%o+sog%sh%o-xmin)/sog%xm%d)+0.5
