module inv_tomo_mod
  use process_types
  use model_types
  use dottest
  use solver_reg_mod
  use solver_tiny_mod
  use solver_prc_mod
  use cgstep_mod
  use helix_chain_mod
  use i_op_mod
  use polydiv
  use sep
  implicit none
  real,allocatable,dimension(:),private     :: mask
  contains
!--------------------------------------------------------------------------------
 subroutine inv_tomo (param,model_param,nelem_model,slow,slow0,eps,oper1,oper2,oper3)

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

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

    interface
      function oper3 (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,reg,prc
    integer                                :: ierr,iter,type,nelem_model,nAop
    type (tomo_param_type)                 :: param
    type (model_param_type)                :: model_param
    real, dimension (nelem_model)          :: slow
    real, dimension (nelem_model)          :: slow0
    real                                   :: eps

    verb=.true.
    nAop=model_param%nz_mod*model_param%nx_mod*model_param%ny_mod
    call from_param("method",reg,.false.)
    call from_param("method",prc,.false.)
    if(not(reg).and.not(prc)) reg=.true.

    if (exist_file("mask")) then
       if (reg) then
          call solver_reg(m=slow,d=slow0,Fop=oper1,Aop=oper2,stepper=cgstep,nAop=nAop,&
            Jop=oper3,niter=param%niter,eps=eps,verb=.true.)
       else if(prc) then
          call solver_prc(m=slow,d=slow0,Fop=oper1,Sop=oper2,stepper=cgstep,nSop=nAop,&
            Jop=oper3,niter=param%niter,eps=eps,verb=.true.)
       end if
    else
       if (reg) then
          call solver_reg(m=slow,d=slow0,Fop=oper1,Aop=oper2,stepper=cgstep,nAop=nAop,&
            niter=param%niter,eps=eps,verb=.true.)
       else if(prc) then
          call solver_prc(m=slow,d=slow0,Fop=oper1,Sop=oper2,stepper=cgstep,nSop=nAop,&
            niter=param%niter,eps=eps,verb=.true.)
       end if
    end if
 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
