module wei_obj_S

  use wei_util
  use wei_process
  use wei_type
  use wei_work
  use wei_xmap
  !---------------
  implicit none

  real,    dimension(:), pointer, private ::  sloline
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  type(weslow),    public :: S   !Slowness
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  type(weslowmap), public :: slo
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!

contains
  !----------------------------------------------------------------
  subroutine weibslow_init(X,tag_S,tag_G)
    type(weslow),       intent(inout) :: X
    character(len=128), intent(in)    :: tag_S,tag_G
#ifdef DEBUG
    character(len=128), parameter :: unit='weibslow_init' 
    if(pro%debug) call in(unit)
#endif
    integer :: er

    if(pro%conv_wave) then
       X%s_s=X%s1
       X%s_g=X%s2
       X%s_s%tag = tag_S
       X%s_g%tag = tag_g
       call slowness_history_init(X%s_s)
       call slowness_history_init(X%s_g)
    else
       X%s_s=X%s1
       X%s_g=X%s1
       X%s_s%tag = tag_S
       call slowness_history_init(X%s_s)
    end if

    !! prestack binary array (PBA)
    allocate(X%pba(pro%amx%n,pro%amy%n, pro%ahx%n,pro%ahy%n, pro%nzs),stat=er)
    if(er/=0) call seperr('cant allocate X%pba')


    call find_max_slow()
   call send_slow(S)
    pro%slow=S

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
    pro%slow=X
  end subroutine weibslow_init


 subroutine find_max_slow()
  integer :: i,ndone,block,j,ierr
  character(len=128) :: tag_1,tag_2
  real,allocatable :: junk(:)
  
  ndone=0;
  allocate(junk(S%s_s%az%n))
  allocate(pro%max_slow(S%s_s%az%n))
#ifdef SEP_MPI
 if((pro%mpi_nthread>1 .and. pro%mpi_ithread==0) .or. pro%mpi_nthread==1 ) then
  tag_1=S%s_s%tag
  tag_2=S%s_g%tag
  S%s_s%tag=s%s_s%tag_0
  S%s_g%tag=s%s_g%tag_0
#endif
  i=0
  do while(ndone< S%s_s%az%n)
    i=i+1
		block=min(s%s_s%az%n-ndone,size(S%s_s%bin,3))
   	call slowness_read(S,i,block)
    if(pro%conv_wave) then
    do j=1,block
      pro%max_slow(j+ndone)=max(maxval(real(S%s_s%bin(:,:,j))),maxval(real(S%s_g%bin(:,:,j))))
    end do
   else
    do j=1,block
      pro%max_slow(j+ndone)=maxval(real(S%s_s%bin(:,:,j)))
    end do
   end if
    ndone=ndone+block
  end do
#ifdef SEP_MPI
  S%s_s%tag=tag_1
  S%s_g%tag=tag_2
  junk=pro%max_slow
 end if 
 if(pro%mpi_nthread>1)  then
  call MPI_BCAST(junk, size(junk),MPI_REAL,0,MPI_COMM_WORLD,ierr)
  pro%max_slow=junk
 end if
 !AAAQ 

#endif
 deallocate(junk)
 end subroutine find_max_slow 


  !----------------------------------------------------------------
  subroutine slowness_history_init(X)
    type(weslow_comp), intent(inout) :: X
#ifdef DEBUG
    character(len=128), parameter :: unit='slowness_history_init' 
    if(pro%debug) call in(unit)
#endif

    call auxinout(X%tag)
    X%from = "history"

    call hinaxis(X%tag,X%amx,1); X%amx%l='mx'
    if(pro%amx%n ==0) pro%amx = X%amx

    call hinaxis(X%tag,X%amy,2); X%amy%l='my'
    if(pro%amy%n ==0) pro%amy = X%amy

    call hinaxis(X%tag,X%az ,3); X%az%l ='z_'
    if(pro%az%n ==0) then
       pro%az = X%az
    else
       if(.not. coax( X%az, pro%az)) then
          call slowness_report(X)
          call seperr('S.az /= pro.az')
       end if
    end if

    call from_param("velocity",X%velocity,.false.)
    call from_aux(trim(X%tag),"esize",X%esize,4)

    call slowness_setup(X)
    
