!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------

subroutine fletcher_reeves (i_loop, i_cyc, finished, u)

  use cesrv_struct
  use cesrv_interface
  use nr

  implicit none

  type (universe_struct), target :: u

  real(rp) small_merit, merit0, merit1, merit_last_init
  real(rp) h_sum, gg, dgg, dm_dv
  real(rp) brent2, xmin, gam, tol
  real(rp) ax, bx, cx, fa, fb, fc

  integer i_loop, i_cyc, i, ix

  logical end_loop, finished, reinit / .false. /, init_needed / .true. /

	interface
		function merit_1dim (x)
		use nrtype
		implicit none
		real(dp), intent(in) :: x
		real(dp) :: merit_1dim
		end function merit_1dim
	end interface

! Init

  if (init_needed) then
    small_merit = 1.0e-3
    init_needed = .false.
  endif

! Set up for first time through (I_LOOP = 1).
! and load gradiants into G and H

  finished = .false.
  end_loop = .false.

  if (i_loop == 1 .or. reinit) then
    reinit = .false.
    call merit_calc(merit0)
    if (i_loop == 1) merit_last_init = merit0
    call dmerit_calc ('normal')

    do ix = lbound(u%var, 1), ubound(u%var, 1)
      if (u%var(ix)%useit) then
        u%var(ix)%G = -u%var(ix)%dmerit * u%var(ix)%step
        u%var(ix)%H = u%var(ix)%G
      endif
    enddo

  endif

! Optimize. This is essentually FRPRMN from Numerica Recipes pg 305
! main loop

  do i = 1, i_cyc

! locate minimum along H
! BRENT is the routine that actually varies the quad k's
! The strange factors of 1.0e10 are for BRENT which does not like to
! have very small numbers.

    H_sum = 0.0

    do ix = lbound(u%var, 1), ubound(u%var, 1)
      if (u%var(ix)%useit) then
        u%var(ix)%old = u%var(ix)%model
        H_sum = H_sum + abs(u%var(ix)%H)
      endif
    enddo

    if (H_sum == 0) then
      print *, 'OPTIMIZER: AT MINIMUM!'
      finished = .true.
      return
    endif

! Go the the minimum in the direction indicated by u%var()%H

    if (abs(h_sum) > 1) then
      logic%h_scale = 1e10
    else
      logic%h_scale = 1
    endif

    ax = 0
    bx = logic%h_scale / h_sum

    call mnbrak (ax, bx, cx, fa, fb, fc, merit_1dim)   ! NR pg 281

    tol = logic%opt_tolerance * logic%h_scale * 1e-10
    merit1 = brent2 (ax, bx, cx, merit_1dim, tol, xmin) ! NR

! Reinitialize if merit not changing too much
! Terminate if merit not changing much from last reinit

    if (merit1 .le. 1.0e-4 * small_merit) then
      print *, 'OPEIMIZER: SUCCESS'
      finished = .true.
      return
    endif

    if (abs(merit0 - merit1) .le. 1.0e-4*(merit0+merit1+small_merit)) then
      i_cyc = i
      reinit = .true.    ! reinitialize next loop
      if (abs(merit_last_init - merit1) .le. &
                    logic%change_min*(merit0+merit1+small_merit)) then
        print *, 'OPTIMIZER: APARENTLY AT MIN'
        finished = .true.
      endif
      merit_last_init = merit1
      return
    endif

!

    call merit_calc(merit0)
    call dmerit_calc ('normal')

    GG = 0
    DGG = 0

!

    do ix = lbound(u%var, 1), ubound(u%var, 1)
      if (u%var(ix)%useit) then
        dm_dv = u%var(ix)%dmerit * u%var(ix)%step
        GG = GG + u%var(ix)%G**2
        DGG = DGG + (dm_dv + u%var(ix)%G) * dm_dv
      endif
    enddo

!

    if (GG == 0) then
      print *, 'OPTIMIZER: SUCCESS, ZERO GRADIANT'
      finished = .true.
      return
    endif

    GAM = DGG / GG

!

    do ix = lbound(u%var, 1), ubound(u%var, 1)
      if (u%var(ix)%useit) then
        u%var(ix)%G = -u%var(ix)%dmerit * u%var(ix)%step
        u%var(ix)%H = u%var(ix)%G + GAM * u%var(ix)%H
      endif
    enddo

  enddo

end subroutine