amun-code/sources/forcing.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: FORCING
!!
!!  This module provides energy injection, e.g. turbulence driving, supernova
!!  explosions, etc.
!!
!!*****************************************************************************
!
module forcing

  implicit none

! interfaces for procedure pointers
!
  abstract interface
    subroutine update_forcing_iface(dt)
      real(kind=8), intent(in)  :: dt
    end subroutine
  end interface

! pointers to the forcing methods
!
  procedure(update_forcing_iface), pointer, save :: update_forcing => null()

! flag indicating if the energy injection is enabled and if it is done in
! the Fourier space
!
  logical, save :: forcing_enabled = .false.
  logical, save :: forcing_fourier = .false.
  logical, save :: energy_profile  = .false.

! implemented injection methods
!
  integer, parameter :: injection_none     = 0
  integer, parameter :: injection_eddy     = 1
  integer, parameter :: injection_oh       = 2
  integer, parameter :: injection_alvelius = 3

! implemented driving force spectral profiles
!
  integer, parameter :: profile_normal    = 1
  integer, parameter :: profile_parabolic = 2
  integer, parameter :: profile_powerlaw  = 3

! driving force method and spectral profile
!
  integer, save :: injection_method  = injection_none
  integer, save :: injection_profile = profile_normal

! forcing parameters
!
  real(kind=8), save :: power = 1.0d+00
  real(kind=8), save :: rate  = 1.0d+03
  real(kind=8), save :: amp   = 1.0d-03
  real(kind=8), save :: del   = 1.0d-01
  real(kind=8), save :: xlo   = 9.9d+99
  real(kind=8), save :: xup   = 9.9d+99
  real(kind=8), save :: xth   = 0.0d+00
  real(kind=8), save :: xfc   = 1.0d+00
  real(kind=8), save :: ylo   = 9.9d+99
  real(kind=8), save :: yup   = 9.9d+99
  real(kind=8), save :: yth   = 0.0d+00
  real(kind=8), save :: yfc   = 1.0d+00
#if NDIMS == 3
  real(kind=8), save :: zlo   = 9.9d+99
  real(kind=8), save :: zup   = 9.9d+99
  real(kind=8), save :: zth   = 0.0d+00
  real(kind=8), save :: zfc   = 1.0d+00
#endif /* NDIMS == 3 */

! Fourier profile parameters
!
  integer     , save :: nmodes = 0
  real(kind=8), save :: fpower = 1.0d+00
  real(kind=8), save :: lscale = 5.0d-01
  real(kind=8), save :: tscale = 1.0d+00
  real(kind=8), save :: vscale = 5.0d-01
  real(kind=8), save :: ascale = 5.0d-01
  real(kind=8), save :: kf     = 2.0d+00
  real(kind=8), save :: kl     = 0.0d+00
  real(kind=8), save :: ku     = 4.0d+00
  real(kind=8), save :: kc     = 5.0d-01
  real(kind=8), save :: pidx   = - 5.0d+00 / 3.0d+00
  real(kind=8), save :: fmin   = 1.0d-08

! module parameters
!
  real(kind=8), save :: dinj = 0.0d+00
  real(kind=8), save :: einj = 0.0d+00
  real(kind=8), save :: rinj = 0.0d+00
  real(kind=8), save :: arms = 0.0d+00

! bound for the driving mode wave number
!
  integer, save :: kmax = 4

! arrays for driving mode positions and amplitudes
!
  integer     , dimension(:,:  ) , allocatable :: kmodes
  real(kind=8), dimension(:)     , allocatable :: fmodes

! vector of driving directions in Alfvelius method
!
  real(kind=8), dimension(:,:)   , allocatable :: e1vecs, e2vecs

! array for driving mode coefficients
!
  complex(kind=8), dimension(:,:), allocatable, target :: fcoefs

! velocity Fourier coefficients in Alfvelius method
!
  complex(kind=8), dimension(:,:), allocatable :: vcoefs

! by default everything is private
!
  private

! declare public subroutines
!
  public :: initialize_forcing, finalize_forcing, print_forcing
  public :: update_forcing
  public :: forcing_enabled, einj, rinj, arms, nmodes, fcoefs

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
  contains
!
!===============================================================================
!!
!!***  PUBLIC SUBROUTINES  *****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine INITIALIZE_FORCING:
! -----------------------------
!
!   Subroutine initializes module FORCING.
!
!   Arguments:
!
!     verbose - flag determining if the subroutine should be verbose;
!     status  - return flag of the procedure execution status;
!
!===============================================================================
!
  subroutine initialize_forcing(verbose, status)

    use constants , only : pi2
    use helpers   , only : print_message
    use parameters, only : get_parameter
    use random    , only : randuni, randnorz

    implicit none

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

    character(len=32) :: wfmt = "(1x,'WARNING:',1x,a)"
    character(len=64) :: injection      = "none"
    character(len=64) :: profile_type   = "gauss"
    character(len=64) :: profile_energy = "off"
    integer           :: i, j, l, k2
#if NDIMS == 3
    integer           :: k = 0
#endif /* NDIMS == 3 */
    real(kind=8)      :: kl2, ku2, kv2, kv
    real(kind=8)      :: fa, fi, uu
#if NDIMS == 3
    real(kind=8)      :: phi
#endif /* NDIMS == 3 */

    real(kind=8), dimension(NDIMS) :: u, v

    integer, dimension(:), allocatable :: kcount

    character(len=*), parameter :: loc = 'FORCING::initialize_forcing()'

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

! obtain the chosen injection method
!
    call get_parameter("driving_method", injection)

! get the injection profile parameters
!
    call get_parameter("driving_x_lower", xlo)
    call get_parameter("driving_x_upper", xup)
    call get_parameter("driving_y_lower", ylo)
    call get_parameter("driving_y_upper", yup)
#if NDIMS == 3
    call get_parameter("driving_z_lower", zlo)
    call get_parameter("driving_z_upper", zup)
#endif /* NDIMS == 3 */
    xth = xup - xlo
    if (xth > 0.0d+00) then
      xfc = 2.5d-01 * pi2 / xth
    else
      xfc = 1.0d+00
    end if
    yth = yup - ylo
    if (yth > 0.0d+00) then
      yfc = 2.5d-01 * pi2 / yth
    else
      yfc = 1.0d+00
    end if
#if NDIMS == 3
    zth = zup - zlo
    if (zth > 0.0d+00) then
      zfc = 2.5d-01 * pi2 / zth
    else
      zfc = 1.0d+00
    end if
#endif /* NDIMS == 3 */

