amun-code/sources/sources.F90

!!******************************************************************************
!!
!!  This file is part of the AMUN source code, a program to perform
!!  Newtonian or relativistic magnetohydrodynamical simulations on uniform or
!!  adaptive mesh.
!!
!!  Copyright (C) 2008-2024 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License as published by
!!  the Free Software Foundation, either version 3 of the License, or
!!  (at your option) any later version.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!
!!******************************************************************************
!!
!! module: SOURCES
!!
!!  This modules adds source terms.
!!
!!******************************************************************************
!
module sources

  implicit none

! interfaces for the extra source terms subroutine
!
  abstract interface
    subroutine update_extra_sources_iface(t, dt, pdata)
      use blocks, only : block_data
      implicit none
      real(kind=8)                    , intent(in)    :: t, dt
      type(block_data), pointer       , intent(inout) :: pdata
    end subroutine
  end interface

! pointer to the extra source terms subroutine
!
  procedure(update_extra_sources_iface), pointer, save :: &
                                                  update_extra_sources => null()

! GLM-MHD source terms type (1 - EGLM, 2 - HEGLM)
!
  integer          , save :: glm_type    = 0
  character(len=32), save :: glm_name = "none"

! viscosity coefficient
!
  real(kind=8)     , save :: viscosity   = 0.0d+00

! resistivity coefficient
!
  real(kind=8)     , save :: resistivity = 0.0d+00
  real(kind=8)     , save :: anomalous   = 0.0d+00
  real(kind=8)     , save :: jcrit       = 1.0d+00

  private

  public :: initialize_sources, finalize_sources, print_sources
  public :: update_sources, update_extra_sources
  public :: viscosity, resistivity

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
  contains
!
!===============================================================================
!!
!!***  PUBLIC SUBROUTINES  *****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine INITIALIZE_SOURCES:
! -----------------------------
!
!   Subroutine initializes module SOURCES.
!
!   Arguments:
!
!     verbose - a logical flag turning the information printing;
!     status  - an integer flag for error return value;
!
!===============================================================================
!
  subroutine initialize_sources(verbose, status)

    use parameters, only : get_parameter

    implicit none

    logical, intent(in)  :: verbose
    integer, intent(out) :: status

    character(len=8) :: tglm = "none"

!-------------------------------------------------------------------------------
!
    status = 0

! get the type of the GLM source terms
!
    call get_parameter("glm_source_terms", tglm)

! set the glm_type variable to correct value
!
    select case(trim(tglm))
    case("eglm", "EGLM")
      glm_type = 1
      glm_name = "EGLM"
    case("heglm", "HEGLM")
      glm_type = 2
      glm_name = "HEGLM"
    case("kepes", "KEPES")
      glm_type = 3
      glm_name = "KEPES"
    case default
      glm_type = 0
      glm_name = "none"
    end select

! get viscosity coefficient
!
      call get_parameter("viscosity", viscosity)

      if (viscosity < 0.0d+00) then
        if (verbose) then
          write(*,*)
          write(*,"(1x,a)") "ERROR!"
          write(*,"(1x,a)") "Negative viscosity coefficient!"
        end if
        status = 1
      end if

! get resistivity coefficients
!
      call get_parameter("resistivity", resistivity)

      if (resistivity < 0.0d+00) then
        if (verbose) then
          write(*,*)
          write(*,"(1x,a)") "ERROR!"
          write(*,"(1x,a)") "Negative resistivity coefficient!"
        end if
        status = 1
      end if

      call get_parameter("anomalous", anomalous)

      if (anomalous < 0.0d+00) then
        if (verbose) then
          write(*,*)
          write(*,"(1x,a)") "ERROR!"
          write(*,"(1x,a)") "Negative anomalous resistivity coefficient!"
        end if
        status = 1
      end if

      call get_parameter("jcrit", jcrit)

      if (jcrit <= 0.0d+00) then
        if (verbose) then
          write(*,*)
          write(*,"(1x,a)") "ERROR!"
          write(*,"(1x,a)") "Non-positive critical current density coefficient!"
        end if
        status = 1
      end if

!-------------------------------------------------------------------------------
!
  end subroutine initialize_sources
!
!===============================================================================
!
! subroutine FINALIZE_SOURCES:
! ---------------------------
!
!   Subroutine releases memory used by the module.
!
!   Arguments:
!
!     status - an integer flag for error return value;
!
!===============================================================================
!
  subroutine finalize_sources(status)

    implicit none

    integer, intent(out) :: status

