module wei_sls
  !! WEI slowness operator
  !! for FK-FX downward continuation
  !! includes salt logic

  use wei_util
  use wei_process
  use wei_obj_S
  !---------------
  implicit none

  real, dimension(:,:,:,:), pointer, private :: sxy

  real,    private :: vsalt, ssalt  
contains
!!$
!!$=head1 NAME
!!$
!!$ weislsN_init - Initialize multiple  slowness operator w/salt
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ierr=weislsN_init()
!!$
!!$=head1 DESCRIPTION
!!$
!!$ Initialize multiple slowness operator w/salt (from_param(vel_salt)
!!$
!!$=head1 SEE ALSO
!!$
!!$L<weisloN_init> L<weislo1> L<weislsN>
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$=cut
!!$

  !----------------------------------------------------------------
  !! N reference slownesses
  integer function weislsN_init() result(st)
    integer :: er
#ifdef DEBUG
    character(len=128), parameter :: unit='weislsN_init'
    if(pro%debug) call in(unit)
#endif

    allocate(sxy(pro%amx%n,pro%amy%n,pro%ahx%n,pro%ahy%n),stat=er)
    if(er/=0) call seperr('cant allocate sxy')

    call from_param("vsalt",vsalt,0.)
    if(vsalt >0.) ssalt=1./vsalt

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
    st=OK
  end function weislsN_init
!!$
!!$=head1 NAME
!!$
!!$ weislsN - Multiple  slowness operator w/salt
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ierr=weislsN()
!!$
!!$=head1 DESCRIPTION
!!$
!!$  multiple slowness operator w/salt 
!!$
!!$=head1 SEE ALSO
!!$
!!$L<weisloN_init> L<weislo1> L<weislsN_init>
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$=cut
!!$
  !----------------------------------------------------------------  
  !! N reference slownesses
  integer function weislsN() result(st)
    integer :: izs,ifk
    real    :: ds, so, sz
    real    :: sl,su
    real    :: smin,smax
    integer :: imx, imy
    real    :: ssalt_true
    real    :: ssalt_mean,ssalt_sum
    real    :: slocal
    integer :: ins,ndo
#ifdef DEBUG
    character(len=128), parameter :: unit='weislsN'
    if(pro%debug) call in(unit)
#endif

    !! spread the slowness in its PBA    
    call weislow_push(S)

    !! min/max slowness in the interval
    izsloop: do izs=1,pro%nzs
       slo%nrs(izs) = pro%nfk
       slo%si(:,:,:,:,izs)=1

       sxy(:,:,:,:) &
       = .5* real(S%pba(:,:,:,:,izs)) &
       + .5*aimag(S%pba(:,:,:,:,izs))

       smin = slowness_min(S,izs)
       smax = slowness_max(S,izs)

 !! salt
       if(vsalt >0. .and. smin <= ssalt .and. pro%nfk>1) then
          slo%nrs(izs) = pro%nfk+1

          !! nfk+1 goes to salt slowness
          !! find mean salt velocity
          ssalt_sum=0.
          ins=0
          do imy=1,size(S%s1%bin,2)
             do imx=1,size(S%s1%bin,1)
                slocal = real(S%s1%bin(imx,imy,izs))
                if(slocal <= ssalt) then
                   ins = ins+1
                   ssalt_sum = ssalt_sum + slocal
                end if
             end do
          end do
             ssalt_mean = ssalt_sum/ins
          slo%ssz(izs,pro%nfk+1) = ssalt_mean
          slo%srz(izs,pro%nfk+1) = ssalt_mean
!          where(sxy <= ssalt)
!          end where

          !! update min slowness in sediments region
          smin=1.e+6
          do imy=1,size(S%s1%bin,2)
             do imx=1,size(S%s1%bin,1)
                slocal = real(S%s1%bin(imx,imy,izs))
                !! am I in the sediment region
                if(   slocal > ssalt) then
                   if(slocal < smin ) smin=slocal
                      end if

             end do
          end do
   else !! no salt
          slo%nrs(izs) = pro%nfk
       end if



       !! equally spaced points in the updated [smin...smax]
       ds = (smax-smin) / pro%nfk
       do ifk=1,pro%nfk
          slo%ssz(izs,ifk) = smin + (ifk-.5)*ds
          slo%srz(izs,ifk) = smin + (ifk-.5)*ds
       end do

         slo%si(:,:,:,:,izs)=pro%nfk+1
       do ifk=1,pro%nfk
          sl = smin+(ifk-1)*ds
          su = smin+(ifk  )*ds
          where( sxy(:,:,:,:)>=sl .and. sxy(:,:,:,:)<=su)
             slo%si(:,:,:,:,izs)=ifk
          end where
       end do
    end do izsloop


!       if(vsalt >0. .and. smin <= ssalt .and. pro%nfk>1) then
!          smin=0.
!          ins=0
!          do imy=1,pro%amy%n
!             do imx=1,pro%amx%n
!                if(real(S%bin(imx,imy,izs)) >  ssalt) then
!                   ins = ins+1
!                   smin = smin + real(S%bin(imx,imy,izs))
!                end if
!             end do
!          end do
!          smin = smin/ins
!
!          !! equally spaced points in the updated [smin...smax]
!          ds = (smax-smin) /(pro%nfk-1)
!          do ifk=1,pro%nfk-1
!             slo%ssz(izs,ifk) = smin + (ifk-.5)*ds
!             slo%srz(izs,ifk) = smin + (ifk-.5)*ds
!          end do
!          slo%si(:,:,:,:,izs) = 1+(sxy-smin)/ds
!          slo%si(:,:,:,:,izs) = max(min(1+nint((sxy-smin)/ds),pro%nfk-1),1)
!
!          !! average slowness in the salt region
!          ssalt_true=0.
!          ins=0
!          do imy=1,pro%amy%n
!             do imx=1,pro%amx%n
!                if(real(S%bin(imx,imy,izs)) <= ssalt) then
!                   ins = ins+1
!                   ssalt_true = ssalt_true + real(S%bin(imx,imy,izs))
!                end if
!             end do
!          end do
!          if(ins/=0) then
!             ssalt_true = ssalt_true/ins 
!          else
!             ssalt_true = ssalt
!          end if
!          
!          slo%ssz(izs,pro%nfk) = ssalt_true
!          slo%srz(izs,pro%nfk) = ssalt_true
!
!          where(sxy <= ssalt)
!             slo%si(:,:,:,:,izs)=pro%nfk
!          end where
!
!       else !! no salt
!          ds = (smax-smin) /pro%nfk
!          do ifk=1,pro%nfk
!             slo%ssz(izs,ifk) = smin + (ifk-.5)*ds
!             slo%srz(izs,ifk) = smin + (ifk-.5)*ds
!          end do
!          slo%si(:,:,:,:,izs) = 1+(sxy-smin)/ds
!          slo%si(:,:,:,:,izs) = max(min(1+nint((sxy-smin)/ds),pro%nfk),1)
!       end if
!
!    end do izsloop

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
    st=OK
  end function weislsN
  !----------------------------------------------------------------
end module wei_sls