! select the energy injection method
!
    select case(trim(injection))
    case("eddies", "eddy")
      forcing_enabled  = .true.
      injection_method = injection_eddy
      update_forcing   => update_forcing_eddies

      call get_parameter("injection_power", power)
      call get_parameter("injection_rate" , rate )
      call get_parameter("eddy_amplitude" , amp  )
      call get_parameter("eddy_width"     , del  )

    case("Ornstein–Uhlenbeck", "ornstein–uhlenbeck", "oh")
      forcing_enabled  = .true.
      forcing_fourier  = .true.
      injection_method = injection_oh
      update_forcing   => update_forcing_oh

! get the driving power
!
      call get_parameter("driving_timescale", tscale)
      if (tscale <= 0.0d+00) then
        if (verbose) then
          write(*,*)
          write(*,wfmt) "'driving_timescale' must be larger than zero!"
          write(*,wfmt) "Resetting it to the default value: 1.0!"
        end if
        tscale = 1.0d+00
      end if

    case("Alvelius", "alvelius", "al")
      forcing_enabled  = .true.
      forcing_fourier  = .true.
      injection_method = injection_alvelius
      update_forcing   => update_forcing_alvelius

! get the driving power
!
      call get_parameter("driving_power", fpower)
      if (fpower <= 0.0d+00) then
        if (verbose) then
          write(*,*)
          write(*,wfmt) "'driving_power' must be larger than zero!"
          write(*,wfmt) "Resetting it to the default value: 1.0!"
        end if
        fpower = 1.0d+00
      end if

    case default
      injection_method = injection_none
      if (verbose .and. injection /= "none") then
        write(*,*)
        write(*,wfmt) "Unknown injection method! Forcing disabled!"
      end if
    end select

! initialization of the Fourier space driving
!
    if (forcing_fourier) then

      call get_parameter("driving_wavenumber"   , kf)
      call get_parameter("driving_profile_width", kc)
      call get_parameter("driving_profile_lower", kl)
      call get_parameter("driving_profile_upper", ku)
      call get_parameter("amplitude_threshold"  , fmin)

      if (kf < 1.0d+00) then
        if (verbose) then
          write(*,*)
          write(*,wfmt) "'driving_wavenumber' must be larger than 1.0!"
          write(*,wfmt) "Resetting it to the default value: 2.0!"
        end if
        kf = 2.0d+00
      end if
      if (kc <= 0.0d+00) then
        if (verbose) then
          write(*,*)
          write(*,wfmt) "'driving_profile_width' must be larger than 0.0!"
          write(*,wfmt) "Resetting it to the default value: 1.0!"
        end if
        kc = 1.0d+00
      end if
      if (kc > kf) then
        if (verbose) then
          write(*,*)
          write(*,wfmt) "'driving_profile_width' must be smaller than " //     &
                        "'driving_wavenumber'!"
          write(*,wfmt) "Resetting it to the half of 'driving_wavenumber'!"
        end if
        kc = 5.0d-01 * kf
      end if
      if (kl >= kf) then
        if (verbose) then
          write(*,*)
          write(*,wfmt) "'driving_profile_lower' must be smaller than " //     &
                        "'driving_wavenumber'!"
          write(*,wfmt) "Setting it to 'driving_wavenumber' " //               &
                        "- 'driving_profile_width'!"
        end if
        kl = kf - kc
      end if
      if (ku <= kf) then
        if (verbose) then
          write(*,*)
          write(*,wfmt) "'driving_profile_upper' must be larger than " //      &
                        "'driving_wavenumber'!"
          write(*,wfmt) "Setting it to 'driving_wavenumber' " //               &
                        "+ 'driving_profile_width'!"
        end if
        ku = kf + kc
      end if
      fmin = max(1.0d-08, fmin)
      kl2  = kl**2
      ku2  = ku**2

! get the upper bound for the wave number
!
      kmax = ceiling(ku)

! allocate arrays to count the modes with the same wave number
!
      allocate(kcount(NDIMS * kmax**2), stat = status)

      if (status == 0) then

! initialize mode counts
!
        kcount(:) =  0

! count the modes with the same k, get the minimum sufficient number of
! driving modes
!
#if NDIMS == 3
        do i = 1, kmax
          do j = -kmax, kmax
            do k = -kmax, kmax
              k2  = i**2 + j**2 + k**2
              kv2 = 1.0d+00 * k2
              if (kl2 <= kv2 .and. kv2 <= ku2) then
                nmodes     = nmodes + 1
                kcount(k2) = kcount(k2) + 2
              end if
            end do
          end do
        end do
        do j = 1, kmax
          do k = -kmax, kmax
            k2  = j**2 + k**2
            kv2 = 1.0d+00 * k2
            if (kl2 <= kv2 .and. kv2 <= ku2) then
              nmodes     = nmodes + 1
              kcount(k2) = kcount(k2) + 2
            end if
          end do
        end do
        do k = 1, kmax
          k2  = k**2
          kv2 = 1.0d+00 * k2
          if (kl2 <= kv2 .and. kv2 <= ku2) then
            nmodes     = nmodes + 1
            kcount(k2) = kcount(k2) + 2
          end if
        end do
#else /* NDIMS == 3 */
        do i = 1, kmax
          do j = -kmax, kmax
            k2  = i**2 + j**2
            kv2 = 1.0d+00 * k2
            if (kl2 <= kv2 .and. kv2 <= ku2) then
              nmodes     = nmodes + 1
              kcount(k2) = kcount(k2) + 2
            end if
          end do
        end do
        do j = 1, kmax
          k2  = j**2
          kv2 = 1.0d+00 * k2
          if (kl2 <= kv2 .and. kv2 <= ku2) then
            nmodes     = nmodes + 1
            kcount(k2) = kcount(k2) + 2
          end if
        end do
#endif /* NDIMS == 3 */

! allocate arrays for driving modes (positions, amplitudes and coefficients)
!
        allocate(e1vecs(nmodes,NDIMS), e2vecs(nmodes,NDIMS),                   &
                 kmodes(nmodes,NDIMS), fmodes(nmodes),                         &
                 fcoefs(nmodes,NDIMS), stat = status)

        if (status == 0) then

! prepare wave vectors
!
          l = 0
