!!$
!!$=head1 NAME
!!$
!!$ Phase - Phaseshift migration (possibly split-step)
!!$
!!$=head1 SYNOPSIS
!!$
!!$Phase par=parfile.P D=input.H  R=output.H S=slow.H
!!$
!!$=head1 COORDINATE SYSTEM FOR DOWNWARD CONTINUATION DESCRIPTION PARAMETERS
!!$
!!$=over 4
!!$
!!$
!!$=item amx__n (integer)
!!$
!!$  Number of CMP in X (axis 1 of D,R and S)
!!$
!!$=item amx_d (integer) 
!!$
!!$  Sampling of CMP in X
!!$
!!$=item amx_o (integer)
!!$
!!$  First CMP  X location 
!!$
!!$=item amy__n (integer)
!!$
!!$  Number of CMP in X (axis 2 of D,R and S)
!!$
!!$=item amy_d (integer) 
!!$
!!$  Sampling of CMP in Y
!!$
!!$=item amy_o (integer)
!!$
!!$  First CMP  Y location 
!!$
!!$=item ahx__n (integer)
!!$
!!$  Number of offsets in X (axis 3 of D and possibly R  )
!!$
!!$=item ahx_d (integer) 
!!$
!!$  Sampling of offsets in X
!!$
!!$=item ahx_o (integer)
!!$
!!$  First offset  X location 
!!$
!!$=item ahy__n (integer)
!!$
!!$  Number of offsets in Y (axis 4 of D and possibly R)
!!$
!!$=item ahy_d (integer) 
!!$
!!$  Sampling of offset in Y
!!$
!!$=item ahy_o (integer)
!!$
!!$  First CMP  Y location 
!!$
!!$
!!$=back
!!$
!!$
!!$=head1 R/S  DESCRIPTION PARAMETERS
!!$
!!$=over 4
!!$
!!$
!!$=item az__n (integer) (axis 3 of S, axis 5 of R)
!!$
!!$  Number of depth steps
!!$
!!$=item az__d (integer)
!!$
!!$  Sampling of depths
!!$
!!$=item az__o (integer)
!!$
!!$  First depth of model
!!$
!!$=item image_real (logcal)
!!$
!!$  Whether or not make the image real
!!$
!!$=back
!!$
!!$
!!$=head1 DATA DESCRIPTION PARAMETERS
!!$
!!$=over 4
!!$
!!$=item aw__n (integer) (axis 5 of D)
!!$
!!$  Number of frequencies
!!$
!!$=item aw__d (integer)
!!$
!!$  Sampling of frequencies
!!$
!!$=item aw__o (integer)
!!$
!!$  First frequency of data
!!$
!!$=back
!!$
!!$=head1 SLOWNESS PARAMETERS
!!$
!!$=over 4
!!$
!!$=item  velocity  (logical)
!!$
!!$ [.false.] whether or not S is velocity
!!$
!!$=item  nfk  (integer)
!!$
!!$ [1] number of velocities for each depth step 
!!$
!!$=back
!!$
!!$=head1 PH GATHERS
!!$
!!$=over 4
!!$ 
!!$=item pCIG (logical)
!!$
!!$ [.false.] whether (.true.) or not (.false.) to create PH gathers
!!$
!!$=item apx__n (integer) 
!!$
!!$  Number of PH in X (axis 3 of R)
!!$
!!$=item apx_d (integer) 
!!$
!!$  Sampling of PH in X
!!$
!!$=item apx_o (integer)
!!$
!!$  First PH X location 
!!$
!!$=item apy__n (integer)
!!$
!!$  Number of PH in Y (axis 4 of R)
!!$
!!$=item apy_d (integer) 
!!$
!!$  Sampling of PH in Y
!!$
!!$=item apy_o (integer)
!!$
!!$  First PH  Y location 
!!$
!!$=item pgiceps (float)
!!$
!!$  [.1] Epsilon parameter for Phgather construction
!!$
!!$=item vsalt (float)
!!$
!!$  [0.] Salt velocity 
!!$
!!$=back
!!$
!!$=head1 MEMORY MANAGEMET
!!$
!!$=over 4
!!$
!!$=item nzs (integer)
!!$
!!$Number of Z steps to do in memory
!!$
!!$=item nws (integer)
!!$
!!$Number of W to hold in memory
!!$
!!$=back
!!$
!!$=head1 DESCRIPTION
!!$
!!$ Perform split step phase shift migration/modeling
!!$
!!$
!!$
!!$
!!$=head1 COMMENTS
!!$
!!$=over 4
!!$
!!$=item  anx_n,ahx_d,ahx_o and axis 1 of D and R  must be the same
!!$
!!$=item  any_n,ahy_d,ahy_o and axis 2 of D and R  must be the same
!!$
!!$=item Slowness must contain entire  span of X and Y
!!$
!!$=item  az__n,az__d,az__o , axis 3 of S and axis 5 of R  must be the same
!!$
!!$=item When doing ph gathers apx_n = ahx_n = n3 of D
!!$
!!$=item When doing ph gathers apy_n = ahy_n = n4 of D
!!$
!!$=item When doing offset gathers  ahx_n,ahx_d,ahx_o = axis 3 of D and R
!!$
!!$=item When doing offset gathers  ahy_n,ahy_d,ahy_o = axis 4 of D and R
!!$
!!$=item nzs must be a multiple of az__n
!!$
!!$=item nws*nmpi_threads*nomp_threads must be a multiple of aw__n 
!!$
!!$=cut
!!$
!!$
!!$Authors:
!!$
!!$ Paul Sava
!!$
!!$ Bob Clapp
!!$
!!$