!-------------------------------------------------------------------------------
!
    status = 0

!-------------------------------------------------------------------------------
!
  end subroutine finalize_sources
!
!===============================================================================
!
! subroutine PRINT_SOURCES:
! ------------------------
!
!   Subroutine prints module parameters.
!
!   Arguments:
!
!     verbose - a logical flag turning the information printing;
!
!===============================================================================
!
  subroutine print_sources(verbose)

    use equations, only : magnetized
    use helpers  , only : print_section, print_parameter

    implicit none

    logical, intent(in) :: verbose

!-------------------------------------------------------------------------------
!
    if (verbose) then

      call print_section(verbose, "Source terms")
      call print_parameter(verbose, "viscosity"       , viscosity  )
      if (magnetized) then
        call print_parameter(verbose, "resistivity"     , resistivity)
        if (anomalous > 0.0d+00) then
          call print_parameter(verbose, "anomalous"       , anomalous  )
          call print_parameter(verbose, "jcrit"           , jcrit      )
        end if
        call print_parameter(verbose, "glm source terms", glm_name   )
      end if

    end if

!-------------------------------------------------------------------------------
!
  end subroutine print_sources
!
!===============================================================================
!
! subroutine UPDATE_SOURCES:
! -------------------------
!
!   Subroutine adds source terms to the array of increments dU.
!
!   Arguments:
!
!     t, dt - the time and time increment;
!     pdata - the pointer to a data block;
!
!===============================================================================
!
  subroutine update_sources(t, dt, pdata)

    use blocks     , only : block_data
    use coordinates, only : nn => bcells
    use coordinates, only : ax, adx, ay, ady, adz
#if NDIMS == 3
    use coordinates, only : az
#endif /* NDIMS == 3 */
    use equations  , only : magnetized
    use equations  , only : inx, iny, inz
    use equations  , only : idn, ivx, ivy, ivz, imx, imy, imz, ien
    use equations  , only : ibx, iby, ibz, ibp
    use gravity    , only : gravity_enabled, gravitational_acceleration
    use helpers    , only : print_message
    use operators  , only : divergence, gradient, laplace, curl
    use workspace  , only : resize_workspace, work, work_in_use

    implicit none

    real(kind=8)                    , intent(in)    :: t, dt
    type(block_data), pointer       , intent(inout) :: pdata

    logical, save :: first = .true.

    integer       :: status
    integer       :: i, j, k
    real(kind=8)  :: fc, gc
    real(kind=8)  :: gx, gy, gz
    real(kind=8)  :: dvxdx, dvxdy, dvxdz, divv
    real(kind=8)  :: dvydx, dvydy, dvydz
    real(kind=8)  :: dvzdx, dvzdy, dvzdz

    integer, save :: nt
!$  integer :: omp_get_thread_num

    real(kind=8), dimension(3)  :: ga, dh
    real(kind=8), dimension(nn) :: x, y
#if NDIMS == 3
    real(kind=8), dimension(nn) :: z
#else /* NDIMS == 3 */
    real(kind=8), dimension( 1) :: z
#endif /* NDIMS == 3 */
    real(kind=8), dimension(:,:,:)    , pointer, save :: db
    real(kind=8), dimension(:,:,:,:,:), pointer, save :: tmp
!$omp threadprivate(first, nt, db, tmp)

    character(len=*), parameter :: loc = 'SOURCES:update_sources()'

!-------------------------------------------------------------------------------
!
    nt = 0
!$  nt = omp_get_thread_num()
    if (first) then
      i = nn**NDIMS
      j = 10 * i

      call resize_workspace(j, status)
      if (status /= 0) then
        call print_message(loc, "Could not resize the workspace!")
        go to 100
      end if

#if NDIMS == 3
        db(1:nn,1:nn,1:nn)          => work(  1:i,nt)
        tmp(1:3,1:3,1:nn,1:nn,1:nn) => work(i+1:j,nt)
#else /* NDIMS == 3 */
        db(1:nn,1:nn,1: 1)          => work(  1:i,nt)
        tmp(1:3,1:3,1:nn,1:nn,1: 1) => work(i+1:j,nt)
#endif /* NDIMS == 3 */

      first = .false.
    end if

