 program OW_resmig

 use sep
 use OW_type
 use OW_parms_areal
 use OW_fftw_areal

 implicit none

 real				:: kz1,kz2,kznew,rho,kmh1,kmh2,kza2,kzb2
 real,allocatable			:: cig(:,:,:)
 complex,allocatable		:: c1(:,:,:),tfld(:,:,:)
 integer				:: ix, iz, ih, ikznew

 call sep_init()
 call sep_begin_prog()

 call owparam()

 open(unit=2,file="resmig.log")

 allocate(c1(sog%z%n,sog%sh%n,sog%xm%n))
 allocate(cig(sog%z%n,sog%sh%n,sog%xm%n),tfld(sog%z%n,sog%sh%n,sog%xm%n))

 call sreed("in",cig,4*sog%z%n*sog%sh%n*sog%xm%n)

 call owfftw3(tfld)
 tfld=cmplx(cig,0.)
 call sfftw_execute(plan3F)
 c1=tfld*sqrt(1./(1.*sog%z%n*sog%sh%n*sog%xm%n));
 tfld=0;

 rho=1.1

 do ix=1,kxm%n
    do ih=1,ksh%n

       kmh1=(ksh%all(ih)+kxm%all(ix))**2
       kmh2=(-ksh%all(ih)+kxm%all(ix))**2

       do iz=1,kz%n
          if (kz%all(iz) .eq. 0.) then
!             if (ksh%all(ih) .eq. 0. .and. kxm%all(ih) .ne. 0.) then
!                kznew=kxm%all(ix)**2*ksh%all(2)/kz%all(2)
!             else if (kxm%all(ih) .eq. 0. .and. ksh%all(ih) .ne. 0.) then
!                kznew=ksh%all(ih)**2*kxm%all(2)/kz%all(2)
!             else if (kxm%all(ih) .eq. 0. .and. ksh%all(ih) .eq. 0.) then
!                kznew=0.
!             end if

          else

             kza2=kz%all(iz)**2+ksh%all(ih)**2
             kzb2=kz%all(iz)**2+kxm%all(ix)**2
             kz1=rho**2*kza2*kzb2/kz%all(iz)**2 - kmh1
             kz2=rho**2*kza2*kzb2/kz%all(iz)**2 - kmh2
             if (kz1 .ge. 0. .and. kz2 .ge. 0.) then
                kznew = 0.5*(sqrt(kz1)+sqrt(kz2))
             end if

             if (kznew .ge. 0. .and. kznew .le. maxval(kz%all)) then
                ikznew=int(kznew/kz%all(2)) + 1.5
                tfld(ikznew,ih,ix) = tfld(ikznew,ih,ix) + c1(iz,ih,ix) 
! write(2,*) kz%all(iz),ksh%all(ih),kxm%all(ix),kznew
             else if (kznew .ge. minval(kz%all) .and. kznew .le. kz%all(kz%n)) then
                ikznew=int((minval(kz%all)-kznew)/kz%all(2))+int(kz%n/2.)+1.5
                tfld(ikznew,ih,ix) = tfld(ikznew,ih,ix) + c1(iz,ih,ix) 
! write(2,*) kz%all(iz),ksh%all(ih),kxm%all(ix),kznew
             end if
               
          end if
       end do
    end do
 end do

 call sfftw_execute(plan3I)
 cig=real(tfld*sqrt(1./(1.*sog%z%n*sog%sh%n*sog%xm%n)))

 call srite("out",cig,4*sog%z%n*sog%sh%n*sog%xm%n)
 deallocate(cig,tfld,c1)

 call sep_end_prog()

 end program
