module OW_dotst

  use OW_type
  use OW_parms_areal
  use OW_wei_areal

  implicit none

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

    real, optional				:: vref(:,:,:)
    integer				:: i,j,k,l,m,ind
    real					:: taper(:,:)
    real					:: vel(kxm%n,kym%n,img%z%n)
    complex				:: recv(:,:,:),source(:,:,:)
    real					:: data1(img%xm%n*nsecdimx*img%ym%n*nsecdimy*img%z%n)
    real					:: data2(img%xm%n*nsecdimx*img%ym%n*nsecdimy*img%z%n)
    real					:: modl1(img%xm%n*img%ym%n*img%z%n)
    real					:: modl2(img%xm%n*img%ym%n*img%z%n)
    double precision, dimension (2)	 	:: dot

    call random_number(modl1)
    call random_number(data2);!data2=2*data2-1
    modl1=modl1/1000000;data2=data2/1000000
    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 i=1,img%z%n
    do j=1,nsecdimy
    do k=1,img%ym%n
    do l=1,nsecdimx
    do m=1,img%xm%n
       ind=(i-1)*img%xm%n*nsecdimx*img%ym%n*nsecdimy+(j-1)*img%xm%n*nsecdimx*img%ym%n+(k-1)*img%xm%n*nsecdimx+(l-1)*img%xm%n+m
       dot(1)=dot(1)+dble(data1(ind))*dble(data2(ind))
    end do
    end do
    end do
    end do
    end do

    do i=1,img%z%n
    do k=1,img%ym%n
    do m=1,img%xm%n
       ind=(i-1)*img%xm%n*img%ym%n+(k-1)*img%xm%n+m
       dot(2)=dot(2)+dble(modl1(ind))*dble(modl2(ind))
    end do
    end do
    end do

    write (0,*) dot

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

    integer				:: i,j,k,l,m,iz,iw
    complex				:: data1a(kxm%n,kym%n)
    complex				:: data1b(kxm%n,kym%n)
    complex				:: data2a(kxm%n,kym%n)
    complex				:: data2b(kxm%n,kym%n)
    real					::   temp(kxm%n,kym%n)
    real					:: modl1(img%xm%n,nsecdimx,img%ym%n,nsecdimy)
    real					:: modl2(img%xm%n,nsecdimx,img%ym%n,nsecdimy)
    double precision, dimension (2)	 	:: dot

     dot=0.
     do iw=1,sou%w%n
     do iz=1,img%z%n
       data1a=0;data1b=0
       call random_number(modl1)
       call random_number(temp)
       modl1=modl1/1000000
       data2a=cmplx(sin(temp),cos(temp))/1000000
       data2b=cmplx(exp(temp),tan(temp))/1000000
       modl2=0.
       call owimage(.false., iw, data2a, data2b, modl2,.false.)
       call owimage(.true. , iw, data1a, data2b, modl1, .true.)
       call owimage(.true. , iw, data2a, data1b, modl1,.false.)

       do k=1,kym%n
       do m=1,kxm%n
          dot(1)=dot(1)+(conjg(data1a(m,k))*data2a(m,k)+conjg(data1b(m,k))*data2b(m,k))/2
       end do
       end do

       do j=1,nsecdimy
       do k=1,img%ym%n
       do l=1,nsecdimx
       do m=1,img%xm%n
          dot(2)=dot(2)+dble(modl1(m,l,k,j))*dble(modl2(m,l,k,j))
       end do
       end do
       end do
       end do
    end do
    end do

    write (0,*) dot

  end subroutine
!!===============================================================
  subroutine owdotst_extrap()

    integer				:: i,j,k,l,m,iz,iw
    complex				:: data1(kxm%n,kym%n,sou%w%n)
    complex				:: data2(kxm%n,kym%n,sou%w%n)
    complex				:: modl1(kxm%n,kym%n,sou%w%n)
    complex				:: modl2(kxm%n,kym%n,sou%w%n)
    real					::  temp(kxm%n,kym%n,sou%w%n)
    double precision, dimension (2)	 	:: dot
    integer, external			:: omp_get_thread_num,omp_get_num_threads

     dot=0.
     call random_number(temp)
     modl1=cmplx(exp(temp),tan(temp))/1000000
     data2=cmplx(sin(temp),cos(temp))/1000000
     modl2=0.;data1=0;
     !$OMP parallel do private(inode,iw,iz) ordered schedule(dynamic,1) num_threads(nodes)
     do iw=1,sou%w%n
        inode=omp_get_thread_num()+1; 
        do iz=1,img%z%n
          if (extrap(1:3).eq."ffd") then
             if (rec%w%all(iwa).ne.0) then
                call owffd( iwa, inode, 1,  data2(:,:,iw), vel(:,:,iz-1), tflds(:,:,inode), "sou")
             end if
          else if (extrap(1:4).eq."pspi") then
             if (present(sref)) then
                call owpspi(  iwa, inode, 1, dws(1,:,:,inode), vel(:,:,iz-1), tflds(:,:,inode), "sou", sref(:,:,iz-1))
             else
                call owpspi(  iwa, inode, 1, dws(1,:,:,inode), vel(:,:,iz-1), tflds(:,:,inode), "sou")
             end if
          end if
          dws(1,:,:,inode)=dws(1,:,:,inode)*taper




       do k=1,kym%n
       do m=1,kxm%n
          dot(1)=dot(1)+(conjg(data1a(m,k))*data2a(m,k)+conjg(data1b(m,k))*data2b(m,k))/2
       end do
       end do

       do j=1,nsecdimy
       do k=1,img%ym%n
       do l=1,nsecdimx
       do m=1,img%xm%n
          dot(2)=dot(2)+dble(modl1(m,l,k,j))*dble(modl2(m,l,k,j))
       end do
       end do
       end do
       end do
    end do
    end do

    write (0,*) dot

  end subroutine
!!===============================================================
end module
