module OW_renorm
  use adj_mod
  use OW_type
  use OW_parms_areal
  use OW_wei_areal
  implicit none
contains
!============================================================
  function owrenorm_lop(adj,add,mm,dd) result(stat)
    logical, intent(in)         :: adj,add
    real              :: mm(:),dd(:)
    integer           :: stat

    call adjnull(adj,add,mm,dd)
    call owrenorm_exe(adj,add,mm,dd)
    stat=0
  end function
!============================================================
subroutine owrenorm_exe(adj,add,mm,dd)
logical             :: adj,add
complex             :: bwr(img%xm%n,nz,sou%w%n),bws(img%xm%n,nz,sou%w%n)
complex             :: dwr(img%xm%n,nz,sou%w%n),dws(img%xm%n,nz,sou%w%n)
real                :: dd(img%xm%n,nz,sou%w%n),mm(img%xm%n,nz)
integer             :: i,j,k,l

do i=1,rec%xd%n
   call owwemva_init(i)
   bwr=0;bws=0;dwr=0;dws=0
   call ow_wvfld_io(bwr,"bkwr","reed")
   call ow_wvfld_io(bws,"bkws","reed")
   call ow_wvfld_io(dwr,"dkwr","reed")
   call ow_wvfld_io(dws,"dkws","reed")
   if (adj) then
      do j=1,rec%w%n
         do k=1,nz
            do l=1,img%xm%n
               mm(l,k)=mm(l,k) + (1+dwr(l,k,j)*bwr(l,k,j)/((conjg(bwr(l,k,j))*bwr(l,k,j))+tiny(1.))+&
                                    dws(l,k,j)*bws(l,k,j)/((conjg(bws(l,k,j))*bws(l,k,j))+tiny(1.)))*dd(l,k,j)
            end do
         end do
      end do
   else
      do j=1,rec%w%n
         do k=1,nz
            do l=1,img%xm%n
               dd(l,k,j)=dd(l,k,j) + (1+dwr(l,k,j)*bwr(l,k,j)/((conjg(bwr(l,k,j))*bwr(l,k,j))+tiny(1.))+&
                                        dws(l,k,j)*bws(l,k,j)/((conjg(bws(l,k,j))*bws(l,k,j))+tiny(1.)))*mm(l,k)
            end do
         end do
      end do
   end if
end do   
end subroutine
end module
