subroutine mrq_func (x, a, y_fit, dy_da)

  use cesrv_struct
  use cesrv_interface
  use super_universe_com

  implicit none

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

  real(rp) :: x(:), a(:)
  real(rp) :: y_fit(:)
  real(rp) :: dy_da(:, :)
  real(rp) merit0, r_step

  integer i, j, k, n, nn, im, iv

  logical limited

! transfer "a" array to model

  k = 0
  u => super%u_(1)
  u%energy_var%v(1)%model => u%orb(0)%vec(6)  ! This should not be needed

  do i = 1, size(u%var)
    u%var(i)%old = u%var(i)%model
    if (cesrv_common%var_useit(i)) then
      k = k + 1
      u%var(i)%model = a(k)
    endif
    u%var(i)%h = u%var(i)%model
  enddo

! if limited then readjust other variables in proportion and reinit mrqmin.

  call limit_calc (limited)

  if (limited) then

    r_step = 1
    do i = 1, size(u%var)
      if (cesrv_common%var_useit(i) .and. u%var(i)%h /= u%var(i)%old) then
        r_step = min (r_step, (u%var(i)%model - u%var(i)%old) / &
                                                   (u%var(i)%h - u%var(i)%old))
      endif
    enddo

    k = 0
    do i = 1, size(u%var)
      if (.not. cesrv_common%var_useit(i)) cycle
      u%var(i)%model = u%var(i)%old + r_step * (u%var(i)%h - u%var(i)%old)
      k = k + 1
      a(k) = u%var(i)%model
      do j = 1, logic%u_num
        uu => super%u_(j)
        uu%var(i)%model = u%var(i)%model
        call var_bookkeeper(uu%var(i), uu%ring, uu%orb)
      enddo
    enddo

    print *, 'Limit retreat factor:', r_step
    cesrv_common%at_limit_flag = .true.

  endif

! calculate derivatives

  k = 0
  do i = 1, size(u%var)
    if (.not. cesrv_common%var_useit(i)) cycle
    k = k + 1
    u%var(i)%model = a(k)
    do j = 1, logic%u_num
      super%u_(j)%var(i)%model = u%var(i)%model
      call var_bookkeeper(super%u_(j)%var(i), super%u_(j)%ring, super%u_(j)%orb)
    enddo
    y_fit(k) = u%var(i)%model
    dy_da(k,:) = 0
    dy_da(k,k) = 1
  enddo

! calculate model fit

  call merit_calc(merit0)

  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_fit(k) = uu%data(i)%delta
      if (.not. uu%ok_status) y_fit(k) = 1e10  ! something large
      im = uu%data(i)%ix_dmeas
      nn = 0
      do n = 1, size(u%var)
        if (cesrv_common%var_useit(n)) then
          nn = nn + 1
          iv = u%var(n)%ix_dvar
          dy_da(k,nn) = uu%dm_dv(im, iv)
        endif
      enddo
    enddo
  enddo

end subroutine