#if NDIMS == 3
          do i = 1, kmax
            do j = -kmax, kmax
              do k = -kmax, kmax
                k2  = i**2 + j**2 + k**2
                kv2 = 1.0d+00 * k2
                if (kl2 <= kv2 .and. kv2 <= ku2) then
                  l = l + 1
                  kmodes(l,1)   = i
                  kmodes(l,2)   = j
                  kmodes(l,3)   = k
                end if
              end do
            end do
          end do
          do j = 1, kmax
            do k = -kmax, kmax
              k2  = j**2 + k**2
              kv2 = 1.0d+00 * k2
              if (kl2 <= kv2 .and. kv2 <= ku2) then
                l = l + 1
                kmodes(l,1)   = 0
                kmodes(l,2)   = j
                kmodes(l,3)   = k
              end if
            end do
          end do
          do k = 1, kmax
            k2  = k**2
            kv2 = 1.0d+00 * k2
            if (kl2 <= kv2 .and. kv2 <= ku2) then
              l = l + 1
              kmodes(l,1)   = 0
              kmodes(l,2)   = 0
              kmodes(l,3)   = k
            end if
          end do
#else /* NDIMS == 3 */
          do i = 1, kmax
            do j = -kmax, kmax
              k2  = i**2 + j**2
              kv2 = 1.0d+00 * k2
              if (kl2 <= kv2 .and. kv2 <= ku2) then
                l = l + 1
                kmodes(l,1) = i
                kmodes(l,2) = j
              end if
            end do
          end do
          do j = 1, kmax
            k2  = j**2
            kv2 = 1.0d+00 * k2
            if (kl2 <= kv2 .and. kv2 <= ku2) then
              l = l + 1
              kmodes(l,1) = 0
              kmodes(l,2) = j
            end if
          end do
#endif /* NDIMS == 3 */

! generate random vectors, perpendicular to the wave vector
!
          do l = 1, nmodes

            u(:) = kmodes(l,:)
#if NDIMS == 3
            if (abs(kmodes(l,3)) < abs(kmodes(l,1))) then
              v(:) = (/ kmodes(l,2), - kmodes(l,1), 0 /)
            else
              v(:) = (/ 0, - kmodes(l,3), kmodes(l,2) /)
            end if
#else /* NDIMS == 3 */
            v(:) = (/ kmodes(l,2), - kmodes(l,1) /)
#endif /* NDIMS == 3 */

            uu = dot_product(v(:), v(:))
            if (uu > 0.0d+00) then

              e1vecs(l,:) = v(:)
#if NDIMS == 3
              e2vecs(l,1) = u(2) * v(3) - u(3) * v(2)
              e2vecs(l,2) = u(3) * v(1) - u(1) * v(3)
              e2vecs(l,3) = u(1) * v(2) - u(2) * v(1)
#endif /* NDIMS == 3 */

              e1vecs(l,:) = e1vecs(l,:) / sqrt(sum(e1vecs(l,:)**2))
#if NDIMS == 3
              e2vecs(l,:) = e2vecs(l,:) / sqrt(sum(e2vecs(l,:)**2))
#else /* NDIMS == 3 */
              e2vecs(l,:) = 0.0d+00
#endif /* NDIMS == 3 */

            else
              call print_message(loc, "v(:) is ZERO!")
            end if

          end do

! determine if the profile is set by energy or amplitude
!
          call get_parameter("profile_by_energy", profile_energy)
          select case(trim(profile_energy))
          case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
            energy_profile = .true.
          case default
            energy_profile = .false.
            kcount(:)      = 1
          end select

! initialize driving mode amplitudes depending on the chosen profile
!
          call get_parameter("driving_profile", profile_type)
          select case(trim(profile_type))

! --- power-law profile ---
!
          case("power-law", "powerlaw", "band")

            injection_profile = profile_powerlaw

! get the power index of the power-law distribution
!
            call get_parameter("powerlaw_index", pidx)

! calculate amplitudes of driving modes
!
            fi = 0.0d+00
            do l = 1, nmodes
              k2        = sum(kmodes(l,:)**2)
              kv        = sqrt(1.0d+00 * k2)
              fa        = kv**pidx / kcount(k2)
              fi        = fi + fa
              fmodes(l) = sqrt(fa)
            end do
            fmodes(:) = fmodes(:) / sqrt(fi)

! --- normal distribution profile ---
!
          case("gauss", "gaussian", "normal")

            injection_profile = profile_normal

! calculate amplitudes of driving modes
!
            fi = 0.0d+00
            do l = 1, nmodes
              k2        = sum(kmodes(l,:)**2)
              kv        = sqrt(1.0d+00 * k2)
              fa        = exp(- 0.5d+00 * ((kv - kf) / kc)**2) / kcount(k2)
              fi        = fi + fa
              fmodes(l) = sqrt(fa)
            end do
            fmodes(:) = fmodes(:) / sqrt(fi)

! --- parabolic profile ---
!
          case("parabolic")

            injection_profile = profile_parabolic

! correct the injection scale and width
!
            kf = 0.5d+00 * (ku + kl)
            kc = 0.5d+00 * (ku - kl)

! calculate amplitudes of driving modes
!
            fi = 0.0d+00
            do l = 1, nmodes
              k2        = sum(kmodes(l,:)**2)
              kv        = sqrt(1.0d+00 * k2)
              fa        = 4.0d+00 * (kv - kl) * (ku - kv) / (ku - kl)**2
              fa        = max(0.0d+00, fa)**2 / kcount(k2)
              fi        = fi + fa
              fmodes(l) = sqrt(fa)
            end do
            fmodes(:) = fmodes(:) / sqrt(fi)

! --- normal distribution profile ---
!
          case default

! warn if the specified profile is not implemented and the normal one is used
!
            if (verbose .and. .not. (profile_type == "normal" .or.             &
                                     profile_type == "gauss"  .or.             &
                                     profile_type == "gaussian")) then
              write(*,*)
              write(*,wfmt) "Unknown spectral profile of driving!"
              write(*,wfmt) "Available profiles: " //                          &
                            "'power-law', 'normal', 'parabolic'."
              write(*,wfmt) "Using default normal distribution profile " //    &
                            "'gaussian'."

            end if

            injection_profile = profile_normal

! calculate amplitudes of driving modes
!
            fi = 0.0d+00
            do l = 1, nmodes
              k2        = sum(kmodes(l,:)**2)
              kv        = sqrt(1.0d+00 * k2)
              fa        = exp(- 0.5d+00 * ((kv - kf) / kc)**2) / kcount(k2)
              fi        = fi + fa
              fmodes(l) = sqrt(fa)
            end do
            do l = 1, nmodes
              k2        = sum(kmodes(l,:)**2)
              fmodes(l) = fmodes(l) / sqrt(fi)
            end do

          end select

