module RWE3D_wemig
  use sep            !! basic SEPlib routine I/O module
  use util_          !! handle axes module
  use fft            !! FFT module - could use others though
  use Kernel3D       !! GOGI's 80 degree code
  use sep_timers_mod !! Timing
  use fft_mod        !! Set up FFTW
  use elliptic       !! Elliptic integrals for Ellipsoidal 
  use RWE3D_geometry !! For calculcating geometry factors
  implicit none
  integer, private :: ne1,ne2,ne3,nw,nxtap,nytap,nx1,nx2,nx3
  integer, private :: ns1,ns2,nr1,nr2,method,nyint,nxint
  real, private :: de1,de2,de3,dw,oe1,oe2,oe3,ow,dr1,dr2,or1,or2
  real, private ::  pi,r0,eps,ox1,ox2,ox3,dx1,dx2,dx3,os1,os2
  real, private :: ds1,ds2,minang,maxang
  real, allocatable,private :: ww(:),tap(:,:),Timg(:,:,:,:)
  real, allocatable,private :: vmin(:),saxmax(:)
  complex, allocatable, private :: rax(:,:,:),sax(:,:,:) 
  logical, private :: verbose,source_norm,timing
  !!
contains
  !----------------------------------------------------------------
  !! . . Initialization routine - parameter passing and array allocation
  subroutine wemig_axes_init(aE1,aE2,aE3,aW,aX1,aX2,aX3,aS1,aS2,aR1,aR2)
    type(myaxis)          :: aE1,aE2,aE3,aW,aX1,aX2,aX3,aS1,aS2,aR1,aR2

    !! . . Cartesian
    nx1=ax1%n; nx2=ax2%n; nx3=ax3%n
    ox1=ax1%o; ox2=ax2%o; ox3=ax3%o
    dx1=ax1%d; dx2=ax2%d; dx3=ax3%d

    !! . . RWE
    ne1=aE1%n; ne2=aE2%n; ne3=aE3%n;
    de1=aE1%d; de2=aE2%d; de3=aE3%d; 
    oe1=aE1%o; oe2=aE2%o; oe3=aE3%o; 

    !! . . Input data
    nr1=aR1%n; nr2=aR2%n; 
    dr1=aR1%d; dr2=aR2%d; 
    or1=aR1%o; or2=aR2%o

    !! . . Shot axis
    ns1=as1%n; ns2=as2%n;
    os1=as1%o; os2=as2%o;
    ds1=as1%d; ds2=as2%d

    !! . . Frequency
    pi=acos(-1.)
    nw=aW%n; dw=aW%d*2.*pi; ow=aW%o*2.*pi

  end subroutine wemig_axes_init
 
 !----------------------------------------------------------------
  !! . . Initialization routine - parameter passing and array allocation
  subroutine wemig_param_init(verbose_in,r0_in,minang_in,maxang_in,eps_in,method_in)
    integer :: method_in
    real    :: r0_in,rbeg_in,rend_in,eps_in,minang_in,maxang_in
    logical :: verbose_in

    !! For OMP
    integer, external :: omp_get_num_threads,omp_get_thread_num
    integer  :: nth,ith,i1,ie3,ie1,ie2

    pi=acos(-1.)
    eps=eps_in    
    method = method_in
    verbose=verbose_in;    
    r0=r0_in
    minang = minang_in
    maxang = maxang_in
    call from_param("nxtap" ,nxtap, 40) !! Width of cosine taper
    call from_param("nytap" ,nytap, 40) !! Width of cosine taper
    call from_param("source_norm",source_norm,.false.)
    call from_param("timing",timing,.false.)

    !$OMP PARALLEL
    nth = omp_get_num_threads()
    !$OMP END PARALLEL

    !! Velocity and fields
    allocate( vmin(ne3),saxmax(nth) )
    allocate( tap(ne1,ne2),sax(ne1,ne2,nth),rax(ne1,ne2,nth),ww(nth))
    allocate( Timg(ne1,ne2,ne3,nth) )

    !! Initalize
    !$OMP PARALLEL DO PRIVATE(ie3,ie2,ie1) 
    do ie3=1,ne3
       do ie2=1,ne2
          do ie1=1,ne1
             Timg(ie1,ie2,ie3,1)=0.
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    if (nth .gt. 1) then
       !$OMP PARALLEL DO PRIVATE(ith,ie3,ie2,ie1) 
       do ith=2,nth
          do ie3=1,ne3
             do ie2=1,ne2
                do ie1=1,ne1
                   Timg(ie1,ie2,ie3,ith)=0.
                end do
             end do
          end do
       end do
       !$OMP END PARALLEL DO
    end if
    call rwetaper_init() !! Initialize the taper program

  end subroutine wemig_param_init
  
  !----------------------------------------------------------------
  !! . . Main wavefield extrapolation and imaging routine
  subroutine wemig(rwf,swf,RImg,ishot,RVel,fields)
    integer :: iw,ie1,ie2,ie3,inn,is,id,ip,ishot
    real    :: vtest,rarg,sarg,ss,w
    real    :: RImg(:,:,:),fields(:,:,:,:),RVel(:,:,:)
    complex :: rwf(:,:,:),swf(:,:,:)


    !! For OMP
    integer, external :: omp_get_num_threads,omp_get_thread_num
    integer  :: nth,ith
    logical :: dofilt

    !$OMP PARALLEL
    nth = omp_get_num_threads()
    !$OMP END PARALLEL


    sarg =-1.; 
    rarg = 1.

    !! Get Minimum values for the Fourier filtering
    do ie3=1,ne3
       !       vmin(ie3) = minval( RVel(:,:,ie3) ) 
       vmin(ie3) = minval( RVel(:,:,ie3) / fields(:,:,ie3,1) ) 
    end do

    if (method .eq. 1) then
       call shot_ker_init(nth, ne1, ne2,    de1, de2,   de3, 0., 0., .11)
    else if  (method .eq. 2) then
       call shot_ker_init(nth, ne1, ne2, de1, r0*de2,r0*de3, 0., 0., .11)
    else if (method .eq. 3)  then
       call shot_ker_init(nth, ne1, ne2, de1,    de2,   de3, 0., 0., .11)
    end if

    !! . . Initialize FFT2D algorithm
    call fftw_init(sax,rax)

    !$OMP PARALLEL DO SCHEDULE(DYNAMIC) PRIVATE(iw,id,ie3,ie2,ie1) &
    SHARED(nw,ow,dw,ne3,ne2,ne1,rarg,sarg,Timg,sax,rax,fields,saxmax,tap)

    Frequency_loop: do iw=1,nw  !! FREQUENCY LOOP         
       id = omp_get_thread_num()+1

       ww(id) = dble(iw-1.)*dw+ow

       !! . . Define the Source wavefield 
       do ie2=1,ne2
          do ie1=1,ne1
             sax(ie1,ie2,id) = swf(ie1,ie2,iw)
          end do
       end do

       !! . . Define Receiver wavefield 
       do ie2=1,ne2
          do ie1=1,ne1
             rax(ie1,ie2,id) = rwf(ie1,ie2,iw)
          end do
       end do

       !! . . Taper wavefields
       call rwetaper(sax,id) 
       call rwetaper(rax,id)

       !! . . Find the max amplitude of the SWF for normalization
       if (source_norm) saxmax(id) = maxval( real(conjg(sax(:,:,id))*sax(:,:,id)))

       Extrap_loop: do ie3=2,ne3  

          !! . . Step Source and Receiver wavefields
        !  call shot_ker_onestep(id,sax,(Rvel(:,:,ie3)+Rvel(:,:,ie3-1))/2.,&
        !  (fields(:,:,ie3,:)+fields(:,:,ie3-1,:))/2.,ww(id),sarg)

          call shot_ker_onestep(id,rax,(Rvel(:,:,ie3)+Rvel(:,:,ie3-1))/2.,&
          (fields(:,:,ie3,:)+fields(:,:,ie3-1,:))/2.,ww(id),rarg)

          !! . . High-angle filter  FFT to kx
          !call SAX_Forward(sax,id)
          call RAX_Forward(rax,id)

          !! . . Phase correction routine
          !call wemig_phs_correction(id,sax,ww(id),vmin(ie3),sarg,&
          !(fields(:,:,ie3,:)+fields(:,:,ie3-1,:))/2.)
          call wemig_phs_correction(id,rax,ww(id),vmin(ie3),rarg,&
          (fields(:,:,ie3,:)+fields(:,:,ie3-1,:))/2.)

          !! . . FFT to x
          !call SAX_Backward(sax,id)
          call RAX_Backward(rax,id)

          !! . . Taper
          !call rwetaper(sax,id) 
          call rwetaper(rax,id) 

          !! . . Imaging condition
          do ie2=1,ne2
             do ie1=1,ne1
                Timg(ie1,ie2,ie3,id)=Timg(ie1,ie2,ie3,id) +&
                     real( rax(ie1,ie2,id) ) 
