module S_mod

  use down_mod
  use xmap_mod
  !---------------
  implicit none

  real,    dimension(:), pointer, private ::  sloline
  type(slowmap), public :: slo
  type(slow),public :: S

contains
  !----------------------------------------------------------------
  !INITIALIZE THE SLOWNESS FUNCTION
  ! GIVEN UP TO TWO TAGS DESCRIBING THE SLOWNESS
  logical function slow_init(X,tag_S,tag_R)
    type(slow),       intent(inout) :: X
    character(len=128), intent(in)    :: tag_S,tag_R
    integer :: er
#ifdef DEBUG
    character(len=128), parameter :: unit='slow_init' 
    if(run%debug) call in(unit)
#endif

    slow_init=F 

    if(down%cwave) then
       X%s_s%tag = tag_S
       X%s_r%tag = tag_R
       if(.not. slowness_history_init(X%s_s)) then
          write(run%io,*) "trouble init slowness history"
          return
       end if
       if(.not. slowness_history_init(X%s_r)) then
          write(run%io,*) "trouble init slowness history"
          return
       end if
    else
       X%s_s%tag = tag_S
       if(.not. slowness_history_init(X%s_s)) then
          write(run%io,*) "trouble init slowness history"
          return
       end if
       X%s_r%tag=X%s_s%tag
       X%s_r%from=X%s_s%from
       X%s_r%bites=X%s_s%bites
       X%s_r%esize=X%s_s%esize
       X%s_r%amx=X%s_s%amx
       X%s_r%amy=X%s_s%amy
       X%s_r%az=X%s_s%az
       X%s_r%bin=>X%s_s%bin
    end if

    call slow_unitcheck(S%s_s)

    down%S=S
    slow_init=T 

#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end function slow_init
!----------------------------------------------------------------
!INITIALIZE SLOWNESS FROM A FILE
!
  logical function slowness_history_init(X)
    type(slow_comp), intent(inout) :: X
    integer  :: stat
    character(len=128) :: temp
#ifdef DEBUG
    character(len=128), parameter :: unit='slowness_history_init' 
    if(run%debug) call in(unit)
#endif

    slowness_history_init=F
!    call auxinout(X%tag)
    X%from = "history"
    call hinaxis(X%tag,X%amx,1); X%amx%l='ix'
    call hinaxis(X%tag,X%amy,2); X%amy%l='iy'
    call hinaxis(X%tag,X%az ,3); X%az%l ='z '

    call from_aux(trim(X%tag),"esize",X%esize,4)
    if(X%esize/=4)call seperr("Velocity esize shoule be real")
    if(.not. slowness_setup(X)) then
       write(run%io,*) "trouble with slowness setup"
       return
    end if
    if(down%watercolumn) then
       if(X%az%o+(X%az%n-1)*X%az%d.lt.down%awtr%o) &
         call seperr("Water column datum depth is greater than velocity model max depth.  Idiot!")
    end if
    slowness_history_init=T

#ifdef DEBUG
    if(run%debug) call slowness_report(X)
    if(run%debug) call out(unit)
#endif
  end function slowness_history_init
!----------------------------------------------------------------
! . . A few sanity checks
subroutine slow_unitcheck(X)
  type(slow_comp)::X
  real:: sloc,loc

  !! . . Sanity checks
  if(down%aix%n.gt.1) then
     if (X%amx%o .gt. down%aix%o )&
       call seperr("Slowness model origin greater than data coverage: x-axis")
     loc=(down%aix%n-1)*down%aix%d+down%aix%o
     sloc=(X%amx%n-1)*X%amx%d+X%amx%o
     if(sloc.lt.loc) call seperr("Slowness model end less than data extent: x-axis")
     if (((down%amx%d/X%amx%d .gt. 100).or.(down%amx%d/X%amx%d .lt. 0.01 )))&
       call seperr("Check units are of same magnitude between slowness data files: x-axis")
     if(X%amx%n.lt.2) call seperr("Slowness needs x-axis.")
     if(X%amx%d.ne.down%amx%d) call seperr("I demand you sample the velocity model appropriately! x-axis")
  end if
  if(down%aiy%n.gt.1) then
     if (X%amy%o .gt.down%aiy%o )&
     call seperr("Slowness model origin greater than data coverage: y-axis")
     loc=(down%aiy%n-1)*down%aiy%d+down%aiy%o
     sloc=(X%amy%n-1)*X%amy%d+X%amy%o
     if(sloc.lt.loc) call seperr("Slowness model end less than data extent: y-axis")
     if (((down%amy%d/X%amy%d .gt. 100).or.(down%amy%d/X%amy%d .lt. 0.01 )).and.(X%amy%n /=1)) &
      call seperr("check units are of same magnitude between slowness data files: y-axis")
     if(X%amy%n.lt.2) call seperr("Slowness needs y-axis.")
     if(X%amy%d.ne.down%amy%d) call seperr("I demand you sample the velocity model appropriately! y-axis")
  end if
