amun-code/sources/user_problem.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) 2017-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: USER_PROBLEM
!!
!!  This module provides subroutines to setup custom problem.
!!
!!*****************************************************************************
!
module user_problem

  implicit none

  integer            :: profile = 1

  real(kind=8), save :: beta = 2.00d+00
  real(kind=8), save :: zeta = 0.00d+00
  real(kind=8), save :: dens = 1.00d+00
  real(kind=8), save :: bamp = 1.00d+00
  real(kind=8), save :: bgui = 0.00d+00
  real(kind=8), save :: dlta = 1.00d-02
  real(kind=8), save :: pres = 5.00d-01
  real(kind=8), save :: pmag = 5.00d-01
  real(kind=8), save :: ptot = 1.00d+00
  real(kind=8), save :: valf = 1.00d+00
  real(kind=8), save :: lund = 1.00d+00
  real(kind=8), save :: eta  = 0.00d+00

  integer(kind=4), save :: runit = 33

  logical, save :: resistive  = .false.

  private
  public :: initialize_user_problem, finalize_user_problem
  public :: user_statistics

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
  contains
!
!===============================================================================
!
! subroutine INITIALIZE_USER_PROBLEM:
! ----------------------------------
!
!   Subroutine initializes user problem. It could read problem parameters which
!   are used in all subroutines defining this specific problem.
!
!   Arguments:
!
!     problem - the problem name
!     rcount  - the run count for restarted jobs
!     verbose - a logical flag turning the information printing;
!     status  - an integer flag for error return value;
!
!===============================================================================
!
  subroutine initialize_user_problem(problem, rcount, verbose, status)

    use boundaries , only : custom_boundary_x, custom_boundary_y
    use equations  , only : adiabatic_index, csnd, csnd2, ipr
    use helpers    , only : print_section, print_parameter, print_message
    use mesh       , only : setup_problem
    use parameters , only : get_parameter

    implicit none

    character(len=64), intent(in)  :: problem
    integer          , intent(in)  :: rcount
    logical          , intent(in)  :: verbose
    integer          , intent(out) :: status

    character(len=64) :: tmp, append = "off", fname
    logical           :: flag

    character(len=*), parameter :: &
                               loc = 'USER_PROBLEM::initialize_user_problem()'
!-------------------------------------------------------------------------------
!
    call get_parameter("profile", tmp)
    select case(tmp)
      case('bhattacharjee', 'sincos', 'trigonometric')
        profile = 2
      case default
        profile = 1
    end select

    call get_parameter("dens", dens)
    call get_parameter("bamp", bamp)
    call get_parameter("bgui", bgui)
    call get_parameter("beta", beta)
    if (beta <= 0.0d+00) then
      if (verbose) &
        call print_message(loc, "Plasma-beta must be positive (beta > 0.0)!")
      status = 1
      return
    end if
    call get_parameter("zeta", zeta)
    if (zeta < 0.0d+00 .or. zeta > 1.0d+00) then
      if (verbose) &
        call print_message(loc, "Parameter 'zeta' must be between 0.0 and 1.0!")
      status = 1
      return
    end if
    if (dens <= 0.0d+00) then
      if (verbose) &
        call print_message(loc, "Density must be positive (dens > 0.0)!")
      status = 1
      return
    end if
    call get_parameter("resistivity", eta)
    if (eta < 0.0d+00) then
      if (verbose) &
        call print_message(loc, "Resistivity cannot be negative!")
      status = 1
      return
    else
      resistive = .true.
      dlta = sqrt(eta)
    end if
    call get_parameter("delta", dlta)
    if (dlta <= 0.0d+00) then
      if (verbose) &
        call print_message(loc, &
                      "The initial thickness must be positive (delta > 0.0)!")
      status = 1
      return
    end if

    pmag = 0.5d+00 * (bamp**2 + bgui**2)
    pres = beta * pmag
    ptot = pres + pmag
    if (ipr > 0) then
      csnd2 = adiabatic_index * pres / dens
    else
      csnd2 = pres / dens
    end if
    csnd  = sqrt(csnd2)
    valf  = bamp / sqrt(dens)
    lund  = valf / max(tiny(eta), eta)

! determine if to append or create another file
!
    call get_parameter("statistics_append", append)
    select case(trim(append))
    case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
      write(fname, "('magnetic-flux.dat')")
      inquire(file = fname, exist = flag)
    case default
      write(fname, "('magnetic-flux_',i2.2,'.dat')") rcount
      flag = .false.
    end select

