module velsimp6  

  use adj_mod

  implicit none
  integer, private :: nt 
  integer, private :: nx 
  integer, private :: ns 
  real, private :: t0 
  real, private :: dt 
  real, private :: x0 
  real, private :: dx 
  real, private :: s0 
  real, private :: ds 
  integer, private :: nthreads

  contains

  subroutine velsimp6_init ( t0_in,dt_in,x0_in,dx_in,s0_in,ds_in,nt_in&
    &,nx_in,ns_in )

!$  integer, external  :: omp_get_num_threads,omp_get_thread_num

    real    :: t0_in 
    real    :: dt_in 
    real    :: x0_in 
    real    :: dx_in 
    real    :: s0_in 
    real    :: ds_in 

    integer    :: nt_in 
    integer    :: nx_in 
    integer    :: ns_in 
    character(len=128) :: desc

    t0 = t0_in 
    dt = dt_in 
    x0 = x0_in 
    dx = dx_in 
    s0 = s0_in 
    ds = ds_in 
    nt = nt_in 
    nx = nx_in 
    ns = ns_in

    nthreads = 1
!!$OMP PARALLEL
!!$      nthreads = omp_get_num_threads()
!!$OMP END PARALLEL
!    write(0,*) 'nthreads = ',nthreads

  end subroutine

  function velsimp6_lop ( adj, add, modl, data) result(stat)

    integer            :: stat 
    logical,intent(in) :: adj,add 
    real,dimension(:)  :: modl,data 
    call adjnull (adj,add,modl,data )
    call velsimp6_lop2(adj,modl,data )

    stat=0

  end function 

  subroutine velsimp6_lop2(adj,modl,data)

    logical,intent(in) :: adj 
    real, dimension (nt,ns)  :: modl 
    real, dimension (nt,nx)  :: data 
    integer it,iz,ix,is,nread,nth,ith
    real x,s,sx,t,z,wt    

if (adj) then

!!$OMP PARALLEL DO PRIVATE(is,s,ix,x,iz,z,sx,t,it,wt)
     do is= 1, ns

      s = s0 + (is-1) * ds               ! slowness

      do ix= 1, nx  
        x = x0 + (ix-1) * dx             ! offset

        do iz= 2, nt  

          z = t0 + (iz-1) * dt
          sx = abs( s * x)               ! normalized offset
          t = sqrt( z * z + sx * sx)
          it = 1.5 + (t-t0)/dt 

           if (it>nt) then
            exit
           end if

          wt = (z/t) / sqrt(t)

            modl(iz,is) =   modl(iz,is) + data(it,ix)*wt

        end do
      end do

    end do 
!!$OMP END PARALLEL DO 

else

!!$OMP PARALLEL DO PRIVATE(ix,x,is,s,iz,z,sx,t,it,wt)
   do ix= 1, nx  

     x = x0 + (ix-1) * dx                ! offset

     do is= 1, ns

      s = s0 + (is-1) * ds               ! slowness

        do iz= 2, nt  

          z = t0 + (iz-1) * dt
          sx = abs( s * x)               ! normalized offset
          t = sqrt( z * z + sx * sx)
          it = 1.5 + (t-t0)/dt 

           if (it>nt) then
            exit
           end if

          wt = (z/t) / sqrt(t)

            data(it,ix) =   data(it,ix) + modl(iz,is)*wt

        end do
      end do

    end do 
!!$OMP END PARALLEL DO 

end if

end subroutine velsimp6_lop2

subroutine velsimp6_close()
end subroutine velsimp6_close

end module