! === Ornstein–Uhlenbeck driving ===
!
          if (injection_method == injection_oh) then

! get the characteristic driving length, velocity, acceleration scales, and
! determine the estimated driving power (since the main injection is through
! the forcing-velocity correlation, it may not represent the actual
! instantenous injection rate)
!
            lscale = 1.0d+00 / kf
            vscale = lscale / tscale
            ascale = sqrt(2.0d+00) * lscale / tscale**2 / sqrt(tscale)
            fpower = vscale**2 / tscale

! normalize the driving amplitude profile to the driving power
!
            fmodes(:) = ascale * fmodes(:)

! generate the initial driving mode coefficients
!
            do l = 1, nmodes

#if NDIMS == 3
              phi         = pi2 * randuni()
              fcoefs(l,:) = fmodes(l) * (e1vecs(l,:) * cos(phi)                &
                                       + e2vecs(l,:) * sin(phi)) * randnorz()
#else /* NDIMS == 3 */
              fcoefs(l,:) = fmodes(l) *  e1vecs(l,:) * randnorz()
#endif /* NDIMS == 3 */

            end do ! l = 1, nmodes

! calculate the rms of acceleration
!
            arms = sqrt(5.0d-01 * sum(real(fcoefs(:,:))**2                     &
                                   + aimag(fcoefs(:,:))**2))

          end if ! injection = 'OH'

! === Alvelius driving ===
!
          if (injection_method == injection_alvelius) then

! get the characteristic driving length, velocity, time, and acceleration scales
!
            lscale = 1.0d+00 / kf
            tscale = sqrt(lscale / sqrt(2.0d+00 * fpower))
            vscale = lscale / tscale
            ascale = vscale / tscale

! normalize the driving amplitude profile to the driving power
!
            fmodes(:) = sqrt(2.0d+00) * ascale * fmodes(:)

! reset the driving coefficients and calculate the rms of acceleration
!
            fcoefs(:,:) = cmplx(0.0d+00, 0.0d+00, kind=8)
            arms        = 0.0d+00

! allocate vectors of velocity Fourier coefficients
!
            allocate(vcoefs(nmodes,NDIMS), stat = status)
            if (status /= 0) &
              call print_message(loc, "Cannot allocate vcoefs(:,:)!")

          end if ! injection = 'Alvelius'

        end if ! allocation of driving mode arrays
      end if ! allocation of kcount

! deallocate mode counter
!
      if (allocated(kcount)) deallocate(kcount)

    end if ! Fourier forcing

!-------------------------------------------------------------------------------
!
  end subroutine initialize_forcing
!
!===============================================================================
!
! subroutine FINALIZE_FORCING:
! ---------------------------
!
!   Subroutine releases memory used by the module variables.
!
!   Arguments:
!
!     status  - return flag of the procedure execution status;
!
!===============================================================================
!
  subroutine finalize_forcing(status)

    implicit none

    integer, intent(out) :: status

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

    if (allocated(kmodes)) deallocate(kmodes)
    if (allocated(fmodes)) deallocate(fmodes)
    if (allocated(fcoefs)) deallocate(fcoefs)
    if (allocated(e1vecs)) deallocate(e1vecs)
    if (allocated(e2vecs)) deallocate(e2vecs)
    if (allocated(vcoefs)) deallocate(vcoefs)

!-------------------------------------------------------------------------------
!
  end subroutine finalize_forcing
!
!===============================================================================
!
! subroutine PRINT_FORCING:
! ------------------------
!
!   Subroutine prints module parameters.
!
!   Arguments:
!
!     verbose - a logical flag turning the information printing;
!
!===============================================================================
!
  subroutine print_forcing(verbose)

    use coordinates, only : periodic
    use helpers    , only : print_section, print_parameter

    implicit none

    logical, intent(in) :: verbose

    character(len=64) :: injection    = "none"
    character(len=64) :: profile_type = "gaussian"

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

      select case(injection_method)
      case(injection_eddy)
        injection = "random eddies"
      case(injection_oh)
        injection = "Ornstein–Uhlenbeck"
      case(injection_alvelius)
        injection = "Alvelius"
      case default
        injection = "none"
      end select

      call print_section(verbose, "Forcing")
      call print_parameter(verbose, "method", injection)

      select case(injection_method)
      case(injection_eddy)
        call print_parameter(verbose, "injection power", power)
        call print_parameter(verbose, "injection rate" , rate)
        call print_parameter(verbose, "eddy amplitude" , amp)
        call print_parameter(verbose, "eddy width"     , del)
      end select

      if (forcing_fourier) then
        if (injection_method == injection_oh) then
          call print_parameter(verbose, "power (estimated)", fpower)
        end if
        if (injection_method == injection_alvelius) then
          call print_parameter(verbose, "power"            , fpower)
        end if
        call print_parameter(verbose, "length scale"      , lscale)
        call print_parameter(verbose, "time scale"        , tscale)
        call print_parameter(verbose, "velocity scale"    , vscale)
        call print_parameter(verbose, "acceleration scale", ascale)
        select case(injection_profile)
        case(profile_powerlaw)
          profile_type = "power-law"
        case(profile_normal)
          profile_type = "gaussian"
        case(profile_parabolic)
          profile_type = "parabolic"
        case default
          profile_type = "unknown"
        end select

        call print_parameter(verbose, "spectrum profile", profile_type)
        if (energy_profile) then
          call print_parameter(verbose, "profile of", "energy")
        else
          call print_parameter(verbose, "profile of", "amplitude")
        end if
        call print_parameter(verbose, "number of driving modes", nmodes)
        call print_parameter(verbose, "driving wavenumber"     , kf)
        call print_parameter(verbose, "lower wavenumber limit" , kl)
        call print_parameter(verbose, "upper wavenumber limit" , ku)
        if (injection_profile == profile_powerlaw) then
          call print_parameter(verbose, "power-law index", pidx)
        end if
        if (injection_profile == profile_normal) then
          call print_parameter(verbose, "spectrum profile width", kc)
        end if
        call print_parameter(verbose, "amplitude threshold", fmin)
      end if

      if (.not. periodic(1)) &
        call print_parameter(verbose, "profile X transition between", xlo, xup)
      if (.not. periodic(2)) &
        call print_parameter(verbose, "profile Y transition between", ylo, yup)
#if NDIMS == 3
      if (.not. periodic(3)) &
        call print_parameter(verbose, "profile Z transition between", zlo, zup)
#endif /* NDIMS == 3 */

    end if

!-------------------------------------------------------------------------------
!
  end subroutine print_forcing