end subroutine slow_unitcheck

  !----------------------------------------------------------------
  ! SETUP SLOWNESS
  ! 
logical function slowness_setup(X)
  type(slow_comp), intent(inout) :: X
  integer :: er,i1,i2,i3
  real::a
  character(len=128) :: temp
#ifdef DEBUG
  character(len=128), parameter :: unit='slowness_setup' 
  if(run%debug) call in(unit)
#endif

  slowness_setup=F
  X%bites = X%esize * product( (/down%amx%n,down%amy%n,down%nzs/) )

  !!--------------------
  allocate(X%bin(down%amx%n,down%amy%n,down%nzs),stat=er)
  if(er/=0) call seperr('cant allocate S.bin')
  do i3=1,down%nzs
     do i2=1,down%amy%n
        do i1=1,down%amx%n
           X%bin(i1,i2,i3)=0.0
        end do
     end do
  end do
  a=4*size(X%bin)
  down%mem=down%mem+a
#ifdef MEMCHECK
  if(run%debug) write(run%io,*) 'S.F90 s%bin: allocate',a/1024./1024.,'Mbytes.  Total=',down%mem/1024./1024.
#endif
  slowness_setup=T
  if(run%verb.gt.2) call slowness_report(X)
#ifdef DEBUG
  if(run%debug) call out(unit)
#endif
end function slowness_setup
!----------------------------------------------------------------
  !WRITE INFO ABOUT SLOWNESS
  subroutine slowness_report(X)
    type(slow_comp), intent(in) :: X
    real::bb
    bb=real(down%az%n)/real(down%nzs)
    write(run%io,*)'        ************************************************'
    write(run%io,*)'        *    ',trim(X%tag),' from ',trim(X%from)
    write(run%io,*)'        *              [n]      [o]      [d]'
    write(run%io,'(a,i4,f10.2,f10.3)       ') "         *         Mx: ",X%amx%n, X%amx%o, X%amx%d
    write(run%io,'(a,i4,f10.2,f10.3)       ') "         *         My: ",X%amy%n, X%amy%o, X%amy%d
    write(run%io,'(a,i4,f10.2,f10.3,a,f5.2,a)') "         *          Z: ",X%az%n,  X%az%o,  X%az%d," ->",bb," blocks"
    if(down%vel) then
       write(run%io,'(a,i1,a)')               "         *       esize=",X%esize,"                 VELOCITY"
    else 
       write(run%io,'(a,i1,a)')               "         *       esize=",X%esize,"                 SLOWNESS"
    end if
    if(associated(X%bin)) then
       write(run%io,'(a,f9.2,a)        ') "         *       block:",4*size(X%bin)/1024./1024.," Mbytes"
    end if

  end subroutine slowness_report
  !--------------------------------------------------------------
  ! READ SLOWNESS
  subroutine read_slowness(k)
    integer, intent(in) :: k
#ifdef DEBUG
    character(len=128), parameter :: unit='slowness_read' 
    if(run%debug) call in(unit)
#endif

    call slowness_location(S%s_s)
    call slowness_read_comp(S%s_s,k)
    if(down%cwave) then
    call slowness_location(S%s_r)
       call slowness_read_comp(S%s_r,k)
    end if
    if(run%verb.gt.2) then
       if (down%vel) then
          write(run%io,*) "        *       Read velocity model"
       else
          write(run%io,*) "        *       Read slowness model"
       end if
    end if
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end subroutine read_slowness
  !----------------------------------------------------------------
  !READ A SLOWNESS COMPONENT
  subroutine slowness_read_comp(X,k)
    type(slow_comp), intent(inout) :: X
    integer, intent(in) :: k
    integer ::i1,i2,i3,i4,i5,ierr,n_dim_cube,n
    integer,dimension(3)::n_cube,n_wind,f_wind,j_wind
#ifdef DEBUG
    character(len=128), parameter :: unit='slowness_read_comp' 
    if(run%debug) call in(unit)
#endif

    n=down%zblock
    n_dim_cube=3
    n_cube=(/   X%amx%n,   X%amy%n,X%az% n /)
    j_wind=(/1         ,1         ,1       /)
    n_wind=(/down%amx%n,down%amy%n,down%nzs/)
    f_wind=(/X%indx -1,X%indy -1,(k-1)*n  /)
    if(run%lock) call sep_open_lock(X%tag)
    ierr=sreed_window(X%tag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,X%esize,X%bin(1,1,1))
    if(run%lock) call sep_close_lock(X%tag)
    if(X%bin(1,1,1).eq.0.) call seperr("Did you Pad your Vmodel without extend=1?")
    if(down%vel)then