! check if the file exists; if not, create a new one, otherwise move to the end
!
    if (flag .and. rcount > 1) then
#ifdef __INTEL_COMPILER
      open(newunit = runit, file = fname, form = 'formatted', status = 'old' &
                            , position = 'append', buffered = 'yes')
#else /* __INTEL_COMPILER */
      open(newunit = runit, file = fname, form = 'formatted', status = 'old' &
                            , position = 'append')
#endif /* __INTEL_COMPILER */
      write(runit,"('#')")
    else
#ifdef __INTEL_COMPILER
      open(newunit = runit, file = fname, form = 'formatted' &
                            , status = 'replace', buffered = 'yes')
#else /* __INTEL_COMPILER */
      open(newunit = runit, file = fname, form = 'formatted' &
                            , status = 'replace')
#endif /* __INTEL_COMPILER */
    end if

! write the integral file header
!
    write(runit,'("#",a24,8a25)') 'time', '|Bx| int'  ,               &
                                          '|Bx| y-adv', '|Bx| z-adv', &
                                          '|Bx| y-shr', '|Bx| z-shr', &
                                          '|Bx| y-dif', '|Bx| z-dif', &
                                          '|Bx| Psi'
    write(runit,"('#')")

    call print_section(verbose, "Parameters (* - set, otherwise calculated)")
    select case(profile)
    case(2)
      call print_parameter(verbose, '(*) flux tubes profile', &
                                                  'Bhattacharjee et al. (2009)')
    case default
      call print_parameter(verbose, '(*) flux tubes profile', &
                                                  'Kowal et al. (2023)')
    end select
    call print_parameter(verbose, '(*) beta (plasma-beta)'   , beta)
    call print_parameter(verbose, '(*) zeta'                 , zeta)
    call print_parameter(verbose, '(*) dens'                 , dens)
    call print_parameter(verbose, '(*) bamp (amplitude)'     , bamp)
    call print_parameter(verbose, '(*) bgui (guide field)'   , bgui)
    call print_parameter(verbose, '(*) delta (thickness)'    , dlta)
    call print_parameter(verbose, '( ) pres (thermal pres.)' , pres)
    call print_parameter(verbose, '( ) pmag (magnetic pres.)', pmag)
    call print_parameter(verbose, '( ) csnd (sound speed)'   , csnd)
    call print_parameter(verbose, '( ) Valf (Alfven speed)'  , valf)
    if (resistive) &
      call print_parameter(verbose, '( ) S (Lundquist number)', lund)

    setup_problem     => setup_user_problem
    custom_boundary_x => user_boundary_x
    custom_boundary_y => user_boundary_y

    status = 0

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

    implicit none

    integer, intent(out) :: status

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

!-------------------------------------------------------------------------------
!
  end subroutine finalize_user_problem
!
!===============================================================================
!
! subroutine SETUP_USER_PROBLEM:
! -----------------------------
!
!   Subroutine sets the initial conditions for the user specific problem.
!
!   Arguments:
!
!     pdata - pointer to the datablock structure of the currently initialized
!             block;
!
!===============================================================================
!
  subroutine setup_user_problem(pdata)

    use blocks     , only : block_data
    use constants  , only : pi, pi2
    use coordinates, only : nn => bcells, xmin, xmax
    use coordinates, only : ax, ay
    use equations  , only : prim2cons, csnd2
    use equations  , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp

    implicit none

    type(block_data), pointer, intent(inout) :: pdata

    integer                     :: i, j, k
    real(kind=8)                :: r1, r2, bx, by, bz
    real(kind=8), dimension(nn) :: x, y
    real(kind=8), dimension(nn) :: snx, csx, sny, csy, thy, ch2
    real(kind=8), dimension(nn,nn) :: az