#if NDIMS == 2
    k = 1
#endif /* NDIMS == 2 */

! proceed only if the gravitational term is enabled
!
    if (gravity_enabled) then

! prepare block coordinates
!
      x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
      y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
#if NDIMS == 3
      z(:) = pdata%meta%zmin + az(pdata%meta%level,:)
#endif /* NDIMS == 3 */

! iterate over all positions in the YZ plane
!
#if NDIMS == 3
      do k = 1, nn
#endif /* NDIMS == 3 */
        do j = 1, nn
          do i = 1, nn

! get gravitational acceleration components
!
            call gravitational_acceleration(t, dt, x(i), y(j), z(k), ga(1:3))

! calculate the gravitational source terms
!
            gx = pdata%q(idn,i,j,k) * ga(1)
            gy = pdata%q(idn,i,j,k) * ga(2)
#if NDIMS == 3
            gz = pdata%q(idn,i,j,k) * ga(3)
#endif /* NDIMS == 3 */

! add source terms to momentum equations
!
            pdata%du(imx,i,j,k) = pdata%du(imx,i,j,k) + gx
            pdata%du(imy,i,j,k) = pdata%du(imy,i,j,k) + gy
#if NDIMS == 3
            pdata%du(imz,i,j,k) = pdata%du(imz,i,j,k) + gz
#endif /* NDIMS == 3 */

! add source terms to total energy equation
!
            if (ien > 0) then

#if NDIMS == 2
              pdata%du(ien,i,j,k) = pdata%du(ien,i,j,k) &
                                                    + gx * pdata%q(ivx,i,j,k) &
                                                    + gy * pdata%q(ivy,i,j,k)
#endif /* NDIMS == 2 */
#if NDIMS == 3
              pdata%du(ien,i,j,k) = pdata%du(ien,i,j,k) &
                                                    + gx * pdata%q(ivx,i,j,k) &
                                                    + gy * pdata%q(ivy,i,j,k) &
                                                    + gz * pdata%q(ivz,i,j,k)
#endif /* NDIMS == 3 */

            end if

          end do ! i = 1, nn
        end do ! j = 1, nn
#if NDIMS == 3
      end do ! k = 1, nn
#endif /* NDIMS == 3 */

    end if ! gravity enabled

! proceed only if the viscosity coefficient is not zero
!
    if (viscosity > 0.0d+00) then

      if (work_in_use(nt)) &
        call print_message(loc, "Workspace is being used right now! " // &
                                "Corruptions can occur!")

      work_in_use(nt) = .true.

! prepare coordinate increments
!
      dh(1) = adx(pdata%meta%level)
      dh(2) = ady(pdata%meta%level)
      dh(3) = adz(pdata%meta%level)

! calculate the velocity Jacobian
!
      call gradient(dh(:), pdata%q(ivx,:,:,:), tmp(inx,inx:inz,:,:,:))
      call gradient(dh(:), pdata%q(ivy,:,:,:), tmp(iny,inx:inz,:,:,:))
      call gradient(dh(:), pdata%q(ivz,:,:,:), tmp(inz,inx:inz,:,:,:))

! iterate over all cells
!
#if NDIMS == 3
      do k = 1, nn
#endif /* NDIMS == 3 */
        do j = 1, nn
          do i = 1, nn

! prepare the νρ factor
!
            gc    = viscosity * pdata%q(idn,i,j,k)
            fc    = 2.0d+00 * gc

! get the velocity Jacobian elements
!
            dvxdx = tmp(inx,inx,i,j,k)
            dvxdy = tmp(inx,iny,i,j,k)
            dvxdz = tmp(inx,inz,i,j,k)
            dvydx = tmp(iny,inx,i,j,k)
            dvydy = tmp(iny,iny,i,j,k)
            dvydz = tmp(iny,inz,i,j,k)
            dvzdx = tmp(inz,inx,i,j,k)
            dvzdy = tmp(inz,iny,i,j,k)
            dvzdz = tmp(inz,inz,i,j,k)
            divv  = (dvxdx + dvydy + dvzdz) / 3.0d+00

