program OWresmig

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 (resmig_lin) then
 allocate(B_image(img%xm%n,nsecdim,img%z%n));B_image=0.
 call sreed("bkgimg",B_image,4*img%xm%n*nsecdim*img%z%n)
 allocate(tfld3(img%xm%n,nsecdim,img%z%n));tfld3=0.
 call owfftw3(tfld3)
    dtrho=1-rho
    tfld3=B_image
    do ix=1,kxm%n
       do ish=1,nsecdim
          tfld3(ix,ish,:)=-cmplx(0.,1.)*tfld3(ix,ish,:)*img%z%all 
       end do
    end do
    call owfft3Fexe(tfld3)
    tfld3=tfld3/sqrt(1.*kz%n*ksh%n*kxm%n)
    call ow_resmig_lin(tfld3, rho)
    call owfft3Iexe(tfld3)
    tfld3=tfld3/sqrt(1.*kz%n*ksh%n*kxm%n)
    B_image=real(tfld3)*dtrho
    call srite("Dimage",B_image,4*img%z%n*nsecdim*img%xm%n)
 else if (resmig) then
 allocate(B_image(img%xm%n,nsecdim,img%z%n));B_image=0.
 call sreed("bkgimg",B_image,4*img%xm%n*nsecdim*img%z%n)
 allocate(tfld3(img%xm%n,nsecdim,img%z%n));tfld3=0.
 call owfftw3(tfld3)
    tfld3=cmplx(B_image,0.)
    call owfft3Fexe(tfld3)
    tfld3=tfld3/sqrt(1.*kz%n*ksh%n*kxm%n)
    call ow_resmig(tfld3, rho)
    call owfft3Iexe(tfld3)
    tfld3=tfld3/sqrt(1.*kz%n*ksh%n*kxm%n)
    B_image=real(tfld3)
    call srite("Dimage",B_image,4*img%z%n*nsecdim*img%xm%n)
 deallocate(B_image)
 deallocate(tfld3)
 else if (resmigx) then
    allocate(B_image(img%xm%n,nsecdim,img%z%n));B_image=0.
    call sreed("bkgimg",B_image,4*img%xm%n*nsecdim*img%z%n)
    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(3,*) "****************"
    write(3,*) "* NUM.THREADS = ",nodes
    write(3,*) "****************"
    allocate(resd3(img%xm%n,nsecdim,img%z%n));resd3=0.
    call owfftw3(resd3)
    resd3=cmplx(B_image,0.)
    B_image=0.
    write(0,*) "Almost there ... wait please"
    call owfft3Fexe(resd3)
    resd3=resd3/sqrt(1.*kz%n*ksh%n*kxm%n)
    call owfftdestroy3()
    allocate(tfld3d(img%xm%n,nsecdim,img%z%n,nodes));tfld3d=0.
    allocate(D_image(img%xm%n,nsecdim,img%z%n,nrho));D_image=0.
    allocate(rhoxz(img%xm%n,img%z%n));rhoxz=0.
    call sreed("rhoxz",rhoxz,4*img%xm%n*img%z%n)
    a=9999.;b=-9999.
    do k=1,img%z%n
       do j=1,img%xm%n
          if (rhoxz(j,k).gt.b) b=rhoxz(j,k)
          if (rhoxz(j,k).lt.a) a=rhoxz(j,k)
       end do
    end do
    orho=a
    drho=(b-a)/(nrho-1)
    write(0,*) a,b,drho
    call owfftw3d(tfld3d,nodes)
    call owexec_time(finish)
    write(0,*) "Finished initialization in", (finish-start)/60," min"
    !$OMP parallel do private(inode,i,rho1,j,k) 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_resmig(tfld3d(:,:,:,inode), rho1)
      call owfft3Iexed(inode)
      D_image(:,:,:,i)=real(tfld3d(:,:,:,inode)/sqrt(1.*kz%n*ksh%n*kxm%n))
      if(mod(i,10).eq.1) write(0,*) i,") rho = ", rho1
    end do
    !$OMP END PARALLEL DO
    if (exist_file("xxx")) then
       call to_history("n1",img%xm%n,"xxx"); call to_history("d1",img%xm%d,"xxx"); call to_history("o1",img%xm%o,"xxx")  
       call to_history("n2",nsecdim,"xxx") ; call to_history("d2",dsecdim,"xxx") ; call to_history("o2",osecdim,"xxx")
       call to_history("n3",img%z%n,"xxx") ; call to_history("d3",img%z%d,"xxx") ; call to_history("o3",img%z%o,"xxx")
       call to_history("n4",nrho,"xxx"); call to_history("d4",drho,"xxx"); call to_history("o4",orho,"xxx")
       call to_history("n5",1,"xxx"); call to_history("d5",1,"xxx"); call to_history("o5",0,"xxx")  
       call to_history("esize",4,"xxx")
       call srite("xxx",D_image,4*img%z%n*nsecdim*img%xm%n*nrho)
    end if
    stat=owlint_init(rhoxz,nrho,orho,drho) 
    stat=owlint_lop(.true.,.false.,B_image,D_image) 
    stat=owlint_close()
    call srite("Dimage",B_image,4*img%z%n*nsecdim*img%xm%n)
    call owfftdestroy3d()
    deallocate(resd3,tfld3d,B_image,D_image,rhoxz)
 else 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(3,*) "****************"
    write(3,*) "* NUM.THREADS = ",nodes
    write(3,*) "****************"
    allocate(resd3(img%xm%n,nsecdim,img%z%n));resd3=0.
    call owfftw3(resd3)
    allocate(B_image(img%xm%n,nsecdim,img%z%n));B_image=0.
    allocate(tfld3d(img%xm%n,nsecdim,img%z%n,nodes));tfld3d=0.
    allocate(D_image(img%xm%n,nsecdim,img%z%n,nrho));D_image=0.
    call owexec_time(finish)
    call owfftw3d(tfld3d,nodes)
    write(0,*) "Finished initialization in", (finish-start)/60," min"
    do iy=1,img%ym%n
       B_image=0;D_image=0;resd3=0;
       call sreed("bkgimg",B_image,4*img%xm%n*nsecdim*img%z%n)
       resd3=cmplx(B_image,0.)
       call owfft3Fexe(resd3)
       resd3=resd3/sqrt(1.*kz%n*ksh%n*kxm%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_resmig(tfld3d(:,:,:,inode), rho1)
         call owfft3Iexed(inode)
         D_image(:,:,:,i)=real(tfld3d(:,:,:,inode)/sqrt(1.*kz%n*ksh%n*kxm%n))
         if(mod(i,10).eq.1) write(0,*) i,") rho = ", rho1, " line = ",(iy-1)*img%ym%d+img%ym%o
       end do
       !$OMP END PARALLEL DO
       call srite("Dimage",D_image,4*img%z%n*nsecdim*img%xm%n*nrho)
    end do
    call owfftdestroy3d()
    call owfftdestroy3()
    deallocate(resd3,tfld3d,B_image,D_image)
 else if (resmigtau) then
 allocate(tfld3(img%xm%n,nsecdim,img%z%n));tfld3=0.
 call owfftw3(tfld3)
    tfld3=cmplx(B_image,0.)
    call owfft3Fexe(tfld3)
    tfld3=tfld3/sqrt(1.*kz%n*ksh%n*kxm%n)
    write(0,*) "FFT"
    call ow_resmig_tau(tfld3, rho)
    write(0,*) "Resmig"
    call owfft3Iexe(tfld3)
    tfld3=tfld3/sqrt(1.*kz%n*ksh%n*kxm%n)
    write(0,*) "IFFT"
    B_image=real(tfld3)
!    B_image=abs(tfld3)
    call srite("Dimage",B_image,4*img%z%n*nsecdim*img%xm%n)
 deallocate(B_image)
 deallocate(tfld3)
 end if  
 call sep_end_prog()
end program
