module solver_cjg_mod
! 0 = F m - d
  use sep
  use OW_parms_areal
  use cdfomel_mod
  use steepest_mod
  implicit none

  integer                       :: nn
  real, private, allocatable    :: vel0(:), taper(:), pick1(:)

  contains
!----------------------------------------------------------------------------
  subroutine slvcjg_init(nn_in,vel0_in,taper_in,pick_in)
    integer          :: nn_in
    real             :: vel0_in(:), taper_in(:), pick_in(:)
    allocate(vel0(size(vel0_in)), taper(size(taper_in)), pick1(size(pick_in)))
    nn=nn_in
    vel0=vel0_in
    taper=taper_in
    pick1=pick_in
  end subroutine
!----------------------------------------------------------------------------
  subroutine solver_cjg(eps1,niter,l,u)
    real, dimension(nn)                :: d      ! data
    real, dimension(size(vel0))        :: m       ! model
    integer,               intent(in)  :: niter   ! number of iterations
    integer                            :: iter    ! iteration number
    real, dimension(size(vel0))        :: g       ! gradient (dm)
    real, dimension(nn)                :: rd      ! data residual (pointer)
    real, dimension(nn)                :: gd      ! conj grad     (pointer)
    double precision                   :: fn, ggg(size(vel0)), mm(size(vel0))
    real                               :: R_image(img%xm%n,nsecdim,img%z%n),eps1
    real,optional                      :: l(:), u(:)! bounds

    m=vel0
    write(0,*) "solver ",eps1 
    call cdplus_init(niter,nn,size(vel0),eps1)
    do iter = 1,niter  
       rd=0.;gd=0.;ggg=0;
       call fgeval_init(m,taper,iter,pick1)
       call feval(fn,R_image)
       write(0,*) " fn = ",fn
       call geval(ggg,R_image)
       g=sngl(ggg)/maxval(abs(ggg))
       call A2g(g,gd)
       rd=reshape(R_image,(/nn/))
       call fgeval_close()
       if (present(l)) then
          call cdplus( iter, m, g, rd, gd, l, u)
       else
          call cdplus( iter, m, g, rd, gd)
       end if
       mm=dble(m)
       call ow_smoovel(mm,vel0,pick1)
       m=sngl(mm)
       call to_history("n4",iter,"Velout")
       call srite("Velout",m,4*size(vel0))
    end do 
    call cdplus_close()
  end subroutine 
!----------------------------------------------------------------------------
  subroutine slvcjg_close()
    deallocate(vel0, taper, pick1)
  end subroutine
!----------------------------------------------------------------------------
end module 
