subroutine dmerit_calc (what)

  use cesrv_struct
  use cesrv_interface
  use super_universe_com

  implicit none

  type (universe_struct), pointer :: u1, uu
  type (lat_struct), pointer :: ring
  type (ele_struct), pointer :: eles(:)
  type (coord_struct), pointer :: orb(:)

  integer i, j, k, im, iv, ix

  real(rp) model0, phase_ave
  real(rp), pointer :: dm_dv(:,:)

  character(*) what  ! 'normal' or 'reinit'

  logical read_digested, ok, version_ok
  logical, allocatable :: msk(:)

! some preliminary stuff.
! Do not try to do anything until the ring has been initialized

  u1 => super%u_(1)
  orb   => u1%orb
  ring  => u1%ring
  eles  => ring%ele

  if (.not. logic%ring_initialized) return

  if (.not. logic%opt_running) then
    do j = 1, logic%u_num
      uu => super%u_(j)
      call set_data_useit_opt (uu%data)
    enddo
  endif

  if (what == 'reinit') then
    do j = 1, logic%u_num
      uu => super%u_(j)
      uu%var(:)%dm_dv_computed = .false.
    enddo
  endif

! reinit %dm_dv array with %ix_dvar and %ix_dmeas pointers if needed

  do j = 1, logic%u_num

    uu => super%u_(j)

    if (.not. associated(uu%dm_dv)) then

      uu%var(:)%ix_dvar = 0
      uu%var(:)%dm_dv_computed = .false.
      uu%data(:)%ix_dmeas = 0

      iv = 0
      do i = 1, size(uu%var)
        if (uu%var(i)%good_opt) then
          iv = iv + 1
          uu%var(i)%ix_dvar = iv
        endif
      enddo

      im = 0
      do i = lbound(uu%data, 1), ubound(uu%data, 1)
        if (uu%data(i)%good_opt) then
          im = im + 1
          uu%data(i)%ix_dmeas = im
        endif
      enddo

      allocate (uu%dm_dv(im, iv))
      uu%dm_dv = 0

      print '(a, 2i5)', ' Dmerit: Allocated dm_dv matrix', im, iv

    endif

  enddo

!---------------------------------------------------------
! if u%dm_dv table needs to be updated...

  if (any(u1%var%useit .and. .not. u1%var%dm_dv_computed)) then

! read from digested dmerit file if not a reinit

   read_digested = (logic%ok_to_read_dmerit_file) .and. &
                   (what /= 'reinit') .and. (logic%opt_vars /= opt_custom$)

   if (read_digested) then
     call read_digested_dmerit (ok, version_ok, u1)
     do j = 2, logic%u_num
       super%u_(j)%dm_dv = u1%dm_dv
     enddo
   endif

! here if a computation is needed to update u%dm_dv

    if (any(u1%var%useit .and. .not. u1%var%dm_dv_computed)) then
      print *, 'Dmerit: Init', count(u1%var%useit .and. .not. u1%var%dm_dv_computed)

      logic%dmerit_calc_on = .true.

      do k = 1, logic%u_num

        uu => super%u_(k)
        dm_dv => uu%dm_dv
        orb   => uu%orb
        ring  => uu%ring
        eles  => ring%ele

        call ring_calc (uu)
        uu%data(:)%old = uu%data(:)%model

        do i = lbound(uu%var, 1), ubound(uu%var, 1)
          if (uu%var(i)%useit .and. .not. uu%var(i)%dm_dv_computed) then
            iv = uu%var(i)%ix_dvar
            model0 = uu%var(i)%model
            uu%var(i)%model = uu%var(i)%model + uu%var(i)%step
            call var_bookkeeper(uu%var(i), uu%ring, uu%orb)
            call ring_calc (uu)
            do j = lbound(uu%data, 1), ubound(uu%data, 1)
              im = uu%data(j)%ix_dmeas
              if (im /= 0) then
                uu%dm_dv(im, iv) = (uu%data(j)%model - uu%data(j)%old) / uu%var(i)%step
              endif
            enddo
            uu%var(i)%model = model0
            uu%var(i)%dm_dv_computed = .true.
            call var_bookkeeper(uu%var(i), uu%ring, uu%orb)
          endif
        enddo

      enddo

      logic%dmerit_calc_on = .false.

    endif

  endif

!---------------------------------------------------------
! renormalize dm_dv for phase data to take into account that the average of
! phase_x and the average of phase_y is adjusted to be zero.

  do k = 1, logic%u_num

    uu => super%u_(k)

    if (allocated(msk)) deallocate(msk)
    allocate (msk(size(uu%dm_dv, 1)))

    ix = count(uu%phase%x%d(:)%useit_opt)

    if (ix /= 0) then
      msk = .false.
      do i = lbound(uu%phase%x%d, 1), ubound(uu%phase%x%d, 1)
        if (uu%phase%x%d(i)%useit_opt) msk(uu%phase%x%d(i)%ix_dmeas) = .true.
      enddo
      do i = 1, size(uu%var)
        if (.not. uu%var(i)%good_opt) cycle
        iv = uu%var(i)%ix_dvar
        dm_dv => uu%dm_dv
        phase_ave = sum(uu%dm_dv(:, iv), mask = msk) / ix
        where (msk) uu%dm_dv(:, iv) = uu%dm_dv(:, iv) - phase_ave
      enddo
    endif

!

    ix = count(uu%phase%y%d(:)%useit_opt)

    if (ix /= 0) then
      msk = .false.
      do i = lbound(uu%phase%y%d, 1), ubound(uu%phase%y%d, 1)
        if (uu%phase%y%d(i)%useit_opt) msk(uu%phase%y%d(i)%ix_dmeas) = .true.
      enddo
      do i = 1, size(uu%var)
        if (.not. uu%var(i)%good_opt) cycle
        iv = uu%var(i)%ix_dvar
        phase_ave = sum(uu%dm_dv(:, iv), mask = msk) / ix
        where (msk) uu%dm_dv(:, iv) = uu%dm_dv(:, iv) - phase_ave
      enddo
    endif

  enddo

!---------------------------------------------------------
! calculate dmerit itself

  u1%var(:)%dmerit = 0

  do k = 1, logic%u_num

    uu => super%u_(k)
    call ring_calc (uu)

    do i = 1, size(u1%var)
      if (u1%var(i)%useit) then
        iv = u1%var(i)%ix_dvar
        do j = 1, size(uu%data)
          if (uu%data(j)%useit_opt) then
            im = uu%data(j)%ix_dmeas
            u1%var(i)%dmerit = u1%var(i)%dmerit + &
                   2 * uu%data(j)%weight * uu%dm_dv(im, iv) * uu%data(j)%delta
          endif
        enddo
        if (k == 1) u1%var(i)%dmerit = u1%var(i)%dmerit + logic%u_num * &
                                       2 * u1%var(i)%weight * u1%var(i)%delta
      endif
    enddo

  enddo

end subroutine