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

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

  real, dimension(:,:,:,:), pointer, private :: sxy
  real, dimension(:), pointer, private :: bound,ref_slow,s_slice

contains
!!$
!!$=head1 NAME
!!$
!!$ weisll_init - Initialize multiple  slowness  lloyd operator
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ierr=weisll_init()
!!$
!!$=head1 DESCRIPTION
!!$
!!$ Initialize multiple slowness operator lloyd
!!$
!!$=head1 SEE ALSO
!!$
!!$L<weisloN_init> L<weislo1> L<weislsN>
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$=cut
!!$

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

    allocate(bound(pro%nfk),ref_slow(pro%nfk))
    allocate(s_slice(size(S%s_s%bin,1)*size(S%s_s%bin,2)))
    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 init_lloyd()


#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
    st=OK
  end function weisll_init
!!$
!!$=head1 NAME
!!$
!!$ weisll - Multiple  slowness operator lloyd vel/find
!!$
!!$=head1 SYNOPSIS
!!$
!!$ ierr=weisll()
!!$
!!$=head1 DESCRIPTION
!!$
!!$  multiple slowness operator lloyd
!!$
!!$=head1 SEE ALSO
!!$
!!$L<weisloN_init> L<weislo1> L<weislsN_init>
!!$
!!$=head1 LIBRARY
!!$
!!$ sepweif90
!!$
!!$=cut
!!$
  !----------------------------------------------------------------  
  !! N reference slownesses
  integer function weisll() result(st)
    integer :: izs,ifk,i4,i3,i2,i1
    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))

       s_slice=reshape(S%s_s%bin(:,:,izs),(/size(S%s_s%bin(:,:,izs))/))
       slo%nrs(izs)=lloyd_go(s_slice,ref_slow,pro%nfk,bound,pro%mpi_ithread==0)
       slo%ssz(izs,:) =ref_slow
       slo%srz(izs,:) = ref_slow
       do  i4=1,size(sxy,4)
       do  i3=1,size(sxy,3)
       do  i2=1,size(sxy,2)
       do  i1=1,size(sxy,1)
         slo%si(i1,i2,i3,i4,izs)=1
         do while(sxy(i1,i2,i3,i4) > bound(slo%si(i1,i2,i3,i4,izs)))
           slo%si(i1,i2,i3,i4,izs)=slo%si(i1,i2,i3,i4,izs)+1
         end do
       end do;end do;end do;end do
    end do izsloop


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