#ifdef SEP_MPI
    if(pro%mpi_nthread>1) then
       X%tag_0=X%tag;
       write(X%tag_local,"(a,a)") trim(X%tag_0),pro%tag_num
       x%tag=X%tag_local
       call auxscr(X%tag)
       call houaxis(X%tag,X%amx,1)
       call houaxis(X%tag,X%amy,2)
       call houaxis(X%tag,X%az, 3)
       call to_history("esize",X%esize,X%tag)
!       if(0/=mpi_sep_tag_bcast(0,trim(X%tag_0),pro%mpi_ithread,trim(X%tag),&
!       X%bites/4,pro%nzb)) call seperr("trouble distributing slowness")
    end if
#endif

    write(X%tag_var,"(a,a,a)") trim(X%tag_0),"_VS",pro%tag_num
    call auxscr(X%tag_var)
    call srite(X%tag_var,0.,4)

#ifdef DEBUG
    if(pro%debug) call slowness_report(X)
    if(pro%debug) call out(unit)
#endif
  end subroutine slowness_history_init
  !----------------------------------------------------------------
  subroutine slowness_param_init(X)
    type(weslow_comp), intent(inout) :: X
#ifdef DEBUG
    character(len=128), parameter :: unit='slowness_param_init' 
    if(pro%debug) call in(unit)
#endif
    
    call auxscr(X%tag)
    X%from = "parameter"
    
    X%amx = pro%amx; 
    X%amy = pro%amy;
    X%az  = pro%az;
    X%esize=8
    X%velocity=.false.
    
    if(pro%mpi_nthread>1) then
       X%tag_0=X%tag;
       write(X%tag,"(a,a)") trim(X%tag_0),pro%tag_num
       call auxscr(X%tag)
       call houaxis(X%tag_0,X%amx,1)
       call houaxis(X%tag_0,X%amy,2)
       call houaxis(X%tag_0,X%az, 3)
       call to_history("esize",8,X%tag_0)
    end if
    
    call houaxis(X%tag,X%amx,1)
    call houaxis(X%tag,X%amy,2)
    call houaxis(X%tag,X%az, 3)
    call to_history("esize",8,X%tag)
    
    call slowness_setup(X)
    call slowness_assign(X,(0.,0.))
    
#ifdef DEBUG
    if(pro%debug) call slowness_report(X)
    if(pro%debug) call out(unit)
#endif
  end subroutine slowness_param_init
  !----------------------------------------------------------------
  subroutine slowness_scratch_init(X)
    type(weslow_comp), intent(inout) :: X
#ifdef DEBUG
    character(len=128), parameter :: unit='slowness_scratch_init' 
    if(pro%debug) call in(unit)
#endif
    
    call auxtmp(X%tag)
    X%from = "parameter"
    
    X%amx = pro%amx; call houaxis(X%tag,X%amx,1)
    X%amy = pro%amy; call houaxis(X%tag,X%amy,2)
    X%az  = pro%az;  call houaxis(X%tag,X%az, 3)
    
    call to_history("esize",8,X%tag)
    
    call slowness_setup(X)
    call slowness_assign(X,(0.,0.))
    
#ifdef DEBUG
    if(pro%debug) call slowness_report(X)
    if(pro%debug) call out(unit)
#endif
  end subroutine slowness_scratch_init
  !----------------------------------------------------------------
  subroutine slowness_setup(X)
    type(weslow_comp), intent(inout) :: X
    integer :: er,i1,i2,i3
#ifdef DEBUG
    character(len=128), parameter :: unit='slowness_setup' 
    if(pro%debug) call in(unit)
