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-2022 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

  real(kind=8), save :: beta = 1.00d+00
  real(kind=8), save :: zeta = 0.00d+00
  real(kind=8), save :: eta  = 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 :: 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 :: dlta = 1.00d-16
  real(kind=8), save :: blim = 1.00d+00

  integer     , save :: pert = 0
  integer     , save :: nper = 10
  real(kind=8), save :: bper = 0.00d+00
  real(kind=8), save :: vper = 0.00d+00
  real(kind=8), save :: kper = 1.00d+00
  real(kind=8), save :: kvec = 1.00d+00
  real(kind=8), save :: xcut = 1.00d+99
  real(kind=8), save :: ycut = 1.00d+99
  real(kind=8), save :: xdec = 1.00d-01
  real(kind=8), save :: ydec = 1.00d-01

  real(kind=8), save :: ylo  = -9.00d+99
  real(kind=8), save :: yup  =  9.00d+99

  real(kind=8), save :: tabs = 1.00d+00
  real(kind=8), save :: adif = 5.00d-01
  real(kind=8), save :: acut = 1.00d+00
  real(kind=8), save :: adec = 1.00d+00
  real(kind=8), save :: yabs = 9.00d+99

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

  logical, save :: absorption = .false.
  logical, save :: resistive  = .false.

  real(kind=8), dimension(:), allocatable :: kx, ky, kz, ux, uy, uz, ph

  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
#if NDIMS == 3
    use constants  , only : pi
#endif /* NDIMS == 3 */
    use constants  , only : pi2
    use coordinates, only : ng => nghosts, ady, xlen, zlen, ymax
    use equations  , only : magnetized, csnd, csnd2
    use helpers    , only : print_section, print_parameter, print_message
    use mesh       , only : setup_problem
    use parameters , only : get_parameter
    use random     , only : randuni, randsym
    use shapes     , only : update_shapes

    implicit none

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

    character(len=64) :: perturbation = "noise", append = "off", fname
    character(len=64) :: enable_absorption = "off"
    logical           :: flag
    integer           :: n, kd
    real(kind=8)      :: thh, fc, ydis = 9.00d+99
#if NDIMS == 3
    real(kind=8)      :: thv, tx, ty, tz, tt
#else /* NDIMS == 3 */
    real(kind=8)      :: ka
#endif /* NDIMS == 3 */

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

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

    if (.not. magnetized) then
      if (verbose) &
        call print_message(loc, &
               "The reconnection problem does not work without magnetic field.")
      status = 1
      return
    end if

! get the reconnection equilibrium parameters
!
    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
    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.
    end if
    call get_parameter("dens", dens)
    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("bamp", bamp)
    call get_parameter("bgui", bgui)

    call get_parameter("blimit", blim)

! calculate the maximum magnetic pressure, thermal pressure from the plasma-β
! parameters, and the sound speed in the case of isothermal equations of state
!
    pmag  = 0.5d+00 * (bamp**2 + bgui**2)
    pres  = beta * pmag
    ptot  = pres + pmag
    csnd2 = pres / dens
    csnd  = sqrt(csnd2)
    valf  = sqrt(2.0d+00 * pmag / dens)
    lund  = valf / max(tiny(eta), eta)
    dlta  = lund**(- 1.0d+00 / 3.0d+00)
    blim  = max(blim, ng * ady(1))

! get the geometry parameters
!
    call get_parameter("delta", dlta)
    if (dlta < 0.0d+00) then
      if (verbose) &
        call print_message(loc, "Density must be positive (dens > 0.0)!")
      status = 1
      return
    end if

! get the perturbation parameters
!
    call get_parameter("perturbation", perturbation)
    call get_parameter("nper"        , nper)
    call get_parameter("bper"        , bper)
    call get_parameter("vper"        , vper)
    call get_parameter("kper"        , kper)
    call get_parameter("xcut"        , xcut)
    call get_parameter("ycut"        , ycut)
    call get_parameter("xdec"        , xdec)
    call get_parameter("ydec"        , ydec)

! choose the perturbation type
!
    select case(perturbation)
    case('noise', 'random')
      pert = 0
    case('mode', 'one mode', 'one-mode', 'one_mode')
      pert = 1
    case('multi-wave', 'random waves', 'random-waves', 'random_waves')
      pert = 2
    case default
      perturbation = 'mode'
      pert = 1
    end select

! prepare the wave vector of the perturbation
!
    kvec = pi2 * kper

! prepare the wave vectors for multi-wave perturbation
!
    if (pert == 2) then

! allocate phase and wave vector components
!
      allocate(kx(nper), ky(nper), kz(nper), ux(nper), uy(nper), uz(nper),     &
               ph(nper), stat = status)

      if (status == 0) then

! choose random wave vector directions
!
        kd = int(xlen / zlen)
        fc = 1.0d+00 / sqrt(1.0d+00 * nper)
        do n = 1, nper
          thh   = pi2 * randuni()
#if NDIMS == 3
          thv   = pi  * randsym()
          tx    = cos(thh) * cos(thv)
          ty    = sin(thh) * cos(thv)
          tz    =            sin(thv)
          kx(n) = pi2 * nint(kper * tx)
          ky(n) = pi2 * nint(kper * ty)
          kz(n) = pi2 * nint(kper * tz / kd) * kd
          tt    = 0.0d+00
          do while(tt < 1.0d-08)
            thh   = pi2 * randuni()
            thv   = pi  * randsym()
            tx    = cos(thh) * cos(thv)
            ty    = sin(thh) * cos(thv)
            tz    =            sin(thv)
            ux(n) = ty * kz(n) - tz * ky(n)
            uy(n) = tz * kx(n) - tx * kz(n)
            uz(n) = tx * ky(n) - ty * kx(n)
            tt    = sqrt(ux(n)**2 + uy(n)**2 + uz(n)**2)
          end do
          ux(n) = fc * ux(n) / tt
          uy(n) = fc * uy(n) / tt
          uz(n) = fc * uz(n) / tt
#else /* NDIMS == 3 */
          kx(n) = pi2 * nint(kper * cos(thh))
          ky(n) = pi2 * nint(kper * sin(thh))
          kz(n) = 0.0d+00
          ka    = sqrt(kx(n)**2 + ky(n)**2)
          ux(n) =   fc * ky(n) / ka
          uy(n) = - fc * kx(n) / ka
          uz(n) = 0.0d+00
#endif /* NDIMS == 3 */
          ph(n) = pi2 * randuni()
        end do


      else
        if (verbose) &
          call print_message(loc, &
                           "Could not allocate space for perturbation vectors!")
        return
      end if
    end if ! pert = 2