! calculate elements of the viscous stress tensor
!
            tmp(inx,inx,i,j,k) = fc * (dvxdx - divv)
            tmp(iny,iny,i,j,k) = fc * (dvydy - divv)
            tmp(inz,inz,i,j,k) = fc * (dvzdz - divv)
            tmp(inx,iny,i,j,k) = gc * (dvxdy + dvydx)
            tmp(inx,inz,i,j,k) = gc * (dvxdz + dvzdx)
            tmp(iny,inz,i,j,k) = gc * (dvydz + dvzdy)
            tmp(iny,inx,i,j,k) = tmp(inx,iny,i,j,k)
            tmp(inz,inx,i,j,k) = tmp(inx,inz,i,j,k)
            tmp(inz,iny,i,j,k) = tmp(iny,inz,i,j,k)

          end do ! i = 1, nn
        end do ! j = 1, nn
#if NDIMS == 3
      end do ! k = 1, nn
#endif /* NDIMS == 3 */

! calculate the divergence of the first tensor row
!
      call divergence(dh(:), tmp(inx,inx:inz,:,:,:), db(:,:,:))

! add viscous source terms to the X momentum equation
!
      pdata%du(imx,:,:,:) = pdata%du(imx,:,:,:) + db(:,:,:)

! calculate the divergence of the second tensor row
!
      call divergence(dh(:), tmp(iny,inx:inz,:,:,:), db(:,:,:))

! add viscous source terms to the Y momentum equation
!
      pdata%du(imy,:,:,:) = pdata%du(imy,:,:,:) + db(:,:,:)

! calculate the divergence of the third tensor row
!
      call divergence(dh(:), tmp(inz,inx:inz,:,:,:), db(:,:,:))

! add viscous source terms to the Z momentum equation
!
      pdata%du(imz,:,:,:) = pdata%du(imz,:,:,:) + db(:,:,:)

! add viscous source term to total energy equation
!
      if (ien > 0) then

! iterate over all cells
!
#if NDIMS == 3
        do k = 1, nn
#endif /* NDIMS == 3 */
          do j = 1, nn
            do i = 1, nn

! calculate scalar product of v and viscous stress tensor τ
!
              gx = pdata%q(ivx,i,j,k) * tmp(inx,inx,i,j,k) &
                 + pdata%q(ivy,i,j,k) * tmp(inx,iny,i,j,k) &
                 + pdata%q(ivz,i,j,k) * tmp(inx,inz,i,j,k)
              gy = pdata%q(ivx,i,j,k) * tmp(iny,inx,i,j,k) &
                 + pdata%q(ivy,i,j,k) * tmp(iny,iny,i,j,k) &
                 + pdata%q(ivz,i,j,k) * tmp(iny,inz,i,j,k)
              gz = pdata%q(ivx,i,j,k) * tmp(inz,inx,i,j,k) &
                 + pdata%q(ivy,i,j,k) * tmp(inz,iny,i,j,k) &
                 + pdata%q(ivz,i,j,k) * tmp(inz,inz,i,j,k)

! update (v.τ), use the first row of the tensor tmp
!
              tmp(inx,inx,i,j,k) = gx
              tmp(inx,iny,i,j,k) = gy
              tmp(inx,inz,i,j,k) = gz

            end do ! i = 1, nn
          end do ! j = 1, nn
#if NDIMS == 3
        end do ! k = 1, nn
#endif /* NDIMS == 3 */

! calculate the divergence of (v.τ)
!
        call divergence(dh(:), tmp(inx,inx:inz,:,:,:), db(:,:,:))

! update the energy increment
!
        pdata%du(ien,:,:,:) = pdata%du(ien,:,:,:) + db(:,:,:)

      end if ! ien > 0

      work_in_use(nt) = .false.

    end if ! viscosity is not zero

!=== add magnetic field related source terms ===
!
    if (magnetized) then

      if (work_in_use(nt)) &
        call print_message(loc, "Workspace is being used right now! " // &
                                "Corruptions can occur!")

      work_in_use(nt) = .true.

! prepare coordinate increments
!
      dh(1) = adx(pdata%meta%level)
      dh(2) = ady(pdata%meta%level)
      dh(3) = adz(pdata%meta%level)

! add the EGLM-MHD and KEPES source terms
!
      if (glm_type > 0) then

! calculate the magnetic field divergence
!
        call divergence(dh(:), pdata%q(ibx:ibz,:,:,:), db(:,:,:))