!-------------------------------------------------------------------------------
!
    x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
    y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)

    select case(profile)

    case(2)
      snx(:) = sin(pi * x(:))
      csx(:) = cos(pi * x(:))
      sny(:) = sin(pi2 * y(:))
      csy(:) = cos(pi2 * y(:))
      thy(:) = tanh(y(:) / dlta)
      ch2(:) = cosh(y(:) / dlta)**2

      do j = 1, nn
        do i = 1, nn
          bx = bamp * csx(i) * (csy(j) * thy(j) + sny(j) / (pi2 * dlta * ch2(j)))
          by = bamp * 5.0d-01 * snx(i) * sny(j) * thy(j)
          bz = sqrt(max(0.0d+00, bgui**2 + zeta * (bamp**2 - bx**2 - by**2)))

          pdata%q(ibx,i,j,:) = bx
          pdata%q(iby,i,j,:) = by
          pdata%q(ibz,i,j,:) = bz
        end do
      end do
      pdata%q(ibp,:,:,:) = 0.0d+00

    case default
      do j = 1, nn
        do i = 1, nn
          r1 = 5.0d-01 - sqrt(x(i)**2 + ((2.0d+00 * y(j) - 5.0d-01))**2)
          r2 = 5.0d-01 - sqrt(x(i)**2 + ((2.0d+00 * y(j) + 5.0d-01))**2)
          az(i,j) = 2.5d-01 * (r1 * (tanh(r1 / dlta) + 1.0d+00) &
                            +  r2 * (tanh(r2 / dlta) + 1.0d+00))
        end do
      end do

      do j = 2, nn - 1
        do i = 2, nn - 1
          bx = (az(i,j+1) - az(i,j-1)) / (y(j+1) - y(j-1))
          by = (az(i-1,j) - az(i+1,j)) / (x(i+1) - x(i-1))
          bz = sqrt(max(0.0d+00, bgui**2 + zeta * (bamp**2 - bx**2 - by**2)))

          pdata%q(ibx,i,j,:) = bx
          pdata%q(iby,i,j,:) = by
          pdata%q(ibz,i,j,:) = bz
        end do
      end do
      pdata%q(ibp,:,:,:) = 0.0d+00

    end select

    if (ipr > 0) then
      pdata%q(idn,:,:,:) = dens
      pdata%q(ipr,:,:,:) =  ptot - 5.0d-01 * sum(pdata%q(ibx:ibz,:,:,:)**2,1)
    else
      pdata%q(idn,:,:,:) = (ptot - 5.0d-01 * sum(pdata%q(ibx:ibz,:,:,:)**2,1)) &
                         / csnd2
    end if

    pdata%q(ivx,:,:,:) = 0.0d+00
    pdata%q(ivy,:,:,:) = 0.0d+00
    pdata%q(ivz,:,:,:) = 0.0d+00

#if NDIMS == 3
    do k = 1, nn
#else /* NDIMS == 3 */
    do k = 1, 1
#endif /* NDIMS == 3 */
      do j = 1, nn
        call prim2cons(pdata%q(:,:,j,k), pdata%u(:,:,j,k), .true.)
      end do
    end do

!-------------------------------------------------------------------------------
!
  end subroutine setup_user_problem