#endif
    
    if( (1.0*X%az%n)/pro%nzs /= X%az%n/pro%nzs ) then
       call slowness_report(X)
       call seperr('S.az.n is not a multiple of nzs')
    end if
    
    pro%nzb = pro%az%n / pro%nzs
    !!--------------------
    allocate(X%bin(X%amx%n,X%amy%n,pro%nzs),stat=er)
    if(er/=0) call seperr('cant allocate S.bin')

    
    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)=(0.0,0.0)
    end do
    end do
    end do
    !!--------------------
    X%bites = X%esize * product( (/X%amx%n,X%amy%n,pro%nzs/) )
    
    
    !! slowness interpolation map
    allocate( sloline(pro%nzs),stat=er)
    if(er/=0) call seperr('cant allocate  sloline')

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine slowness_setup
  !----------------------------------------------------------------
  subroutine slowness_assign(X,v)
    type(weslow_comp), intent(inout) :: X
    complex, intent(in) :: v
    integer :: izb
#ifdef DEBUG
    character(len=128), parameter :: unit='slowness_assign' 
    if(pro%debug) call in(unit)
#endif

    X%bin=v
    do izb=1,pro%nzb
       call slowness_write(X,izb)
    end do

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine slowness_assign
  !----------------------------------------------------------------
  subroutine slowness_report(X)
    type(weslow_comp), intent(in) :: X

    write(0,*)'        ************************************************'
    write(0,*)'        *       ',trim(X%tag),' from ',trim(X%from)
    write(0,'(a,i3,f10.2,f10.3)       ') "         *    Mx: ",X%amx%n, X%amx%o, X%amx%d
    write(0,'(a,i3,f10.2,f10.3)       ') "         *    My: ",X%amy%n, X%amy%o, X%amy%d
    write(0,'(a,i3,f10.2,f10.3,a,i3,a)') "         *     Z: ",X%az%n,  X%az%o,  X%az%d," ->",pro%nzb," blocks"
    if(X%velocity) then
       write(0,'(a,i1,a)')               "         *    esize=",X%esize," velocity"
    else 
       write(0,'(a,i1,a)')               "         *    esize=",X%esize," slowness"
    end if
    if(associated(X%bin)) then
       write(0,'(a,f6.2,a,f6.2,a)        ') "         * block: ",X%bites/1024.," / ",8*size(X%bin)/1024.," Kbytes"
       write(0,*)'        ************************************************'
    end if
    
  end subroutine slowness_report
  !--------------------------------------------------------------
  subroutine slowness_read(X,k,nzs)
    type(weslow), intent(inout) :: X
    integer, intent(in) :: k
    integer, optional :: nzs
    integer :: stat,nz,i
    real,allocatable  :: temp_buf(:,:,:)
#ifdef DEBUG
    character(len=128), parameter :: unit='slowness_read' 
    if(pro%debug) call in(unit)
#endif

   
     call slowness_read_comp(X%s_s,k,nzs)
    if(pro%conv_wave) then
     call slowness_read_comp(X%s_g,k,nzs)
    end if
    
#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine slowness_read
  subroutine slowness_read_comp(X,k,nzs)
    type(weslow_comp), intent(inout) :: X
    integer, intent(in) :: k
    integer, optional :: nzs
    integer :: stat,nz,i
    real,allocatable  :: temp_buf(:,:,:)
#ifdef DEBUG
    character(len=128), parameter :: unit='slowness_read_comp' 
    if(pro%debug) call in(unit)