! update the momentum component increments, i.e.
!     d/dt (ρv) + ∇.F = - (∇.B)B
!
        pdata%du(imx,:,:,:) = pdata%du(imx,:,:,:) &
                                               - db(:,:,:) * pdata%q(ibx,:,:,:)
        pdata%du(imy,:,:,:) = pdata%du(imy,:,:,:) &
                                               - db(:,:,:) * pdata%q(iby,:,:,:)
        pdata%du(imz,:,:,:) = pdata%du(imz,:,:,:) &
                                               - db(:,:,:) * pdata%q(ibz,:,:,:)

! add the HEGLM-MHD and KEPES source terms
!
        if (glm_type > 1) then

! update magnetic field component increments, i.e.
!     d/dt B + ∇.F = - (∇.B)v
!
          pdata%du(ibx,:,:,:) = pdata%du(ibx,:,:,:) &
                                               - db(:,:,:) * pdata%q(ivx,:,:,:)
          pdata%du(iby,:,:,:) = pdata%du(iby,:,:,:) &
                                               - db(:,:,:) * pdata%q(ivy,:,:,:)
          pdata%du(ibz,:,:,:) = pdata%du(ibz,:,:,:) &
                                               - db(:,:,:) * pdata%q(ivz,:,:,:)

! update the energy equation
!
          if (ien > 0) then

! calculate scalar product of velocity and magnetic field
!
            tmp(inx,inx,:,:,:) = sum(pdata%q(ivx:ivz,:,:,:) &
                                                   * pdata%q(ibx:ibz,:,:,:), 1)

! add the divergence potential source term to the energy equation, i.e.
!     d/dt E + ∇.F = - (∇.B) (v.B)
!
            pdata%du(ien,:,:,:) = pdata%du(ien,:,:,:) &
                                               - db(:,:,:) * tmp(inx,inx,:,:,:)

          end if ! ien > 0

        end if ! glm_type > 1

        if  (ibp > 0) then
          if  (glm_type == 3) then

! calculate the divergence of velocity
!
            call divergence(dh(:), pdata%q(ivx:ivz,:,:,:), db(:,:,:))

! add the divergence potential source term to the energy equation, i.e.
!     d/dt ψ + ∇.F = ½ψ(∇.v)
!
            pdata%du(ibp,:,:,:) = pdata%du(ibp,:,:,:) &
                          + 5.0d-01 * pdata%q(ibp,:,:,:) * db(:,:,:)

          else if (ien > 0) then

