





  do iter=1, niter

!    perform line search to update the model
!    ---------------------------------------
     stat = bwi_linsec_oper(ftag_%csou_, ftag_%crec_, ftag_%bvel_, ftag_%bimg_, ftag_%hess_, ftag_%srch_, ftag_%imsk_, &
                            vgrid, nsrc, alph, fobj(iter-1), hcut, dsotyp, objtyp, objeps, wntreg, epsreg)

!    migrate with current model
!    --------------------------
     stat = bwi_kerimg_oper(T, F, ftag_%csou_, ftag_%crec_, ftag_%bvel_, ftag_%bimg_)

!    if Hessian is requested
!    -----------------------
     if (hgrid%wnthes) then
!       compute the imaging hessian
!       ---------------------------
        stat = bwi_kerhes_oper(F, ftag_%sfnc_, ftag_%gsou_, ftag_%grec_, ftag_%bvel_, ftag_%hess_)
!       apply the imaging Hessian
!       -------------------------
        stat = bwi_hesnor_oper(ftag_%hess_, igrid, hgrid)
     end if

!    compute the residual and objective function with current model
!    --------------------------------------------------------------
     stat = bwi_restom_oper(ftag_%bimg_, ftag_%hess_, ftag_%imsk_, igrid, hgrid, fobj(iter), hcut, dsotyp, objtyp, objep
s)

!    compute the gradient with the current model
!    -------------------------------------------
     stat = bwi_kertom_oper(T, F, ftag_%csou_, ftag_%crec_, ftag_%bvel_, ftag_%grdp_, ftag_%bimg_)

!    apply regularization
     if (wntreg) then
        stat = bwi_objreg_oper(ftag_%bvel_, vgrid, epsreg, fobj(iter))
        stat = bwi_grdreg_oper(ftag_%grdp_, ftag_%bvel_, vgrid, epsreg)
     end if

     if (smtgrd) then
        stat = bwi_grdsmt_oper(ftag_%grdp_, vgrid, rect1, rect2, rect3, nsmt)
     end if

!    scale the gradient
!    ------------------
     stat = bwi_gscale_oper(ftag_%grdp_, vgrid)

     if (vgrid%grdmsk) then
        stat = bwi_grdmsk_oper(ftag_%grdp_, ftag_%gmsk_, vgrid) 
     end if 
!    get the search direction
!    ------------------------
     select case (method)
       case ("sdstep")
            stat = bwi_sdstep_oper(ftag_%srch_, ftag_%grdp_, vgrid)
       case ("cgstep")
            stat = bwi_cgstep_oper(ftag_%srch_, ftag_%grdo_, ftag_%grdp_, vgrid)
     end select

     if (modmov) then
        call sseek(ftag_%bvel_, 0, 0)
        call sep_read(vmod, ftag_%bvel_)
        call sep_write(vmod, mmov_)
     end if 

     write(0,*) "iter = ", iter, fobj(iter), alph

  end do

  if (objmov) call sep_write(real(fobj), fobj_)


  write(0,*)"***********************"
  write(0,*)"      Finished!!!      "
  write(0,*)"***********************"





























  integer function bwi_linsec_oper(csou_, crec_, bvel_, bimg_, hess_, srch_, imsk_, vgrid, nsrc, alph, fobj, &
                                   hcut, dsotyp, objtyp, objeps, wntreg, epsreg) result(stat)
    character(len=*) :: csou_, crec_, bvel_, bimg_, hess_, srch_, imsk_
    type(vel_grid)   :: vgrid
    integer          :: nsrc
    real             :: alph
    double precision :: fobj
    real             :: hcut
    integer          :: dsotyp
    character(len=6) :: objtyp
    real             :: objeps 
    logical          :: wntreg
    real             :: epsreg
 
    double precision :: f(nsrc), fopt, norm
    real             :: sfac
    integer          :: iopt, isrc

    sfac = 0.5*(nsrc+1)

!   read the background model to be updated
!   ---------------------------------------
    call bwi_fio_vmod_updt(grdo, vgrid, bvel_, "reed")

!   read the search direction
!   -------------------------
    call bwi_fio_vmod_updt(srch, vgrid, srch_, "reed")

!   normalize
!   ---------
    norm = maxval(abs(srch))
    srch = srch/norm

!   perturb the model
!   -----------------
    do isrc=1, nsrc
!      trial velocity model
!      -------------------- 
       grdp = grdo + alph*isrc*srch
       call bwi_fio_vmod_updt(grdp, vgrid, bvel_, "rite")
!      compute the image using the trial velocity
!      ------------------------------------------
       stat = bwi_kerimg_oper(.true., .false., csou_, crec_, bvel_, bimg_)
!      compute the objective function
!      ------------------------------
       stat = bwi_restom_oper(bimg_, hess_, imsk_, igrid, hgrid, f(isrc), hcut, dsotyp, objtyp, objeps)

       if (wntreg) stat = bwi_objreg_oper(bvel_, vgrid, epsreg, fobj)
!       write(0,*)"alph = ", alph*isrc, "f = ", f(isrc)
    end do

!   find the minimum one
!   --------------------
    fopt = f(1)
    iopt = 1
    do isrc=1, nsrc
       if (f(isrc)<=fopt) then
          iopt = isrc
          fopt = f(isrc)
       end if
    end do

!   compare with the current objective function value
!   -------------------------------------------------
    if (f(iopt)<=fobj) then
       grdp = grdo + alph*iopt*srch
       alph = alph*iopt/sfac
       fobj = f(iopt)
    else
       grdp = grdo
       alph = alph * 0.5
    end if

    call bwi_fio_vmod_updt(grdp, vgrid, bvel_, "rite")

    srch = 0.; grdo = 0.; grdp = 0.

    stat = 0
  end function