!
!===============================================================================
!
! subroutine UPDATE_FORCING_EDDIES:
! --------------------------------
!
!   Subroutine adds the energy injection terms.
!
!   Arguments:
!
!     dt - time increment;
!
!===============================================================================
!
  subroutine update_forcing_eddies(dt)

    use coordinates, only : xmin, ymin, zmin, xlen, ylen, zlen
    use random     , only : randuni, randsym

    implicit none

    real(kind=8), intent(in) :: dt

    integer                    :: ni, n
#if NDIMS == 3
    real(kind=8)               :: tmp
#endif /* NDIMS == 3 */
    real(kind=8), dimension(3) :: xp, ap

!-------------------------------------------------------------------------------
!
! calculate the number of eddies to be injected during this interval
!
    dinj = dinj + rate * dt
    ni   = int(floor(dinj))
    dinj = dinj - 1.0d+00 * ni

! inject the required number of eddies
!
    if (ni > 0) then

      do n = 1, ni

! get random position
!
        xp(1) = xmin + xlen * randuni()
        xp(2) = ymin + ylen * randuni()
        xp(3) = zmin + zlen * randuni()

! get random orientation
!
#if NDIMS == 3
        tmp = 0.0d+00
        do while(tmp <= tiny(0.0d+00))
          ap(1) = randsym()
          ap(2) = randsym()
          ap(3) = randsym()
          tmp   = sqrt(ap(1)**2 + ap(2)**2 + ap(3)**2)
        end do
        ap(:) = amp * (ap(:) / tmp) / del
#else /* NDIMS == 3 */
        ap(:) = sign(1.0d+00, randsym()) * (/ 0.0d+00, 0.0d+00, amp / del /)
#endif /* NDIMS == 3 */
        ap(:) = ap(:) * dt

! iterate over data blocks and add forcing components
!
        call inject_eddy(xp(:), ap(:))

      end do

    end if

!-------------------------------------------------------------------------------
!
  end subroutine update_forcing_eddies
!
!===============================================================================
!
! subroutine UPDATE_FORCING_OH:
! ----------------------------
!
!   Subroutine adds the energy injection terms for Ornstein–Uhlenbeck driving.
!
!   Arguments:
!
!     dt - time increment;
!
!===============================================================================
!
  subroutine update_forcing_oh(dt)

#if NDIMS == 3
    use constants, only : pi2
#endif /* NDIMS == 3 */
    use random   , only : randuni, randnorz

    implicit none

    real(kind=8), intent(in) :: dt

    integer      :: l
    real(kind=8) :: acoeff, dcoeff
#if NDIMS == 3
    real(kind=8) :: phi
#endif /* NDIMS == 3 */
    real(kind=8) :: dinj

    complex(kind=8), dimension(NDIMS) :: finj

!-------------------------------------------------------------------------------
!
! calculate drift and diffusion coefficients
!
    acoeff = exp(- dt / tscale)
    dcoeff = sqrt(1.0d+00 - acoeff**2)

! iterate over driving modes
!
    do l = 1, nmodes

! generate gaussian random vector along the driving modes
!
#if NDIMS == 3
      phi         = pi2 * randuni()
      finj(:)     = fmodes(l) * (e1vecs(l,:) * cos(phi)                        &
                               + e2vecs(l,:) * sin(phi)) * randnorz()
#else /* NDIMS == 3 */
      finj(:)     = fmodes(l) *  e1vecs(l,:) * randnorz()
#endif /* NDIMS == 3 */

! advance the driving force modes
!
      fcoefs(l,:) = acoeff * fcoefs(l,:) + dcoeff * finj(:)

    end do

! calculate the rms of acceleration
!
    arms = sqrt(5.0d-01 * sum(real(fcoefs(:,:))**2 + aimag(fcoefs(:,:))**2))

! store previously injected energy
!
    dinj = einj

! inject driving modes
!
    call inject_fmodes(dt)

! calculate the energy injected during this step
!
    dinj = einj - dinj

! calculate the injection rate
!
    rinj = dinj / dt

!-------------------------------------------------------------------------------
!
  end subroutine update_forcing_oh
!
!===============================================================================
!
! subroutine UPDATE_FORCING_ALVELIUS:
! ----------------------------------
!
!   Subroutine adds the energy injection terms for Alvelius driving.
!
!   Arguments:
!
!     dt - time increment;
!
!===============================================================================
!
  subroutine update_forcing_alvelius(dt)

    use constants, only : pi2
    use random   , only : randuni

    implicit none

    real(kind=8), intent(in) :: dt

    integer         :: l
    real(kind=8)    :: th1, dinj, sqdt
#if NDIMS == 3
    real(kind=8)    :: th2, phi, psi, ga, gb
#endif /* NDIMS == 3 */
    complex(kind=8) :: aran
#if NDIMS == 3
    complex(kind=8) :: bran
    complex(kind=8) :: xi1, xi2
#endif /* NDIMS == 3 */

!-------------------------------------------------------------------------------
!
! determine velocify coefficients in Fourier space
!
    call get_vcoefs()

! square root of dt
!
    sqdt = sqrt(dt)

! iterate over driving modes
!
    do l = 1, nmodes

#if NDIMS == 3
! get xi1 and xi2
!
      xi1 = dot_product(vcoefs(l,:), e1vecs(l,:))
      xi2 = dot_product(vcoefs(l,:), e2vecs(l,:))

! get random angles phi and psi and determine th1 and th2 to reduce
! velocity-force correlation
!
      phi = pi2 * randuni()
      ga  = sin(phi)
      gb  = cos(phi)
      psi = pi2 * randuni()
      th1 = - ga * aimag(xi1) + gb * (sin(psi) * real(xi2)                     &
                                    - cos(psi) * aimag(xi2))
      if (abs(th1) > epsilon(0.0d+00)) then
        th1 = atan((ga * real(xi1) + gb * (sin(psi) * aimag(xi2)               &
                                         + cos(psi) * real(xi2))) / th1)
      end if
      th2 = psi + th1

! calculate amplitudes
!
      aran = cmplx(cos(th1), sin(th1), kind=8) * ga
      bran = cmplx(cos(th2), sin(th2), kind=8) * gb

! calculate driving mode coefficient
!
      fcoefs(l,:) = fmodes(l) * (aran * e1vecs(l,:)                            &
                               + bran * e2vecs(l,:)) / sqdt
#else /* NDIMS == 3 */
! get random phase
!
      th1 = pi2 * randuni()
      aran = cmplx(cos(th1), sin(th1), kind=8)

