module OW_rho

  use OW_type
  use OW_parms_rho
  use adj_mod

  implicit none
  real, allocatable, private  :: dips(:,:)

  contains

!===============================================================================
  function ow_rho_init(dips_in) result(stat)
    integer                     :: stat
    real                        :: dips_in(img%z%n,img%xm%n)
    allocate(dips(img%z%n,img%xm%n))
    dips=dips_in
    stat=0
  end function 
!===============================================================================
  function ow_rho_lop(adj, add, modl, data) result(stat)
    integer                     :: stat 
    logical,intent(in)          :: adj,add 
    real                        :: modl(:)
    real                        :: data(:)

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

!! internal variables
    integer		:: inode,ix,iah,irho,izm,izd,idx,nodes
    real			:: aah,tgh,rhox,zm,zd,rdx,fx,gx
    integer, external         :: omp_get_thread_num,omp_get_num_threads

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

!$OMP PARALLEL DO PRIVATE  (inode,ix,iah,aah,tgh,irho,rhox,izm,zm,zd,rdx,idx,izd,fx,gx) schedule(dynamic,1)  num_threads(nodes) 
do ix=1, img%xm%n
   inode = omp_get_thread_num()+1
   do iah=1,img%ah%n
      aah = img%ah%o + (iah-1)*img%ah%d
      tgh = tand(aah)*tand(aah)
      do irho=1,img%rho%n
         rhox=(irho-1)*img%rho%d+img%rho%o
         do izm=1,img%z%n
            zm  = img%z%o + (izm-1) * img%z%d
            zd  = zm - (rhox -1)*tgh*zm
            rdx = (zd - img%z%o) / img%z%d
            idx = rdx
            izd = idx + 1
            fx = rdx - idx
            gx = 1.0 - fx
            if ( izd >= img%z%n .or. izd < 1) cycle
            if (adj) then
               modl(izm,irho,ix) = modl(izm,irho,ix) +  (gx*data(izd,iah,ix) +fx*data(izd+1,iah,ix))
            else
               data(izd+1,iah,ix)  = data(izd+1,iah,ix)  + modl(izm,irho,ix)*fx
               data(izd,iah,ix  )  = data(izd,iah,ix  )  + modl(izm,irho,ix)*gx
            end if 
         end do
      end do
   end do
   if (mod(ix,20).eq.1) write(0,*) "DONE CIG", ix, " FROM ", img%xm%n
end do
!$OMP END PARALLEL DO
 end subroutine  
!===============================================================================
  function ow_rho_close() result(stat)
    integer             :: stat
    deallocate(dips)
    stat=0
  end function 
!=================================================================================  
end module 