!$OMP PARALLEL DO   
       do i3=1,size(X%bin,3)
       do i2=1,size(X%bin,2)
       do i1=1,size(X%bin,1)          
          X%bin(i1,i2,i3)=1./X%bin(i1,i2,i3)
       end do;end do; end do
!$OMP END PARALLEL DO
    end if


#ifdef DEBUG 
    if(run%debug) call out(unit)
#endif
  end subroutine slowness_read_comp

!--------------------------------------------------------------
  subroutine slowness_location(X)
    type(slow_comp), intent(inout)   ::X
    integer:: ierr,i,j,k,l,m,n,p,o
#ifdef DEBUG
    character(len=128), parameter :: unit='slowness_location'
    if(run%debug) call in(unit)
#endif
    X%indx= (1.5+ (down%amx%o-X%amx%o)/down%amx%d)
    if (X%indx .lt. 1) call seperr("Your velocity model is for the wrong place (x).")

    X%indy=1
    if (X%amy%n.gt.1) then
       X%indy= (1.5+ (down%amy%o-X%amy%o)/down%amy%d)
       if (X%indy .lt. 1) call seperr("Your velocity model is for the wrong place (y).")
    end if

    X%finx=X%indx + down%amx%n -1
    X%finy=X%indy + down%amy%n -1
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end subroutine slowness_location
!----------------------------------------------------------------
! WRITE OUT SLOWNESS
  subroutine slowness_write(X,k)
    type(slow_comp), intent(inout) :: X
    integer :: k
    integer :: stat
    integer :: iz
#ifdef DEBUG
    character(len=128), parameter :: unit='slowness_write'
    if(run%debug) call in(unit)
#endif

    stat=sseek_block(X%tag,k-1,X%bites,0)
    stat= srite(X%tag, X%bin, X%bites)

#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end subroutine slowness_write
!----------------------------------------------------------------
!ALLOCATE SPACE FOR SLOWNESS
  subroutine slow_allocate()
    integer :: er
    real::a
#ifdef DEBUG
    character(len=128), parameter :: unit='slow_allocate' 
    if(run%debug) call in(unit)
#endif
    !! FKop slowness arrays
    allocate(slo%sz_s(down%nzs,down%nref+1),stat=er)
    if(er/=0) call seperr('cant allocate slo%sz_s')
    allocate(slo%nrs_s(down%nzs),stat=er)
    if(er/=0) call seperr('cant allocate slo%nrs_s')
    allocate(slo%bound_s(down%nzs,down%nref),stat=er)
    if(er/=0) call seperr('cant allocate slo%bound_s')
    !! FXop interpolation array 
    allocate(slo%si_s(down%amx%n,down%amy%n,down%nzs),stat=er)   
    if(er/=0) call seperr('cant allocate slo%si_s')

    if(down%cwave)then
       allocate(slo%sz_r(down%nzs,down%nref+1),stat=er)
       if(er/=0) call seperr('cant allocate slo%sz_r')
       allocate(slo%nrs_r(down%nzs),stat=er)
       if(er/=0) call seperr('cant allocate slo%nrs_r')
       allocate(slo%bound_r(down%nzs,down%nref),stat=er)
       if(er/=0) call seperr('cant allocate slo%bound_r')  
       allocate(slo%si_r(down%amx%n,down%amy%n,down%nzs),stat=er)   
       if(er/=0) call seperr('cant allocate slo%si_r')
    else
       slo%sz_r=>slo%sz_s
       slo%nrs_r=>slo%nrs_r
       slo%bound_r=>slo%bound_s
       slo%si_r=>slo%si_s
    end if

    slo%si=>slo%si_s
    slo%bound=>slo%bound_s
    slo%nrs=>slo%nrs_s
    slo%sz=>slo%sz_s

    a=4*(size(slo%sz_s)+size(slo%nrs_s)+size(slo%bound_s)+size(slo%si_s))
    if(down%cwave) a=a*2.
    down%mem=down%mem+a
#ifdef MEMCHECK
    if(run%debug) write(run%io,*) 'S.F90 working arrays: allocate',a/1024./1024.,'Mbytes.  Total=',down%mem/1024./1024.
#endif

#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end subroutine slow_allocate

!----------------------------------------------------------------
  subroutine sellect_slow_comp(source)
    logical :: source  
    slo%source=source
    if(source) then
       slo%si=>slo%si_s
       slo%nrs=>slo%nrs_s
       slo%bound=>slo%bound_s
       slo%sz=>slo%sz_s
       S%bin=>S%s_s%bin
    else
       slo%si=>slo%si_r
       slo%nrs=>slo%nrs_r
       slo%bound=>slo%bound_r
       slo%sz=>slo%sz_r
       S%bin=>S%s_r%bin
    end if
  end subroutine sellect_slow_comp