!
!===============================================================================
!
! subroutine UPDATE_USER_SOURCES:
! ------------------------------
!
!   Subroutine adds the user defined source terms.
!
!   Arguments:
!
!     t, dt - the time and time increment;
!     pdata - the pointer to a data block;
!
!===============================================================================
!
  subroutine update_user_sources(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 update_user_sources
!
!===============================================================================
!
! subroutine UPDATE_USER_SHAPES:
! -----------------------------
!
!   Subroutine defines the regions updated by user.
!
!   Arguments:
!
!     pdata - pointer to the data block structure of the currently initialized
!             block;
!     time  - time at the moment of update;
!     dt    - time step since the last update;
!
!===============================================================================
!
  subroutine update_user_shapes(pdata, time, dt)

    use blocks, only : block_data

    implicit none

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

!-------------------------------------------------------------------------------
!

!-------------------------------------------------------------------------------
!
  end subroutine update_user_shapes
!
!===============================================================================
!
! subroutine USER_BOUNDARY_X:
! --------------------------
!
!   Subroutine updates ghost zones within the specific region along
!   the X direction.
!
!   Arguments:
!
!     is      - the block side along the X direction for the ghost zone update;
!     jl, ju  - the cell index limits for the Y direction;
!     kl, ku  - the cell index limits for the Z direction;
!     t, dt   - time and time increment;
!     x, y, z - the block coordinates;
!     qn      - the array of variables to update;
!
!===============================================================================
!
  subroutine user_boundary_x(is, jl, ju, kl, ku, t, dt, x, y, z, qn)

    use coordinates, only : nn => bcells, nb, ne, nbl, neu
    use equations  , only : ivx, ibx

    implicit none

    integer                         , intent(in)    :: is, jl, ju, kl, ku
    real(kind=8)                    , intent(in)    :: t, dt
    real(kind=8), dimension(:)      , intent(in)    :: x, y, z
    real(kind=8), dimension(:,:,:,:), intent(inout) :: qn

    integer :: i, j, ir

!-------------------------------------------------------------------------------
!
    if (is == 1) then
      do i = nbl, 1, -1
        ir = nb + (nbl - i)
        do j = jl, ju
          qn( : ,i,j,:) =   qn( : ,ir,j,:)

          qn(ivx,i,j,:) = - qn(ivx,ir,j,:)
          qn(ibx,i,j,:) = - qn(ibx,ir,j,:)
        end do
      end do
    else ! is == 1
      do i = neu, nn
        ir = ne - (i - neu)
        do j = jl, ju
          qn( : ,i,j,:) =   qn( : ,ir,j,:)

          qn(ivx,i,j,:) = - qn(ivx,ir,j,:)
          qn(ibx,i,j,:) = - qn(ibx,ir,j,:)
        end do
      end do
    end if

!-------------------------------------------------------------------------------
!
  end subroutine user_boundary_x
!
!===============================================================================
!
! subroutine USER_BOUNDARY_Y:
! --------------------------
!
!   Subroutine updates ghost zones within the specific region along
!   the Y direction.
!
!   Arguments:
!
!     js      - the block side along the Y direction for the ghost zone update;
!     il, iu  - the cell index limits for the X direction;
!     kl, ku  - the cell index limits for the Z direction;
!     t, dt   - time and time increment;
!     x, y, z - the block coordinates;
!     qn      - the array of variables to update;
!
!===============================================================================
!
  subroutine user_boundary_y(js, il, iu, kl, ku, t, dt, x, y, z, qn)

    use coordinates, only : nn => bcells, nb, ne, nbl, neu
    use equations  , only : ivy, iby

    implicit none

    integer                         , intent(in)    :: js, il, iu, kl, ku
    real(kind=8)                    , intent(in)    :: t, dt
    real(kind=8), dimension(:)      , intent(in)    :: x, y, z
    real(kind=8), dimension(:,:,:,:), intent(inout) :: qn

    integer :: i, j, jr

!-------------------------------------------------------------------------------
!
    if (js == 1) then
      do j = nbl, 1, -1
        jr = nb + (nbl - j)
        do i = il, iu
          qn( : ,i,j,:) =   qn( : ,i,jr,:)

          qn(ivy,i,j,:) = - qn(ivy,i,jr,:)
          qn(iby,i,j,:) = - qn(iby,i,jr,:)
        end do
      end do
    else ! js == 1
      do j = neu, nn
        jr = ne - (j - neu)
        do i = il, iu
          qn( : ,i,j,:) =   qn( : ,i,jr,:)

          qn(ivy,i,j,:) = - qn(ivy,i,jr,:)
          qn(iby,i,j,:) = - qn(iby,i,jr,:)
        end do
      end do
    end if

!-------------------------------------------------------------------------------
!
  end subroutine user_boundary_y
!
!===============================================================================
!
! subroutine USER_BOUNDARY_Z:
! --------------------------
!
!   Subroutine updates ghost zones within the specific region along
!   the Z direction.
!
!   Arguments:
!
!     ks      - the block side along the Z direction for the ghost zone update;
!     il, iu  - the cell index limits for the X direction;
!     jl, ju  - the cell index limits for the Y direction;
!     t, dt   - time and time increment;
!     x, y, z - the block coordinates;
!     qn      - the array of variables to update;
!
!===============================================================================
!
  subroutine user_boundary_z(ks, il, iu, jl, ju, t, dt, x, y, z, qn)

    implicit none

    integer                         , intent(in)    :: ks, il, iu,jl, ju
    real(kind=8)                    , intent(in)    :: t, dt
    real(kind=8), dimension(:)      , intent(in)    :: x, y, z
    real(kind=8), dimension(:,:,:,:), intent(inout) :: qn

!-------------------------------------------------------------------------------
!

!-------------------------------------------------------------------------------
!
  end subroutine user_boundary_z
!
!===============================================================================
!
! subroutine USER_STATISTICS:
! -------------------------------
!
!   Subroutine can be use to store user defined time statistics. The file to
!   store these statistics should be properly created in subroutine
!   initialize_user_problem() and closed in finalize_user_problem().
!
!   Arguments:
!
!     time - the current simulation time;
!
!===============================================================================
!
  subroutine user_statistics(time)

    use blocks     , only : block_meta, block_data, list_data
    use coordinates, only : nc => ncells, nn => bcells, nb, ne, nbm, nep
#if NDIMS == 3
    use coordinates, only : periodic, ymin, ymax, zmin, zmax
#else /* NDIMS == 3 */
    use coordinates, only : periodic, ymin, ymax
#endif /* NDIMS == 3 */
    use coordinates, only : ax, adx, ady, adz
#if NDIMS == 3
    use equations  , only : ivx, ivy, ivz, ibx, iby, ibz, ibp
#else /* NDIMS == 3 */
    use equations  , only : ivx, ivy, ibx, iby, ibz, ibp
#endif /* NDIMS == 3 */
    use helpers    , only : print_message, flush_and_sync
#ifdef MPI
    use mpitools   , only : reduce_sum
#endif /* MPI */
    use operators  , only : curl, derivative_1st
    use workspace  , only : resize_workspace, work, work_in_use

    implicit none
    real(kind=8), intent(in) :: time
    logical, save :: first = .true.

    type(block_meta), pointer :: pmeta
    type(block_data), pointer :: pdata

    integer, save :: nt, i
!$  integer :: omp_get_thread_num
    integer                     :: nl, nu, status
    real(kind=8)                :: dyz
    real(kind=8), dimension(3)  :: dh
    real(kind=8), dimension(16) :: rterms
    real(kind=8), dimension(nn) :: x

    real(kind=8), dimension(:,:,:,:), pointer, save :: cur
    real(kind=8), dimension(:,:)    , pointer, save :: vv, bb, jj
!$omp threadprivate(first, nt, i, cur, vv, bb, jj)

    character(len=*), parameter :: loc = 'USER_PROBLEM:user_time_statistics()'

!-------------------------------------------------------------------------------
!
    nt = 0
!$  nt = omp_get_thread_num()
    if (first) then
      if (resistive) then
        nu = 3 * nn**NDIMS + 9 * nc**(NDIMS - 2)
      else
        nu = 3 * nn**NDIMS + 6 * nc**(NDIMS - 2)
      end if

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

      nl = 1
      nu = nl - 1 + 3 * nn**NDIMS
#if NDIMS == 3
      cur(1:3,1:nn,1:nn,1:nn) => work(nl:nu,nt)
#else /* NDIMS == 3 */
      cur(1:3,1:nn,1:nn,1: 1) => work(nl:nu,nt)
#endif /* NDIMS == 3 */
      nl = nu + 1
      nu = nl - 1 + 3 * nc**(NDIMS - 2)
#if NDIMS == 3
      vv(1:3,1:nc) => work(nl:nu,nt)
#else /* NDIMS == 3 */
      vv(1:3,1: 1) => work(nl:nu,nt)
#endif /* NDIMS == 3 */
      nl = nu + 1
      nu = nl - 1 + 3 * nc**(NDIMS - 2)
#if NDIMS == 3
      bb(1:3,1:nc) => work(nl:nu,nt)
#else /* NDIMS == 3 */
      bb(1:3,1: 1) => work(nl:nu,nt)
#endif /* NDIMS == 3 */
      if (resistive) then
        nl = nu + 1
        nu = nl - 1 + 3 * nc**(NDIMS - 2)
#if NDIMS == 3
        jj(1:3,1:nc) => work(nl:nu,nt)
#else /* NDIMS == 3 */
        jj(1:3,1: 1) => work(nl:nu,nt)
#endif /* NDIMS == 3 */
      end if

      first = .false.
    end if
    rterms(:) = 0.0d+00

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

    work_in_use(nt) = .true.

    pdata => list_data
    do while(associated(pdata))
      pmeta => pdata%meta

      if (pmeta%xmin <= 0.0d+00 .and. pmeta%xmax >= 0.0d+00) then

        dh(1) = adx(pmeta%level)
        dh(2) = ady(pmeta%level)
        dh(3) = adz(pmeta%level)

        x(:) = pmeta%xmin + ax(pmeta%level,:)

        i = nb
        do while(abs(x(i)) > dh(1) .and. i <= ne)
          i = i + 1
        end do

        if (nb < i .and. i <= ne) then

          dyz = ady(pmeta%level) * adz(pmeta%level)

! the integral of |Bx|
!
          rterms(1) = rterms(1) &
#if NDIMS == 3
                   + sum(abs(pdata%u(ibx,i,nb:ne,nb:ne))) * dyz
#else /* NDIMS == 3 */
                   + sum(abs(pdata%u(ibx,i,nb:ne,  :  ))) * dyz
#endif /* NDIMS == 3 */

! calculate current density (J = ∇xB)
!
          if (resistive) &
            call curl(dh(:), pdata%q(ibx:ibz,:,:,:), cur(1:3,:,:,:))

          if (.not. periodic(2)) then

            if (pmeta%ymin - dh(2) < ymin) then

#if NDIMS == 3
              vv(1,:) = 5.0d-01 * sum(pdata%q(ivx,i,nbm:nb,nb:ne), 1)
              vv(2,:) = 5.0d-01 * sum(pdata%q(ivy,i,nbm:nb,nb:ne), 1)
              bb(1,:) = 5.0d-01 * sum(pdata%q(ibx,i,nbm:nb,nb:ne), 1)
              bb(2,:) = 5.0d-01 * sum(pdata%q(iby,i,nbm:nb,nb:ne), 1)
#else /* NDIMS == 3 */
              vv(1,:) = 5.0d-01 * sum(pdata%q(ivx,i,nbm:nb,  :  ), 1)
              vv(2,:) = 5.0d-01 * sum(pdata%q(ivy,i,nbm:nb,  :  ), 1)
              bb(1,:) = 5.0d-01 * sum(pdata%q(ibx,i,nbm:nb,  :  ), 1)
              bb(2,:) = 5.0d-01 * sum(pdata%q(iby,i,nbm:nb,  :  ), 1)
#endif /* NDIMS == 3 */

! advection of Bx across the Y-boundary
!
              rterms(2) = rterms(2) &
                        + sum(sign(bb(1,:) * vv(2,:), bb(1,:))) * dh(3)

! shear of By along the Y-boundary
!
              rterms(4) = rterms(4) &
                        - sum(sign(bb(2,:) * vv(1,:), bb(1,:))) * dh(3)

              if (resistive) then

#if NDIMS == 3
                jj(3,:) = 5.0d-01 * sum(cur(3,i,nbm:nb,nb:ne), 1)
#else /* NDIMS == 3 */
                jj(3,:) = 5.0d-01 * sum(cur(3,i,nbm:nb,  :  ), 1)
#endif /* NDIMS == 3 */

! diffusion of Bx at the Y-boundary
!
                rterms(6) = rterms(6) &
                          + sum(sign(jj(3,:), bb(1,:))) * dh(3)

              end if ! resistivity
            end if

            if (pmeta%ymax + dh(2) > ymax) then

#if NDIMS == 3
              vv(1,:) = 5.0d-01 * sum(pdata%q(ivx,i,ne:nep,nb:ne), 1)
              vv(2,:) = 5.0d-01 * sum(pdata%q(ivy,i,ne:nep,nb:ne), 1)
              bb(1,:) = 5.0d-01 * sum(pdata%q(ibx,i,ne:nep,nb:ne), 1)
              bb(2,:) = 5.0d-01 * sum(pdata%q(iby,i,ne:nep,nb:ne), 1)
#else /* NDIMS == 3 */
              vv(1,:) = 5.0d-01 * sum(pdata%q(ivx,i,ne:nep,  :  ), 1)
              vv(2,:) = 5.0d-01 * sum(pdata%q(ivy,i,ne:nep,  :  ), 1)
              bb(1,:) = 5.0d-01 * sum(pdata%q(ibx,i,ne:nep,  :  ), 1)
              bb(2,:) = 5.0d-01 * sum(pdata%q(iby,i,ne:nep,  :  ), 1)
#endif /* NDIMS == 3 */

! advection of Bx across the Y-boundary
!
              rterms(2) = rterms(2) &
                        - sum(sign(bb(1,:) * vv(2,:), bb(1,:))) * dh(3)

! shear of By along the Y-boundary
!
              rterms(4) = rterms(4) &
                        + sum(sign(bb(2,:) * vv(1,:), bb(1,:))) * dh(3)

              if (resistive) then

#if NDIMS == 3
                jj(3,:) = 5.0d-01 * sum(cur(3,i,ne:nep,nb:ne), 1)
#else /* NDIMS == 3 */
                jj(3,:) = 5.0d-01 * sum(cur(3,i,ne:nep,  :  ), 1)
#endif /* NDIMS == 3 */

! diffusion of Bx at the Y-boundary
!
                rterms(6) = rterms(6) &
                          - sum(sign(jj(3,:), bb(1,:))) * dh(3)

              end if ! resistivity
            end if
          end if ! non-periodic along Y

#if NDIMS == 3
          if (.not. periodic(3)) then

            if (pmeta%zmin - dh(3) < zmin) then

              vv(1,:) = 5.0d-01 * sum(pdata%q(ivx,i,nb:ne,nbm:nb), 1)
              vv(3,:) = 5.0d-01 * sum(pdata%q(ivz,i,nb:ne,nbm:nb), 1)
              bb(1,:) = 5.0d-01 * sum(pdata%q(ibx,i,nb:ne,nbm:nb), 1)
              bb(3,:) = 5.0d-01 * sum(pdata%q(ibz,i,nb:ne,nbm:nb), 1)

! advection of Bx across the Z-boundary
!
              rterms(3) = rterms(3) &
                        + sum(sign(bb(1,:) * vv(3,:), bb(1,:))) * dh(2)

! shear of Bz along the Z-boundary
!
              rterms(5) = rterms(5) &
                        - sum(sign(bb(3,:) * vv(1,:), bb(1,:))) * dh(2)

              if (resistive) then
                jj(2,:) = 5.0d-01 * sum(cur(2,i,nb:ne,nbm:nb), 1)

! diffusion of Bx at the Z-boundary
!
                rterms(7) = rterms(7) &
                          + sum(sign(jj(2,:), bb(1,:))) * dh(2)

              end if ! resistivity
            end if

            if (pmeta%zmax + dh(3) > zmax) then

              vv(1,:) = 5.0d-01 * sum(pdata%q(ivx,i,nb:ne,ne:nep), 1)
              vv(3,:) = 5.0d-01 * sum(pdata%q(ivz,i,nb:ne,ne:nep), 1)
              bb(1,:) = 5.0d-01 * sum(pdata%q(ibx,i,nb:ne,ne:nep), 1)
              bb(3,:) = 5.0d-01 * sum(pdata%q(ibz,i,nb:ne,ne:nep), 1)

! advection of Bx across the Z-boundary
!
              rterms(3) = rterms(3) &
                        - sum(sign(bb(1,:) * vv(3,:), bb(1,:))) * dh(2)

! shear of Bz along the Z-boundary
!
              rterms(5) = rterms(5) &
                        + sum(sign(bb(3,:) * vv(1,:), bb(1,:))) * dh(2)

              if (resistive) then

                jj(2,:) = 5.0d-01 * sum(cur(2,1,nb:ne,ne:nep), 1)

! diffusion of Bx at the Z-boundary
!
                rterms(7) = rterms(7) &
                          - sum(sign(jj(2,:), bb(1,:))) * dh(2)

              end if ! resistivity
            end if
          end if ! non-periodic along Z
#endif /* NDIMS == 3 */

! calculate X-derivative of ψ
!
          call derivative_1st(1, dh(1), pdata%q(ibp,:,:,:), cur(1,:,:,:))

          rterms(8) = rterms(8) &
#if NDIMS == 3
                     - sum(sign(cur(1,i,nb:ne,nb:ne), &
                          pdata%q(ibx,i,nb:ne,nb:ne))) * dyz
#else /* NDIMS == 3 */
                     - sum(sign(cur(1,i,nb:ne,  :  ), &
                          pdata%q(ibx,i,nb:ne,  :  ))) * dyz
#endif /* NDIMS == 3 */

        end if
      end if

      pdata => pdata%next
    end do

    work_in_use(nt) = .false.

#ifdef MPI
    call reduce_sum(rterms(:))
#endif /* MPI */

    rterms(6) = eta * rterms(6)
    rterms(7) = eta * rterms(7)

    write(runit,"(9es25.15e3)") time, rterms(1:8)
    call flush_and_sync(runit)

    100 continue

!-------------------------------------------------------------------------------
!
  end subroutine user_statistics

!===============================================================================
!
end module user_problem