module OW_dotst

  use OW_type
  use OW_exec_time
  use OW_parms_areal
  use OW_taper_areal
  use OW_weicg_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(kxm%n,img%z%n)
    complex				:: recv(:,:),source(:,:)
    real					:: data1(img%xm%n,nsecdim,img%z%n),data2(img%xm%n,nsecdim,img%z%n)
    real					:: modl1(img%xm%n,img%z%n),modl2(img%xm%n,img%z%n)
    real, dimension (2)		 	:: dot

    call random_number(modl1)
    call random_number(data2)
    modl1=modl1/1000;data2=data2/1000
    modl2=0.;data1=0.;dot=0.
    call owwemva_init(1)
    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 k=1,img%z%n
       dot(2)=dot(2)+dot_product(modl1(:,k),modl2(:,k))
       do j=1,nsecdim
          dot(1)=dot(1)+dot_product(data1(:,j,k),data2(:,j,k))
       end do
    end do

    write (0,*) dot

  end subroutine
!!===============================================================
  subroutine owdotst_wemva1()

    integer				:: i,j,k
    real					:: data1(img%xm%n,nsecdim,img%z%n),data2(img%xm%n,nsecdim,img%z%n)
    real					:: modl1(img%xm%n,img%z%n),modl2(img%xm%n,img%z%n)
    real, dimension (2)		 	:: dot

    call random_number(modl1)
    call random_number(data2)
    modl1=modl1/1000;data2=data2/1000
    modl2=0.;data1=0.;dot=0.
    call owwemva_lop2(.false.,.false.,modl1,data1)
    call owwemva_lop2(.true., .false.,modl2,data2)
    do k=1,img%z%n
       dot(2)=dot(2)+dot_product(modl1(:,k),modl2(:,k))
       do j=1,nsecdim
          dot(1)=dot(1)+dot_product(data1(:,j,k),data2(:,j,k))
       end do
    end do

    write (0,*) dot

  end subroutine
!!===============================================================
  subroutine owdotst_dso()

    integer				:: i,j,k
    real					:: modl1(img%xm%n,nsecdim,img%z%n)
    real					:: data1(img%xm%n,nsecdim,img%z%n)
    real					:: modl2(img%xm%n,nsecdim,img%z%n)
    real					:: data2(img%xm%n,nsecdim,img%z%n)
    real					:: temp(img%xm%n,nsecdim,img%z%n)
    real, dimension (2)		 	:: dot

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

    temp=modl1
    call owdso(.false.,modl1)
    data1=modl1
    modl1=temp
    temp=data2
    call owdso(.true.,data2)
    modl2=data2
    data2=temp

    do k=1,img%z%n
       do j=1,nsecdim
          dot(1)=dot(1)+dot_product(data1(:,j,k),data2(:,j,k))
          dot(2)=dot(2)+dot_product(modl1(:,j,k),modl2(:,j,k))
       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                             :: vel(:,:)
    real, allocatable			:: modl1(:,:,:),modl2(:,:,:),tempd(:,:),tempm(:,:,:)
    real, dimension (2)		 	:: dot

    allocate(tempd(img%xm%n,img%z%n-50));tempd=0
    allocate(data1a(img%xm%n,img%z%n));data1a=0
    allocate(data1b(img%xm%n,img%z%n));data1b=0
    allocate(data2a(img%xm%n,img%z%n));data2a=0
    allocate(data2b(img%xm%n,img%z%n));data2b=0
    allocate(tempm(img%xm%n,nsecdim,img%z%n-50));tempm=0
    allocate(modl1(img%xm%n,nsecdim,img%z%n));modl1=0
    allocate(modl2(img%xm%n,nsecdim,img%z%n));modl2=0
    allocate(vel(img%xm%n,img%z%n));vel=0

    call random_number(tempm)
    modl1(:,:,26:img%z%n-25)=sin(3.14159*tempm)
    call random_number(tempd)
    data2a(:,26:img%z%n-25)=cmplx(tempd,sin(3.14159*tempd))
    call random_number(tempd)
    data2b(:,26:img%z%n-25)=cmplx(tempd,cos(3.14159*tempd))
    vel=tempd
    
    modl2=0.;data1a=0.;data1b=0.;dot=0.

    iwa=10

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

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

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

    write (0,*) dot

    deallocate(tempd,tempm)
    deallocate(data1a,data1b)
    deallocate(data2a,data2b)
    deallocate(modl1,modl2,vel)

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

    integer				:: i,j,k
    real					:: vel(kxm%n,img%z%n)
    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
    complex, 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%xm%n,img%z%n,rec%w%n))
    allocate(modl1(img%xm%n,img%z%n,rec%w%n))
    allocate(modl2(img%xm%n,img%z%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/1000;data2=data2/1000
    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
allocate(refl(kxm%n,nodes));refl=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(:,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
      ws(:,inode)=0
      do iz=img%z%n,2,-1
         ws(:,inode)=ws(:,inode)+modl1(:,iz,iwa)
         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
      data1(:,iwa)=ws(:,inode)
      data1(:,iwa)=data1(:,iwa)+modl1(:,1,iwa)
   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)+conjg(data1(i,j))*data2(i,j)
          do k=1,img%z%n
             dot(2)=dot(2)+conjg(modl1(i,k,j))*modl2(i,k,j)
          end do
       end do
    end do

    write (0,*) dot

  end subroutine

end module