!                Timg(ie1,ie2,ie3,id)=Timg(ie1,ie2,ie3,id)+&
!                     real((rax(ie1,ie2,id))*conjg(sax(ie1,ie2,id)) ) 
             end do
          end do

       end do Extrap_loop

       if (verbose .and. id .eq. 1 ) &
       write(0,'(a,i5,i5,i5,2F8.2,2F11.4)')  &
       '         * ST: F%,is,ns,w,im: ',floor(100*dble(iw)/dble(nw)),ishot,ns1*ns2,ww(id)/2./pi,maxval(Timg(:,:,:,id))

    end do Frequency_loop
    !$OMP END PARALLEL DO
    write(0,*) ' '

    !! SUMMATION OVER OUTSIDE IMAGE CONDITION
    do ith=1,nth
       !$OMP PARALLEL DO PRIVATE(ie3,ie2,ie1) 
       do ie3=1,ne3
          do ie2=1,ne2
             do ie1=1,ne1
                RImg(ie1,ie2,ie3)=RImg(ie1,ie2,ie3) + Timg(ie1,ie2,ie3,ith)
             end do
          end do
       end do
       !$OMP END PARALLEL DO
    end do

   !$OMP PARALLEL DO PRIVATE(ie3,ie2,ie1,ith) 
    do ith=1,nth
       do ie3=1,ne3
          do ie2=1,ne2
             do ie1=1,ne1
                Timg(ie1,ie2,ie3,ith) = 0.
             end do
          end do
       end do
    end do
   !$OMP END PARALLEL DO


    call shot_ker_release()

  end subroutine wemig

  !----------------------------------------------------------------
  !! . . Apply Taper to wavefield
  subroutine rwetaper(dax,id)
    complex :: dax(:,:,:)
    integer :: ie1,ie2,id
    do ie2=1,ne2
       do ie1=1,ne1
          dax(ie1,ie2,id) = dax(ie1,ie2,id) * tap(ie1,ie2)
       end do
    end do
  end subroutine rwetaper

  !----------------------------------------------------------------
  !! . . Read in the parameters required for setting up geometry...etc
  subroutine Read_In_Parameters(aE3,verbose,norm,nsx,nsy,nsz,eps,minang,&
                                maxang,method,newimage,nxint,nyint)
    type(myaxis) :: aE3
    integer:: norm,nsx,nsy,nsz,method,nxint,nyint
    real :: eps,minang,maxang,amax,b,c
    logical :: verbose,newimage
    call from_param("norm",norm,1)  
    call from_param("nsx",nsx,2)
    call from_param("nsy",nsy,2)
    call from_param("nsz",nsz,2)
    call from_param("ne3",ne3,500)
    call from_param("oe3",oe3,0.) 
    call from_param("de3",de3,0.003)
    call from_param("eps",eps,0.001)
    call from_param("minang",minang, 10.)
    call from_param("maxang",maxang,170.)
    call from_param("method",method,3)
    call from_param("verbose",verbose,.false.) 
    call from_param("newimage",newimage,.false.)
    call from_param("nxint",nxint,3)
    call from_param("nyint",nyint,3)

    aE3%n=ne3
    aE3%o=oe3
    aE3%d=de3
  end subroutine Read_In_Parameters

  !----------------------------------------------------------------
  !! . . Read in the header information describing wavefield geometry
  subroutine RWF_Info(aE1,aE2,aW,aS1,aS2,aR1,aR2)
    type(myaxis) ::   aE1,aE2,aW,aS1,aS2,aR1,aR2
    character(len=128) :: name
    name = 'rwf'
    call ainaxis(aE1,1,name)
    call ainaxis(aE2,2,name)
    call ainaxis(aW ,3,name)
    call ainaxis(aS1,4,name)
    call ainaxis(aS2,5,name)
    call ainaxis(aR1,1,name)
    call ainaxis(aR2,2,name)

  end subroutine RWF_Info

  !----------------------------------------------------------------
  !! . . Read in the header information describing wavefield geometry
  subroutine IMAGE_Out_Info(aX1,aX2,aX3)
    type(myaxis) ::     aX1,aX2,aX3
    character(len=128) :: name
    name = 'image'
    call aouaxis(aX3,1,name)
    call aouaxis(aX1,2,name)
    call aouaxis(aX2,3,name)
    call auxputch("esize","i",4,"image")
  end subroutine IMAGE_Out_Info

  !----------------------------------------------------------------
  !! . . Read in the header information describing wavefield geometry
  subroutine IMAGE_Info(aX1,aX2,aX3)
    type(myaxis) ::     aX1,aX2,aX3
    character(len=128) :: name
    name = 'image'
    call ainaxis(aX3,1,name)
    call ainaxis(aX1,2,name)
    call ainaxis(aX2,3,name)
  end subroutine IMAGE_Info

 !----------------------------------------------------------------
 !! . . Read in the header information describing Cartesian velocity/image space
  subroutine VEL_Info(aX1,aX2,aX3)
    type(myaxis) :: aX1,aX2,aX3
    character(len=128) :: name
    name = 'vel'
    call ainaxis(aX3,1,name)
    call ainaxis(aX1,2,name)
    call ainaxis(aX2,3,name)
  end subroutine VEL_Info

 !----------------------------------------------------------------
 !! . . Write out header information on the image
  subroutine write_out(aX1,aX2,aX3,aE1,aE2,aE3,aS1,aS2,aR1,aR2)
    type(myaxis) ::    aX1,aX2,aX3,aE1,aE2,aE3,aS1,aS2,aR1,aR2
    character (len=128 ) :: name

    !! . . Create output header file associated with "image" tag
    name= "RImage.H"
    call aouaxis(aE1,1,name)
    call aouaxis(aE2,2,name)
    call aouaxis(aE3,3,name)
    call auxputch("esize","i",4,"image")

    !! . . Create output header file associated with the rays (QC)
    name= "rays"
    call aouaxis(aE1,1,name)
    call aouaxis(aE2,2,name)
    call aouaxis(aE3,3,name)
    call auxputch("n4","i",3,name)
    call auxputch("o4","r",1.,name)
    call auxputch("d4","r",1.,name)
    call auxputch("esize","i",4,name)

    !! . . Write(0,*) out taper
    name="taper.H"
    call aouaxis(aE1,1,name)
    call aouaxis(aE2,2,name)
    call auxputch("esize","i",4,name)

  end subroutine write_out

 !----------------------------------------------------------------
 !! . . Zero a bunch of important arrays
  subroutine Zero_arrays(RVel,CVel,RImg,CImg,ncz,ncx,ncy,aE1,aE2,aE3)
    type(myaxis) :: aE1,aE2,aE3
    integer :: ix,iy,iz,ie1,ie2,ie3,ncz,ncx,ncy
    real    :: RVel(:,:,:),CVel(:,:,:),CImg(:,:,:),RImg(:,:,:)

    !$OMP PARALLEL DO PRIVATE(ie3,ie2,ie1)
    do ie3=1,aE3%n
       do ie2=1,aE2%n
          do ie1=1,aE1%n
             RVel(ie1,ie2,ie3)=0.
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO PRIVATE(ie3,ie2,ie1)
    do ie3=1,aE3%n
       do ie2=1,aE2%n
          do ie1=1,aE1%n
             RImg(ie1,ie2,ie3)=0.
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO PRIVATE(ix,iy,iz)
    do iy=1,ncy
       do ix=1,ncx
          do iz=1,ncz
             CVel(iz,ix,iy)=0.
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO PRIVATE(ix,iy,iz)
    do iy=1,ncy
       do ix=1,ncx
          do iz=1,ncz
             CImg(iz,ix,iy)=0.
          end do
       end do
    end do
    !$OMP END PARALLEL DO


  end subroutine Zero_arrays
  
 !----------------------------------------------------------------
 !! . . Set up the Cartesian geometry for migration
  subroutine CART_Setup(raysin,fields,gnorm,n3,o3,d3)
    integer :: iE1,iE2,iE3,n3
    real:: o3,d3
    real :: raysin(:,:,:,:),fields(:,:,:,:),gnorm(:,:,:)


    !! Do Cartesian coordinate setup
    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
    do iE3=1,ne3
       do iE2=1,ne2
          do iE1=1,ne1
             raysin(iE1,iE2,iE3,1) =  dble(iE1-1.)*de1 + oe1 + os1 
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
    do iE3=1,ne3
       do iE2=1,ne2
          do iE1=1,ne1
             raysin(iE1,iE2,iE3,2) =  dble(iE2-1.)*de2 + oe2 + os2 
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
    do iE3=1,ne3
       do iE2=1,ne2
          do iE1=1,ne1
             raysin(iE1,iE2,iE3,3) =  dble(iE3-1.)*de3 + oe3 
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
    do iE3=1,ne3
       do iE2=1,ne2
          do iE1=1,ne1
             fields(iE1,iE2,iE3,:)= 1.
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
    do iE3=1,ne3
       do iE2=1,ne2
          do iE1=1,ne1
             gnorm(iE1,iE2,iE3)= 1.
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    write(0,*) '        ************************************************'
    write(0,*) '        * '
    write(0,*) '        *      3D  Cartesian Migration'
    write(0,*) '        * '
    write(0,*) '        *            written by Jeff Shragge ' 
    write(0,*) '        *  '
    write(0,*) '        ************************************************'
    write(0,*) ' ' 
    write(0,*) '        ************************************************'
    write(0,*) '        * CARTESIAN GRID CHECK '
    write(0,'(a,2f11.2,2f11.2)') '         * X-dir',minval(raysin(:,:,:,1)),maxval(raysin(:,:,:,1))
    write(0,'(a,2f11.2,2f11.2)') '         * Y-dir',minval(raysin(:,:,:,2)),maxval(raysin(:,:,:,2))
    write(0,'(a,2f11.2,2f11.2)') '         * Z-dir',minval(raysin(:,:,:,3)),maxval(raysin(:,:,:,3))
    write(0,*) '        ************************************************'
    write(0,*) ' '

  end subroutine CART_Setup
  
 !----------------------------------------------------------------
 !! . .Set up the cylindrical elliptic  coordinate system for migration
  subroutine CYNEL_Setup(raysin,fields,gnorm,n2,o2,d2)
    integer              :: ii,iE1,iE2,iE3,n2
    real                 :: o2,d2
    real                 :: fields(:,:,:,:),raysin(:,:,:,:),gnorm(:,:,:)
    real,    allocatable :: rline(:),tline(:)
    allocate( tline(ne3),rline(ne2) )

    !! . . Define the extrapolation axis
    if (ne2 .ne. 1) then
       d2 = (maxang-minang)/180.*pi/dble(ne2-1.)
       o2 = pi*minang/180.
       r0 = dble(ne2-1.)*de2/2.
    else
       d2 = (maxang-minang)/180.*pi/dble(ne3-1.)
       o2 = pi*minang/180.
       r0 = dble(ne3-1.)*de1/2.
    end if

    tline=(/ (iE3,iE3=0,ne3-1) /)*de3+oe3
    if (ne2 .eq. 1) then
       rline(1) = pi/2.
    else
       rline=(/ (iE2,iE2=0,ne2-1) /)*d2+o2
    end if


    !! -------------------------------------
    !! . . COORDINATES

    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
    do iE3=1,ne3
       do iE2=1,ne2
          do iE1=1,ne1
             raysin(iE1,iE2,iE3,1)=dble(iE1-1.)*de1+oe1+os1 
          end do
       end do
    end do
    !$OMP END PARALLEL DO


    if (ne2 .eq. 1) then
       !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
       do iE3=1,ne3
          do iE2=1,ne2
             do iE1=1,ne1
                raysin(iE1,iE2,iE3,2)=-abs(r0)*cosh(tline(iE3))*cos(rline(iE2))+oe2+os2
             end do
          end do
       end do
       !$OMP END PARALLEL DO
    else
       !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
       do iE3=1,ne3
          do iE2=1,ne2
             do iE1=1,ne1
                raysin(iE1,iE2,iE3,2)=-abs(r0)*cosh(tline(iE3))*cos(rline(iE2))+oe2+os2+r0
             end do
          end do
       end do
       !$OMP END PARALLEL DO
    end if
    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
    do iE3=1,ne3
       do iE2=1,ne2
          do iE1=1,ne1
             raysin(iE1,iE2,iE3,3)=abs(r0)*sinh(tline(iE3))*sin(rline(iE2))
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    !! -------------------------------------
    !! . . FIELDS 

    !! . . Slowness stretch factor
    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
    do iE3=1,ne3
       do iE2=1,ne2
          do iE1=1,ne1
             fields(iE1,iE2,iE3,1)= sqrt(sinh(tline(iE3))**2+sin(rline(iE2))**2)
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    !! . . Cartesian Inline (x-direction)
    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
    do iE3=1,ne3
       do iE2=1,ne2
          do iE1=1,ne1
             fields(iE1,iE2,iE3,2)= sqrt(sinh(tline(iE3))**2+sin(rline(iE2))**2)
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    !! . . Elliptic Crossline (y-direction)
    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
    do iE3=1,ne3
       do iE2=1,ne2
          do iE1=1,ne1
             fields(iE1,iE2,iE3,3)= 1.
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    !! . . GNORM
    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
    do iE3=1,ne3
       do iE2=1,ne2
          do iE1=1,ne1
             gnorm(iE1,iE2,iE3)= sqrt(sinh(tline(iE3))**2+sin(rline(iE2))**2)
          end do
       end do
    end do
    !$OMP END PARALLEL DO


    write(0,*) '        ************************************************'
    write(0,*) '        * 3D RWE migration in CYLINDRICAL ELLIPIC coordinates '
    write(0,*) '        * '
    write(0,*) '        *       written by Jeff Shragge ' 
    write(0,*) '        *  '
    write(0,*) '        ************************************************'
    write(0,*) ' ' 
    write(0,*) '        ************************************************'
    write(0,*) '        * ELLIPTIC GRID CHECK:'
    write(0,'(a,2f11.2,2f11.2)') '         * MIN/MAX Angles  :',minval(rline)*180./pi,maxval(rline)*180./pi
    write(0,'(a,2f11.2,2f11.2)') '         * MIN/MAX Ex step :',minval(tline),maxval(tline)
    write(0,'(a,2f11.2,2f11.2)') '         * Rline min/max   :',minval(rline)*180./pi,maxval(rline)*180./pi
    write(0,*) '        ************************************************'
    write(0,*) ' '
    deallocate( tline, rline )

  end subroutine CYNEL_Setup
  
 !----------------------------------------------------------------
 !! . .Set up the cylindrical elliptic  coordinate system for migration
  subroutine ELLIP_Setup(raysin,fields,gnorm,aE1,aE2,aE3,aR1,aR2)
    type(myaxis)      :: aE1,aE2,aE3,aR1,aR2
    integer           :: ii,iE1,iE2,iE3,n2,chunk
    real              :: o2,d2,F,E,Z,aa,bb,gg,m,mp,tol,Kp,K,b0
    real              :: b,c,de1,de2,de3,sn,cn,dn,am,bmax,gmax,amax
    real              :: fields(:,:,:,:),raysin(:,:,:,:),denom,gnorm(:,:,:)
    real, allocatable :: e1(:),e2(:),e3(:),se1(:),se2(:)
    character (len=128 ) :: name
    allocate( e3(ne3),e1(ne1),e2(ne2) )

    call from_param("b",b,0.99)
    call from_param("c",c,1.)

    call from_param("bmax",bmax,0.95)
    call from_param("gmax",gmax,0.5)
    call from_param("amax",amax,0.8)

    write(0,*) 'USING b,c,gmax,bmax,amax: ',b,c,gmax,bmax,amax

    !! Coefficients
    m=b/c
    mp = sqrt(1.-b**2/c**2)
    tol = 2.2204*10.**(-16)
    call from_param("chunk",chunk,400)


    call elliptic_setup(chunk,tol)
    call elliptic12(pi/2.,m, K ,E,Z)
    call elliptic12(pi/2.,mp,Kp,E,Z)

    !! -----------------------------------------------------------------
    !!
    !! . . Define the beta (e1), gamma (e2) and alpha(e3) directions

    !! . . This is actually the crossline

    if (ne1 .eq. 1) then
       call ellipj(0.,mp,sn,cn,dn,am)
       e1(1) = b/dn
       de1=1.
       aE1%d=1.
    else 

       de1 = 2.*bmax*sqrt(mp)/dble(ne1-1.)

       do iE1 = 1,ne1/2
          call ellipj(bmax*sqrt(mp)-(iE1-1.)*de1,mp,sn,cn,dn,am)
          e1(iE1) = b/dn
       end do

       !! . . Set end
       call ellipj(bmax*sqrt(mp),mp,sn,cn,dn,am)
       e1(ne1) = b/dn
       aE1%d = abs(2.*bmax*sqrt(mp)/dble(ne1-1.))
    end if

    do iE1 = 1,ne1/2
       e1(ne1-iE1+1) = e1(iE1)
    end do
