module sll_mod
  !! Lloyds Algorithm slowness operator
  !! for FK-FX downward continuation
  !! includes salt logic ??
  use type_mod

  use util_mod
  use S_mod
  use lloyd_mod
  !---------------
  implicit none

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

contains
  !----------------------------------------------------------------
  !! N reference slownesses
  integer function sll_init() result(st)
    integer :: er
    real::a
#ifdef DEBUG
    character(len=128), parameter :: unit='slsN_init'
    if(run%debug) call in(unit)
#endif
    allocate(bound(down%nref),ref_slow(down%nref))
    allocate(s_slice(size(S%s_s%bin,1)*size(S%s_s%bin,2)))
    allocate(sxy(down%amx%n,down%amy%n),stat=er)
    if(er/=0) call seperr('cant allocate sxy')

    a=4*(size(bound)+size(s_slice)+size(sxy))
    down%mem=down%mem+a
#ifdef MEMCHECK
    if(run%debug) write(run%io,*) 'sll.F90: allocate',a/1024./1024.,'Mbytes.  Total=',down%mem/1024./1024.
#endif
    call init_lloyd()
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
    st=OK
  end function sll_init


  !----------------------------------------------------------------  
  !! N reference slownesses
  integer function sll() result(st)
    integer :: izs,ifk
    real    :: ds, so, sz
    real    :: sl,su
    real    :: smin,smax
    integer :: imx, imy, i1,i2
    real    :: ssalt_true
    real    :: ssalt_mean,ssalt_sum
    real    :: slocal
    integer :: ins,ndo,ierr

#ifdef DEBUG
    character(len=128), parameter :: unit='slsN'
    if(run%debug) call in(unit)
#endif

    call sellect_slow_comp(.true.)

    do izs=1,down%nzs
       slo%nrs(izs) = down%nref
       do i2=1,size(sxy,2)
       do i1=1,size(sxy,1)
          sxy(i1,i2) =S%bin(i1,i2,izs)
       end do; end do
       s_slice=reshape(sxy,(/size(sxy)/))
       call lloyd_setup(izs)
    end do

    if (down%cwave) then
       call sellect_slow_comp(.false.)
       do izs=1,down%nzs
          slo%nrs(izs) = down%nref
          do i2=1,size(sxy,2)
          do i1=1,size(sxy,1)
             sxy(i1,i2) =S%bin(i1,i2,izs)
          end do; end do
          s_slice=reshape(sxy,(/size(sxy)/))
          call lloyd_setup(izs) 
       end do
    else 
       slo%nrs_r=>slo%nrs_s
       slo%sz_r =>slo%sz_s
       slo%si_r =>slo%si_s
       slo%bound_r=>slo%bound_s
    end if

#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
    st=OK
  end function sll
  !----------------------------------------------------------------
  subroutine lloyd_setup(izs)
    integer :: i1,i2,izs

    slo%nrs(izs)=lloyd_go(s_slice,ref_slow,down%nref,bound)
    slo%sz(izs,1:size(ref_slow)) =ref_slow
    slo%bound(izs,1:size(bound))= bound

    do  i2=1,size(sxy,2)
    do  i1=1,size(sxy,1)
       slo%si(i1,i2,izs)=1
       do while(sxy(i1,i2) > bound(slo%si(i1,i2,izs)))
          slo%si(i1,i2,izs)=slo%si(i1,i2,izs)+1
       end do
    end do
    end do
  end subroutine lloyd_setup
!----------------------------------------------------------------
end module sll_mod