! calculate driving mode coefficient
!
      fcoefs(l,:) = fmodes(l) * aran * e1vecs(l,:) / sqdt
#endif /* NDIMS == 3 */

    end do

! calculate the rms of acceleration
!
    arms = sqrt(5.0d-01 * sum(real(fcoefs(:,:))**2 + aimag(fcoefs(:,:))**2))

! store previously injected energy
!
    dinj = einj

! inject driving modes
!
    call inject_fmodes(dt)

! calculate the energy injected during this step
!
    dinj = einj - dinj

! calculate the injection rate
!
    rinj = dinj / dt

!-------------------------------------------------------------------------------
!
  end subroutine update_forcing_alvelius
!
!===============================================================================
!!
!!***  PRIVATE SUBROUTINES  ****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine INJECT_EDDY:
! ----------------------
!
!   Subroutine injects one eddy to the domain at given position.
!
!   Arguments:
!
!     xp, ap - eddy position and amplitude vectors;
!
!===============================================================================
!
  subroutine inject_eddy(xp, ap)

    use blocks, only : block_data, list_data

    implicit none

    real(kind=8), dimension(3), intent(in) :: xp, ap

    type(block_data), pointer :: pdata

!-------------------------------------------------------------------------------
!
! assign pdata with the first block on the data block list
!
    pdata => list_data

! iterate over all data blocks
!
    do while (associated(pdata))

! inject eddy into the current block
!
      call inject_eddy_block(pdata, xp(:), ap(:))

! assign pdata to the next block
!
      pdata => pdata%next

    end do ! over data blocks

!-------------------------------------------------------------------------------
!
  end subroutine inject_eddy
!
!===============================================================================
!
! subroutine INJECT_EDDY_BLOCK:
! ----------------------------
!
!   Subroutine injects one eddy to the local block.
!
!   Arguments:
!
!     xp, ap - eddy position and amplitude vectors;
!
!===============================================================================
!
  subroutine inject_eddy_block(pdata, xp, ap)

    use blocks     , only : block_data
    use coordinates, only : nm => bcells
    use coordinates, only : ax, ay, xlen, ylen
#if NDIMS == 3
    use coordinates, only : az, zlen
#endif /* NDIMS == 3 */
    use coordinates, only : periodic
    use equations  , only : idn, imx, imy, imz, ien

    implicit none

    type(block_data), pointer , intent(inout) :: pdata
    real(kind=8), dimension(3), intent(in)    :: xp, ap

    integer      :: i, j, k
    real(kind=8) :: x2, y2, r2
#if NDIMS == 3
    real(kind=8) :: z2
#endif /* NDIMS == 3 */
    real(kind=8) :: fx, fy, fz, fp, e1, e2

    real(kind=8), dimension(nm)    :: x, y, z

!-------------------------------------------------------------------------------
!
! prepare block coordinates
!
    if (periodic(1)) then
      fx   = 0.5d+00 * xlen
      x(:) = pdata%meta%xmin + ax(pdata%meta%level,:) - xp(1)
      x(:) = abs(x(:) / fx)
      x(:) = (abs(1.0d+00 - abs(1.0d+00 - abs(x(:))))) * fx / del
    else
      x(:) = (pdata%meta%xmin + ax(pdata%meta%level,:) - xp(1)) / del
    end if
    if (periodic(2)) then
      fy   = 0.5d+00 * ylen
      y(:) = pdata%meta%ymin + ay(pdata%meta%level,:) - xp(2)
      y(:) = abs(y(:) / fy)
      y(:) = (abs(1.0d+00 - abs(1.0d+00 - abs(y(:))))) * fy / del
    else
      y(:) = (pdata%meta%ymin + ay(pdata%meta%level,:) - xp(2)) / del
    end if
#if NDIMS == 3
    if (periodic(3)) then
      fz   = 0.5d+00 * zlen
      z(:) = pdata%meta%zmin + az(pdata%meta%level,:) - xp(3)
      z(:) = abs(z(:) / fz)
      z(:) = (abs(1.0d+00 - abs(1.0d+00 - abs(z(:))))) * fz / del
    else
      z(:) = (pdata%meta%zmin + az(pdata%meta%level,:) - xp(3)) / del
    end if
#else /* NDIMS == 3 */
    z(:) = 0.0d+00
#endif /* NDIMS == 3 */

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

! iterate over the block coordinates
!
    if (ien > 0) then
#if NDIMS == 3
      do k = 1, nm
        z2 = z(k)**2
#endif /* NDIMS == 3 */
        do j = 1, nm
          y2 = y(j)**2
          do i = 1, nm
            x2 = x(i)**2
#if NDIMS == 3
            r2 = x2 + y2 + z2
#else /* NDIMS == 3 */
            r2 = x2 + y2
#endif /* NDIMS == 3 */

            fp = pdata%u(idn,i,j,k) * exp(- 0.5d+00 * r2)

#if NDIMS == 3
            fx = (- ap(3) * y(j) + ap(2) * z(k)) * fp
            fy = (- ap(1) * z(k) + ap(3) * x(i)) * fp
            fz = (- ap(2) * x(i) + ap(1) * y(j)) * fp
#else /* NDIMS == 3 */
            fx = (- ap(3) * y(j)) * fp
            fy = (+ ap(3) * x(i)) * fp
            fz = 0.0d+00
#endif /* NDIMS == 3 */

            e1 = 0.5d+00 * sum(pdata%u(imx:imx,i,j,k)**2) / pdata%u(idn,i,j,k)

            pdata%u(imx,i,j,k) = pdata%u(imx,i,j,k) + fx
            pdata%u(imy,i,j,k) = pdata%u(imy,i,j,k) + fy
            pdata%u(imz,i,j,k) = pdata%u(imz,i,j,k) + fz

            e2 = 0.5d+00 * sum(pdata%u(imx:imz,i,j,k)**2) / pdata%u(idn,i,j,k)

            pdata%u(ien,i,j,k) = pdata%u(ien,i,j,k) + (e2 - e1)

          end do
        end do
#if NDIMS == 3
      end do
#endif /* NDIMS == 3 */
    else
#if NDIMS == 3
      do k = 1, nm
        z2 = z(k)**2
#endif /* NDIMS == 3 */
        do j = 1, nm
          y2 = y(j)**2
          do i = 1, nm
            x2 = x(i)**2
#if NDIMS == 3
            r2 = x2 + y2 + z2
#else /* NDIMS == 3 */
            r2 = x2 + y2
#endif /* NDIMS == 3 */

            fp = pdata%u(idn,i,j,k) * exp(- 0.5d+00 * r2)

