module inv_tomo_mod
  use process_types
  use model_types
  use dottest
  use solver_reg_mod
  use solver_tiny_mod
  use cgstep_mod
  use helix_chain_mod
  use i_op_mod
  use polydiv
  use OW_steer_parms
  use sep
 
  implicit none
  real,allocatable,dimension(:),private     :: mask
  contains
!--------------------------------------------------------------------------------
 subroutine inv_tomo (slow,slow0,nelem_model)!,weight)

!    interface
!      function oper (adj, add, x, dat) result (stat)
!        integer :: stat
!        logical, intent (in) :: adj, add
!        real, dimension (:) :: x, dat
!      end function
!    end interface

    logical                                :: add,adj,verb
    integer                                :: ierr,iter,type,nAop,nelem_model
!    type (model_param_type)                :: model_param
    real, dimension (nelem_model)          :: slow
    real, dimension (nelem_model)          :: slow0
    real                                   :: eps

    nAop=model_param%nz_mod*model_param%nx_mod*model_param%ny_mod
    if (allocated(mask)) then
!       call solver_reg(m=slow,d=slow0,Fop=polydiv_lop,Aop=i_op,stepper=cgstep,nAop=nAop,&
       call solver_reg(m=slow,d=slow0,Fop=i_op,Aop=polydiv_lop,stepper=cgstep,nAop=nAop,&
            Jop=mask_op,niter=steerniter,eps1=epspoly,verb=.true.)
    else
       call solver_reg(m=slow,d=slow0,Fop=polydiv_lop,Aop=i_op,stepper=cgstep,nAop=nAop,&
            niter=steerniter,eps1=epspoly,verb=.true.)
    end if
    call cgstep_close()
 end subroutine
!--------------------------------------------------------------------------------
 subroutine mask_init(mask_in)
   real, dimension(:)        :: mask_in
   allocate(mask(size(mask_in)))
   mask=mask_in
 end subroutine
!--------------------------------------------------------------------------------
function mask_op(adj,add,m,d) result(stat)
  real,dimension(:)          :: m,d
  logical                    :: adj,add
  integer                    :: stat
  call adjnull (adj,add,m,d)
  call mask_op2(adj,add,m,d)
  stat=0
end function
!--------------------------------------------------------------------------------
subroutine mask_op2(adj,add,m,d)
  real,dimension(:)          :: m,d
  logical                    :: adj,add
  if (adj) then
     m=m+d*mask
  else
     d=d+m*mask
  end if
end subroutine
!--------------------------------------------------------------------------------
subroutine mask_close()
  deallocate(mask)
end subroutine
!--------------------------------------------------------------------------------
end module
