module OW_dotst

  use OW_type
  use OW_exec_time
  use OW_parms_areal
  use OW_taper_areal
  use OW_wei_areal
  use OW_image_areal
  use OW_fftw_areal

  implicit none

  contains
!!===============================================================
  subroutine owdotst_wemva(taper,recv,source,vel,vref)

    real, optional				:: vref(:,:,:)
    integer				:: i,j,k
    real					:: taper(:)
    real					:: vel(:,:)
    complex				:: recv(:,:),source(:,:)
    real, allocatable			:: data1(:,:,:),modl1(:,:)
    real, allocatable			:: data2(:,:,:),modl2(:,:)
    real, dimension (2)		 	:: dot

    allocate(modl1(img%z%n,img%xm%n))
    allocate(modl2(img%z%n,img%xm%n))
    allocate(data1(img%z%n,img%sh%n,img%xm%n))
    allocate(data2(img%z%n,img%sh%n,img%xm%n))

    call random_number(modl1)
    call random_number(data2)
    modl1=modl1;data2=data2
    modl2=0.;data1=0.;dot=0.

    if (present(vref)) then
       call owwemva(.false.,taper,recv,source,vel,modl1,data1,vref)
       call owwemva(.true., taper,recv,source,vel,modl2,data2,vref)
    else
       call owwemva(.false.,taper,recv,source,vel,modl1,data1)
       call owwemva(.true., taper,recv,source,vel,modl2,data2)
    end if
    do i=1,img%xm%n
       do k=1,img%z%n
          do j=1,img%sh%n
             dot(1)=dot(1)+data1(k,j,i)*data2(k,j,i)
          end do
          dot(2)=dot(2)+modl1(k,i)*modl2(k,i)
       end do
    end do

    write (0,*) dot

  end subroutine
!!===============================================================
  subroutine owdotst_image()

    integer				:: i,j,k,iwa
    complex, allocatable			:: data1a(:,:),data1b(:,:)
    complex, allocatable			:: data2a(:,:),data2b(:,:)
    real, allocatable			:: modl1(:,:,:),modl2(:,:,:),temp(:,:)
    real, dimension (2)		 	:: dot

    allocate(temp(img%z%n,img%xm%n))
    allocate(data1a(img%z%n,img%xm%n))
    allocate(data1b(img%z%n,img%xm%n))
    allocate(data2a(img%z%n,img%xm%n))
    allocate(data2b(img%z%n,img%xm%n))
    allocate(modl1(img%z%n,img%sh%n,img%xm%n))
    allocate(modl2(img%z%n,img%sh%n,img%xm%n))

    call random_number(modl1)
    modl1=cos(3.14159*modl1)
    call random_number(temp)
    data2a=cmplx(temp,sin(3.14159*temp))
    call random_number(temp)
    data2b=cmplx(temp,cos(3.14159*temp))

    modl2=0.;data1a=0.;data1b=0.;dot=0.

    iwa=10

    data1a=data1a/10;data1b=data1b/10
    data2a=data2a/10;data2b=data2b/10
    modl1=modl1/10;  modl1=modl1/10

    do k=1,img%z%n
       call owimage(.false., iwa, data2a(k,:), data2b(k,:), modl2(k,:,:), .false.)
       call owimage(.true.,  iwa, data1a(k,:), data2a(k,:), modl1(k,:,:),  .true.)
       call owimage(.true.,  iwa, data2b(k,:), data1b(k,:), modl1(k,:,:), .false.)
    end do

    do i=1,img%xm%n
       do k=1,img%z%n
          do j=1,img%sh%n
             dot(2)=dot(2)+modl1(k,j,i)*modl2(k,j,i)
          end do
             dot(1)=dot(1)+(conjg(data1a(k,i))*data2b(k,i)+conjg(data2a(k,i))*data1b(k,i))/2
!             dot(1)=dot(1)+real(data1a(k,i))*real(data1b(k,i))+real(data2a(k,i))*real(data2b(k,i))
!             dot(1)=dot(1)+data2a(k,i)*data2b(k,i)
       end do
    end do

    write (0,*) dot

    deallocate(temp)
    deallocate(data1a,data1b)
    deallocate(data2a,data2b)
    deallocate(modl1,modl2)

