amun-code/sources/forcing.F90
Grzegorz Kowal 586df0bfa6 FORCING: Fix generation of vectors e1 and e2 in Alvelius method.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
2020-04-28 14:25:15 -03:00

1899 lines
50 KiB
Fortran
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

!!******************************************************************************
!!
!! 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-2020 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
#ifdef PROFILE
! import external subroutines
!
use timers, only : set_timer, start_timer, stop_timer
#endif /* PROFILE */
! module variables are not implicit by default
!
implicit none
#ifdef PROFILE
! timer indices
!
integer, save :: ifi, ifu
#endif /* PROFILE */
! interfaces for procedure pointers
!
abstract interface
subroutine update_forcing_iface(t, dt)
real(kind=8), intent(in) :: t, 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_powerlaw = 0
integer, parameter :: profile_normal = 1
integer, parameter :: profile_parabolic = 2
! driving force method and spectral profile
!
integer, save :: injection_method = injection_none
integer, save :: injection_profile = profile_powerlaw
! 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
! 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 = 1.0d+00
real(kind=8), save :: ascale = 1.0d+00
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 = 1.0d+00
real(kind=8), save :: pidx = - 5.0d+00 / 3.0d+00
real(kind=8), save :: fmin = 1.0d-08
real(kind=8), save :: sole = 1.0d+00
real(kind=8), save :: norm = 3.0d+00 / sqrt(2.0d+00)
! 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
! solenoidal projection tensor
!
real(kind=8), dimension(:,:,:), allocatable :: kprjct
! array for driving mode coefficients
!
complex(kind=8), dimension(:,:), allocatable :: 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
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
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)
! import external procedures and variables
!
use iso_fortran_env, only : error_unit
use parameters , only : get_parameter
use random , only : randsym, randnorz
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
logical, intent(in) :: verbose
integer, intent(out) :: status
! local variables
!
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"
logical :: test
integer :: i, j, k = 0, l, k2
real(kind=8) :: kl2, ku2, kv2, kv
real(kind=8) :: fa, fi, uu
! local vectors
!
real(kind=8), dimension(NDIMS) :: u, v
! allocatable arrays
!
integer, dimension(:), allocatable :: kcount
! local parameters
!
character(len=*), parameter :: loc = 'FORCING::initialize_forcing()'
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! set timer descriptions
!
call set_timer('forcing:: initialization', ifi)
call set_timer('forcing:: update' , ifu)
! start accounting time for initialization/finalization
!
call start_timer(ifi)
#endif /* PROFILE */
! reset the status flag
!
status = 0
! obtain the chosen injection method
!
call get_parameter("injection_method", injection)
! 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("OrnsteinUhlenbeck", "ornsteinuhlenbeck", "oh")
forcing_enabled = .true.
forcing_fourier = .true.
injection_method = injection_oh
update_forcing => update_forcing_oh
! get the solenoidal factor
!
call get_parameter("solenoidal_factor", sole)
if (sole < 0.0d+00 .or. sole > 1.0d+00) then
if (verbose) then
write(*,*)
write(*,wfmt) "'solenoidal_factor' must be between 0.0 and 1.0!"
write(*,wfmt) "Resetting it to the default value: 1.0!"
end if
sole = 1.0d+00
end if
case("Alvelius", "alvelius", "al")
forcing_enabled = .true.
forcing_fourier = .true.
injection_method = injection_alvelius
update_forcing => update_forcing_alvelius
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
! 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
! get the autocorrelation/turnover time
!
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
! get the characteristic velocity scale
!
call get_parameter("driving_velocity", vscale)
if (vscale <= 0.0d+00) then
if (verbose) then
write(*,*)
write(*,wfmt) "'driving_velocity' must be larger than zero!"
write(*,wfmt) "Resetting it to the default value: 1.0!"
end if
vscale = 1.0d+00
end if
call get_parameter("injection_scale" , kf)
call get_parameter("injection_width" , kc)
call get_parameter("lower_scale" , kl)
call get_parameter("upper_scale" , ku)
call get_parameter("amplitude_threshold", fmin)
if (kf < 1.0d+00) then
if (verbose) then
write(*,*)
write(*,wfmt) "'injection_scale' 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) "'injection_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) "'injection_width' must be smaller than " // &
"'injection_scale'!"
write(*,wfmt) "Resetting it to half of 'injection_scale'!"
end if
kc = 5.0d-01 * kf
end if
if (kl >= kf) then
if (verbose) then
write(*,*)
write(*,wfmt) "'lower_scale' must be smaller than 'injection_scale'!"
write(*,wfmt) "Setting it to 'injection_scale' - 'injection_width'!"
end if
kl = kf - kc
end if
if (ku <= kf) then
if (verbose) then
write(*,*)
write(*,wfmt) "'upper_scale' must be larger than 'injection_scale'!"
write(*,wfmt) "Setting it to 'injection_scale' + 'injection_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(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 */
! determine if the profile is set by energy or amplitude
!
call get_parameter("energy_profile", 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("injection_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("power_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
do l = 1, nmodes
k2 = sum(kmodes(l,:)**2)
fmodes(l) = fmodes(l) / sqrt(fi)
end do
! --- 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
do l = 1, nmodes
k2 = sum(kmodes(l,:)**2)
fmodes(l) = fmodes(l) / sqrt(fi)
end do
! --- 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
do l = 1, nmodes
k2 = sum(kmodes(l,:)**2)
fmodes(l) = fmodes(l) / sqrt(fi)
end do
! --- 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
! === OrnsteinUhlenbeck driving ===
!
if (injection_method == injection_oh) then
! get the characteristic driving length, velocity, time, and acceleration scales
!
lscale = 1.0d+00 / kf
ascale = vscale / tscale * sqrt(2.0d+00 / tscale)
fpower = vscale**2 / tscale
! normalize the driving amplitude profile to the driving power
!
fmodes(:) = ascale * fmodes(:)
! allocate arrays for projection tensors
!
allocate(kprjct(nmodes,NDIMS,NDIMS), stat = status)
if (status == 0) then
! generate the driving mode coefficients with random phases and directions
!
do l = 1, nmodes
! prepare the projection tensors
!
uu = sum(kmodes(l,:)**2)
fa = (1.0d+00 - 2.0d+00 * sole) / uu
do j = 1, NDIMS
kprjct(l,j,j) = sole
do i = 1, NDIMS
kprjct(l,i,j) = kprjct(l,i,j) &
+ fa * kmodes(l,i) * kmodes(l,j)
end do
end do
! generate a random vector, project it and normalize
!
test = .true.
do while (test)
do i = 1, NDIMS
v(i) = randsym()
end do
do i = 1, NDIMS
u(i) = dot_product(v(:), kprjct(l,:,i))
end do
uu = dot_product(u(:), u(:))
if (uu > 0.0d+00) then
u(:) = u(:) / sqrt(uu)
test = .false.
end if
end do ! while
! calculate the initial driving mode coefficient
!
fcoefs(l,:) = fmodes(l) * u(:) * randnorz()
end do ! l = 1, nmodes
! calculate the rms of acceleration
!
arms = sqrt(abs(sum(fcoefs(:,:) * conjg(fcoefs(:,:)))) / NDIMS)
end if ! allocation of kprjct
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(5.0d-01 * lscale / sqrt(fpower))
vscale = lscale / tscale
ascale = vscale / tscale
! normalize the driving amplitude profile to the driving power
!
fmodes(:) = ascale * fmodes(:)
! allocate perpendicular vectors
!
allocate(e1vecs(nmodes,NDIMS), e2vecs(nmodes,NDIMS), &
vcoefs(nmodes,NDIMS), stat = status)
if (status == 0) then
! 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
write(error_unit,"('[', a, ']: ', a)") &
trim(loc), "v(:) is ZERO!"
end if
end do
fcoefs(:,:) = cmplx(0.0d+00, 0.0d+00)
! calculate the rms of acceleration
!
arms = sqrt(abs(sum(fcoefs(:,:) * conjg(fcoefs(:,:)))) / NDIMS)
end if ! allocation of e1 and e2 vectors
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
#ifdef PROFILE
! stop accounting time
!
call stop_timer(ifi)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
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)
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
integer, intent(out) :: status
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for initialization/finalization
!
call start_timer(ifi)
#endif /* PROFILE */
if (allocated(kmodes)) deallocate(kmodes)
if (allocated(fmodes)) deallocate(fmodes)
if (allocated(fcoefs)) deallocate(fcoefs)
if (allocated(kprjct)) deallocate(kprjct)
if (allocated(e1vecs)) deallocate(e1vecs)
if (allocated(e2vecs)) deallocate(e2vecs)
if (allocated(vcoefs)) deallocate(vcoefs)
#ifdef PROFILE
! stop accounting time
!
call stop_timer(ifi)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine finalize_forcing
!
!===============================================================================
!
! subroutine PRINT_FORCING:
! ------------------------
!
! Subroutine prints module parameters.
!
! Arguments:
!
! verbose - a logical flag turning the information printing;
!
!===============================================================================
!
subroutine print_forcing(verbose)
! import external procedures and variables
!
use helpers , only : print_section, print_parameter
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
logical, intent(in) :: verbose
! local variables
!
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 = "OrnsteinUhlenbeck"
case(injection_alvelius)
injection = "Alvelius"
case default
injection = "none"
end select
call print_section(verbose, "Forcing")
call print_parameter(verbose, "injection 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
call print_parameter(verbose, "driving power" , fpower)
call print_parameter(verbose, "driving length scale" , lscale)
call print_parameter(verbose, "driving time scale" , tscale)
call print_parameter(verbose, "driving velocity scale", vscale)
call print_parameter(verbose, "driving acceleration scale", ascale)
if (injection_method == injection_oh) then
call print_parameter(verbose, "solenoidal factor" , sole)
end if
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, "injection 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, "injection scale" , kf)
call print_parameter(verbose, "lower scale limit", kl)
call print_parameter(verbose, "upper scale 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, "gaussian profile width", kc)
end if
call print_parameter(verbose, "amplitude threshold", fmin)
call print_parameter(verbose, "number of modes" , nmodes)
end if
end if
!-------------------------------------------------------------------------------
!
end subroutine print_forcing
!
!===============================================================================
!
! subroutine UPDATE_FORCING_EDDIES:
! --------------------------------
!
! Subroutine adds the energy injection terms.
!
! Arguments:
!
! t, dt - time and its increment;
!
!===============================================================================
!
subroutine update_forcing_eddies(t, dt)
! import external procedures and variables
!
use coordinates, only : xmin, ymin, zmin, xlen, ylen, zlen
use random , only : randuni, randsym
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), intent(in) :: t, dt
! local variables
!
integer :: ni, n
real(kind=8) :: tmp
real(kind=8), dimension(3) :: xp, ap
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for forcing term update
!
call start_timer(ifu)
#endif /* PROFILE */
! proceed only if forcing is enabled
!
if (forcing_enabled) then
! 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 == 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 if ! forcing enabled
#ifdef PROFILE
! stop accounting time
!
call stop_timer(ifu)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine update_forcing_eddies
!
!===============================================================================
!
! subroutine UPDATE_FORCING_OH:
! ----------------------------
!
! Subroutine adds the energy injection terms for OrnsteinUhlenbeck driving.
!
! Arguments:
!
! t, dt - time and its increment;
!
!===============================================================================
!
subroutine update_forcing_oh(t, dt)
! import external procedures and variables
!
use random, only : randsym, randnorz
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), intent(in) :: t, dt
! local variables
!
logical :: test
integer :: i, j, k, l
real(kind=8) :: acoeff, dcoeff
real(kind=8) :: dinj, uu
! local vectors
!
complex(kind=8), dimension(NDIMS) :: finj
real(kind=8) , dimension(NDIMS) :: u, v
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for forcing term update
!
call start_timer(ifu)
#endif /* PROFILE */
! proceed only if forcing is enabled
!
if (forcing_enabled) then
! 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 a random vector, project it and normalize
!
test = .true.
do while (test)
do i = 1, NDIMS
v(i) = randsym()
end do
do i = 1, NDIMS
u(i) = dot_product(v(:), kprjct(l,:,i))
end do
uu = dot_product(u(:), u(:))
if (uu > 0.0d+00) then
u(:) = u(:) / sqrt(uu)
test = .false.
end if
end do ! while
! generate gaussian random vector along the driving modes
!
finj(:) = fmodes(l) * u(:) * randnorz()
! advance the driving force modes
!
fcoefs(l,:) = acoeff * fcoefs(l,:) + dcoeff * finj(:)
end do
! calculate the rms of acceleration
!
arms = sqrt(abs(sum(fcoefs(:,:) * conjg(fcoefs(:,:)))) / NDIMS)
! 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 if ! forcing enabled
#ifdef PROFILE
! stop accounting time
!
call stop_timer(ifu)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine update_forcing_oh
!
!===============================================================================
!
! subroutine UPDATE_FORCING_ALVELIUS:
! ----------------------------------
!
! Subroutine adds the energy injection terms for Alvelius driving.
!
! Arguments:
!
! t, dt - time and its increment;
!
!===============================================================================
!
subroutine update_forcing_alvelius(t, dt)
! import external procedures and variables
!
use constants, only : pi2
use random , only : randuni
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), intent(in) :: t, dt
! local variables
!
logical :: test
integer :: i, j, k, l
real(kind=8) :: th1, th2, phi, psi, ga, gb, dinj, sqdt
complex(kind=8) :: aran, bran, xi1, xi2
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for forcing term update
!
call start_timer(ifu)
#endif /* PROFILE */
! proceed only if forcing is enabled
!
if (forcing_enabled) then
! determine velocify coefficients in Fourier space
!
call get_vcoefs()
! square root of dt
!
sqdt = sqrt(dt)
! iterate over driving modes
!
do l = 1, nmodes
! get xi1 and xi2
!
xi1 = dot_product(vcoefs(l,:), e1vecs(l,:))
xi2 = dot_product(vcoefs(l,:), e2vecs(l,:))
! get random angles theta1, theta2, and phi
!
phi = pi2 * randuni()
ga = sin(phi)
gb = cos(phi)
psi = pi2 * randuni()
th1 = - ga * dimag(xi1) + gb * (sin(psi) * dreal(xi2) &
- cos(psi) * dimag(xi2))
if (th1 /= 0.0d+00) then
th1 = atan((ga * dreal(xi1) + gb * (sin(psi) * dimag(xi2) &
+ cos(psi) * dreal(xi2))) / th1)
end if
th2 = psi + th1
! calculate amplitude
!
aran = cmplx(cos(th1), sin(th1)) * ga
bran = cmplx(cos(th2), sin(th2)) * gb
! advance the driving force modes
!
fcoefs(l,:) = fmodes(l) * (aran * e1vecs(l,:) &
+ bran * e2vecs(l,:)) / sqdt
end do
! calculate the rms of acceleration
!
arms = sqrt(abs(sum(fcoefs(:,:) * conjg(fcoefs(:,:)))) / NDIMS)
! 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 if ! forcing enabled
#ifdef PROFILE
! stop accounting time
!
call stop_timer(ifu)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
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)
! include external variables
!
use blocks, only : block_data, list_data
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), dimension(3), intent(in) :: xp, ap
! local pointers
!
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)
! include external variables
!
use blocks , only : block_data
use coordinates, only : nm => bcells
use coordinates, only : ax, ay, az
use coordinates, only : xlen, ylen, zlen
use coordinates, only : periodic
use equations , only : nv
use equations , only : idn, imx, imy, imz, ien
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
type(block_data), pointer , intent(inout) :: pdata
real(kind=8), dimension(3), intent(in) :: xp, ap
! local variables
!
integer :: i, j, k = 1
real(kind=8) :: x2, y2, z2, r2
real(kind=8) :: fx, fy, fz, fp, e1, e2
! local arrays
!
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 */
! 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)
! include external variables
!
use blocks, only : block_data, list_data
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), intent(in) :: dt
! local pointers
!
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)
! include external variables
!
use blocks , only : block_data
use constants , only : pi2
use coordinates, only : nm => bcells, nb, ne
use coordinates, only : ax, ay, az, advol
use equations , only : idn, imx, imy, imz, ien, ivx, ivy, ivz
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
type(block_data), pointer, intent(inout) :: pdata
real(kind=8) , intent(in) :: dt
! local variables
!
integer :: i, j, k = 1, l, n
real(kind=8) :: cs, sn, tt
real(kind=8) :: dvol
! local arrays
!
real(kind=8), dimension(nm):: x, y, z, kx, ky, kz
real(kind=8), dimension(nm):: snkx, snky, snkz
real(kind=8), dimension(nm):: cskx, csky, cskz
#if NDIMS == 3
real(kind=8), dimension(3,nm,nm,nm) :: acc
real(kind=8), dimension( nm,nm,nm) :: den
#else /* NDIMS == 3 */
real(kind=8), dimension(2,nm,nm, 1) :: acc
real(kind=8), dimension( nm,nm, 1) :: den
#endif /* NDIMS == 3 */
!
!-------------------------------------------------------------------------------
!
! 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)
! reset the acceleration
!
acc(:,:,:,:) = 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, nm
#endif /* NDIMS == 3 */
do j = 1, nm
do i = 1, nm
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(:,i,j,k) = acc(:,i,j,k) + (dreal(fcoefs(l,:)) * cs &
+ dimag(fcoefs(l,:)) * sn)
end do
end do
#if NDIMS == 3
end do
#endif /* NDIMS == 3 */
end if
end do ! l = 1, nmodes
! 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
!-------------------------------------------------------------------------------
!
end subroutine inject_fmodes_block
!
!===============================================================================
!
! subroutine GET_VCOEFS:
! ---------------------
!
! Subroutine determines velocity Fourier coefficients.
!
!
!===============================================================================
!
subroutine get_vcoefs()
! include external variables
!
use blocks , only : block_data, list_data
#ifdef MPI
use mpitools, only : reduce_sum_complex_array
#endif /* MPI */
! local variables are not implicit by default
!
implicit none
! local variables
!
integer :: status
! local pointers
!
type(block_data), pointer :: pdata
!
!-------------------------------------------------------------------------------
!
! reset vcoefs
!
vcoefs(:,:) = cmplx(0.0d+00, 0.0d+00)
! 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_complex_array(size(vcoefs), vcoefs, status)
#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)
! include external variables
!
use blocks , only : block_data
use constants , only : pi2
use coordinates, only : nm => bcells, nb, ne
use coordinates, only : ax, ay, az, advol
use equations , only : ivx, ivy, ivz
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
type(block_data), pointer, intent(inout) :: pdata
! local variables
!
integer :: i, j, k = 1, l
real(kind=8) :: cs, sn, tt, dvol
complex(kind=8) :: cf
! local arrays
!
real(kind=8), dimension(nm):: x, y, z, kx, ky, kz
real(kind=8), dimension(nm):: snkx, snky, snkz
real(kind=8), dimension(nm):: cskx, csky, cskz
!
!-------------------------------------------------------------------------------
!
! 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) * 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