program Phase   
  use sep 
  use wei_util 
  use wei_process 
  use wei_mig 
  use wei_cig 
  use wei_kxmig 
  use wei_sll 
  use wei_slo 
  use wei_sls 
  use wei_ssf 
  use wei_wem 

  implicit none
  integer :: stat 

  integer :: i,n,ierr,iverb
  logical :: verb
  real :: vsalt
  integer :: nfk
  logical :: pCIG,sll


  !----------------------------------------------------------------
#ifdef SEP_MPI
  call MPI_INIT(stat)
  call  set_no_putch()
#endif
  call sep_init(SOURCE)
  !----------------------------------------------------------------
#ifdef SEP_MPI
  call MPI_COMM_SIZE(MPI_COMM_WORLD, n,ierr)
  call MPI_COMM_RANK(MPI_COMM_WORLD, i,ierr)
  if(i==0) then
     call from_param("verb",verb,.false.)
     iverb=0;if(verb) iverb=1
     call MPI_SEP_SEND_ARGS(n,10,iverb)
  else
     call MPI_SEP_RECEIVE_ARGS()
  end if
#endif


  call from_param("nfk",nfk,1)
  call from_param("vsalt",vsalt,0.)
  call from_param("pCIG",pCIG,.false.)
  call from_param("sll",sll,.false.)


  !----------------------------------------------------------------
  pro%operator=".R_.D"
  !----------------------------------------------------------------
if(sll) then

    call weimig_init(   SLin=weisll_init	&
    ,                WCin=weikxmign_init	&
    ,                   FKin=weiwem_init	&
    ,                   FXin=weissf_init	&
    ,                   IGin=weihcig_init	)
     stat = weimig(.true.,pro%add,pro%R,pro%D &
     ,              SLop=weisll	&
     ,           WCop=weikxmign	&
     ,              FKop=weiwem_fast	&
     ,              FXop=weissf	&
     ,              IGop=weihcig)

else

if(nfk==1)  then
  if(.not. pCIG) then
    write(0,*) "IN THIS ONb"

    call weimig_init(   SLin=weislo1_init	&
    ,                WCin=weikxmig1_init	&
    ,                   FKin=weiwem_init	&
    ,                   FXin=weissf_init	&
    ,                   IGin=weihcig_init	)
     stat = weimig(.true.,pro%add,pro%R,pro%D &
     ,              SLop=weislo1	&
     ,           WCop=weikxmig1	&
     ,              FKop=weiwem_fast	&
     ,              FXop=weissf	&
     ,              IGop=weihcig)
   else !1 velocity pCIG 
    write(0,*) "IN THIS ONa"
    call weimig_init(   SLin=weislo1_init	&
    ,                WCin=weikxmig1_init	&
    ,                   FKin=weiwem_init	&
    ,                   FXin=weissf_init	&
    ,                   IGin=weipcig_init	)
     stat = weimig(.true.,pro%add,pro%R,pro%D &
     ,              SLop=weislo1	&
     ,           WCop=weikxmig1	&
     ,              FKop=weiwem_fast	&
     ,              FXop=weissf	&
     ,              IGop=weipcig)
   end if
else if(vsalt==0.) then  !multiple velocities
  if(.not. pCIG) then
    call weimig_init(   SLin=weisloN_init	&
    ,                WCin=weikxmign_init	&
    ,                   FKin=weiwem_init	&
    ,                   FXin=weissf_init	&
    ,                   IGin=weihcig_init	)
     stat = weimig(.true.,pro%add,pro%R,pro%D &
     ,              SLop=weisloN	&
     ,           WCop=weikxmign	&
     ,              FKop=weiwem_fast	&
     ,              FXop=weissf	&
     ,              IGop=weihcig)
   else !1 velocity pCIG 
    call weimig_init(   SLin=weisloN_init	&
    ,                WCin=weikxmign_init	&
    ,                   FKin=weiwem_init	&
    ,                   FXin=weissf_init	&
    ,                   IGin=weipcig_init	)
     stat = weimig(.true.,pro%add,pro%R,pro%D &
     ,              SLop=weisloN	&
     ,           WCop=weikxmign	&
     ,              FKop=weiwem_fast	&
     ,              FXop=weissf	&
     ,              IGop=weipcig)
   end if
else  !multiple velocities

  if(.not. pCIG) then
    call weimig_init(   SLin=weisloN_init	&
    ,                WCin=weikxmign_init	&
    ,                   FKin=weiwem_init	&
    ,                   FXin=weissf_init	&
    ,                   IGin=weihcig_init	)
    write(0,*) "IN THIS ON1"
     stat = weimig(.true.,pro%add,pro%R,pro%D &
     ,              SLop=weisloN	&
     ,           WCop=weikxmign	&
     ,              FKop=weiwem_fast	&
     ,              FXop=weissf	&
     ,              IGop=weihcig)
   else !1 velocity pCIG 
    call weimig_init(   SLin=weisloN_init	&
    ,                WCin=weikxmign_init	&
    ,                   FKin=weiwem_init	&
    ,                   FXin=weissf_init	&
    ,                   IGin=weipcig_init	)
     stat = weimig(.true.,pro%add,pro%R,pro%D &
     ,              SLop=weisloN	&
     ,           WCop=weikxmign	&
     ,              FKop=weiwem_fast	&
     ,              FXop=weissf	&
     ,              IGop=weipcig)
   end if

end if


endif

#ifdef SEP_MPI
  call MPI_FINALIZE(ierr) 
#endif
  call sep_close()
  !----------------------------------------------------------------
  call exit(0)
end program Phase
