program OW3dresmig

use sep 
use OW_type
use OW_exec_time
use OW_parms_areal
use OW_residmig_mod
use OW_fftw_areal
use OW_linterp
implicit none

real, allocatable					:: B_image(:,:,:,:),D_image(:,:,:,:,:),rhoxz(:,:)
real                                                        :: dtrho,rho1,start,finish,a,b
integer                                                     :: ix,ish,i,inode,nodes,j,k,stat,iy
complex, allocatable				:: tfld3(:,:,:,:),tfld3d(:,:,:,:,:),resd3(:,:,:,:)
integer, external                                           :: omp_get_thread_num,omp_get_num_threads

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 call sep_init()
 call sep_begin_prog()
 call owparam()

 write(0,*) "version 1.0"
 if (scan) then
    call owexec_time(start)
    write(0,*) "Initialization phase ..."
    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(resd3(img%xm%n,img%hx%n,img%z%n,img%ym%n));resd3=0.
    allocate(B_image(img%xm%n,img%hx%n,img%z%n,img%ym%n));B_image=0.
    allocate(tfld3d(img%xm%n,img%hx%n,img%z%n,img%ym%n,nodes));tfld3d=0.
    call owfftw3(resd3)
    call owfftw3d(tfld3d,nodes)
    call owexec_time(finish)
    write(0,*) "Finished initialization in", (finish-start)/60," min"
    call sreed("bkgimg",B_image,4*img%xm%n*img%hx%n*img%z%n*img%ym%n)
    resd3=cmplx(B_image,0.)
    deallocate(B_image)
    call owfft3Fexe(resd3)
    resd3=resd3/sqrt(1.*kz%n*kshx%n*kxm%n*kym%n)
    !$OMP parallel do private(inode,i,rho1) ordered schedule(dynamic,1) num_threads(nodes)
    do i=1,nrho
      inode = omp_get_thread_num()+1
      rho1 = (i-1)*drho+orho
      tfld3d(:,:,:,:,inode)=resd3
      call ow_CAMresmig(tfld3d(:,:,:,:,inode), rho1)
      call owfft3Iexed(inode)
      tfld3d(:,:,:,:,inode)=tfld3d(:,:,:,:,inode)/sqrt(1.*kz%n*kshx%n*kxm%n*kym%n)
      if(mod(i,10).eq.1) write(0,*) i,") rho = ", rho1
      !$OMP ORDERED
        call srite("Dimage",real(tfld3d(:,:,:,:,inode)),4*img%z%n*img%hx%n*img%xm%n*img%ym%n)
      !$OMP END ORDERED
    end do
    !$OMP END PARALLEL DO
    call owfftdestroy3d()
    call owfftdestroy3()
    deallocate(resd3,tfld3d)
 end if  
 call sep_end_prog()
end program
