module OW_combine

use OW_type
use OW_parms_areal
use OW_exec_time

implicit none

contains

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

integer,optional					:: igap,nodes
real,optional					:: rand1(:,:,:)
integer						:: iw,iz,inode
integer						:: jxs,ixs
real						:: start,finish
complex 						:: inset(:,:,:)
real						:: intmp(:,:,:)
real, allocatable					:: weight(:,:,:), temp(:,:)
integer, external 					:: omp_get_thread_num,omp_get_num_threads

if (exist_file("weight")) then
   allocate(weight(sog%xm%n,sog%z%n,sou%w%n-subf),temp(sog%xm%n,sog%z%n))
   call sreed("weight",temp,4*sog%xm%n*sog%z%n)
   do iw=1,sou%w%n-subf
      weight(:,:,iw)=temp      
   end do
   deallocate(temp)
end if

if (encode.and.(rand.or.exist_file("InCode"))) then
   call owexec_time(start)
   !$OMP parallel do private (   iw,&
                                 iz,&
                                ixs,&
                                jxs,&
                               inode) schedule(dynamic,1) num_threads(nodes)
   do iw = 1,sou%w%n-subf
!      inode=omp_get_thread_num+1
      do ixs=1,sog%xm%n
         jxs=ixs+int((sog%xm%o+sog%sh%o-xmin)/sog%xm%d)+0.5
         do iz=nzmax,nzmin,-1
!            if (allocated(weight)) then
!               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*(weight(ixs,iz,iw))))!!*rand1(iw,ixs,iz)
!            else
!               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,iz)))
!            end if
            if ((iz-1)*sog%z%d > 0.8) then
               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)))
            else
               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 if
         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)
   if (allocated(weight)) deallocate(weight)
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+int((sog%xm%o+sog%sh%o-xmin)/sog%xm%d)+0.5
            do iz=nzmax,nzmin,-1
               if ((iz-1)*sog%z%d > 0.8) then
                  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)))
               else
                  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 if
            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