end subroutine
!!===============================================================
  subroutine owdotst_pspi(vel,sref)

    integer				:: i,j,k
    real					:: vel(:,:)
    complex, allocatable			:: modl1(:,:,:),modl2(:,:,:)
    real, allocatable                             :: temp1(:,:),temp2(:,:,:)
    complex, allocatable			:: data1(:,:),data2(:,:)
    real, dimension (2)		 	:: dot
    real, optional				:: sref(:,:,:)
    logical				:: adj
    integer				:: iwa,iz,inode,nodes
    real, allocatable                             :: refl(:,:)
    complex, allocatable			:: ws(:,:)
    complex,allocatable			:: tflds(:,:)
    integer, external 			:: omp_get_thread_num,omp_get_num_threads

    allocate(temp1(img%xm%n,rec%w%n))
    allocate(temp2(img%z%n,img%xm%n,rec%w%n))
    allocate(modl1(img%z%n,img%xm%n,rec%w%n))
    allocate(modl2(img%z%n,img%xm%n,rec%w%n))
    allocate(data1(img%xm%n,rec%w%n))
    allocate(data2(img%xm%n,rec%w%n))

    call random_number(temp1)
    call random_number(temp2)
    data2=cmplx(temp1,sin(3.14159*temp1))
    modl1=cmplx(temp2,cos(3.14159*temp2))
    deallocate(temp1,temp2)
    modl1=modl1/10;data2=data2/10
    modl2=0.;data1=0.;dot=0.

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,*) "****************"

allocate(ws(kxm%n,nodes));ws=0
allocate(tflds(kxm%n,nodes));tflds=0
call owfftws(tflds,nodes)

   !$OMP parallel do private(inode,iwa,iz) ordered schedule(dynamic,1) num_threads(nodes)
   do iwa = 1,rec%w%n
      inode=omp_get_thread_num()+1
      ws(tapertr+padtraces+1:rec%h%n+padtraces+tapertr,inode)=data2(:,iwa)
      modl2(:,:,iwa)=0.
      modl2(1,:,iwa)=ws(:,inode)
      do iz=2,img%z%n
         if (present(sref)) then
            call owpspi(  iwa, inode, 1,  ws(:,inode), vel(iz-1,:), tflds(:,inode), "sou", sref(:,:,iz-1))
         else
            call owpspi(  iwa, inode, 1,  ws(:,inode), vel(iz-1,:), tflds(:,inode), "sou")
         end if
         modl2(iz,:,iwa)=ws(:,inode)
      end do
      if (mod(iwa,50)==1) write(0,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%xd%n, &
                                     " FROM NODE: ",inode
   end do
   !$OMP END PARALLEL DO
   ws=0;tflds=0
   allocate(refl(kxm%n,nodes));refl=0
   !$OMP parallel do private(inode,iwa,iz) ordered schedule(dynamic,1) num_threads(nodes)
   do iwa = 1,rec%w%n
      inode=omp_get_thread_num()+1
      ws(:,inode)=0
      do iz=img%z%n,2,-1
         tflds(:,inode)=modl1(iz,:,iwa)
         call owfftexes(tflds(:,inode),inode)
         refl(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)
         tflds(:,inode)=0.
         tflds(:,inode)=ws(:,inode)
         call owfftexes(tflds(:,inode),inode)
         ws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)
         ws(:,inode)=ws(:,inode)+refl(:,inode)
         if (present(sref)) then
            call owpspi(  iwa, inode,-1,  ws(:,inode), vel(iz-1,:), tflds(:,inode), "sou", sref(:,:,iz-1))
         else
            call owpspi(  iwa, inode,-1,  ws(:,inode), vel(iz-1,:), tflds(:,inode), "sou")
         end if
      end do
      if (mod(iwa,50)==1) write(0,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%xd%n, &
				" FROM NODE: ",inode
      data1(:,iwa)=ws(:,inode)
   end do
   !$OMP END PARALLEL DO
!   call owfftdestroy()
   deallocate(tflds,ws,refl)

6 format(A13,i3,A8,i3,A14,i3,A12,i3)

    do i=1,img%xm%n
       do j=1,rec%w%n
          dot(1)=dot(1)+(data1(i,j))*data2(i,j)
!          dot(1)=dot(1)+data1(i,j)*data2(i,j)
          do k=1,img%z%n
             dot(2)=dot(2)+(modl1(k,i,j))*modl2(k,i,j)
          end do
       end do
    end do

    write (0,*) dot

  end subroutine

end module