!!$


!!$    if (ne1 .eq. 1) then
!!$       call ellipj(0.,mp,sn,cn,dn,am)
!!$       e1(1)=b/dn
!!$       de1=1.
!!$       aE1%d=1.
!!$    else
!!$       !! Alternate
!!$       !! . . zero
!!$       call ellipj( bmax*sqrt(mp) ,mp,sn,cn,dn,am)
!!$       e1(  1) =-b/dn+b
!!$       e1(ne1) = b/dn-b
!!$       de1 = 2.*(b/dn-b)/ dble(ne1-1.)
!!$       do iE1=2,ne1-1
!!$          e1(iE1) = (iE1-1.)*de1+e1(1)
!!$       end do
!!$       e1 = abs(e1)+b
!!$       aE1%d = de1
!!$    end if
!!$
!!$    do iE1 = 1,ne1/2
!!$       e1(ne1-iE1+1) = e1(iE1)
!!$    end do


!!$    !! --------------------------------------------------------------
    !! . . This is actually the inline

!!$    if (ne2 .eq. 1) then
!!$       call ellipj(0.,m,sn,cn,dn,am)
!!$       e2(1) = b*sn
!!$       de2=1.
!!$       aE2%d=1.
!!$    else
!!$       !! . . Gamma 
!!$       de2 = 2.*gmax*K/dble(ne2)
!!$       
!!$       !! . . Set zero
!!$       call ellipj(0.,m,sn,cn,dn,am)
!!$       e2(ne2/2) = b*sn
!!$       do iE2 = 1,ne2/2
!!$          bb = (iE2-1.)*de2
!!$          call ellipj(bb,m,sn,cn,dn,am)
!!$          e2(ne2/2+iE2-1) = b*sn
!!$          e2(ne2/2-iE2+1) = b*sn
!!$       end do
!!$       !! . . Set end
!!$       call ellipj(gmax*K,m,sn,cn,dn,am)
!!$       e2(ne2) = b*sn
!!$       aE2%d=de2
!!$    end if

    !! Alternate
    !! . . zero
    if (ne2 .eq. 1) then
       call ellipj(0.,m,sn,cn,dn,am)
       e2(1) = b*sn
       de2=1.
       aE2%d=1.
    else

       e2(1) = -gmax*b
       de2 = 2.*b*gmax / dble(ne2-1.)

       do iE2=2,ne2/2
          e2(iE2) = (iE2-1.)*de2+e2(1)
       end do

       do iE2 = 1,ne2/2
          e2(ne2-iE2+1) = e2(iE2)
       end do
    end if
    e2 = abs(e2)
    aE2%d =abs( de2 )

    !! --------------------------------------------------------------

