module vrms2int_mod 
! Transform from RMS to interval velocity
  use causint
  use igrad1
  use cgstep_mod
  use solver_mod
  use gradchain
implicit none
  contains
! --------------------------------------------------------------------------------
  subroutine vrms2int( niter, eps, weight, refl, vrms, vint, prec, err, res )
    integer,            intent( in)     :: niter     ! iterations
    real,               intent( in)     :: eps       
    ! scaling parameter
    real, dimension(:), intent( in out) :: vrms      
    ! RMS velocity
    real, dimension(:), intent( out)    :: vint      
    ! interval velocity
    real, dimension(:), intent( in out) :: err       
    ! residual error per iteration
    real, dimension(:), intent( in out) :: res       
    ! final residual error
    real, dimension(:), pointer         :: weight    ! data weighting
    real, dimension(:), pointer         :: refl      
    ! picked strong reflections
    integer                             :: st,it,nt  
    logical                             :: prec
    logical, dimension( size( vint))    :: mask
    real,    dimension( size( vrms))    :: dat 
    nt = size( vrms)
    do it= 1, nt  
      dat( it) = vrms( it) * vrms( it) * it
      weight( it) = weight( it)*(1./it)           
      ! decrease weight with time
    end do 
    mask = .false.
    mask( 1) = .true.            ! constrain first point
    vint = 0.     
    vint( 1) = dat( 1)
    if ( prec) then
!			call gradchain_init( weight )
!
!      call solver_prec( causint_lop, cgstep, x=vint, dat=dat, niter=niter, known=mask,
!                        x0=vint, wt=weight, prec=weighted_causint_lop, eps=eps, nprec=nt,
!                        err=err, res=res )
      call solver_prec( causint_lop, cgstep, x=vint, dat=dat, niter=niter&
        &, known=mask,x0=vint, wt=weight, prec=causint_lop, eps=eps,&
        & nprec=nt,err=err, res=res )
      st = causint_lop( .false., .false., vint, dat)
!      st = weighted_causint_lop( .false., .false., vint, dat)
    else
!      call solver_reg( causint_lop, cgstep, igrad1_lop, nt, vint, dat, niter, eps,
!                       known=mask, x0=vint, wt=weight, err=err, res=res )
      if ( sum(refl) < nt) then
        refl = 1.0-refl
      end if
      call gradchain_init( refl )
      call solver_reg( causint_lop, cgstep, weighted_grad_lop, nt, vint&
        &, dat, niter, eps,known=mask, x0=vint, wt=weight, err=err,&
        & res=res,verb=.true. )
    end if 
    res = res*weight
    call cgstep_close()
    do it= 1, nt  
      vrms( it) = sqrt( dat( it)/it)
    end do 
    vint = sqrt( vint)
  end subroutine
! --------------------------------------------------------------------------------
end module 