#if NDIMS == 3
            fx = (- ap(3) * y(j) + ap(2) * z(k)) * fp
            fy = (- ap(1) * z(k) + ap(3) * x(i)) * fp
            fz = (- ap(2) * x(i) + ap(1) * y(j)) * fp
#else /* NDIMS == 3 */
            fx = (- ap(3) * y(j)) * fp
            fy = (+ ap(3) * x(i)) * fp
            fz = 0.0d+00
#endif /* NDIMS == 3 */

            pdata%u(imx,i,j,k) = pdata%u(imx,i,j,k) + fx
            pdata%u(imy,i,j,k) = pdata%u(imy,i,j,k) + fy
            pdata%u(imz,i,j,k) = pdata%u(imz,i,j,k) + fz

          end do
        end do
#if NDIMS == 3
      end do
#endif /* NDIMS == 3 */
    end if

!-------------------------------------------------------------------------------
!
  end subroutine inject_eddy_block
!
!===============================================================================
!
! subroutine INJECT_FMODES:
! ------------------------
!
!   Subroutine injects Fourier driving modes.
!
!   Arguments:
!
!     dt - the time increment;
!
!===============================================================================
!
  subroutine inject_fmodes(dt)

    use blocks, only : block_data, list_data

    implicit none

    real(kind=8), intent(in) :: dt

    type(block_data), pointer :: pdata

!-------------------------------------------------------------------------------
!
! assign pdata with the first block on the data block list
!
    pdata => list_data

! iterate over all data blocks
!
    do while (associated(pdata))

! inject eddy into the current block
!
      call inject_fmodes_block(pdata, dt)

! assign pdata to the next block
!
      pdata => pdata%next

    end do ! over data blocks

!-------------------------------------------------------------------------------
!
  end subroutine inject_fmodes
!
!===============================================================================
!
! subroutine INJECT_FMODES_BLOCK:
! ------------------------------
!
!   Subroutine injects Fourier driving modes to the local block.
!
!   Arguments:
!
!     pdata - a pointer to the data block;
!     dt    - the time increment;
!
!===============================================================================
!
  subroutine inject_fmodes_block(pdata, dt)

    use blocks     , only : block_data
    use constants  , only : pi2
    use coordinates, only : nn => bcells, nb, ne
    use coordinates, only : ax, ay, advol
#if NDIMS == 3
    use coordinates, only : az
#endif /* NDIMS == 3 */
    use coordinates, only : periodic
    use equations  , only : idn, imx, imy, imz, ien
    use helpers    , only : print_message
    use workspace  , only : resize_workspace, work, work_in_use

    implicit none

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

    logical, save :: first = .true.

    integer      :: i, j, k, l, n, status
    real(kind=8) :: cs, sn
#if NDIMS == 3
    real(kind=8) :: tt
#endif /* NDIMS == 3 */
    real(kind=8) :: dvol

    real(kind=8), dimension(nn):: x, y, z
    real(kind=8), dimension(nn):: kx, ky, snkx, snky, cskx, csky
#if NDIMS == 3
    real(kind=8), dimension(nn):: kz, snkz, cskz
#endif /* NDIMS == 3 */

    real(kind=8), dimension(:,:,:,:), pointer, save :: acc
    real(kind=8), dimension(:,:,:)  , pointer, save :: den

    integer, save :: nt
!$  integer :: omp_get_thread_num
!$omp threadprivate(first, nt, acc, den)

    character(len=*), parameter :: loc = 'FORCING:inject_fmodes_block()'

!-------------------------------------------------------------------------------
!
    nt = 0
!$  nt = omp_get_thread_num()
    if (first) then
      i = 3 * nn**NDIMS
      j = 4 * 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
      acc(1:3,1:nn,1:nn,1:nn) => work(  1:i,nt)
      den(1:nn,1:nn,1:nn)     => work(i+1:j,nt)