!!$    !! . . Alpha
!!$    de3 = amax*K/dble(ne3-1.)
!!$    aE3%d = de3
!!$    do iE3 = 1,ne3
!!$       aa=(iE3-1.)*de3
!!$       call ellipj(aa,m,sn,cn,dn,am)
!!$       e3(iE3) = c*dn/cn-c
!!$    end do
!!$ 
    !! . . Get first
    call ellipj(0.,m,sn,cn,dn,am)
    e3(  1) = c *dn/cn -c

    !! . . Get last
    call ellipj(amax*K,m,sn,cn,dn,am)
    e3(ne3) = c *dn/cn -c

    !! . . Linear interpolation between
    de3 = (e3(ne3)-e3(1))/dble(ne3-1.)
    do iE3 = 2,ne3
       e3(iE3) = ((iE3-1.)*de3+e3(1))**2
    end do

    aE3%d = de3


    !! . . X-Direction
    denom = (b*c)**2
    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
    do iE3=1,ne3
       do iE2=1,ne2
          do iE1=1,ne1
             raysin(iE1,iE2,iE3,1) = sign(sqrt(abs((e3(iE3)+c)**2*&
             e1(iE1)**2*e2(iE2)**2/denom)),dble(iE2-ne2/2.-1.))  
          end do
       end do
    end do
    !$OMP END PARALLEL DO



    !! . . Y-Direction
    denom = (b**2*(c**2-b**2))
    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
    do iE3=1,ne3
       do iE2=1,ne2
          do iE1=1,ne1
             raysin(iE1,iE2,iE3,2) = sign(sqrt(abs(((e3(iE3)+c)**2-b**2)*&
             (e1(iE1)**2-b**2)*(b**2-e2(iE2)**2)/denom)),dble(iE1-ne1/2-1.))

          end do
       end do
    end do
    !$OMP END PARALLEL DO

    !! . . Z-Direction
    denom = (c**2*(c**2-b**2))
    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
    do iE3=1,ne3
       do iE2=1,ne2
          do iE1=1,ne1
             raysin(iE1,iE2,iE3,3) = sqrt(abs(((e3(iE3)+c)**2-c**2)*(c**2-e1(iE1)**2)*&
             (c**2-e2(iE2)**2)/denom))
          end do
       end do
    end do
    !$OMP END PARALLEL DO

    !! --------------------------------------------------------------
    !! . . Define r0
    !! . . Fix scaling to be equal to input grid
    if      (ne2 .eq. 1 ) then
       r0 = (aR1%n-1.)*aR1%d / (raysin(ne1,1,1,1)-raysin(1,1,1,1))
       aE1%d =    aR1%d
       aE2%d =    aR1%d
       aE3%d = r0*aE3%d
    else if (ne1 .eq. 1) then
       r0 = (aR2%n-1.)*aR2%d / (raysin(1,ne2,1,2)-raysin(1,1,1,2))
       aE2%d =    aR2%d
       aE1%d =    aR2%d
       aE3%d = r0*aE3%d
    else
       r0 = (aR2%n-1.)*aR2%d /( maxval(raysin(:,:,1,1))-minval(raysin(:,:,1,1)))
       aE1%d = r0*aE1%d
       aE2%d = r0*aE2%d