! determine if to append or create another file
!
    call get_parameter("reconnection_append", append)
    select case(trim(append))
    case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
      write(fname, "('reconnection.dat')")
      inquire(file = fname, exist = flag)
    case default
      write(fname, "('reconnection_',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,26a25)')                                           &
                'time', '|Bx| int', '|Bx| inf'                               &
                  , '|Bx| y-adv', '|Bx| z-adv', '|Bx| y-shr', '|Bx| z-shr'   &
                  , '|Bx| y-dif', '|Bx| z-dif', '|Bx| Psi'                   &
                  , 'Vin lower' , 'Vin upper'                                &
                  , 'Emag', 'Emag inf'                                       &
                  , 'Emag x-adv', 'Emag y-adv', 'Emag z-adv'                 &
                  , 'Emag x-v.b', 'Emag y-v.b', 'Emag z-v.b'                 &
                  , 'Emag x-dif', 'Emag y-dif', 'Emag z-dif'                 &
                  , 'Emag-Ekin', 'Emag-Eint', 'Emag-Psi'
    write(runit,"('#')")

    call get_parameter("ydistance", ydis)
    ydis = min(ydis, ymax)
    ylo  = - ydis
    yup  =   ydis

! absorption parameters
!
    call get_parameter("enable_shapes", enable_absorption)
    select case(trim(enable_absorption))
    case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
      absorption = .true.
    case default
      absorption = .false.
    end select
    call get_parameter("tabs", tabs)
    call get_parameter("adif", adif)
    call get_parameter("acut", acut)
    call get_parameter("adec", adec)
    yabs = ymax - abs(acut)

! print information about the user problem setup
!
    call print_section(verbose, "Parameters (* - set, otherwise calculated)")
    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'                 , bamp)
    call print_parameter(verbose, '(*) bgui (guide field)'   , bgui)
    call print_parameter(verbose, '( ) pres (thermal pres.)' , pres)
    call print_parameter(verbose, '( ) pmag (magnetic pres.)', pmag)
    call print_parameter(verbose, '( ) ptot (total pressure)', ptot)
    call print_parameter(verbose, '( ) csnd (sound speed)'   , csnd)
    call print_parameter(verbose, '( ) Valf (Alfven speed)'  , valf)
    if (resistive) then
      call print_parameter(verbose, '( ) S (Lundquist number)', lund)
    end if
    call print_parameter(verbose, '(*) delta (thickness)', dlta)
    call print_parameter(verbose, '(*) blim'             , blim)
    call print_parameter(verbose, '(*) perturbation', perturbation)
    call print_parameter(verbose, '(*) bper'        , bper)
    call print_parameter(verbose, '(*) vper'        , vper)
    if (pert >= 1) then
      call print_parameter(verbose, '(*) kper'      , kper)
    end if
    if (pert == 2) then
      call print_parameter(verbose, '(*) nper'      , nper)
    end if
    call print_parameter(verbose, '(*) xcut'        , xcut)
    call print_parameter(verbose, '(*) ycut'        , ycut)
    call print_parameter(verbose, '(*) xdec'        , xdec)
    call print_parameter(verbose, '(*) ydec'        , ydec)
    call print_parameter(verbose, '(*) ydistance (Vrec)', ydis)
    if (absorption) then
      call print_parameter(verbose, '(*) tabs (absorption)', tabs)
      call print_parameter(verbose, '(*) adif (diffusion)' , adif)
      call print_parameter(verbose, '(*) acut'             , acut)
      call print_parameter(verbose, '(*) adec'             , adec)
      call print_parameter(verbose, '( ) yabs'             , yabs)
    end if

! set procedure pointers for the problem setup subroutine, the shapes update
! for the absorption layer, and boundary conditions;
!
    setup_problem     => setup_user_problem
    update_shapes     => update_user_shapes
    custom_boundary_x => user_boundary_x
    custom_boundary_y => user_boundary_y

!-------------------------------------------------------------------------------
!
  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)

    use helpers, only : print_message

    implicit none

    integer, intent(out) :: status

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

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

    close(runit)

    if (pert == 2) then
      deallocate(kx, ky, kz, ux, uy, uz, ph, stat = status)
      if (status /= 0) &
        call print_message(loc, &
                         "Could not deallocate space for perturbation vectors!")
    end if

!-------------------------------------------------------------------------------
!
  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
    use coordinates, only : ax, ay, adx, ady, ylen
#if NDIMS == 3
    use coordinates, only : az, adz
#endif /* NDIMS == 3 */
    use equations  , only : prim2cons
    use equations  , only : nv, ns
    use equations  , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp, isl
    use equations  , only : csnd2
    use helpers    , only : print_message
    use operators  , only : curl
    use random     , only : randsym
    use workspace  , only : resize_workspace, work, work_in_use

    implicit none

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

    logical, save :: first = .true.

    integer      :: i, j, k, n, status
    real(kind=8) :: xpl, xpu, ypl, ypu
    real(kind=8) :: xp, yp
    real(kind=8) :: yrat
    real(kind=8) :: kv, fv
    real(kind=8) :: vx = 0.0d+00, vy = 0.0d+00, vv, va
#if NDIMS == 3
    real(kind=8) :: vz = 0.0d+00
#endif /* NDIMS == 3 */
    real(kind=8) :: pz = 0.0d+00, pp, ba
#if NDIMS == 3
    real(kind=8) :: px = 0.0d+00, py = 0.0d+00
#endif /* NDIMS == 3 */

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

    real(kind=8), dimension(:,:,:,:), pointer, save :: qpert, pot
!$omp threadprivate(first, nt, qpert, pot)

    real(kind=8), dimension(nv,nn) :: q, u, qprof
    real(kind=8), dimension(nn)    :: xl, xu, xc, fx
    real(kind=8), dimension(nn)    :: yl, yu, yc, fy
#if NDIMS == 3
    real(kind=8), dimension(nn)    :: zc
#endif /* NDIMS == 3 */
    real(kind=8), dimension(3)     :: dh

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

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

    nt = 0
!$  nt = omp_get_thread_num()
    if (first) then
      n = (nv + 3) * nn**NDIMS

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

      i = nv * nn**NDIMS
#if NDIMS == 3
      qpert(1:nv,1:nn,1:nn,1:nn) => work(  1:i,nt)
      pot   (1:3,1:nn,1:nn,1:nn) => work(i+1:n,nt)
#else /* NDIMS == 3 */
      qpert(1:nv,1:nn,1:nn,1: 1) => work(  1:i,nt)
      pot   (1:3,1:nn,1:nn,1: 1) => work(i+1:n,nt)
#endif /* NDIMS == 3 */

      first = .false.
    end if

    dh(1) = adx(pdata%meta%level)
    dh(2) = ady(pdata%meta%level)
#if NDIMS == 3
    dh(3) = adz(pdata%meta%level)