#else /* NDIMS == 3 */
      acc(1:2,1:nn,1:nn,1: 1) => work(  1:i,nt)
      den(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 */

    x(:) = - pi2 * (pdata%meta%xmin + ax(pdata%meta%level,:))
    y(:) = - pi2 * (pdata%meta%ymin + ay(pdata%meta%level,:))
#if NDIMS == 3
    z(:) = - pi2 * (pdata%meta%zmin + az(pdata%meta%level,:))
#else /* NDIMS == 3 */
    z(:) = 0.0d+00
#endif /* NDIMS == 3 */
    dvol = advol(pdata%meta%level)

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

    work_in_use(nt) = .true.

    acc(1:NDIMS,:,:,:) = 0.0d+00

! iterate over driving modes and compute the acceleration in the real space
!
    do l = 1, nmodes

      if (fmodes(l) > fmin) then

        kx(:) = kmodes(l,1) * x(:)
        ky(:) = kmodes(l,2) * y(:)
#if NDIMS == 3
        kz(:) = kmodes(l,3) * z(:)
#endif /* NDIMS == 3 */
        cskx(:) = cos(kx(:))
        snkx(:) = sin(kx(:))
        csky(:) = cos(ky(:))
        snky(:) = sin(ky(:))
#if NDIMS == 3
        cskz(:) = cos(kz(:))
        snkz(:) = sin(kz(:))
#endif /* NDIMS == 3 */

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

              cs = cskx(i) * csky(j) - snkx(i) * snky(j)
              sn = cskx(i) * snky(j) + snkx(i) * csky(j)
#if NDIMS == 3
              tt = cs
              cs = tt * cskz(k) - sn * snkz(k)
              sn = tt * snkz(k) + sn * cskz(k)
#endif /* NDIMS == 3 */

              acc(1:NDIMS,i,j,k) = acc(1:NDIMS,i,j,k)                          &
                                 + (real(fcoefs(l,:)) * cs                     &
                                 + aimag(fcoefs(l,:)) * sn)
            end do
          end do
#if NDIMS == 3
        end do
#endif /* NDIMS == 3 */

      end if

    end do ! l = 1, nmodes

! apply the injection profiles
!
    if (.not. periodic(1)) then
      do i = 1, nn
        if (abs(x(i)) > xup) then
          acc(1:NDIMS,i,:,:) = 0.0d+00
        else if (abs(x(i)) >= xlo) then
          acc(1:NDIMS,i,:,:) = acc(1:NDIMS,i,:,:) &
                             * cos(xfc * (abs(x(i)) - xlo))**2
        end if
      end do
    end if
    if (.not. periodic(2)) then
      do j = 1, nn
        if (abs(y(j)) > yup) then
          acc(1:NDIMS,:,j,:) = 0.0d+00
        else if (abs(y(j)) >= ylo) then
          acc(1:NDIMS,:,j,:) = acc(1:NDIMS,:,j,:) &
                             * cos(yfc * (abs(y(j)) - ylo))**2
        end if
      end do
    end if
#if NDIMS == 3
    if (.not. periodic(3)) then
      do k = 1, nn
        if (abs(z(k)) > zup) then
          acc(1:NDIMS,:,:,k) = 0.0d+00
        else if (abs(z(k)) >= zlo) then
          acc(1:NDIMS,:,:,k) = acc(1:NDIMS,:,:,k) &
                             * cos(zfc * (abs(z(k)) - zlo))**2
        end if
      end do
    end if
#endif /* NDIMS == 3 */

! multiply the acceleration by density and time step
!
    do n = 1, NDIMS
      acc(n,:,:,:) = pdata%u(idn,:,:,:) * acc(n,:,:,:) * dt
    end do

! calculate the kinetic energy before the injection
!
    den(:,:,:) = sum(pdata%u(imx:imz,:,:,:)**2, 1)

! add the momentum increment
!
    pdata%u(imx,:,:,:) = pdata%u(imx,:,:,:) + acc(1,:,:,:)
    pdata%u(imy,:,:,:) = pdata%u(imy,:,:,:) + acc(2,:,:,:)
#if NDIMS == 3
    pdata%u(imz,:,:,:) = pdata%u(imz,:,:,:) + acc(3,:,:,:)
#endif /* NDIMS == 3 */

! calculate the injected energy
!
    den(:,:,:) = 5.0d-01 * (sum(pdata%u(imx:imz,:,:,:)**2, 1)                &
               - den(:,:,:)) / pdata%u(idn,:,:,:)

! calculate the injected energy
!
#if NDIMS == 3
    einj = einj + sum(den(nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */
    einj = einj + sum(den(nb:ne,nb:ne, 1   )) * dvol
#endif /* NDIMS == 3 */

! add the energy increment
!
    if (ien > 0) then
      pdata%u(ien,:,:,:) = pdata%u(ien,:,:,:) + den(:,:,:)
    end if

    work_in_use(nt) = .false.

    100 continue

!-------------------------------------------------------------------------------
!
  end subroutine inject_fmodes_block
!
!===============================================================================
!
! subroutine GET_VCOEFS:
! ---------------------
!
!   Subroutine determines velocity Fourier coefficients.
!
!
!===============================================================================
!
  subroutine get_vcoefs()

    use blocks  , only : block_data, list_data
#ifdef MPI
    use mpitools, only : reduce_sum
#endif /* MPI */

    implicit none

    type(block_data), pointer :: pdata

!-------------------------------------------------------------------------------
!
! reset vcoefs
!
    vcoefs(:,:) = cmplx(0.0d+00, 0.0d+00, kind=8)

! assign pdata with the first block on the data block list
!
    pdata => list_data

! iterate over all data blocks
!
    do while (associated(pdata))

! get contribution of velocity coefficients from the current block
!
      call get_vcoefs_block(pdata)

! assign pdata to the next block
!
      pdata => pdata%next

    end do ! over data blocks

#ifdef MPI
! reduce velocity coefficients over all processes
!
    call reduce_sum(vcoefs)
#endif /* MPI */

!-------------------------------------------------------------------------------
!
  end subroutine get_vcoefs
!
!===============================================================================
!
! subroutine GET_VCOEFS_BLOCK:
! ---------------------------
!
!   Subroutine determines Fourier coefficients of velocity.
!
!   Arguments:
!
!     pdata - a pointer to the data block;
!
!===============================================================================
!
  subroutine get_vcoefs_block(pdata)

    use blocks     , only : block_data
    use constants  , only : pi2
    use coordinates, only : nm => bcells, nb, ne
    use coordinates, only : ax, ay, advol
#if NDIMS == 3
    use coordinates, only : az
#endif /* NDIMS == 3 */
    use equations  , only : ivx, ivy
#if NDIMS == 3
    use equations  , only : ivz
#endif /* NDIMS == 3 */

    implicit none

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

    integer         :: i, j, k, l
    real(kind=8)    :: cs, sn, dvol
#if NDIMS == 3
    real(kind=8)    :: tt
#endif /* NDIMS == 3 */
    complex(kind=8) :: cf

    real(kind=8), dimension(nm):: x, y, z
    real(kind=8), dimension(nm):: kx, ky, snkx, snky, cskx, csky
#if NDIMS == 3
    real(kind=8), dimension(nm):: kz, snkz, cskz
#endif /* NDIMS == 3 */

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

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

! iterate over driving modes and compute the velocity Fourier coefficients
!
    do l = 1, nmodes

      kx(:) = kmodes(l,1) * x(:)
      ky(:) = kmodes(l,2) * y(:)
#if NDIMS == 3
      kz(:) = kmodes(l,3) * z(:)
#endif /* NDIMS == 3 */
      cskx(:) = cos(kx(:))
      snkx(:) = sin(kx(:))
      csky(:) = cos(ky(:))
      snky(:) = sin(ky(:))
#if NDIMS == 3
      cskz(:) = cos(kz(:))
      snkz(:) = sin(kz(:))
#endif /* NDIMS == 3 */

#if NDIMS == 3
      do k = nb, ne
#endif /* NDIMS == 3 */
        do j = nb, ne
          do i = nb, ne

            cs = cskx(i) * csky(j) - snkx(i) * snky(j)
            sn = cskx(i) * snky(j) + snkx(i) * csky(j)
#if NDIMS == 3
            tt = cs
            cs = tt * cskz(k) - sn * snkz(k)
            sn = tt * snkz(k) + sn * cskz(k)
#endif /* NDIMS == 3 */

            cf = cmplx(cs, sn, kind=8) * dvol

            vcoefs(l,1) = vcoefs(l,1) + pdata%q(ivx,i,j,k) * cf
            vcoefs(l,2) = vcoefs(l,2) + pdata%q(ivy,i,j,k) * cf
#if NDIMS == 3
            vcoefs(l,3) = vcoefs(l,3) + pdata%q(ivz,i,j,k) * cf
#endif /* NDIMS == 3 */

          end do
        end do
#if NDIMS == 3
      end do
#endif /* NDIMS == 3 */

    end do ! l = 1, nmodes

!-------------------------------------------------------------------------------
!
  end subroutine get_vcoefs_block

!===============================================================================
!
end module forcing