subroutine marquardt (i_loop, a_lambda, finished, u, at_limit)

use cesrv_struct
use cesrv_interface
use super_universe_com
use nr

implicit none

type (universe_struct), target :: u
type (universe_struct), pointer :: uu

real(rp), allocatable, save :: x(:), y(:), sig(:), a(:)
real(rp), allocatable, save :: covar(:,:), alpha(:,:)
real(rp), allocatable, save :: y_fit(:)
real(rp), allocatable, save :: dy_da(:, :)
real(rp) alambda, chi_sq
real(rp) merit0, a_lambda

integer i, j, k
integer i_loop, m_a, n_data

logical, allocatable, save :: mask_a(:)
logical :: finished, init_needed = .true.
logical :: at_limit

!

call merit_calc (merit0)

finished = .false.
at_limit = .false.
 
if (i_loop == 1 .or. init_needed) then
  a_lambda = -1
  m_a = count(u%var(:)%useit)
  n_data = m_a
  do j = 1, logic%u_num
    n_data = n_data + count(super%u_(j)%data(:)%useit_opt .and. &
                                  (super%u_(j)%data(:)%weight /= 0))
  enddo
  if (allocated(x)) deallocate(x, y, sig, a, covar, alpha, mask_a, &
                                                              y_fit, dy_da)
  allocate (x(n_data), y(n_data), sig(n_data), y_fit(n_data))
  allocate (a(m_a), covar(m_a,m_a), alpha(m_a,m_a), mask_a(m_a))
  allocate (dy_da(n_data, m_a))
  mask_a = .true.
  init_needed = .false.
endif

if (a_lambda > 0) a_lambda = max(a_lambda, 1e-28) ! so no underflow

if (a_lambda > 1e10) then
  print *, 'MARQUARDT: AT MINIMUM'
  finished = .true.
  init_needed = .true.   ! shouldn't need this but...
  a_lambda = 0  ! tell mrqmin we are finished
  alambda = a_lambda
  call mrqmin (x, y, sig, a, mask_a, covar, alpha, chi_sq, mrq_func, alambda)
  a_lambda = alambda
  return
endif

! transfer model to "a" and "y" arrays

k = 0
do i = 1, size(u%var)
  if (.not. u%var(i)%useit) cycle
  k = k + 1
  a(k) = u%var(i)%model
  y(k) = u%var(i)%model - u%var(i)%delta ! = %saved when no base or ref.
  if (u%var(i)%weight == 0) then
    sig(k) = 1e10
  else
    sig(k) = sqrt(1/u%var(i)%weight)
  endif
enddo

do j = 1, logic%u_num
  uu => super%u_(j)
  do i = 1, size(uu%data)
    if (.not. uu%data(i)%useit_opt) cycle
    if (uu%data(i)%weight == 0) cycle
    k = k + 1
    y(k) = 0
    sig(k) = sqrt(1/uu%data(i)%weight)
  enddo
enddo

! call mrqmin.
! save var%useit in case we hit a limit.

cesrv_common%at_limit_flag = .false.
cesrv_common%var_useit = u%var(:)%useit

alambda = a_lambda
call mrqmin (x, y, sig, a, mask_a, covar, alpha, chi_sq, mrq_func, alambda)
a_lambda = alambda
call mrq_func (x, a, y_fit, dy_da)  ! put a -> model

at_limit = cesrv_common%at_limit_flag

if (cesrv_common%at_limit_flag .or. i_loop == logic%opt_loops) then
  a_lambda = 0   ! tell mrqmin we are finished
  alambda = a_lambda
  call mrqmin (x, y, sig, a, mask_a, covar, alpha, chi_sq, mrq_func, alambda)
  a_lambda = alambda
  init_needed = .true.
endif

end subroutine