!       aE1%d = aR1%d
!       aE2%d = aR2%d
       aE3%d = r0*aE3%d
    end if

    raysin = r0 * raysin
    write(0,*) 'USING de1,de2,de3 of:',aE1%d,aE2%d,aE3%d
    call RWE3D_geometry_init(aE1,aE2,aE3)
    call RWE3D_geometry_calc(raysin,fields,gnorm)
    call RWE3D_geometry_close()

!!$    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
!!$    do iE3=1,ne3
!!$       do iE2=1,ne2
!!$          do iE1=1,ne1
!!$             fields(iE1,iE2,iE3,1)= 1./(sqrt(e1(iE1)**2-e2(iE2)**2))
!!$          end do
!!$       end do
!!$    end do
!!$    !$OMP END PARALLEL DO
!!$
!!$    !! . . Inline propagation
!!$    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
!!$    do iE3=1,ne3
!!$       do iE2=1,ne2
!!$          do iE1=1,ne1
!!$             fields(iE1,iE2,iE3,3)= sqrt((e3(iE3)+c)**2-e2(iE2)**2)!/&
!!$!                        (sqrt(e1(iE1)**2-e2(iE2)**2))
!!$          end do
!!$       end do
!!$    end do
!!$    !$OMP END PARALLEL DO
!!$    
!!$    !! . . Crossline propagation
!!$    !$OMP PARALLEL DO PRIVATE(iE3,iE2,iE1)
!!$    do iE3=1,ne3
!!$       do iE2=1,ne2
!!$          do iE1=1,ne1
!!$             fields(iE1,iE2,iE3,2)= sqrt((e3(iE3)+c)**2-e1(iE1)**2)!/&
!!$ !                       (sqrt(e1(iE1)**2-e2(iE2)**2))
!!$          end do
!!$       end do
!!$    end do
!!$    !$OMP END PARALLEL DO


    !! . . Create output header file associated with the rays (QC)
    name= "fields.H"
    call aouaxis(aE1,1,name)
    call aouaxis(aE2,2,name)
    call aouaxis(aE3,3,name)
    call auxputch("n4","i",3,name)
    call auxputch("o4","r",1.,name)
    call auxputch("d4","r",1.,name)
    call auxputch("esize","i",4,name)
    call srite(name,fields,4*size(fields))

    write(0,*) '        ************************************************'
    write(0,*) '        * 3D RWE migration in ELLIPSOIDAL coordinates '
    write(0,*) '        * '
    write(0,*) '        *       written by Jeff Shragge ' 
    write(0,*) '        *  '
    write(0,*) '        ************************************************'
    write(0,*) ' ' 
    write(0,*) '        ************************************************'
    write(0,*) '        * ELLIPSOIDAL GRID CHECK:'
    write(0,'(a,2f11.4,2f11.4)') '         * MIN/MAX E1:',minval(e1),maxval(e1)
    write(0,'(a,2f11.4,2f11.4)') '         * MIN/MAX E2:',minval(e2),maxval(e2)
    write(0,'(a,2f11.4,2f11.4)') '         * MIN/MAX E3:',minval(e3),maxval(e3)
    write(0,*) '        ************************************************'
    write(0,*) ' '

    call deallocate_elliptic()
    deallocate (e1,e2,e3)
  end subroutine ELLIP_Setup

 !----------------------------------------------------------------
 !! . . Write out information for visual verification
  subroutine Geometry_QC(raysin,aE1,aE2,aE3,aS1,aS2,aW,aR1,aR2)
    type(myaxis)             :: aE1,aE2,aE3,aS1,aS2,aW,aR1,aR2
    real :: raysin(:,:,:,:)
    character(len=128) :: name
    !! . . Screen output for QC
    write(0,*) '        ************************************************'
    write(0,*) '        * GEOMETRY CHECK'
    write(0,'(a,2f11.2,2f11.2)')    "         * E1-coord min/max         :", minval(raysin(:,:,:,1)),maxval(raysin(:,:,:,1))
    write(0,'(a,2f11.2,2f11.2)')    "         * E2-coord min/max         :", minval(raysin(:,:,:,2)),maxval(raysin(:,:,:,2))
    write(0,'(a,2f11.2,2f11.2)')    "         * E3-coord min/max         :", minval(raysin(:,:,:,3)),maxval(raysin(:,:,:,3))
    write(0,'(a,i4,2f11.2,2f11.2)') '         * RECEIVERS inline    n,o,d:',nr1,or1,dr1
    write(0,'(a,i4,2f11.2,2f11.2)') '         * RECEIVERS Xline     n,o,d:',nr2,or2,dr2
    write(0,'(a,i4,2f11.2,2f11.2)') '         * EXTRAPOLATION STEPS n,o,d:',ne3,oe3,de3
    write(0,'(a,i4,2f11.2,2f11.2)') '         * SHOTS inline        n,o,d:',aS1%n,aS1%o,aS1%d
    write(0,'(a,i4,2f11.2,2f11.2)') '         * SHOTS Xline         n,o,d:',aS2%n,aS2%o,aS2%d
    write(0,'(a,i4,2f11.2,2f11.2)') '         * FREQUENCIES         n,o,d:',aW%n,aW%o,aW%d
    write(0,*) '        ************************************************'
    write(0,*) ' '
  end subroutine Geometry_QC

 !----------------------------------------------------------------
 !! . . Write out velocity information for visual verification
  subroutine Velocity_Check(CVel,RVel)
    real :: RVel(:,:,:),CVel(:,:,:)
    !! . . Screen output for QC
    write(0,*) '        *************************************************'
    write(0,*) '        * VELOCITY CHECK'
    write(0,'(a,i4,i4,i4,2f11.2,2f11.2)')      '         * Cart VEL n1,n2,n3,min,max:',size(CVel,2),size(CVel,3),size(CVel,1),minval(CVel),maxval(CVel)
    write(0,'(a,i4,i4,i4,2f11.2,2f11.2)')      '         * RWE  VEL n1,n2,n3,min,max:',size(RVel,1),size(RVel,2),size(RVel,3),minval(RVel),maxval(RVel)
    write(0,*) '        ************************************************'
    write(0,*) ' '
  end subroutine Velocity_Check

 !----------------------------------------------------------------
 !! . . Reposition the wavefield on elliptic coordinates
  subroutine CYNEL_Position_XWF(xwfin,xwfout,raysin)
    integer :: iE1,iX1,iE2,iX2
    real :: raysin(:,:,:,:),y,dx1,dx2,dd1,dd2
    complex  :: xwfin(:,:,:),xwfout(:,:,:)

    do iE2 = 1,ne2
       !! . . Generate interpolation coefficients
       !! . . independent of y-direction
       y = raysin(1,iE2,1,2) - os2
       iX2 = floor( ( y - or2 ) / dr2) + 1
       dx1 = y - (or2 + dble(iX2-1.)*dr2 )
       dx2 = (or2 + dble(iX2)*dr2 )-y
       dd1 = dx2/(dx1+dx2)
       dd2 = dx1/(dx1+dx2)

       !! . . Boundary Check
       if (iX2 .ge. ne2) iX2 = ne2-1
       if (iX2 .le. 1 )  iX2 = 1

       if (ne2 .eq. 1) then
          xwfout(:,1,:) = xwfin(:,1,:)
       else
          do iE1 = 1,ne1
             xwfout(iE1,iX2,:) = dd1*xwfin(iE1,iX2,:)+dd2*xwfin(iE1,iX2+1,:)
          end do
       end if

    end do
    write(0,*) '        ************************************************'
    write(0,*) '        * WAVEFIELD CHECK'
    write(0,*) '        * MIN/MAX Wavefield: ',minval(real(xwfout)),maxval(real(xwfout))
    write(0,*) '        ************************************************'
    write(0,*) ' '
 end subroutine CYNEL_Position_XWF

 !----------------------------------------------------------------
 !! . . Reposition the wavefield on elliptic coordinates
 subroutine ELLIP_Position_XWF(xwfin,xwfout,raysin,ar1,ar2)
   type (myaxis) :: ar1,ar2
   integer :: iE1,iE2,iX1,iX2
   real :: raysin(:,:,:,:),x,y,dx1,dx2,ddx1,ddx2,dy1,dy2,ddy1,ddy2
   complex  :: xwfin(:,:,:),xwfout(:,:,:)

   xwfout = 0.
   if (ne2 .eq. 1) then

      do iE1 = 1,ne1

         !! . . Generate interpolation coefficients
         x = raysin(iE1,1,1,1) - os1
         iX1 = floor( ( x - or1 ) / dr1) + 1
         dx1 = x - (ar1%o + dble(iX1-1.)*dr1 )
         dx2 = (ar1%o + dble(iX1)*dr1 )-x
         ddx1 = dx2/(dx1+dx2)
         ddx2 = dx1/(dx1+dx2)
         !! . . Generate interpolation coefficients

         iX2 = 1

         !! . . Boundary Check
         if (iX1 .ge. ne1) iX1 = ne1-1
         if (iX1 .le.   1) iX1 = 1
         xwfout(iX1,1,:) =  ddx1*xwfin(iX1,1,:)+ddx2*xwfin(iX1+1,1,:)

      end do

   else if (ne1 .eq. 1) then

      do iE2 = 1,ne2
         !! . . Generate interpolation coefficients
         y = raysin(1,iE2,1,2) - os2
         iX2 = floor( ( y - or2 ) / dr2) + 1
         dx1 = y - (ar2%o + dble(iX2-1.)*dr2 )
         dx2 = (ar2%o + dble(iX2)*dr2 )-y
         ddx1 = dx2/(dx1+dx2)
         ddx2 = dx1/(dx1+dx2)

         ! write(0,*) iE2,iX2,raysin(1,iE2,1,2),dr2,or2
         !! . . Boundary Check
         if (iX2 .ge. ne2) iX2 = ne2-1
         if (iX2 .le.   1) iX2 = 1

         xwfout(1,iX2,:) =  ddx1*xwfin(1,iX2,:)+ddx2*xwfin(1,iX2+1,:)
      end do

   else 

      do iE1 = 1,ne1
         do iE2 = 1,ne2

            !! . . Generate interpolation coefficients
            x = raysin(iE1,iE2,1,1) - os1
            iX1 = floor( ( x - or1 ) / dr1) + 1
            dx1 = x - (ar1%o + dble(iX1-1.)*dr1 )
            dx2 = (ar1%o + dble(iX1)*dr1 )-x
            ddx1 = dx2/(dx1+dx2)
            ddx2 = dx1/(dx1+dx2)

            !! . . Generate interpolation coefficients
            y = raysin(iE1,iE2,1,2) - os2
            iX2 = floor( ( y - or2 ) / dr2) + 1
            dy1 = y - (ar2%o + dble(iX2-1.)*dr2 )
            dy2 = (ar2%o + dble(iX2)*dr2 ) - y
            ddy1 = dy2/(dy1+dy2)
            ddy2 = dy1/(dy1+dy2)

            !! . . Boundary Check
            if (iX1 .ge. ne1) iX1 = ne1-1
            if (iX1 .le.     1)  iX1 = 1
            !! . . Boundary Check
            if (iX2 .ge. ne2) iX2 = ne2-1
            if (iX2 .le.     1)  iX2 = 1

            xwfout(iX1,iX2,:) =  xwfout(iX1,iX2,:) &
            +ddx1*ddy1*xwfin(iX1  ,iX2    ,:) &
            +ddx1*ddy2*xwfin(iX1  ,iX2+1  ,:) &
            +ddx2*ddy1*xwfin(iX1+1,iX2    ,:) &
            +ddx2*ddy2*xwfin(iX1+1,iX2+1  ,:)
         end do
      end do

   end if

   write(0,*) '        ************************************************'
   write(0,*) '        * WAVEFIELD CHECK'
   write(0,*) '        * MIN/MAX Wavefield: ',minval(real(xwfout)),maxval(real(xwfout))
   write(0,*) '        ************************************************'
   write(0,*) ' '
 end subroutine ELLIP_Position_XWF
 !----------------------------------------------------------------
 !! . . Find the corresponding Cartesian velocity indicies
 subroutine Get_Vel_Indicies(raysin,iix1,iix2,oo1,oo2)
   integer :: iix1(:),iix2(:)
   real :: oo1(:),oo2(:)
   real :: ee,raysin(:,:,:,:)
   real :: x_tmp,y_tmp

   !! X1 index
   oo1 (1) = minval(raysin(:,:,:,1))
   oo1 (2) = maxval(raysin(:,:,:,1))
   iix1(1) = floor   ((oo1(1) - ox1)/dx1) + 1
   iix1(2) = ceiling ((oo1(2) - ox1)/dx1) 
   iix1(1) = max(1,min(iix1(1),nx1))
   iix1(2) = max(1,min(iix1(2),nx1))

   !! X2 index
   oo2 (1) = minval(raysin(:,:,:,2))
   oo2 (2) = maxval(raysin(:,:,:,2))
   iix2(1) = floor   ((oo2(1) - ox2)/dx2) + 1
   iix2(2) = ceiling ((oo2(2) - ox2)/dx2) 
   iix2(1) = max(1,min(iix2(1),nx2))
   iix2(2) = max(1,min(iix2(2),nx2))

 end subroutine Get_Vel_Indicies

 !----------------------------------------------------------------

 subroutine Get_Img_Indicies(raysin,iix1,iix2,oo1,oo2,vox1,vox2)
   integer :: iix1(:),iix2(:)
   real :: oo1(:),oo2(:),vox1(:),vox2(:)
   real :: raysin(:,:,:,:)