#endif
    
    nz=size(X%bin,3); if(present(nzs)) nz=nzs

    stat=sseek_block(X%tag,k-1,X%bites,0)
    if(X%esize==4) then
       allocate(temp_buf(size(X%bin,1),size(X%bin,2),size(X%bin,3)),stat=stat)
       if(stat/=0) call erexit("trouble allocating real temp slow buf")
       call sreed(X%tag, temp_buf, X%bites/size(X%bin,3)*nz)
       if(X%velocity) temp_buf=1./temp_buf
       X%bin=cmplx(temp_buf,0.)
       deallocate(temp_buf)
    else
       call sreed(X%tag, X%bin, X%bites)
    end if

    if(nz <  size(X%bin,3)) then
      do i=nz,size(X%bin,3)
        X%bin(:,:,i)=X%bin(:,:,nz)
     end do
    end if
    
#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine slowness_read_comp

  !--------------------------------------------------------------
  subroutine distribute_slow(X)
    type(weslow), intent(inout) :: X
  end subroutine distribute_slow

  subroutine send_slow(X)
    type(weslow), intent(inout) :: X

#ifdef SEP_MPI
    if(pro%mpi_nthread>1 ) then
        call sseek(trim(X%s_s%tag),0,0)
       if(pro%mpi_ithread==0) then
          call sseek(trim(X%s_s%tag_0),0,0)
        end if
       if(0/=mpi_sep_tag_bcast(0,X%s_s%tag_0,pro%mpi_ithread,X%s_s%tag,X%s_s%bites/4,&
       pro%nzb)) call seperr("trouble distributing slow \n")
    if(pro%conv_wave) then
        call sseek(trim(X%s_g%tag),0,0)
       if(pro%mpi_ithread==0) call sseek(trim(X%s_g%tag_0),0,0)
       if(0/=mpi_sep_tag_bcast(0,X%s_g%tag_0,pro%mpi_ithread,X%s_g%tag,X%s_g%bites/4,&
       pro%nzb)) call seperr("trouble distributing slow \n")
    end if
   end if
#endif
  end subroutine send_slow
  !--------------------------------------------------------------
  subroutine combine_slow(X)
    type(weslow), intent(inout) :: X
#ifdef SEP_MPI
    if(pro%mpi_nthread>1) call aux_unlink(X%s_s%tag_local)
    call aux_unlink(X%s_s%tag_var)
    if(pro%conv_wave) then
      if(pro%mpi_nthread>1) call aux_unlink(X%s_s%tag_local)
      call aux_unlink(X%s_s%tag_var)
    end if
#endif   
  end subroutine combine_slow
  !----------------------------------------------------------------
  !----------------------------------------------------------------
  subroutine slowness_write(X,k)
    type(weslow_comp), intent(inout) :: X
    integer :: k
    integer :: stat
    integer :: iz
#ifdef DEBUG
    character(len=128), parameter :: unit='slowness_write'
    if(pro%debug) call in(unit)
#endif

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

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine slowness_write
  !----------------------------------------------------------------



  subroutine weislow_allocate()
    integer :: er
#ifdef DEBUG
    character(len=128), parameter :: unit='weislow_allocate' 
    if(pro%debug) call in(unit)
