module OW_linterp

use OW_type
use OW_parms_areal

  implicit none

    real,   private, allocatable        :: rhoxz1(:,:)
    real,   private                     :: orho1,drho1
    integer,private                     :: nrho1,nodes
  
  contains

!===============================================================================
  function owlint_init(rhoxz_in,nrho_in,orho_in,drho_in) result(stat)
    real                        :: rhoxz_in(:,:),orho_in,drho_in
    integer                     :: stat,nrho_in,nodes_in
    integer, external         :: omp_get_thread_num,omp_get_num_threads

    allocate(rhoxz1(img%xm%n,img%z%n))
    stat=1
    rhoxz1=rhoxz_in
    nrho1=nrho_in
    orho1=orho_in
    drho1=drho_in
    if (nodes_in.eq.9999) then
    !$OMP PARALLEL
      nodes = omp_get_num_threads()
    !$OMP END PARALLEL
    else
      nodes=node
    end if
    write(3,*) "****************"
    write(3,*) "* NUM.THREADS = ",nodes
    write(3,*) "****************"
    stat=0
  end function
!===============================================================================
  function owlint_lop(adj, add, modl, data) result(stat)
    integer                     :: stat 
    logical,intent(in)          :: adj,add 
    real		  	:: modl(img%xm%n,nsecdim,img%z%n) 
    real		  	:: data(img%xm%n,nsecdim,img%z%n,nrho)

    if (adj) then
       modl=0.
    else
       data=0.
    end if
    call owlint_exe(adj,add,modl,data )
    stat=0
  end function 
!=================================================================================  
  subroutine owlint_exe(adj,add,modl,data)
    logical,intent(in)        :: adj,add
    real		  	:: modl(img%xm%n,nsecdim,img%z%n) 
    real		  	:: data(img%xm%n,nsecdim,img%z%n,nrho)

!! internal variables
    integer		:: inode
    integer 		:: iax,ix,ixm,iah,ih,izm,izd
    integer		:: minx,maxx,indd,indm,idx
    real			:: ax,tgx,xx,ah,tgh,hx,zm,zd,rdx,fx,gx

    integer, external         :: omp_get_thread_num,omp_get_num_threads

if (nodes.eq.9999) then
!$OMP PARALLEL
    nodes = omp_get_num_threads()
!$OMP END PARALLEL
end if

!$OMP PARALLEL DO PRIVATE  (ix,ih,izm,izd,idx,rdx,fx,gx) schedule(dynamic,1)  num_threads(nodes) 
                               
    do izm=1, img%z%n
       do ih=1,nsecdim
          do ix=1, img%xm%n
             rdx = (rhoxz1(ix,izm) - orho1) / drho1 + 1
             idx = rdx
             izd = idx + 1
             fx = rdx - idx
             gx = 1.0 - fx
             if ( izd > nrho .or. izd < 1) cycle
             if (adj) then
                modl(ix,ih,izm) = modl(ix,ih,izm) +  gx*data(ix,ih,izm,idx) +fx*data(ix,ih,izm,izd)
             else
                data(ix,ih,izm,izd) = data(ix,ih,izm,izd) + modl(ix,ih,izm)*fx
                data(ix,ih,izm,idx) = data(ix,ih,izm,idx) + modl(ix,ih,izm)*gx
             end if 
          end do
       end do
       if (mod(izm,20).eq.1) write(0,*) "DONE DEPTH", izm, " FROM ", img%z%n
    end do
!$OMP END PARALLEL DO
 end subroutine  
!===============================================================================
  function owlint_close() result(stat)
    integer                     :: stat
    if (allocated(rhoxz1)) deallocate(rhoxz1)
    stat=0
  end function
!===============================================================================

end module 
