module lapfillx  
  use laplac2
  use cgstep_mod
  use adj_mod
  use random_mod
  use chain0_mod
  use solver_reg_mod
!    use solver_reg1_mod
  implicit none
  integer, private, save                     :: m1, m2, seed
  real,    private, save                     :: scale
  logical, private, save                     :: verb
  real, dimension(:), private, save, pointer :: resid_wt, model_wt,&
    & mask
  contains
!---------------------------------------------------------------------------------
  subroutine lapfill2_init( m1_in, m2_in, verb_in, mask_in, scale_in,&
    & seed_in,   resid_wt_in, model_wt_in )
    optional                    :: resid_wt_in, model_wt_in
    integer                     :: m1_in, m2_in, seed_in
    logical                     :: verb_in
    real                        :: scale_in
    real, dimension(:), pointer :: resid_wt_in, model_wt_in
    real, dimension(:)          :: mask_in
    m1=m1_in
    m2=m2_in
    verb = verb_in
    scale = scale_in
    seed = seed_in
    allocate( resid_wt(m1*m2) )
    resid_wt = 1
    if ( present(resid_wt_in)) then
      resid_wt = resid_wt_in
    end if
    allocate( model_wt(m1*m2) )
    model_wt = 1
    if ( present(model_wt_in)) then
      model_wt = model_wt_in
    end if
    allocate( mask(m1*m2) )
    mask=1
    where(mask_in.eq.0.0)
      mask = 0
    end where
  end subroutine
!---------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
  subroutine lapfill2x( niter, model, eps )
    integer,                    intent (in)     :: niter
    integer                                     :: stat,i
    real                                        :: mrm,mmo,srm,smo,varv
    real,                       intent (in)     :: eps
    real,    dimension (m1*m2), intent (in out) :: model          
    ! model
    real,    dimension (m1*m2)			:: dd          
    real,    dimension (:), pointer             :: xnull,resm          
    ! random
    allocate( xnull(m1*m2),resm(m1*m2) )
    call random_init(seed)
    call rand_real_array_1d( xnull)

    where( mask.eq.1.0)
      xnull=0.0
    end where
!
!   Computes the variances
!
    varv=dot_product(xnull-sum(xnull)/size(xnull),xnull-sum(xnull)/size(xnull))
    scale=sqrt(sqrt(dot_product(model-sum(model)/size(model),model-sum(model)/size(model))/varv))
!    scale=20.
    xnull = scale*(xnull-0.5)
    write(0,*) scale

    write(0,*) 'min/max/size(xnull)',minval(xnull),maxval(xnull),size&
      &(xnull)
    write(0,*) 'min/max/size(model)',minval(model),maxval(model),size&
      &(model)

    call solver_reg(m=model, d=model, &
    Fop=identity_lop,Aop=lapwght_lop,stepper=cgstep,nAop=m1*m2,niter=niter,eps=eps,&
    rm0=xnull, resm=resm, verb=.false.)

    deallocate(xnull)
    call cgstep_close()
  end subroutine
!---------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
  integer function identity_lop( adj, add, model, data)
    real, dimension(:)            :: model, data
    logical                       :: adj,add
    call adjnull( adj, add, model, data)
    if (adj) then
      model = model + data*mask
    else
      data  = data + model*mask
    end if 
    identity_lop = 0
  end function
!---------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
  integer function lapwght_lop( adj, add, model, data)
    integer                       :: stat
    real, dimension(:)            :: model, data
    real, dimension(:), pointer   :: temp
    logical                       :: adj,add
    allocate( temp(m1*m2) )
    temp=0
    call laplac2_init( m1, m2)
    call chain0( modelwt_lop, laplac2_lop, adj, add, model, data, temp&
      & )
    call laplac2_close()
    deallocate( temp )
    lapwght_lop = 0
  end function
!---------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
  integer function modelwt_lop( adj, add, model, data)
    real, dimension(:)            :: model, data
    logical                       :: adj,add
    call adjnull( adj, add, model, data)
    if (adj) then
      model = model + data*model_wt
    else
      data  = data + model*model_wt
    end if 
    modelwt_lop = 0
  end function
!---------------------------------------------------------------------------------
end module 