#endif /* NDIMS == 3 */

    xc(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
    yc(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
#if NDIMS == 3
    zc(:) = pdata%meta%zmin + az(pdata%meta%level,:)
#endif /* NDIMS == 3 */

    xl(:) = pdata%meta%xmin + ax(pdata%meta%level,:) - 5.0d-01 * dh(1)
    xu(:) = xl(:) + dh(1)
    yl(:) = pdata%meta%ymin + ay(pdata%meta%level,:) - 5.0d-01 * dh(2)
    yu(:) = yl(:) + dh(2)

! set the equilibrium profile
!
    do j = 1, nn
      qprof(ibx,j) = bamp * max(-1.0d+00, min(1.0d+00,                         &
                     dlta * (log_cosh(yu(j) / dlta)                            &
                           - log_cosh(yl(j) / dlta)) / dh(2)))
    end do
    qprof(iby,:) = 0.0d+00
    qprof(ibz,:) = zeta * sqrt(bamp**2 - qprof(ibx,:)**2) + bgui
    if (ipr > 0) then
      qprof(idn,:) =  dens
      qprof(ipr,:) =  ptot - 0.5d+00 * sum(qprof(ibx:ibz,:)**2, 1)
    else
      qprof(idn,:) = (ptot - 0.5d+00 * sum(qprof(ibx:ibz,:)**2, 1)) / csnd2
    end if
    qprof(ivx,:) = 0.0d+00
    qprof(ivy,:) = 0.0d+00
    qprof(ivz,:) = 0.0d+00
    qprof(ibp,:) = 0.0d+00
    if (ns > 0) then
      qprof(isl,:) = sign(1.0d+00, yc(:))
    end if

! ratio of the current sheet thickness to the cell size
!
    yrat  = dlta / dh(2)

! prepare decaying factors
!
    fv = 0.5d+00 * pi
    do i = 1, nn
      xp    = fv * min(1.0d+00, max(0.0d+00, abs(xc(i)) - xcut) / xdec)
      fx(i) = cos(xp)**2
    end do ! i = 1, nn
    do j = 1, nn
      yp    = fv * min(1.0d+00, max(0.0d+00, abs(yc(j)) - ycut) / ydec)
      fy(j) = cos(yp)**2
    end do ! i = 1, nn

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

    work_in_use(nt) = .true.

    qpert(:,:,:,:) = 0.0d+00

    if (pert == 0) then

! of velocity
!
      if (abs(vper) > 0.0d+00) then

! initiate the random velocity components
!
#if NDIMS == 3
        do k = 1, nn
#endif /* NDIMS == 3 */
          do j = 1, nn
            if (fy(j) > 0.0d+00) then
              do i = 1, nn
                if (fx(i) > 0.0d+00) then

! calculate the profile of perturbation amplitude
!
                  va = vper * fx(i) * fy(j)

! get the random direction
!
                  vv = 0.0d+00
                  do while(vv < 1.0d-08)
                    vx = randsym()
                    vy = randsym()
#if NDIMS == 3
                    vz = randsym()
#endif /* NDIMS == 3 */

#if NDIMS == 3
                    vv = sqrt(vx * vx + vy * vy + vz * vz)
#else /* NDIMS == 3 */
                    vv = sqrt(vx * vx + vy * vy)
#endif /* NDIMS == 3 */
                  end do

                  qpert(ivx,i,j,k) = va * (vx / vv)
                  qpert(ivy,i,j,k) = va * (vy / vv)
#if NDIMS == 3
                  qpert(ivz,i,j,k) = va * (vz / vv)
#endif /* NDIMS == 3 */

                end if ! |x| < xcut
              end do ! i = 1, nn
            end if ! |y| < ycut
          end do ! j = 1, nn
#if NDIMS == 3
        end do ! k = 1, nn
#endif /* NDIMS == 3 */

      end if ! vper /= 0.0

! of magnetic field
!
      if (abs(bper) > 0.0d+00) then

! reset the potential
!
        pot(:,:,:,:) = 0.0d+00

! initiate the random magnetic field components
!
#if NDIMS == 3
        do k = 1, nn
#endif /* NDIMS == 3 */
          do j = 1, nn
            if (fy(j) > 0.0d+00) then
              do i = 1, nn
                if (fx(i) > 0.0d+00) then

! calculate the profile of perturbation amplitude
!
                  ba = dh(1) * bper * fx(i) * fy(j)

! get the random direction
!
                  pp = 0.0d+00
                  do while(pp < 1.0d-08)
#if NDIMS == 3
                    px = randsym()
                    py = randsym()
#endif /* NDIMS == 3 */
                    pz = randsym()

#if NDIMS == 3
                    pp = sqrt(px * px + py * py + pz * pz)
#else /* NDIMS == 3 */
                    pp = abs(pz)
#endif /* NDIMS == 3 */
                  end do

#if NDIMS == 3
                  pot(1,i,j,k) = ba * (px / pp)
                  pot(2,i,j,k) = ba * (py / pp)
#endif /* NDIMS == 3 */
                  pot(3,i,j,k) = ba * (pz / pp)

                end if ! |x| < xcut
              end do ! i = 1, nn
            end if ! |y| < ycut
          end do ! j = 1, nn
#if NDIMS == 3
        end do ! k = 1, nn
#endif /* NDIMS == 3 */

! calculate magnetic field perturbation components from vector potential
!
        call curl(dh(1:3), pot(:,:,:,:), qpert(ibx:ibz,:,:,:))

      end if ! bper /= 0.0

    end if ! pert == 0

! one-mode perturbation
!
    if (pert == 1) then

! of velocity
!
      if (abs(vper) > 0.0d+00) then

! reset the potential
!
        pot(:,:,:,:) = 0.0d+00

! calculate the perturbation factor
!
        fv = vper * ylen / (pi2 * pi2 * pi * dh(1) * dh(2) * kper)

! calculate the perturbation profile
!
#if NDIMS == 3
        do k = 1, nn
#endif /* NDIMS == 3 */
          do j = 1, nn
            if (fy(j) > 0.0d+00) then
              ypl = pi * yl(j) / ylen
              ypu = pi * yu(j) / ylen
              do i = 1, nn
                if (fx(i) > 0.0d+00) then
                  xpl = pi2 * kper * xl(i)
                  xpu = pi2 * kper * xu(i)

                  pot(3,i,j,k) = fv * fx(i) * fy(j) * (cos(xpl) - cos(xpu))    &
                                                    * (cos(ypl) - cos(ypu))

                end if ! |x| < xcut
              end do ! i = 1, nn
            end if ! |y| < ycut
          end do ! j = 1, nn
#if NDIMS == 3
        end do ! k = 1, nn
#endif /* NDIMS == 3 */

! calculate magnetic field components from vector potential
!
        call curl(dh(1:3), pot(:,:,:,:), qpert(ivx:ivz,:,:,:))

      end if ! vper /= 0.0

! of magnetic field
!
      if (abs(bper) > 0.0d+00) then

! reset the potential
!
        pot(:,:,:,:) = 0.0d+00

! calculate the perturbation factor
!
        fv = bper * ylen / (pi2 * pi2 * pi * dh(1) * dh(2) * kper)

! calculate the perturbation profile
!
#if NDIMS == 3
        do k = 1, nn
#endif /* NDIMS == 3 */
          do j = 1, nn
            if (fy(j) > 0.0d+00) then
              ypl = pi * yl(j) / ylen
              ypu = pi * yu(j) / ylen
              do i = 1, nn
                if (fx(i) > 0.0d+00) then
                  xpl = pi2 * kper * xl(i)
                  xpu = pi2 * kper * xu(i)

                  pot(3,i,j,k) = fv * fx(i) * fy(j) * (sin(xpl) - sin(xpu))    &
                                                    * (sin(ypl) - sin(ypu))

                end if ! |x| < xcut
              end do ! i = 1, nn
            end if ! |y| < ycut
          end do ! j = 1, nn
#if NDIMS == 3
        end do ! k = 1, nn
#endif /* NDIMS == 3 */

! calculate magnetic field components from vector potential
!
        call curl(dh(1:3), pot(:,:,:,:), qpert(ibx:ibz,:,:,:))

      end if ! bper /= 0.0

    end if ! pert == 1

! prepare the random perturbation of velocity
!
    if (pert == 2) then

      if (abs(vper) > 0.0d+00) then

! iterate over the block position and initiate the velocity perturbation
!
#if NDIMS == 3
        do k = 1, nn
#endif /* NDIMS == 3 */
          do j = 1, nn
            if (fy(j) > 0.0d+00) then
              do i = 1, nn
                if (fx(i) > 0.0d+00) then

! calculate the velocity amplitude profile
!
                  fv = vper * fx(i) * fy(j)

! add perturbation components
!
                  do n = 1, nper
#if NDIMS == 3
                    kv = kx(n) * xc(i) + ky(n) * yc(j) + kz(n) * zc(k) + ph(n)
#else /* NDIMS == 3 */
                    kv = kx(n) * xc(i) + ky(n) * yc(j)                 + ph(n)
#endif /* NDIMS == 3 */
                    va = fv * sin(kv)

                    qpert(ivx,i,j,k) = qpert(ivx,i,j,k) + va * ux(n)
                    qpert(ivy,i,j,k) = qpert(ivy,i,j,k) + va * uy(n)
#if NDIMS == 3
                    qpert(ivz,i,j,k) = qpert(ivz,i,j,k) + va * uz(n)
#endif /* NDIMS == 3 */
                  end do

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

      end if ! vper /= 0.0

    end if ! pert == 2

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

! set the variable profiles with perturbations
!
        q(:,:) = qprof(:,:) + qpert(:,i,:,k)

! convert the primitive variables to the conservative ones
!
        call prim2cons(q(:,:), u(:,:), .true.)

! copy the primitive variables to the current block
!
        pdata%q(:,i,:,k) = q(:,:)

! copy the conserved variables to the current block
!
        pdata%u(:,i,:,k) = u(:,:)

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

    work_in_use(nt) = .false.

    100 continue

!-------------------------------------------------------------------------------
!
  end subroutine setup_user_problem
!
!===============================================================================
!
! subroutine UPDATE_USER_SOURCES:
! ------------------------------
!
!   Subroutine adds the user defined source terms.
!
!   Arguments:
!
!     pdata - the pointer to a data block;
!     t, dt - the time and time increment;
!     du    - the array of variable increment;
!
!===============================================================================
!
  subroutine update_user_sources(pdata, t, dt, du)

    use blocks, only : block_data

    implicit none

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

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

!-------------------------------------------------------------------------------
!
  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
    use constants  , only : pi
    use coordinates, only : nn => bcells
    use coordinates, only : ay, adx, ady
#if NDIMS == 3
    use coordinates, only : adz
#endif /* NDIMS == 3 */
    use equations  , only : nv
    use equations  , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
    use equations  , only : prim2cons
    use helpers    , only : print_message
    use operators  , only : laplace
    use workspace  , only : resize_workspace, work, work_in_use

    implicit none

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

    logical, save :: first = .true.

    integer       :: j, k, status
    real(kind=8)  :: fl, fr, fd, fa, fb

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

    real(kind=8), dimension(3)     :: dh = 1.0d+00
    real(kind=8), dimension(nn)    :: yc, fc
    real(kind=8), dimension(nv,nn) :: q, u

    real(kind=8), dimension(:,:,:), pointer, save :: q2
!$omp threadprivate(first, nt, q2)

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

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

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

      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
      q2(1:nn,1:nn,1:nn) => work(1:j,nt)
#else /* NDIMS == 3 */
      q2(1:nn,1:nn,1: 1) => work(1:j,nt)
#endif /* NDIMS == 3 */

      first = .false.
    end if

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

    work_in_use(nt) = .true.

    yc(:) = pdata%meta%ymin + ay(pdata%meta%level,:)

    if (yc(1) < -yabs .or. yc(nn) > yabs) then

      dh(1) = adx(pdata%meta%level)
      dh(2) = ady(pdata%meta%level)
#if NDIMS == 3
      dh(3) = adz(pdata%meta%level)
#endif /* NDIMS == 3 */

      fa = 1.0d+00 - exp(- dt / tabs)
      fb = 5.0d-01 * adif * minval(dh(1:NDIMS))**2

      call laplace(dh, pdata%q(ivy,:,:,:), q2)

      fc(:) = (abs(yc(:)) - yabs) / adec

      do j = 1, nn
        if (fc(j) > 0.0d+00) then
          if (fc(j) < 1.0d+00) then
            fr = fa * sin(5.0d-01 * pi * fc(j))**2
          else
            fr = fa
          end if
          fl = 1.0d+00 - fr
          fd = fr * fb

#if NDIMS == 3
          do k = 1, nn
#endif /* NDIMS == 3 */
            q(idn,:) = fl * pdata%q(idn,:,j,k) + fr * dens
            q(ivx,:) = fl * pdata%q(ivx,:,j,k)
            q(ivy,:) =      pdata%q(ivy,:,j,k) + fd * q2(:,j,k)
            q(ivz,:) = fl * pdata%q(ivz,:,j,k)
            q(ibx,:) = fl * pdata%q(ibx,:,j,k) + fr * sign(bamp, yc(j))
            q(iby,:) = fl * pdata%q(iby,:,j,k)
            q(ibz,:) = fl * pdata%q(ibz,:,j,k) + fr * bgui
            q(ibp,:) = fl * pdata%q(ibp,:,j,k)
            if (ipr > 0) q(ipr,:) = fl * pdata%q(ipr,:,j,k) + fr * pres

            call prim2cons(q(:,:), u(:,:))

            pdata%q(:,:,j,k) = q(:,:)
            pdata%u(:,:,j,k) = u(:,:)
#if NDIMS == 3
          end do
#endif /* NDIMS == 3 */
        end if
      end do
    end if

    work_in_use(nt) = .false.

    100 continue

!-------------------------------------------------------------------------------
!
  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 : magnetized, ivx, ibx, iby, ibp
#if NDIMS == 3
    use equations      , only : ibz
#endif /* NDIMS == 3 */

    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      :: im2, im1, i    , ip1, ip2
    integer      :: jm2, jm1, j    , jp1, jp2
    integer      :: k
#if NDIMS == 3
    integer      :: km2, km1, kp1, kp2
#endif /* NDIMS == 3 */
    real(kind=8) :: dx, dy, dxy
#if NDIMS == 3
    real(kind=8) :: dz, dxz
#endif /* NDIMS == 3 */

!-------------------------------------------------------------------------------
!
    if (magnetized) then

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

      dx  = x(2) - x(1)
      dy  = y(2) - y(1)
#if NDIMS == 3
      dz  = z(2) - z(1)
#endif /* NDIMS == 3 */
      dxy = dx / dy
#if NDIMS == 3
      dxz = dx / dz
#endif /* NDIMS == 3 */

      if (is == 1) then

! iterate over left-side ghost layers
!
        do i = nbl, 1, -1

! calculate neighbor cell indices
!
          ip1 = min(nn, i + 1)
          ip2 = min(nn, i + 2)

! iterate over boundary layer
!
#if NDIMS == 3
          do k = kl, ku
            km2 = max( 1, k - 2)
            km1 = max( 1, k - 1)
            kp1 = min(nn, k + 1)
            kp2 = min(nn, k + 2)
#endif /* NDIMS == 3 */
            do j = jl, ju
              jm2 = max( 1, j - 2)
              jm1 = max( 1, j - 1)
              jp1 = min(nn, j + 1)
              jp2 = min(nn, j + 2)

! make the normal derivative zero
!
              qn(:,i,j,k) = qn(:,nb,j,k)

! prevent the inflow
!
              qn(ivx,i,j,k) = min(0.0d+00, qn(ivx,nb,j,k))

! update the normal component of magnetic field from divergence-free condition
!
              qn(ibx,i,j,k) = qn(ibx,ip2,j,k)                                  &
                            + (qn(iby,ip1,jp1,k) - qn(iby,ip1,jm1,k)) * dxy
#if NDIMS == 3
              qn(ibx,i,j,k) = qn(ibx,i  ,j,k)                                  &
                            + (qn(ibz,ip1,j,kp1) - qn(ibz,ip1,j,km1)) * dxz
#endif /* NDIMS == 3 */
              qn(ibp,i,j,k) = 0.0d+00
            end do ! j = jl, ju
#if NDIMS == 3
          end do ! k = kl, ku
#endif /* NDIMS == 3 */
        end do ! i = nbl, 1, -1

      else ! is == 1

! iterate over right-side ghost layers
!
        do i = neu, nn

! calculate neighbor cell indices
!
          im1 = max( 1, i - 1)
          im2 = max( 1, i - 2)

! iterate over boundary layer
!
#if NDIMS == 3
          do k = kl, ku
            km1 = max( 1, k - 1)
            kp1 = min(nn, k + 1)
            km2 = max( 1, k - 2)
            kp2 = min(nn, k + 2)
#endif /* NDIMS == 3 */
            do j = jl, ju
              jm1 = max( 1, j - 1)
              jp1 = min(nn, j + 1)
              jm2 = max( 1, j - 2)
              jp2 = min(nn, j + 2)

! make the normal derivative zero
!
              qn(:,i,j,k) = qn(:,ne,j,k)

! prevent the inflow
!
              qn(ivx,i,j,k) = max(0.0d+00, qn(ivx,ne,j,k))

! update the normal component of magnetic field from divergence-free condition
!
              qn(ibx,i,j,k) = qn(ibx,im2,j,k)                              &
                            + (qn(iby,im1,jm1,k) - qn(iby,im1,jp1,k)) * dxy
#if NDIMS == 3
              qn(ibx,i,j,k) = qn(ibx,i  ,j,k)                              &
                            + (qn(ibz,im1,j,km1) - qn(ibz,im1,j,kp1)) * dxz
#endif /* NDIMS == 3 */
              qn(ibp,i,j,k) = 0.0d+00
            end do ! j = jl, ju
#if NDIMS == 3
          end do ! k = kl, ku
#endif /* NDIMS == 3 */
        end do ! i = neu, nn
      end if ! is == 1
    else ! ibx > 0
      if (is == 1) then
        do i = nbl, 1, -1
#if NDIMS == 3
          qn( : ,i,jl:ju,kl:ku) =     qn( : ,nb,jl:ju,kl:ku)
          qn(ivx,i,jl:ju,kl:ku) = min(qn(ivx,nb,jl:ju,kl:ku), 0.0d+00)
#else /* NDIMS == 3 */
          qn( : ,i,jl:ju,  :  ) =     qn( : ,nb,jl:ju,  :  )
          qn(ivx,i,jl:ju,  :  ) = min(qn(ivx,nb,jl:ju,  :  ), 0.0d+00)
#endif /* NDIMS == 3 */
        end do ! i = nbl, 1, -1
      else
        do i = neu, nn
#if NDIMS == 3
          qn( : ,i,jl:ju,kl:ku) =     qn( : ,ne,jl:ju,kl:ku)
          qn(ivx,i,jl:ju,kl:ku) = max(qn(ivx,ne,jl:ju,kl:ku), 0.0d+00)
#else /* NDIMS == 3 */
          qn( : ,i,jl:ju,  :  ) =     qn( : ,ne,jl:ju,  :  )
          qn(ivx,i,jl:ju,  :  ) = max(qn(ivx,ne,jl:ju,  :  ), 0.0d+00)
#endif /* NDIMS == 3 */
        end do ! i = neu, nn
      end if
    end if ! magnetized

!-------------------------------------------------------------------------------
!
  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 : magnetized, nv
    use equations  , only : idn, ivy, ipr, ibx, iby, ibz, ibp

    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, im1, ip1
    integer      :: j, jm1, jp1, jm2, jp2
    integer      :: k
#if NDIMS == 3
    integer      :: km1, kp1
#endif /* NDIMS == 3 */
    real(kind=8) :: dx, dy, dyx
#if NDIMS == 3
    real(kind=8) :: dz, dyz
#endif /* NDIMS == 3 */
    real(kind=8) :: fl, fr

!-------------------------------------------------------------------------------
!
    if (magnetized) then

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

      dx  = x(2) - x(1)
      dy  = y(2) - y(1)
#if NDIMS == 3
      dz  = z(2) - z(1)
#endif /* NDIMS == 3 */
      dyx = dy / dx
#if NDIMS == 3
      dyz = dy / dz
#endif /* NDIMS == 3 */

! process left and right side boundary separatelly
!
      if (js == 1) then

! iterate over left-side ghost layers
!
        do j = nbl, 1, -1

! calculate neighbor cell indices
!
          jp1 = min(nn, j + 1)
          jp2 = min(nn, j + 2)

! calculate variable decay coefficients
!
          fr  = (dy * (nb - j - 5.0d-01)) / blim
          fl  = 1.0d+00 - fr

! iterate over boundary layer
!
#if NDIMS == 3
          do k = kl, ku
            km1 = max( 1, k - 1)
            kp1 = min(nn, k + 1)
#endif /* NDIMS == 3 */
            do i = il, iu
              im1 = max( 1, i - 1)
              ip1 = min(nn, i + 1)

! make normal derivatives zero
!
              qn(1:nv,i,j,k) = qn(1:nv,i,nb,k)

! decay density and pressure to their limits
!
              qn(idn,i,j,k) = fl * qn(idn,i,nb,k) + fr * dens
              if (ipr > 0) qn(ipr,i,j,k) = fl * qn(ipr,i,nb,k) + fr * pres

! decay magnetic field to its limit
!
              qn(ibx,i,j,k) = fl * qn(ibx,i,nb,k) - fr * bamp
              qn(ibz,i,j,k) = fl * qn(ibz,i,nb,k) + fr * bgui

! update By from div(B)=0
!
              qn(iby,i,j,k) = qn(iby,i,jp2,k)                              &
                           + (qn(ibx,ip1,jp1,k) - qn(ibx,im1,jp1,k)) * dyx
#if NDIMS == 3
              qn(iby,i,j,k) = qn(iby,i,j  ,k)                              &
                           + (qn(ibz,i,jp1,kp1) - qn(ibz,i,jp1,km1)) * dyz
#endif /* NDIMS == 3 */
              qn(ibp,i,j,k) = 0.0d+00
            end do ! i = il, iu
#if NDIMS == 3
          end do ! k = kl, ku
#endif /* NDIMS == 3 */
        end do ! j = nbl, 1, -1
      else ! js = 1

! iterate over right-side ghost layers
!
        do j = neu, nn

! calculate neighbor cell indices
!
          jm1 = max( 1, j - 1)
          jm2 = max( 1, j - 2)

! calculate variable decay coefficients
!
          fr  = (dy * (j - ne - 5.0d-01)) / blim
          fl  = 1.0d+00 - fr

! iterate over boundary layer
!
#if NDIMS == 3
          do k = kl, ku
            km1 = max( 1, k - 1)
            kp1 = min(nn, k + 1)
#endif /* NDIMS == 3 */
            do i = il, iu
              im1 = max( 1, i - 1)
              ip1 = min(nn, i + 1)

! make normal derivatives zero
!
              qn(1:nv,i,j,k) = qn(1:nv,i,ne,k)

! decay density and pressure to their limits
!
              qn(idn,i,j,k) = fl * qn(idn,i,ne,k) + fr * dens
              if (ipr > 0) qn(ipr,i,j,k) = fl * qn(ipr,i,ne,k) + fr * pres

! decay magnetic field to its limit
!
              qn(ibx,i,j,k) = fl * qn(ibx,i,ne,k) + fr * bamp
              qn(ibz,i,j,k) = fl * qn(ibz,i,ne,k) + fr * bgui

! update By from div(B)=0
!
              qn(iby,i,j,k) = qn(iby,i,jm2,k)                              &
                           + (qn(ibx,im1,jm1,k) - qn(ibx,ip1,jm1,k)) * dyx
#if NDIMS == 3
              qn(iby,i,j,k) = qn(iby,i,j  ,k)                              &
                           + (qn(ibz,i,jm1,km1) - qn(ibz,i,jm1,kp1)) * dyz
#endif /* NDIMS == 3 */
              qn(ibp,i,j,k) = 0.0d+00
            end do ! i = il, iu
#if NDIMS == 3
          end do ! k = kl, ku
#endif /* NDIMS == 3 */
        end do ! j = neu, nn
      end if ! js = 1
    else ! ibx > 0
      if (js == 1) then
        do j = nbl, 1, -1
#if NDIMS == 3
          qn(1:nv,il:iu,j,kl:ku) = qn(1:nv,il:iu,nb,kl:ku)
          qn(ivy ,il:iu,j,kl:ku) = min(0.0d+00, qn(ivy,il:iu,nb,kl:ku))
#else /* NDIMS == 3 */
          qn(1:nv,il:iu,j,  :  ) = qn(1:nv,il:iu,nb,  :  )
          qn(ivy ,il:iu,j,  :  ) = min(0.0d+00, qn(ivy,il:iu,nb,  :  ))
#endif /* NDIMS == 3 */
        end do ! j = nbl, 1, -1
      else
        do j = neu, nn
#if NDIMS == 3
          qn(1:nv,il:iu,j,kl:ku) = qn(1:nv,il:iu,ne,kl:ku)
          qn(ivy ,il:iu,j,kl:ku) = max(0.0d+00, qn(ivy,il:iu,ne,kl:ku))
#else /* NDIMS == 3 */
          qn(1:nv,il:iu,j,  :  ) = qn(1:nv,il:iu,ne,  :  )
          qn(ivy ,il:iu,j,  :  ) = max(0.0d+00, qn(ivy,il:iu,ne,  :  ))
#endif /* NDIMS == 3 */
        end do ! j = neu, nn
      end if
    end if ! ibx > 0

!-------------------------------------------------------------------------------
!
  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 : nn => bcells, nb, ne
    use coordinates, only : adx, ady, adz, advol, ay, yarea
    use equations  , only : ivx, ivy, ivz, ibx, iby, ibz, ibp
    use helpers    , only : print_message, flush_and_sync
#ifdef MPI
    use mpitools   , only : reduce_sum
#endif /* MPI */
    use operators  , only : curl, gradient
    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
!$  integer :: omp_get_thread_num

    integer                     :: ni, nl, nu, nm, np, status
    real(kind=8)                :: dvol, dxz
    real(kind=8), dimension(3)  :: dh
    real(kind=8), dimension(nn) :: yc
    real(kind=8), dimension(32) :: rterms

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

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

!-------------------------------------------------------------------------------
!
    nt = 0
!$  nt = omp_get_thread_num()
    if (first) then
      ni = 3 * nn**NDIMS

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

#if NDIMS == 3
      cur(1:3,1:nn,1:nn,1:nn) => work(1:ni,nt)
#else /* NDIMS == 3 */
      cur(1:3,1:nn,1:nn,1: 1) => work(1:ni,nt)
#endif /* NDIMS == 3 */

      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

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

      dvol  = advol(pmeta%level)

      dxz   = adx(pmeta%level) * adz(pmeta%level)

      yc(:) = pmeta%ymin + ay(pmeta%level,:)
      ni   = -1
      nl   = nb
      nu   = ne
      do while (yc(nl) < ylo .and. nl < ne)
        nl = nl + 1
      end do
      do while (yc(nu) > yup .and. nu > nb)
        nu = nu - 1
      end do

      if (nl < nu) then

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

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

! the integral of magnetic energy
!
        rterms(12) = rterms(12) &
#if NDIMS == 3
                   + sum(pdata%u(ibx:ibz,nb:ne,nl:nu,nb:ne)**2) * dvol
#else /* NDIMS == 3 */
                   + sum(pdata%u(ibx:ibz,nb:ne,nl:nu,  :  )**2) * dvol
#endif /* NDIMS == 3 */

        if (pmeta%ymin <= ylo .and. ylo < pmeta%ymax) then
          ni = nl
          nm = ni - 1

! get the inflow speed
!
          rterms(10) = rterms(10) &
#if NDIMS == 3
                     + sum(pdata%q(ivy,nb:ne,ni,nb:ne)) * dxz
#else /* NDIMS == 3 */
                     + sum(pdata%q(ivy,nb:ne,ni,  :  )) * dxz
#endif /* NDIMS == 3 */

! mean Bx at boundary
!
          rterms(2) = rterms(2) &
#if NDIMS == 3
                    + sum(abs(pdata%q(ibx,nb:ne,nm:ni,nb:ne))) * dxz
#else /* NDIMS == 3 */
                    + sum(abs(pdata%q(ibx,nb:ne,nm:ni,  :  ))) * dxz
#endif /* NDIMS == 3 */

! advection of Bx along Y
!
          rterms(3) = rterms(3) &
#if NDIMS == 3
                    + sum(abs(pdata%q(ibx,nb:ne,nm:ni,nb:ne)) &
                            * pdata%q(ivy,nb:ne,nm:ni,nb:ne)) * dxz
#else /* NDIMS == 3 */
                    + sum(abs(pdata%q(ibx,nb:ne,nm:ni,  :  )) &
                            * pdata%q(ivy,nb:ne,nm:ni,  :  )) * dxz
#endif /* NDIMS == 3 */

! shear of By along X
!
          rterms(5) = rterms(5) &
#if NDIMS == 3
                    - sum(sign(pdata%q(iby,nb:ne,nm:ni,nb:ne) &
                             * pdata%q(ivx,nb:ne,nm:ni,nb:ne), &
                               pdata%q(ibx,nb:ne,nm:ni,nb:ne))) * dxz
#else /* NDIMS == 3 */
                    - sum(sign(pdata%q(iby,nb:ne,nm:ni,  :  ) &
                             * pdata%q(ivx,nb:ne,nm:ni,  :  ), &
                               pdata%q(ibx,nb:ne,nm:ni,  :  ))) * dxz
#endif /* NDIMS == 3 */

! mean magnetic energy
!
          rterms(13) = rterms(13) &
#if NDIMS == 3
                     + sum(pdata%q(ibx:ibz,nb:ne,nm:ni,nb:ne)**2) * dxz
#else /* NDIMS == 3 */
                     + sum(pdata%q(ibx:ibz,nb:ne,nm:ni,  :  )**2) * dxz
#endif /* NDIMS == 3 */

! advection of magnetic energy
!
          rterms(15) = rterms(15) &
#if NDIMS == 3
                     + sum(sum(pdata%q(ibx:ibz,nb:ne,nm:ni,nb:ne)**2,1) &
                             * pdata%q(ivy    ,nb:ne,nm:ni,nb:ne)) * dxz
#else /* NDIMS == 3 */
                     + sum(sum(pdata%q(ibx:ibz,nb:ne,nm:ni,  :  )**2,1) &
                             * pdata%q(ivy    ,nb:ne,nm:ni,  :  )) * dxz
#endif /* NDIMS == 3 */

! advection of magnetic energy
!
          rterms(18) = rterms(18) &
#if NDIMS == 3
                     - sum(sum(pdata%q(ivx:ivz,nb:ne,nm:ni,nb:ne) &
                             * pdata%q(ibx:ibz,nb:ne,nm:ni,nb:ne),1) &
                             * pdata%q(iby    ,nb:ne,nm:ni,nb:ne)) * dxz
#else /* NDIMS == 3 */
                     - sum(sum(pdata%q(ivx:ivz,nb:ne,nm:ni,  :  ) &
                             * pdata%q(ibx:ibz,nb:ne,nm:ni,  :  ),1) &
                             * pdata%q(iby    ,nb:ne,nm:ni,  :  )) * dxz
#endif /* NDIMS == 3 */

          if (resistive) then

! diffusion of Bx through
!
            rterms(7) = rterms(7) &
#if NDIMS == 3
                      + sum(sign(cur(3,nb:ne,nm:ni,nb:ne), &
                           pdata%q(ibx,nb:ne,nm:ni,nb:ne))) * dxz
#else /* NDIMS == 3 */
                      + sum(sign(cur(3,nb:ne,nm:ni,  :  ), &
                           pdata%q(ibx,nb:ne,nm:ni,  :  ))) * dxz
#endif /* NDIMS == 3 */

! diffusion of magnetic energy through the lower Y boundary
!
            rterms(21) = rterms(21) &
#if NDIMS == 3
                         + sum(cur(3,nb:ne,nm:ni,nb:ne) &
                       * pdata%q(ibx,nb:ne,nm:ni,nb:ne) &
                             - cur(1,nb:ne,nm:ni,nb:ne) &
                       * pdata%q(ibz,nb:ne,nm:ni,nb:ne)) * dxz
#else /* NDIMS == 3 */
                         + sum(cur(3,nb:ne,nm:ni,  :  ) &
                       * pdata%q(ibx,nb:ne,nm:ni,  :  ) &
                             - cur(1,nb:ne,nm:ni,  :  ) &
                       * pdata%q(ibz,nb:ne,nm:ni,  :  )) * dxz
#endif /* NDIMS == 3 */

          end if ! resistivity

        end if

        if (pmeta%ymin < yup .and. yup <= pmeta%ymax) then
          ni = nu
          np = ni + 1

! get the inflow speed
!
          rterms(11) = rterms(11) &
#if NDIMS == 3
                     - sum(pdata%q(ivy,nb:ne,ni,nb:ne)) * dxz
#else /* NDIMS == 3 */
                     - sum(pdata%q(ivy,nb:ne,ni,  :  )) * dxz
#endif /* NDIMS == 3 */

! mean Bx at boundary
!
          rterms(2) = rterms(2) &
#if NDIMS == 3
                    + sum(abs(pdata%q(ibx,nb:ne,ni:np,nb:ne))) * dxz
#else /* NDIMS == 3 */
                    + sum(abs(pdata%q(ibx,nb:ne,ni:np,  :  ))) * dxz
#endif /* NDIMS == 3 */

! advection of Bx along Y
!
          rterms(3) = rterms(3) &
#if NDIMS == 3
                    - sum(abs(pdata%q(ibx,nb:ne,ni:np,nb:ne)) &
                            * pdata%q(ivy,nb:ne,ni:np,nb:ne)) * dxz
#else /* NDIMS == 3 */
                    - sum(abs(pdata%q(ibx,nb:ne,ni:np,  :  )) &
                            * pdata%q(ivy,nb:ne,ni:np,  :  )) * dxz
#endif /* NDIMS == 3 */

! shear of By along X
!
          rterms(5) = rterms(5) &
#if NDIMS == 3
                    + sum(sign(pdata%q(iby,nb:ne,ni:np,nb:ne) &
                             * pdata%q(ivx,nb:ne,ni:np,nb:ne), &
                               pdata%q(ibx,nb:ne,ni:np,nb:ne))) * dxz
#else /* NDIMS == 3 */
                    + sum(sign(pdata%q(iby,nb:ne,ni:np,  :  ) &
                             * pdata%q(ivx,nb:ne,ni:np,  :  ), &
                               pdata%q(ibx,nb:ne,ni:np,  :  ))) * dxz
#endif /* NDIMS == 3 */

! mean magnetic energy
!
          rterms(13) = rterms(13) &
#if NDIMS == 3
                     + sum(pdata%q(ibx:ibz,nb:ne,ni:np,nb:ne)**2) * dxz
#else /* NDIMS == 3 */
                     + sum(pdata%q(ibx:ibz,nb:ne,ni:np,  :  )**2) * dxz
#endif /* NDIMS == 3 */

! advection of magnetic energy
!
          rterms(15) = rterms(15) &
#if NDIMS == 3
                     - sum(sum(pdata%q(ibx:ibz,nb:ne,ni:np,nb:ne)**2,1) &
                             * pdata%q(ivy    ,nb:ne,ni:np,nb:ne)) * dxz
#else /* NDIMS == 3 */
                     - sum(sum(pdata%q(ibx:ibz,nb:ne,ni:np,  :  )**2,1) &
                             * pdata%q(ivy    ,nb:ne,ni:np,  :  )) * dxz
#endif /* NDIMS == 3 */
          rterms(18) = rterms(18) &
#if NDIMS == 3
                     + sum(sum(pdata%q(ivx:ivz,nb:ne,ni:np,nb:ne) &
                             * pdata%q(ibx:ibz,nb:ne,ni:np,nb:ne),1) &
                             * pdata%q(iby    ,nb:ne,ni:np,nb:ne)) * dxz
#else /* NDIMS == 3 */
                     + sum(sum(pdata%q(ivx:ivz,nb:ne,ni:np,  :  ) &
                             * pdata%q(ibx:ibz,nb:ne,ni:np,  :  ),1) &
                             * pdata%q(iby    ,nb:ne,ni:np,  :  )) * dxz
#endif /* NDIMS == 3 */

          if (resistive) then

! diffusion of Bx
!
            rterms(7) = rterms(7) &
#if NDIMS == 3
                      - sum(sign(cur(3,nb:ne,ni:np,nb:ne), &
                           pdata%q(ibx,nb:ne,ni:np,nb:ne))) * dxz
#else /* NDIMS == 3 */
                      - sum(sign(cur(3,nb:ne,ni:np,  :  ), &
                           pdata%q(ibx,nb:ne,ni:np,  :  ))) * dxz
#endif /* NDIMS == 3 */

! diffusion of magnetic energy
!
            rterms(21) = rterms(21) &
#if NDIMS == 3
                         - sum(cur(3,nb:ne,ni:np,nb:ne) &
                       * pdata%q(ibx,nb:ne,ni:np,nb:ne) &
                             - cur(1,nb:ne,ni:np,nb:ne) &
                       * pdata%q(ibz,nb:ne,ni:np,nb:ne)) * dxz
#else /* NDIMS == 3 */
                         - sum(cur(3,nb:ne,ni:np,  :  ) &
                       * pdata%q(ibx,nb:ne,ni:np,  :  ) &
                             - cur(1,nb:ne,ni:np,  :  ) &
                       * pdata%q(ibz,nb:ne,ni:np,  :  )) * dxz
#endif /* NDIMS == 3 */

            end if ! resistivity

        end if

! get the conversion between the magnetic energy and kinetic and
! internal energies
!
        rterms(23) = rterms(23) &
#if NDIMS == 3
                   + sum((pdata%q(ivy,nb:ne,nl:nu,nb:ne) &
                        * pdata%q(ibz,nb:ne,nl:nu,nb:ne) &
                        - pdata%q(ivz,nb:ne,nl:nu,nb:ne) &
                        * pdata%q(iby,nb:ne,nl:nu,nb:ne)) &
                              * cur(1,nb:ne,nl:nu,nb:ne)) * dvol
#else /* NDIMS == 3 */
                   + sum((pdata%q(ivy,nb:ne,nl:nu,  :  ) &
                        * pdata%q(ibz,nb:ne,nl:nu,  :  ) &
                        - pdata%q(ivz,nb:ne,nl:nu,  :  ) &
                        * pdata%q(iby,nb:ne,nl:nu,  :  )) &
                              * cur(1,nb:ne,nl:nu,  :  )) * dvol
#endif /* NDIMS == 3 */

        rterms(23) = rterms(23) &
#if NDIMS == 3
                   + sum((pdata%q(ivz,nb:ne,nl:nu,nb:ne) &
                        * pdata%q(ibx,nb:ne,nl:nu,nb:ne) &
                        - pdata%q(ivx,nb:ne,nl:nu,nb:ne) &
                        * pdata%q(ibz,nb:ne,nl:nu,nb:ne)) &
                              * cur(2,nb:ne,nl:nu,nb:ne)) * dvol
#else /* NDIMS == 3 */
                   + sum((pdata%q(ivz,nb:ne,nl:nu,  :  ) &
                        * pdata%q(ibx,nb:ne,nl:nu,  :  ) &
                        - pdata%q(ivx,nb:ne,nl:nu,  :  ) &
                        * pdata%q(ibz,nb:ne,nl:nu,  :  )) &
                              * cur(2,nb:ne,nl:nu,  :  )) * dvol
#endif /* NDIMS == 3 */

        rterms(23) = rterms(23) &
#if NDIMS == 3
                   + sum((pdata%q(ivx,nb:ne,nl:nu,nb:ne) &
                        * pdata%q(iby,nb:ne,nl:nu,nb:ne) &
                        - pdata%q(ivy,nb:ne,nl:nu,nb:ne) &
                        * pdata%q(ibx,nb:ne,nl:nu,nb:ne)) &
                              * cur(3,nb:ne,nl:nu,nb:ne)) * dvol
#else /* NDIMS == 3 */
                   + sum((pdata%q(ivx,nb:ne,nl:nu,  :  ) &
                        * pdata%q(iby,nb:ne,nl:nu,  :  ) &
                        - pdata%q(ivy,nb:ne,nl:nu,  :  ) &
                        * pdata%q(ibx,nb:ne,nl:nu,  :  )) &
                              * cur(3,nb:ne,nl:nu,  :  )) * dvol
#endif /* NDIMS == 3 */

        if (resistive) then
          rterms(24) = rterms(24) &
#if NDIMS == 3
                     - sum(cur(1:3,nb:ne,nl:nu,nb:ne)**2) * dvol
#else /* NDIMS == 3 */
                     - sum(cur(1:3,nb:ne,nl:nu,  :  )**2) * dvol
#endif /* NDIMS == 3 */
        end if

! calculate gradient (∇ψ)
!
        call gradient(dh(:), pdata%q(ibp,:,:,:), cur(1:3,:,:,:))

        rterms(25) = rterms(25) &
#if NDIMS == 3
                   - sum(pdata%q(ibx:ibz,nb:ne,nl:nu,nb:ne) &
                               * cur(1:3,nb:ne,nl:nu,nb:ne)) * dvol
#else /* NDIMS == 3 */
                   - sum(pdata%q(ibx:ibz,nb:ne,nl:nu,  :  ) &
                               * cur(1:3,nb:ne,nl:nu,  :  )) * dvol
#endif /* NDIMS == 3 */

! contribution to |Bx| from ∇ψ
!
        rterms(9) = rterms(9) &
#if NDIMS == 3
                  - sum(sign(cur(1,nb:ne,nl:nu,nb:ne), &
                       pdata%q(ibx,nb:ne,nl:nu,nb:ne))) * dvol
#else /* NDIMS == 3 */
                  - sum(sign(cur(1,nb:ne,nl:nu,  :  ), &
                       pdata%q(ibx,nb:ne,nl:nu,  :  ))) * dvol
#endif /* NDIMS == 3 */

      end if ! nl < nu

      pdata => pdata%next
    end do

    work_in_use(nt) = .false.

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

    rterms(2)     = 2.50d-01 * rterms(2) / yarea
    rterms(12)    = 5.00d-01 * rterms(12)
    rterms(13)    = 1.25d-01 * rterms(13) / yarea
    rterms(3:6)   = 5.00d-01 * rterms(3:6)
    rterms(14:19) = 5.00d-01 * rterms(14:19)
    rterms(7:8)   = 5.00d-01 * eta * rterms(7:8)
    rterms(20:22) = 5.00d-01 * eta * rterms(20:22)
    rterms(24)    = eta * rterms(24)

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

    100 continue

!-------------------------------------------------------------------------------
!
  end subroutine user_statistics
!
!===============================================================================
!
! subroutine LOG_COSH:
! -------------------
!
!   Function calculates the logarithm of the hyperbolic cosine, which is
!   the result of the integration of tanh(x). Direct calculation using
!   Fortran intrinsic subroutines fails for large values of x, therefore
!   the logarithm of cosh is approximated as |x| + log(1/2) for
!   |x| > threshold.
!
!   Arguments:
!
!     x - function argument;
!
!===============================================================================
!
  function log_cosh(x) result(y)

    implicit none

    real(kind=8), intent(in) :: x
    real(kind=8)             :: y

    real(kind=8), parameter :: th = acosh(huge(x)), lh = log(0.5d+00)

!-------------------------------------------------------------------------------
!
    if (abs(x) < th) then
      y = log(cosh(x))
    else
      y = abs(x) + lh
    end if

!-------------------------------------------------------------------------------
!
  end function log_cosh

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