write(0,*) vox1(1),vox1(2),vox2(1),vox2(2)
   !! X1 index
   oo1 (1) = max(minval(raysin(:,:,:,1)),ox1)
   oo1 (2) = min(maxval(raysin(:,:,:,1)),(nx1-1.)*dx1+ox1)
!   oo1 (1) = minval(raysin(:,:,:,1))
!   oo1 (2) = maxval(raysin(:,:,:,1))
   iix1(1) = floor  ( (oo1(1)-vox1(1))/dx1) + 1 
   iix1(2) = ceiling( (oo1(2)-vox1(1))/dx1) 
!   iix1(1) = floor  ( oo1(1)/dx1) + 1 
!   iix1(2) = ceiling( oo1(2)/dx1) 
   iix1(1) = max(1,min(iix1(1),nx1))
   iix1(2) = max(1,min(iix1(2),nx1))

   !! X2 index
   oo2 (1) = max(minval(raysin(:,:,:,2)),ox2)
   oo2 (2) = min(maxval(raysin(:,:,:,2)),(nx2-1.)*dx2+ox2)
!   oo2 (1) = minval(raysin(:,:,:,2))
!   oo2 (2) = maxval(raysin(:,:,:,2))
   iix2(1) = floor  ( (oo2(1)-vox2(1))/dx2) + 1
   iix2(2) = ceiling( (oo2(2)-vox2(1))/dx2) 

