subroutine merit_calc (merit)

use cesrv_struct
use cesrv_interface
use super_universe_com

implicit none

type (universe_struct), pointer :: u

integer i, j
real(rp) merit, ix, psum

! sum over all universes

merit = 0

do j = 1, logic%u_num

  u => super%u_(j)
  u%data(:)%merit = 0

! merit from the data.
! Phase renormalization: The average delta is, by definition, zero.

  if (.not. logic%opt_running) call set_data_useit_opt (u%data)
  call ring_calc (u) ! recalc the twiss parameters

  if (.not. u%ok_status) return

  if (plot_type_has(plot_ref$, logic%opt_base) .and. &
                   plot_type_has(plot_model_and_base$, logic%opt_base)) then
    where (u%data(:)%useit_opt) u%data(:)%delta = u%data(:)%model - u%data(:)%meas + u%data(:)%ref - u%data(:)%base_model

  elseif (plot_type_has(plot_ref$, logic%opt_base)) then
    where (u%data(:)%useit_opt) u%data(:)%delta = u%data(:)%model - u%data(:)%meas + u%data(:)%ref - u%data(:)%design

  elseif (plot_type_has(plot_model_and_base$, logic%opt_base)) then
    where (u%data(:)%useit_opt) u%data(:)%delta = u%data(:)%model - u%data(:)%meas - u%data(:)%base_model

  else
    where (u%data(:)%useit_opt) u%data(:)%delta = u%data(:)%model - u%data(:)%meas
  endif

  psum = sum(u%phase%x%d(:)%delta, mask=u%phase%x%d(:)%useit_opt)
  ix = count(u%phase%x%d(:)%useit_opt)
  if (ix /= 0) where (u%phase%x%d(:)%useit_opt) u%phase%x%d(:)%delta = u%phase%x%d(:)%delta - psum / ix

  psum = sum(u%phase%y%d(:)%delta, mask=u%phase%y%d(:)%useit_opt)
  ix = count(u%phase%y%d(:)%useit_opt)
  if (ix /= 0) where (u%phase%y%d(:)%useit_opt) u%phase%y%d(:)%delta = u%phase%y%d(:)%delta - psum / ix

  where (u%data(:)%useit_opt) u%data(:)%merit = u%data(:)%weight * u%data(:)%delta**2
  merit = merit + sum(u%data(:)%merit, mask = u%data(:)%useit_opt)

enddo

! merit from the variables

u => super%u_(1)
u%var(:)%merit = 0

do i = 1, size(u%var)
  if (.not. u%var(i)%useit) cycle

  if (plot_type_has(plot_ref$, logic%opt_base) .and. &
                   plot_type_has(plot_model_and_base$, logic%opt_base)) then
    u%var(i)%delta = (u%var(i)%model - u%var(i)%base_model) - (u%var(i)%saved - u%var(i)%saved_ref)
  elseif (plot_type_has(plot_ref$, logic%opt_base)) then
    u%var(i)%delta = (u%var(i)%model - u%var(i)%design) - (u%var(i)%saved - u%var(i)%saved_ref)
  elseif (plot_type_has(plot_model_and_base$, logic%opt_base)) then
    u%var(i)%delta = u%var(i)%model - u%var(i)%saved - u%var(i)%base_model
  else
    u%var(i)%delta = u%var(i)%model - u%var(i)%saved
  endif
  u%var(i)%merit = logic%u_num * u%var(i)%weight * u%var(i)%delta**2
enddo

merit = merit + sum(u%var(:)%merit, mask = u%var(:)%useit)

end subroutine