#endif

    !! FKop slowness arrays
    allocate(slo%ssz(pro%nzs,pro%nfk+1),stat=er)
    if(er/=0) call seperr('cant allocate slo%ssz')
    allocate(slo%srz(pro%nzs,pro%nfk+1),stat=er)
    if(er/=0) call seperr('cant allocate slo%srz')
    allocate(slo%nrs(pro%nzs),stat=er)
    if(er/=0) call seperr('cant allocate slo%nrs')


    !! FXop interpolation array
    allocate(slo%si(pro%amx%n,pro%amy%n, pro%ahx%n,pro%ahy%n, pro%nzs),stat=er)   
    if(er/=0) call seperr('cant allocate slo%si')

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine weislow_allocate
  !----------------------------------------------------------------
  real function slowness_mean(X,k) result(smean)
    type(weslow), intent(in) :: X
    integer :: k

   if(pro%conv_wave) then
      smean = (sum(real(X%s_s%bin(:,:,k)))+sum(real(X%s_g%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
  !----------------------------------------------------------------
  real function slowness_max(X,k) result(smax)
    type(weslow), intent(in) :: X
    integer :: k

   if(pro%conv_wave) then
 smax =max(maxval(real(X%s_s%bin(:,:,k))),maxval(real(X%s_g%bin(:,:,k))))
else
    smax = maxval(real(X%s_s%bin(:,:,k)))
 end if
  end function slowness_max
  !----------------------------------------------------------------
  real function slowness_min(X,k) result(smin)
    type(weslow), intent(in) :: X
    integer :: k

   if(pro%conv_wave) then
    smin =min(minval(real(X%s_s%bin(:,:,k))),minval(real(X%s_g%bin(:,:,k))))
  else
    smin = minval(real(X%s_s%bin(:,:,k)))
  end if
  end function slowness_min
  !----------------------------------------------------------------
  !! push slowness from bin to ss and sr
  subroutine weislow_push(X)
    type(weslow), intent(inout) :: X
    integer :: imx,imy,ihx,ihy
    real    :: xc,yc
#ifdef DEBUG
    character(len=128), parameter :: unit='weislow_push' 
    if(pro%debug) call in(unit)
#endif

    X%pba=(0.,0.)

    do ihy=1,pro%ahy%n
       do ihx=1,pro%ahx%n
          do imy=1,pro%amy%n
             do imx=1,pro%amx%n

                !! X%ss = f(X%bin,x,y)
                xc=xmap%sx(imx,ihx)
                yc=xmap%sy(imy,ihy)
                call weislow_interp(.false.,X%s_s,sloline,xc,yc)

                X%pba(imx,imy,ihx,ihy,:) = &
                X%pba(imx,imy,ihx,ihy,:) + &
                (1,0)*sloline

                !! X%sr = f(X%bin)
                xc=xmap%rx(imx,ihx)
                yc=xmap%ry(imy,ihy)
if(pro%conv_wave) then
                call weislow_interp(.false.,X%s_g,sloline,xc,yc)
else
                call weislow_interp(.false.,X%s_s,sloline,xc,yc)
end if

                X%pba(imx,imy,ihx,ihy,:) = &
                X%pba(imx,imy,ihx,ihy,:) + &
                (0,1)*sloline

             end do
          end do
       end do
    end do

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine weislow_push
  !----------------------------------------------------------------
  subroutine weislow_interp(adj,X,line,xc,yc)
    type(weslow_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
    real    :: fx,fy
    real    :: w00,w01,w10,w11, sumw
#ifdef DEBUG
    character(len=128), parameter :: unit='weislow_interp'
    if(pro%debug) call in(unit)
#endif

    cube=>X%bin

    dx = xc
    dy = yc

    !! don't 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 = floor( 1+ (dx-X%amx%o) / X%amx%d )
    iy = floor( 1+ (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
       cube(ix,iy,:) = cube(ix,iy,:) + w00 * line
       cube(jx,iy,:) = cube(jx,iy,:) + w10 * line
       cube(ix,jy,:) = cube(ix,jy,:) + w01 * line
       cube(jx,jy,:) = cube(jx,jy,:) + w11 * line
    else         !! pull
       line    = &
       w00*real(cube(ix,iy,:)) + &
       w10*real(cube(jx,iy,:)) + &
       w01*real(cube(ix,jy,:)) + &
       w11*real(cube(jx,jy,:))
    end if

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine weislow_interp
  subroutine var_slow(iwb)
    integer :: iwb
    call create_modified_slowness(S%s_s,pro%z_samp(iwb)%dz,pro%z_samp(iwb)%zloc) 
    s%s_s%tag=s%s_s%tag_var
    if(pro%conv_wave) then
      call create_modified_slowness(S%s_g,pro%z_samp(iwb)%dz,pro%z_samp(iwb)%zloc) 
      s%s_g%tag=s%s_g%tag_var
    end if
  end subroutine var_slow
  !----------------------------------------------------------------


  !----------------------------------------------------------------
end module wei_obj_S
