submodule (mod_boundary_manager:smod_natural_boundaries) smod_natural_bounds_resistive
  implicit none

contains

  module procedure add_natural_resistive_terms
    real(dp)  :: eps, deps
    real(dp)  :: eta
    real(dp)  :: B02, dB02, drB02
    real(dp)  :: B03, dB03
    real(dp) :: gamma_1
    type(matrix_elements_t) :: elements

    if (.not. settings%physics%resistivity%is_enabled()) return

    gamma_1 = settings%physics%get_gamma_1()

    eps = grid%get_eps(x)
    deps = grid%get_deps()
    eta = physics%resistivity%eta(x)
    B02 = background%magnetic%B02(x)
    dB02 = background%magnetic%dB02(x)
    B03 = background%magnetic%B03(x)
    dB03 = background%magnetic%dB03(x)

    drB02 = deps * B02 + eps * dB02
    elements = new_matrix_elements(state_vector=settings%get_state_vector())

    ! ==================== Quadratic * Quadratic ====================
    call elements%add( &
      2.0d0 * ic * gamma_1 * eta * (k3 * drB02 - k2 * dB03), &
      "T", &
      "a1", &
      spline1=h_quad, &
      spline2=h_quad &
    )
    ! ==================== Quadratic * dCubic ====================
    call elements%add( &
      2.0d0 * ic * gamma_1 * eta * dB03, "T", "a2", spline1=h_quad, spline2=dh_cubic &
    )
    call elements%add( &
      -2.0d0 * ic * gamma_1 * eta * drB02, "T", "a3", spline1=h_quad, spline2=dh_cubic &
    )
    ! ==================== Cubic * Quadratic ====================
    call elements%add(-ic * eta * k2, "a2", "a1", spline1=h_cubic, spline2=h_quad)
    call elements%add(-ic * eta * eps * k3, "a3", "a1", spline1=h_cubic, spline2=h_quad)
    ! ==================== Cubic * dCubic ====================
    call elements%add(ic * eta, "a2", "a2", spline1=h_cubic, spline2=dh_cubic)
    call elements%add(ic * eta * eps, "a3", "a3", spline1=h_cubic, spline2=dh_cubic)

    call add_to_quadblock(quadblock, elements, weight, settings%dims)
    call elements%delete()
  end procedure add_natural_resistive_terms

end submodule smod_natural_bounds_resistive
