module OW_combine

use OW_type
use OW_parms_areal
use OW_exec_time

implicit none

contains

!------------------------------------------------------
subroutine owcombine(intmp,inset,nodes,rand1,igap)

character*80					:: hor
real,allocatable					:: pick1(:,:)
integer						:: par1(2)
integer						:: igap
integer,optional					:: nodes
real,optional					:: rand1(:,:,:)
integer						:: iw,iz,inode
integer						:: jxs,ixs,n
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
      if (exist_file("horsplit")) then
         call from_param('horsplit',hor)
         open(51,file=hor,status='old')
         allocate(pick1(sog%xm%n,sog%sh%n));pick1=9999.
         read(51,*,end=1) par1
         do n=1,par1(2)
            read(51,*,end=1) pick1(par1(1)+n-1,:)
         end do
         1 close(51)
      end if
      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+int((sog%xm%o+sog%sh%o-xmin)/sog%xm%d)+0.5
            do iz=nzmax,nzmin,-1
               where (iz > pick1(ixs,:)) 
                  inset(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)=inset(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)+ &
                                         cmplx(intmp(iz,:,ixs),0)*exp(cmplx(0,eps*rand1(iw,ixs,1)))
               elsewhere
                  inset(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)=inset(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)+ &
                                         cmplx(intmp(iz,:,ixs),0)*exp(cmplx(0,eps*rand1(iw,ixs,2)))
               end where
            end do
         end do
      end do
      !$OMP end parallel do
   else
      !$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(iz,:,ixs),0)*exp(cmplx(0,eps*rand1(iw,ixs,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 (exist_file("horsplit")) then
      call from_param('horsplit',hor)
      open(51,file=hor,status='old')
      allocate(pick1(sog%xm%n,sog%sh%n));pick1=9999.
      read(51,*,end=2) par1
      do n=1,par1(2)
         read(51,*,end=2) pick1(par1(1)+n-1,:)
      end do
      2 close(51)
   end if
   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+int((sog%xm%o+sog%sh%o-xmin)/sog%xm%d)+0.5
            do iz=nzmax,nzmin,-1
               where (iz > pick1(jxs,:)) 
                  inset(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)=inset(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)+ &
!                                                    cmplx(intmp(iz,:,ixs),0)*exp(cmplx(0,eps*rand1(iw,1,1)))
                                                    cmplx(intmp(iz,:,ixs),0)*exp(cmplx(0,eps*rand1(iw,ixs,1)))
               elsewhere
                  inset(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)=inset(iz-nzmin+1,jxs:jxs+sog%sh%n-1,iw)+ &
!                                                    cmplx(intmp(iz,:,ixs),0)*exp(cmplx(0,eps*rand1(iw,1,2)))
                                                    cmplx(intmp(iz,:,ixs),0)*exp(cmplx(0,eps*rand1(iw,ixs,2)))
               end where
            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
            inset(:,jxs:jxs+sog%sh%n-1,iw)=inset(:,jxs:jxs+sog%sh%n-1,iw)+ & 
                                  cmplx(intmp(nzmin:nzmax,:,ixs),0)*exp(cmplx(0,eps*rand1(iw,ixs,1)))
         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 
end module