!!----------------------------------------------------------------
!!FIND THE MEAN VALUE OF THE SLOWNESS
!  real function slowness_mean(X,k) result(smean)
!    type(slow), intent(in) :: X
!    integer :: k
!
!    if(down%cwave) then
!       smean = (sum(real(X%s_s%bin(:,:,k)))+&
!       sum(real(X%s_r%bin(:,:,k))))/(2* X%s_s%amx%n * X%s_s%amy%n)
!    else
!       smean = sum(real(X%s_s%bin(:,:,k)))/(X%s_s%amx%n * X%s_s%amy%n)
!    end if
!  end function slowness_mean
!!----------------------------------------------------------------
!!FIND THE MAX VALUE OF THE SLOWNESS
!  real function slowness_max(X,k) result(smax)
!    type(slow), intent(in) :: X
!    integer :: k
!
!    if(down%cwave) then
!       smax =max(maxval(real(X%s_s%bin(:,:,k))),maxval(real(X%s_r%bin(:,:,k))))
!    else
!       smax = maxval(real(X%s_s%bin(:,:,k)))
!    end if
!  end function slowness_max
!!----------------------------------------------------------------
!!FIND THE MIN VALUE OF THE SLOWNESS
!  real function slowness_min(X,k) result(smin)
!    type(slow), intent(in) :: X
!    integer :: k
!
!    if(down%cwave) then
!       smin =min(minval(real(X%s_s%bin(:,:,k))),minval(real(X%s_r%bin(:,:,k))))
!    else
!       smin = minval(real(X%s_s%bin(:,:,k)))
!    end if
!  end function slowness_min
!----------------------------------------------------------------
!  subroutine slow_interp(adj,X,line,xc,yc)
!    type(slow_comp), intent(inout) :: X
!    logical,                   intent(in) :: adj
!    complex, dimension(:,:,:), pointer    :: cube
!    real,    dimension(:),     pointer    :: line
!    real,                      intent(in) :: xc,yc
!    real    :: dx,dy
!    integer :: ix,iy, jx,jy , i3
!    real    :: fx,fy
!    real    :: w00,w01,w10,w11, sumw
!#ifdef DEBUG
!    character(len=128), parameter :: unit='slow_interp'
!    if(run%debug) call in(unit)
!#endif
!
!    cube=>X%bin
!
!    dx = xc
!    dy = yc
!
!    !! dont allow coordinates outside the slowness grid
!    dx = min(dx,X%amx%o + (X%amx%n-1)* X%amx%d)
!    dx = max(dx,X%amx%o)
!
!    dy = min(dy,X%amy%o + (X%amy%n-1)* X%amy%d)
!    dy = max(dy,X%amy%o)
!
!    !! lower index
!    ix = ( 1.5+ (dx-X%amx%o) / X%amx%d )
!    iy = ( 1.5+ (dy-X%amy%o) / X%amy%d )
!
!    !! keep indices in bounds
!    ix = min(ix,X%amx%n)
!    iy = min(iy,X%amy%n)
!    ix = max(ix,        1  ); jx=min(ix+1,X%amx%n)
!    iy = max(iy,        1  ); jy=min(iy+1,X%amy%n)
!
!    !! linear weights
!    fx = (dx - (X%amx%o + (ix-1)*X%amx%d) ) / X%amx%d
!    fy = (dy - (X%amy%o + (iy-1)*X%amy%d) ) / X%amy%d
!    fx = max(0., fx)
!    fy = max(0., fy)
!
!    !! 00   10       --> x
!    !!    .         |
!    !! 01   11      v
!    !!               y
!    w00=(1.-fx)*(1.-fy)
!    w10=(   fx)*(1.-fy)
!    w01=(1.-fx)*(   fy)
!    w11=(   fy)*(   fy)
!
!    !! normalize the sum of the weights
!    sumw = w00+w10+w01+w11
!    w00 = w00 / sumw
!    w10 = w10 / sumw
!    w01 = w01 / sumw
!    w11 = w11 / sumw
!
!    if(adj) then !! push
!
!       do i3=1,size(cube,3)
!          cube(ix,iy,i3) = cube(ix,iy,i3) + w00 * line(i3)
!          cube(jx,iy,i3) = cube(jx,iy,i3) + w10 * line(i3)
!          cube(ix,jy,i3) = cube(ix,jy,i3) + w01 * line(i3)
!          cube(jx,jy,i3) = cube(jx,jy,i3) + w11 * line(i3)
!       end do
!    else         !! pull
!       line(i3)    = &
!       w00*real(cube(ix,iy,i3)) + &
!       w10*real(cube(jx,iy,i3)) + &
!       w01*real(cube(ix,jy,i3)) + &
!       w11*real(cube(jx,jy,i3))
!    end if
!
!#ifdef DEBUG
!    if(run%debug) call out(unit)
!#endif
!  end subroutine slow_interp

end module S_mod
