Refactoring: Refactor Broyden mixing such that the routine doesn't contain if(kerker) statements for every term

Refactor mix_inc.F90 to have cleaner top-level abstractions. For example:

Example 1.

  if (this%kerker) then
    call X(kerker_preconditioner)(this%der, d1, d2, d3, coeff, this%kerker_factor, f, precond_residual)
    vnew(1:d1, 1:d2, 1:d3) = vin(1:d1, 1:d2, 1:d3) + precond_residual(1:d1, 1:d2, 1:d3)
  else
    vnew(1:d1, 1:d2, 1:d3) = vin(1:d1, 1:d2, 1:d3) + coeff * f(1:d1, 1:d2, 1:d3)
  end if

goes to:

call X(precondition_residual)(this%der, d1, d2, d3, coeff, this%kerker_factor, f, residual)
vnew(1:d1, 1:d2, 1:d3) = vin(1:d1, 1:d2, 1:d3) + residual(1:d1, 1:d2, 1:d3)

where

subroutine precondition_residual(args, ..., residual)
  if (this%kerker) then
    call X(kerker_preconditioner)(this%der, d1, d2, d3, coeff, this%kerker_factor, f, residual)
  else
    residual = coeff * f(1:d1, 1:d2, 1:d3)
  end if
end subroutine

Example 2

  ! other terms
  if (this%kerker) then
    do i = 1, iter_used
      gamma = ww*sum(beta(:, i)*work(:))
      do l = 1, d3
        do k = 1, d2
          call X(kerker_preconditioner_vector)(this%der, coeff, this%kerker_factor, &
            df(1:d1, k, l, i), precond_residual(1:d1, k, l))
          vnew(1:d1, k, l) = vnew(1:d1, k, l) - ww*gamma*(precond_residual(1:d1, k, l) + dv(1:d1, k, l, i))
        end do
      end do
    end do
  else
    do i = 1, iter_used
      gamma = ww*sum(beta(:, i)*work(:))
      vnew(1:d1, 1:d2, 1:d3) = vnew(1:d1, 1:d2, 1:d3) - ww*gamma*(coeff*df(1:d1, 1:d2, 1:d3, i) + dv(1:d1, 1:d2, 1:d3, i))
    end do
  end if

which would become

    do i = 1, iter_used
      call X(precondition_residual_other_terms) (args, ...,beta, work, i, name_me)
      vnew(1:d1, 1:d2, 1:d3) = vnew(1:d1, 1:d2, 1:d3) - ww * name_me(1:d1, 1:d2, 1:d3)
    end do

with

subroutine X(precondition_residual_other_terms) (args,..., beta, work, i, term)

gamma = ww*sum(beta(:, i)*work(:))

if(kerker)then
      do l = 1, d3
        do k = 1, d2
          call X(kerker_preconditioner_vector)(this%der, coeff, this%kerker_factor, &
            df(1:d1, k, l, i), precond_residual(1:d1, k, l))
          term(1:d1, k, l) = gamma * (precond_residual(1:d1, k, l) + dv(1:d1, k, l, i))
        end do
else
      term(1:d1, 1:d2, 1:d3) = gamma *( coeff * df(1:d1, 1:d2, 1:d3, i) + dv(1:d1, 1:d2, 1:d3, i))
endif

end subroutine

Or, one could write two versions of precondition_residual_other_terms with the same signature, then assign the appropriate one (kerker or no kerker) with a pointer.