! calculate the gradient of the divergence correcting field
!
            call gradient(dh(:), pdata%q(ibp,:,:,:), tmp(inx:inz,inx,:,:,:))

            db(:,:,:) = sum(pdata%q(ivx:ivz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)

! update the divergence correcting field
!     d/dt ψ + ∇.F = - (v.∇)ψ
!
            pdata%du(ibp,:,:,:) = pdata%du(ibp,:,:,:) - db(:,:,:)

! add the divergence potential source term to the energy equation, i.e.
!     d/dt E + ∇.F = - B.(∇ψ)
!
            pdata%du(ien,:,:,:) = pdata%du(ien,:,:,:) &
                     - sum(pdata%q(ibx:ibz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)

          end if ! glm == 3
        end if ! ibp
      end if ! glm_type > 0

! if anomalous resistivity is enabled
!
      if (anomalous > 0.0d+00) then

! calculate current density (J = ∇xB)
!
        call curl(dh(:), pdata%q(ibx:ibz,:,:,:), tmp(inx:inz,inx,:,:,:))

! calculate the normalized absolute value of current density (|J|/Jcrit)
!
        tmp(inx,iny,:,:,:) = sqrt(sum(tmp(inx:inz,inx,:,:,:)**2, 1)) / jcrit

! calculate the local resistivity [ηu + ηa (|J|/Jcrit - 1) H(|J|/Jcrit)]
!
        tmp(iny,iny,:,:,:) = resistivity + &
                      anomalous * max(0.0d+00, (tmp(inx,iny,:,:,:) - 1.0d+00))

! multiply the current density vector by the local resistivity (ηJ)
!
        tmp(inx,inz,:,:,:) = tmp(iny,iny,:,:,:) * tmp(inx,inx,:,:,:)
        tmp(iny,inz,:,:,:) = tmp(iny,iny,:,:,:) * tmp(iny,inx,:,:,:)
        tmp(inz,inz,:,:,:) = tmp(iny,iny,:,:,:) * tmp(inz,inx,:,:,:)

! calculate the curl of (ηJ)
!
        call curl(dh(:), tmp(inx:inz,inz,:,:,:), tmp(inx:inz,iny,:,:,:))

! update magnetic field component increments
!
        pdata%du(ibx,:,:,:) = pdata%du(ibx,:,:,:) - tmp(inx,iny,:,:,:)
        pdata%du(iby,:,:,:) = pdata%du(iby,:,:,:) - tmp(iny,iny,:,:,:)
        pdata%du(ibz,:,:,:) = pdata%du(ibz,:,:,:) - tmp(inz,iny,:,:,:)

! update energy equation
!
        if (ien > 0) then

! calculate the vector product Bx(η ∇xB)
!
          tmp(inx,iny,:,:,:) = pdata%q(iby,:,:,:) * tmp(inz,inz,:,:,:) &
                             - pdata%q(ibz,:,:,:) * tmp(iny,inz,:,:,:)
          tmp(iny,iny,:,:,:) = pdata%q(ibz,:,:,:) * tmp(inx,inz,:,:,:) &
                             - pdata%q(ibx,:,:,:) * tmp(inz,inz,:,:,:)
          tmp(inz,iny,:,:,:) = pdata%q(ibx,:,:,:) * tmp(iny,inz,:,:,:) &
                             - pdata%q(iby,:,:,:) * tmp(inx,inz,:,:,:)

! calculate the divergence ∇.[Bx(η ∇xB)]
!
          call divergence(dh(:), tmp(inx:inz,iny,:,:,:), db(:,:,:))

! add the second resistive source term to the energy equation, i.e.
!     d/dt E + ∇.F = η J²
!
          pdata%du(ien,:,:,:) = pdata%du(ien,:,:,:) + db(:,:,:)

        end if ! energy equation present

      else if (resistivity > 0.0d+00) then

! calculate the Laplace operator of B, i.e. Δ(B)
!
        call laplace(dh(:), pdata%q(ibx,:,:,:), tmp(inx,inx,:,:,:))
        call laplace(dh(:), pdata%q(iby,:,:,:), tmp(inx,iny,:,:,:))
        call laplace(dh(:), pdata%q(ibz,:,:,:), tmp(inx,inz,:,:,:))

! multiply by the resistivity coefficient
!
        tmp(iny,inx:inz,:,:,:) = resistivity * tmp(inx,inx:inz,:,:,:)

! update magnetic field component increments
!
        pdata%du(ibx,:,:,:) = pdata%du(ibx,:,:,:) + tmp(iny,inx,:,:,:)
        pdata%du(iby,:,:,:) = pdata%du(iby,:,:,:) + tmp(iny,iny,:,:,:)
        pdata%du(ibz,:,:,:) = pdata%du(ibz,:,:,:) + tmp(iny,inz,:,:,:)

! update energy equation
!
        if (ien > 0) then

! add the first resistive source term to the energy equation, i.e.
!     d/dt E + ∇.F = η B.[Δ(B)]
!
          pdata%du(ien,:,:,:) = pdata%du(ien,:,:,:) &
                             + (pdata%q(ibx,:,:,:) * tmp(iny,inx,:,:,:) &
                             +  pdata%q(iby,:,:,:) * tmp(iny,iny,:,:,:) &
                             +  pdata%q(ibz,:,:,:) * tmp(iny,inz,:,:,:))

! calculate current density J = ∇xB
!
          call curl(dh(:), pdata%q(ibx:ibz,:,:,:), tmp(inz,inx:inz,:,:,:))

! calculate J²
!
          db(:,:,:) = tmp(inz,inx,:,:,:)**2 &
                    + tmp(inz,iny,:,:,:)**2 &
                    + tmp(inz,inz,:,:,:)**2

! add the second resistive source term to the energy equation, i.e.
!     d/dt E + ∇.F = η J²
!
          pdata%du(ien,:,:,:) = pdata%du(ien,:,:,:) + resistivity * db(:,:,:)

        end if ! energy equation present

      end if ! resistivity is not zero

      work_in_use(nt) = .false.

    end if ! magnetized

! add extra source terms
!
    if (associated(update_extra_sources)) &
                                        call update_extra_sources(t, dt, pdata)

    100 continue

!-------------------------------------------------------------------------------
!
  end subroutine update_sources

!===============================================================================
!
end module sources