!   iix2(1) = floor  ( oo2(1)/dx2) + 1
!   iix2(2) = ceiling( oo2(2)/dx2) 

   iix2(1) = max(1,min(iix2(1),nx2))
   iix2(2) = max(1,min(iix2(2),nx2))

 end subroutine Get_Img_Indicies

  !----------------------------------------------------------------
  !! . . Taper Initialization Program
 subroutine rwetaper_init()
   integer :: i1,i2,j1,j2,ie1,ie2,i,jtap

   do ie2=1,ne2
      do ie1=1,ne1
         tap(ie1,ie2)=1.
      end do
   end do

   if (nxtap .gt. 1 .and. nytap .gt. 1) then

      do i2=1,nytap
         do i1=1,nxtap
            j1 = abs(nxtap-i1)
            j2 = abs(nytap-i2)
            tap(   i1 , i2 ) = cos(pi/2.* j1/(nxtap-1))*cos(pi/2.* j2/(nytap-1))
         end do
      end do

      do i2=1,nytap
         do i1=1,nxtap
            j1 = abs(nxtap-i1)
            j2 = abs(nytap-i2)
            tap(ne1-i1+1,i2) = cos(pi/2.* j1/(nxtap-1))*cos(pi/2.* j2/(nytap-1))
         end do
      end do

      do i2=1,nytap
         do i1=1,nxtap
            j1 = abs(nxtap-i1)
            j2 = abs(nytap-i2)
            tap(   i1 , ne2-i2 +1) = cos(pi/2.* j1/(nxtap-1))*cos(pi/2.* j2/(nytap-1))
         end do
      end do

      do i2=1,nytap
         do i1=1,nxtap
            j1 = abs(nxtap-i1)
            j2 = abs(nytap-i2)
            tap(ne1-i1+1,ne2-i2+1) = cos(pi/2.* j1/(nxtap-1))*cos(pi/2.* j2/(nytap-1))
         end do
      end do

      do i2=1,nytap
         do i1=nxtap+1,ne1-nytap
            j2 = abs(nytap-i2)
            tap(i1,i2)         = cos(pi/2.* j2/(nytap-1))
         end do
      end do

      do i2=1,nytap
         do i1=nxtap+1,ne1-nytap
            j2 = abs(nytap-i2)
            tap(i1,ne2-i2+1) = cos(pi/2.* j2/(nytap-1))
         end do
      end do

      do i2=nytap+1,ne2-nytap
         do i1=1,nxtap
            j1 = abs(nytap-i1)
            tap(i1,i2)         = cos(pi/2.* j1/(nxtap-1))
         end do
      end do

      do i2=nytap+1,ne2-nytap
         do i1=1,nxtap
            j1 = abs(nytap-i1)
            tap(ne1-i1+1 ,i2) = cos(pi/2.* j1/(nxtap-1))
         end do
      end do


   else if (nxtap .gt. 1 .and. nytap .le. 1) then

      do i1=1,nxtap
         j1 = abs(nxtap-i1)
         tap(   i1 , : ) = cos(pi/2.* j1/(nxtap-1))
      end do

      do i1=1,nxtap
         j1 = abs(nxtap-i1)
         tap(ne1-i1+1,:) = cos(pi/2.* j1/(nxtap-1))
      end do

   else if (nytap .gt. 1 .and. nxtap .le. 1) then

      do i2=1,nytap
         j2 = abs(nytap-i2)
         tap(:,  i2    ) = cos(pi/2.* j2/(nytap-1))
      end do

      do i2=1,nytap
         j2 = abs(nytap-i2)
         tap(:,ne2-i2+1) = cos(pi/2.* j2/(nytap-1))
      end do

   end if
   call srite("taper.H",tap,4*size(tap) )
 end subroutine rwetaper_init
 !----------------------------------------------------------------

end module  RWE3D_wemig
