m_helper.fpp.f90

Go to the documentation of this file.

# 1 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
# 1 "/home/runner/work/MFC/MFC/src/common/include/case.fpp" 1
! This file exists so that Fypp can be run without generating case.fpp files for
! each target. This is useful when generating documentation, for example. This
! should also let MFC be built with CMake directly, without invoking mfc.sh.
 
! For pre-process.
# 9 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
 
! For moving immersed boundaries in simulation
# 14 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
# 2 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp" 2
# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
 
# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 76 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 151 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 192 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 206 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 231 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 242 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 244 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 284 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 294 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 304 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 313 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 340 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 347 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 353 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 359 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 365 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 371 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
 
# 377 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
! New line at end of file is required for FYPP
# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
 
# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
 
# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
 
# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
! New line at end of file is required for FYPP
# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
 
! GPU parallel region (scalar reductions, maxval/minval)
# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! GPU parallel loop over threads (most common GPU macro)
# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Required closing for GPU_PARALLEL_LOOP
# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Mark routine for device compilation
# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Declare device-resident data
# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Inner loop within a GPU parallel region
# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Scoped GPU data region
# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Host code with device pointers (for MPI with GPU buffers)
# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Allocate device memory (unscoped)
# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Free device memory
# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Atomic operation on device
# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! End atomic capture block
# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Copy data between host and device
# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Synchronization barrier
# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Import GPU library module (openacc or omp_lib)
# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Emit code only for AMD compiler
# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Emit code for non-Cray compilers
# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Emit code only for Cray compiler
# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
! Emit code for non-NVIDIA compilers
# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
 
# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
! New line at end of file is required for FYPP
# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
 
# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
! example see misc/nvidia_uvm/bind.sh. NVIDIA unified memory page placement hint
# 57 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Allocate and create GPU device memory
# 77 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Free GPU device memory and deallocate
# 85 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Cray-specific GPU pointer setup for vector fields
# 109 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Cray-specific GPU pointer setup for scalar fields
# 125 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
! Cray-specific GPU pointer setup for acoustic source spatials
# 150 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 156 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
 
# 163 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
! New line at end of file is required for FYPP
# 3 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp" 2
 
!>
!! @file
!! @brief Contains module m_helper
 
!> @brief Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions
module m_helper
 
    use m_derived_types
    use m_global_parameters
    use ieee_arithmetic  !< for checking nan
 
    implicit none
 
    private
    public :: s_comp_n_from_prim, s_comp_n_from_cons, s_initialize_bubbles_model, s_initialize_nonpoly, s_simpson, s_transcoeff, &
        & s_int_to_str, s_transform_vec, s_transform_triangle, s_transform_model, s_swap, f_cross, f_create_transform_matrix, &
        & f_create_bbox, s_print_2d_array, f_xor, f_logical_to_int, associated_legendre, real_ylm, double_factorial, factorial, &
        & f_cut_on, f_cut_off, s_downsample_data, s_upsample_data, s_cross_product
 
contains
 
    !> Computes the bubble number density n from the primitive variables
    subroutine s_comp_n_from_prim(vftmp, Rtmp, ntmp, weights)
 
 
# 28 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#if MFC_OpenACC
# 28 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
!$acc routine seq 
# 28 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#elif MFC_OpenMP
# 28 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 28 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 28 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
!$omp declare target device_type(any) 
# 28 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#endif
        real(wp), intent(in)                :: vftmp
        real(wp), dimension(nb), intent(in) :: rtmp
        real(wp), intent(out)               :: ntmp
        real(wp), dimension(nb), intent(in) :: weights
        real(wp)                            :: r3
 
        r3 = dot_product(weights, rtmp**3._wp)
        ntmp = (3._wp/(4._wp*pi))*vftmp/r3
 
    end subroutine s_comp_n_from_prim
 
    !> Compute the bubble number density from the conservative void fraction and weighted bubble radii.
    subroutine s_comp_n_from_cons(vftmp, nRtmp, ntmp, weights)
 
 
# 43 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#if MFC_OpenACC
# 43 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
!$acc routine seq 
# 43 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#elif MFC_OpenMP
# 43 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 43 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 43 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
!$omp declare target device_type(any) 
# 43 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#endif
        real(wp), intent(in)                :: vftmp
        real(wp), dimension(nb), intent(in) :: nrtmp
        real(wp), intent(out)               :: ntmp
        real(wp), dimension(nb), intent(in) :: weights
        real(wp)                            :: nr3
 
        nr3 = dot_product(weights, nrtmp**3._wp)
        ntmp = sqrt((4._wp*pi/3._wp)*nr3/vftmp)
 
    end subroutine s_comp_n_from_cons
 
    !> Print a 2D real array to standard output, optionally dividing each element by a given scalar.
    impure subroutine s_print_2d_array(A, div)
 
        real(wp), dimension(:,:), intent(in) :: a
        real(wp), optional, intent(in)       :: div
        integer                              :: i, j
        integer                              :: local_m, local_n
        real(wp)                             :: c
 
        local_m = size(a, 1)
        local_n = size(a, 2)
 
        if (present(div)) then
            c = div
        else
            c = 1._wp
        end if
 
        print *, local_m, local_n
 
        do i = 1, local_m
            do j = 1, local_n
                write (*, fmt="(F12.4)", advance="no") a(i, j)/c
            end do
            write (*, fmt="(A1)") " "
        end do
        write (*, fmt="(A1)") " "
 
    end subroutine s_print_2d_array
 
    !> Initialize bubble model arrays for Euler or Lagrangian bubbles with polytropic or non-polytropic gas.
    impure subroutine s_initialize_bubbles_model()
 
        ! Allocate memory
        if (bubbles_euler) then
#ifdef MFC_DEBUG
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
    block
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
        use iso_fortran_env, only: output_unit
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
        print *, 'm_helper.fpp:90: ', '@:ALLOCATE(weight(nb), R0(nb))'
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
        call flush (output_unit)
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
    end block
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#endif
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
    allocate (weight(nb), r0(nb))
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#if defined(MFC_OpenACC)
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
!$acc enter data create(weight, R0) 
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#elif defined(MFC_OpenMP)
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
!$omp target enter data map(always,alloc:weight, R0) 
# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#endif
            if (.not. polytropic) then
#ifdef MFC_DEBUG
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
    block
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
        use iso_fortran_env, only: output_unit
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
        print *, 'm_helper.fpp:92: ', '@:ALLOCATE(pb0(nb), Pe_T(nb), k_g(nb), k_v(nb), mass_g0(nb), mass_v0(nb))'
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
        call flush (output_unit)
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
    end block
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#endif
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
    allocate (pb0(nb), pe_t(nb), k_g(nb), k_v(nb), mass_g0(nb), mass_v0(nb))
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#if defined(MFC_OpenACC)
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
!$acc enter data create(pb0, Pe_T, k_g, k_v, mass_g0, mass_v0) 
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#elif defined(MFC_OpenMP)
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
!$omp target enter data map(always,alloc:pb0, Pe_T, k_g, k_v, mass_g0, mass_v0) 
# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#endif
#ifdef MFC_DEBUG
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
    block
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
        use iso_fortran_env, only: output_unit
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
        print *, 'm_helper.fpp:93: ', '@:ALLOCATE(Re_trans_T(nb), Re_trans_c(nb), Im_trans_T(nb), Im_trans_c(nb))'
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
        call flush (output_unit)
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
    end block
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#endif
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
    allocate (re_trans_t(nb), re_trans_c(nb), im_trans_t(nb), im_trans_c(nb))
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#if defined(MFC_OpenACC)
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
!$acc enter data create(Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c) 
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#elif defined(MFC_OpenMP)
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
!$omp target enter data map(always,alloc:Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c) 
# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#endif
            else if (qbmm) then
#ifdef MFC_DEBUG
# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
    block
# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
        use iso_fortran_env, only: output_unit
# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
        print *, 'm_helper.fpp:95: ', '@:ALLOCATE(pb0(nb))'
# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
        call flush (output_unit)
# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
    end block
# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#endif
# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
    allocate (pb0(nb))
# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#if defined(MFC_OpenACC)
# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
!$acc enter data create(pb0) 
# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#elif defined(MFC_OpenMP)
# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
!$omp target enter data map(always,alloc:pb0) 
# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#endif
            end if
 
            ! Compute quadrature weights and nodes for polydisperse simulations
            if (nb > 1) then
                call s_simpson(weight, r0)
            else if (nb == 1) then
                r0 = 1._wp
                weight = 1._wp
            else
                stop 'Invalid value of nb'
            end if
            r0 = r0*bub_pp%R0ref
        end if
 
        ! Initialize bubble variables
        call s_initialize_bubble_vars()
 
    end subroutine s_initialize_bubbles_model
 
    !> Set bubble physical parameters and nondimensional numbers from the input configuration.
    impure subroutine s_initialize_bubble_vars()
 
        r0ref = bub_pp%R0ref; p0ref = bub_pp%p0ref
        rho0ref = bub_pp%rho0ref
        ss = bub_pp%ss; pv = bub_pp%pv; vd = bub_pp%vd
        mu_l = bub_pp%mu_l; mu_v = bub_pp%mu_v; mu_g = bub_pp%mu_g
        gam_v = bub_pp%gam_v; gam_g = bub_pp%gam_g
        if (.not. polytropic) then
            if (bubbles_euler) then
                m_v = bub_pp%M_v; m_g = bub_pp%M_g
                k_v = bub_pp%k_v; k_g = bub_pp%k_g
            end if
            r_v = bub_pp%R_v; r_g = bub_pp%R_g
            tw = bub_pp%T0ref
        end if
        if (bubbles_lagrange) then
            cp_v = bub_pp%cp_v; cp_g = bub_pp%cp_g
            k_vl = bub_pp%k_v; k_gl = bub_pp%k_g
        end if
 
        ! Input quantities
        if (bubbles_euler .and. (.not. polytropic)) then
            if (thermal == 2) then
                gam_m = 1._wp
            else
                gam_m = gam_g
            end if
        end if
 
        ! Nondimensional numbers
        eu = p0ref
        ca = eu - pv
        if (.not. f_is_default(bub_pp%ss)) web = 1._wp/ss
        if (.not. f_is_default(bub_pp%mu_l)) re_inv = mu_l
        if (.not. polytropic) pe_c = 1._wp/vd
 
        if (bubbles_euler) then
            ! Initialize variables for non-polytropic (Preston) model
            if (.not. polytropic) then
                call s_initialize_nonpoly()
            end if
            ! Initialize pb based on surface tension for qbmm (polytropic)
            if (qbmm .and. polytropic) then
                pb0 = eu
                if (.not. f_is_default(web)) then
                    pb0 = pb0 + 2._wp/web/r0
                end if
            end if
        end if
 
    end subroutine s_initialize_bubble_vars
 
    !> Initializes non-polydisperse bubble modeling
    impure subroutine s_initialize_nonpoly()
 
        integer                 :: ir
        real(wp), dimension(nb) :: chi_vw0, cp_m0, k_m0, rho_m0, x_vw, omegan
        real(wp), parameter     :: k_poly = 1._wp  !< polytropic index used to compute isothermal natural frequency
        ! Chapman-Enskog transport coefficients for vapor-gas mixture, Ando JAS (2010)
 
        phi_vg = (1._wp + sqrt(mu_v/mu_g)*(m_g/m_v)**(0.25_wp))**2/(sqrt(8._wp)*sqrt(1._wp + m_v/m_g))
        phi_gv = (1._wp + sqrt(mu_g/mu_v)*(m_v/m_g)**(0.25_wp))**2/(sqrt(8._wp)*sqrt(1._wp + m_g/m_v))
 
        ! Initial internal bubble pressure (Euler number + Laplace pressure)
        pb0 = eu + 2._wp/web/r0
 
        ! Vapor mass fraction at bubble wall, Ando JAS (2010)
        chi_vw0 = 1._wp/(1._wp + r_v/r_g*(pb0/pv - 1._wp))
 
        ! Mixture specific heat from mass-weighted vapor/gas contributions
        cp_m0 = chi_vw0*r_v*gam_v/(gam_v - 1._wp) + (1._wp - chi_vw0)*r_g*gam_g/(gam_g - 1._wp)
 
        ! mole fraction of vapor (Eq. 2.23 in Ando 2010)
        x_vw = m_g*chi_vw0/(m_v + (m_g - m_v)*chi_vw0)
 
        ! thermal conductivity for gas/vapor mixture (Eq. 2.21 in Ando 2010)
        k_m0 = x_vw*k_v/(x_vw + (1._wp - x_vw)*phi_vg) + (1._wp - x_vw)*k_g/(x_vw*phi_gv + 1._wp - x_vw)
        k_g(:) = k_g(:)/k_m0(:)
        k_v(:) = k_v(:)/k_m0(:)
 
        ! mixture density (Eq. 2.20 in Ando 2010)
        rho_m0 = pv/(chi_vw0*r_v*tw)
 
        ! mass of gas/vapor
        mass_g0(:) = (4._wp*pi/3._wp)*(pb0(:) - pv)/(r_g*tw)*r0(:)**3
        mass_v0(:) = (4._wp*pi/3._wp)*pv/(r_v*tw)*r0(:)**3
 
        ! Peclet numbers
        pe_t(:) = rho_m0*cp_m0(:)/k_m0(:)
 
        ! Bubble natural frequency, Ando JAS (2010)
        omegan(:) = sqrt(3._wp*k_poly*ca + 2._wp*(3._wp*k_poly - 1._wp)/(web*r0))/r0/sqrt(rho0ref)
        do ir = 1, nb
            call s_transcoeff(omegan(ir)*r0(ir), pe_t(ir)*r0(ir), re_trans_t(ir), im_trans_t(ir))
            call s_transcoeff(omegan(ir)*r0(ir), pe_c*r0(ir), re_trans_c(ir), im_trans_c(ir))
        end do
        im_trans_t = 0._wp
 
    end subroutine s_initialize_nonpoly
 
    !> Computes the transfer coefficient for the non-polytropic bubble compression process
    elemental subroutine s_transcoeff(omega, peclet, Re_trans, Im_trans)
 
        real(wp), intent(in)  :: omega, peclet
        real(wp), intent(out) :: re_trans, im_trans
        complex(wp)           :: imag, trans, c1, c2, c3
 
        imag = (0._wp, 1._wp)
 
        c1 = imag*omega*peclet
        c2 = sqrt(c1)
        c3 = (exp(c2) - exp(-c2))/(exp(c2) + exp(-c2))  ! TANH(c2)
        trans = ((c2/c3 - 1._wp)**(-1) - 3._wp/c1)**(-1)  ! transfer function
 
        re_trans = trans
        im_trans = aimag(trans)
 
    end subroutine s_transcoeff
 
    !> Convert an integer to its trimmed string representation.
    elemental subroutine s_int_to_str(i, res)
 
        integer, intent(in)             :: i
        character(len=*), intent(inout) :: res
 
        write (res, '(I0)') i
        res = trim(res)
 
    end subroutine s_int_to_str
 
    !> Computes the Simpson weights for quadrature
    subroutine s_simpson(local_weight, local_R0)
 
        real(wp), dimension(:), intent(inout) :: local_weight
        real(wp), dimension(:), intent(inout) :: local_r0
        integer                               :: ir
        real(wp)                              :: r0mn, r0mx, dphi, tmp, sd
        real(wp), dimension(nb)               :: phi
 
        sd = poly_sigma
        r0mn = 0.8_wp*exp(-2.8_wp*sd)
        r0mx = 0.2_wp*exp(9.5_wp*sd) + 1._wp
 
        ! phi = ln( R0 ) & return R0
        do ir = 1, nb
            phi(ir) = log(r0mn) + (ir - 1._wp)*log(r0mx/r0mn)/(nb - 1._wp)
            local_r0(ir) = exp(phi(ir))
        end do
 
# 266 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
            dphi = phi(2) - phi(1)
# 268 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
        ! weights for quadrature using Simpson's rule
        do ir = 2, nb - 1
            ! Gaussian
            tmp = exp(-0.5_wp*(phi(ir)/sd)**2)/sqrt(2._wp*pi)/sd
            if (mod(ir, 2) == 0) then
                local_weight(ir) = tmp*4._wp*dphi/3._wp
            else
                local_weight(ir) = tmp*2._wp*dphi/3._wp
            end if
        end do
        tmp = exp(-0.5_wp*(phi(1)/sd)**2)/sqrt(2._wp*pi)/sd
        local_weight(1) = tmp*dphi/3._wp
        tmp = exp(-0.5_wp*(phi(nb)/sd)**2)/sqrt(2._wp*pi)/sd
        local_weight(nb) = tmp*dphi/3._wp
 
    end subroutine s_simpson
 
    !> Compute the cross product of two vectors.
    pure function f_cross(a, b) result(c)
 
 
# 289 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#if MFC_OpenACC
# 289 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
!$acc routine seq 
# 289 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#elif MFC_OpenMP
# 289 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 289 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 289 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
!$omp declare target device_type(any) 
# 289 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#endif
 
        real(wp), dimension(3), intent(in) :: a, b
        real(wp), dimension(3)             :: c
 
        c(1) = a(2)*b(3) - a(3)*b(2)
        c(2) = a(3)*b(1) - a(1)*b(3)
        c(3) = a(1)*b(2) - a(2)*b(1)
 
    end function f_cross
 
    !> @brief Computes the cross product c = a x b of two 3D vectors.
    subroutine s_cross_product(a, b, c)
 
 
# 303 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#if MFC_OpenACC
# 303 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
!$acc routine seq 
# 303 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#elif MFC_OpenMP
# 303 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 303 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
 
# 303 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
!$omp declare target device_type(any) 
# 303 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
#endif
        real(wp), intent(in)  :: a(3), b(3)
        real(wp), intent(out) :: c(3)
 
        c(1) = a(2)*b(3) - a(3)*b(2)
        c(2) = a(3)*b(1) - a(1)*b(3)
        c(3) = a(1)*b(2) - a(2)*b(1)
 
    end subroutine s_cross_product
 
    !> This procedure swaps two real numbers.
    !! @param lhs Left-hand side.
    !! @param rhs Right-hand side.
    elemental subroutine s_swap(lhs, rhs)
 
        real(wp), intent(inout) :: lhs, rhs
        real(wp)                :: ltemp
 
        ltemp = lhs
        lhs = rhs
        rhs = ltemp
 
    end subroutine s_swap
 
    !> Create a transformation matrix.
    function f_create_transform_matrix(param, center) result(out_matrix)
 
        type(ic_model_parameters), intent(in)          :: param
        real(wp), dimension(1:3), optional, intent(in) :: center
        real(wp), dimension(1:4,1:4)                   :: sc, rz, rx, ry, tr, t_back, t_to_origin, out_matrix
 
        sc = transpose(reshape([param%scale(1), 0._wp, 0._wp, 0._wp, 0._wp, param%scale(2), 0._wp, 0._wp, 0._wp, 0._wp, &
                       & param%scale(3), 0._wp, 0._wp, 0._wp, 0._wp, 1._wp], shape(sc)))
 
        rz = transpose(reshape([cos(param%rotate(3)), -sin(param%rotate(3)), 0._wp, 0._wp, sin(param%rotate(3)), &
                       & cos(param%rotate(3)), 0._wp, 0._wp, 0._wp, 0._wp, 1._wp, 0._wp, 0._wp, 0._wp, 0._wp, 1._wp], shape(rz)))
 
        rx = transpose(reshape([1._wp, 0._wp, 0._wp, 0._wp, 0._wp, cos(param%rotate(1)), -sin(param%rotate(1)), 0._wp, 0._wp, &
                       & sin(param%rotate(1)), cos(param%rotate(1)), 0._wp, 0._wp, 0._wp, 0._wp, 1._wp], shape(rx)))
 
        ry = transpose(reshape([cos(param%rotate(2)), 0._wp, sin(param%rotate(2)), 0._wp, 0._wp, 1._wp, 0._wp, 0._wp, &
                       & -sin(param%rotate(2)), 0._wp, cos(param%rotate(2)), 0._wp, 0._wp, 0._wp, 0._wp, 1._wp], shape(ry)))
 
        tr = transpose(reshape([1._wp, 0._wp, 0._wp, param%translate(1), 0._wp, 1._wp, 0._wp, param%translate(2), 0._wp, 0._wp, &
                       & 1._wp, param%translate(3), 0._wp, 0._wp, 0._wp, 1._wp], shape(tr)))
 
        if (present(center)) then
            ! Translation matrix to move center to the origin
            t_to_origin = transpose(reshape([1._wp, 0._wp, 0._wp, -center(1), 0._wp, 1._wp, 0._wp, -center(2), 0._wp, 0._wp, &
                                    & 1._wp, -center(3), 0._wp, 0._wp, 0._wp, 1._wp], shape(tr)))
 
            ! Translation matrix to move center back to original position
            t_back = transpose(reshape([1._wp, 0._wp, 0._wp, center(1), 0._wp, 1._wp, 0._wp, center(2), 0._wp, 0._wp, 1._wp, &
                               & center(3), 0._wp, 0._wp, 0._wp, 1._wp], shape(tr)))
 
            out_matrix = matmul(tr, matmul(t_back, matmul(ry, matmul(rx, matmul(rz, matmul(sc, t_to_origin))))))
        else
            out_matrix = matmul(ry, matmul(rx, rz))
        end if
 
    end function f_create_transform_matrix
 
    !> Transform a vector by a matrix.
    subroutine s_transform_vec(vec, matrix)
 
        real(wp), dimension(1:3), intent(inout)  :: vec
        real(wp), dimension(1:4,1:4), intent(in) :: matrix
        real(wp), dimension(1:4)                 :: tmp
 
        tmp = matmul(matrix, [vec(1), vec(2), vec(3), 1._wp])
        vec = tmp(1:3)
 
    end subroutine s_transform_vec
 
    !> Transform a triangle by a matrix, one vertex at a time.
    subroutine s_transform_triangle(triangle, matrix, matrix_n)
 
        type(t_triangle), intent(inout)          :: triangle
        real(wp), dimension(1:4,1:4), intent(in) :: matrix, matrix_n
        integer                                  :: i
 
        do i = 1, 3
            call s_transform_vec(triangle%v(i,:), matrix)
        end do
 
        call s_transform_vec(triangle%n(1:3), matrix_n)
 
    end subroutine s_transform_triangle
 
    !> Transform a model by a matrix, one triangle at a time.
    subroutine s_transform_model(model, matrix, matrix_n)
 
        type(t_model), intent(inout)             :: model
        real(wp), dimension(1:4,1:4), intent(in) :: matrix, matrix_n
        integer                                  :: i
 
        do i = 1, size(model%trs)
            call s_transform_triangle(model%trs(i), matrix, matrix_n)
        end do
 
    end subroutine s_transform_model
 
    !> Create a bounding box for a model.
    function f_create_bbox(model) result(bbox)
 
        type(t_model), intent(in) :: model
        type(t_bbox)              :: bbox
        integer                   :: i, j
 
        if (size(model%trs) == 0) then
            bbox%min = 0._wp
            bbox%max = 0._wp
            return
        end if
 
        bbox%min = model%trs(1)%v(1,:)
        bbox%max = model%trs(1)%v(1,:)
 
        do i = 1, size(model%trs)
            do j = 1, 3
                bbox%min = min(bbox%min, model%trs(i)%v(j,:))
                bbox%max = max(bbox%max, model%trs(i)%v(j,:))
            end do
        end do
 
    end function f_create_bbox
 
    !> Perform XOR on lhs and rhs.
    elemental function f_xor(lhs, rhs) result(res)
 
        logical, intent(in) :: lhs, rhs
        logical             :: res
 
        res = (lhs .and. .not. rhs) .or. (.not. lhs .and. rhs)
 
    end function f_xor
 
    !> Convert a logical to 1 or 0.
    elemental function f_logical_to_int(predicate) result(int)
 
        logical, intent(in) :: predicate
        integer             :: int
 
        if (predicate) then
            int = 1
        else
            int = 0
        end if
 
    end function f_logical_to_int
 
    !> Real spherical harmonic Y_lm(theta, phi). theta = polar angle from +z (acos(z/r)), phi = atan2(y,x). Uses associated Legendre
    !! P_l^|m|(cos theta). Standard normalisation.
    function real_ylm(theta, phi, l, m) result(Y)
 
        integer, intent(in)  :: l, m
        real(wp), intent(in) :: theta, phi
        real(wp)             :: y, x, prefac
        integer              :: m_abs
 
        m_abs = abs(m)
        if (m_abs > l) then
            y = 0._wp
            return
        end if
        x = cos(theta)
        prefac = sqrt((2*l + 1)*real(factorial(l - m_abs), wp)/real(factorial(l + m_abs), wp)/(4._wp*pi))
        if (m == 0) then
            y = prefac*associated_legendre(x, l, 0)
        else if (m > 0) then
            y = prefac*sqrt(2._wp)*associated_legendre(x, l, m_abs)*cos(m*phi)
        else
            y = prefac*sqrt(2._wp)*associated_legendre(x, l, m_abs)*sin(m_abs*phi)
        end if
 
    end function real_ylm
 
    !> Associated Legendre polynomial P_l^m(x) (Ferrers function, Condon-Shortley phase). Valid for integer l >= 0, 0 <= m <= l, and
    !! x in [-1,1]. Returns 0 for |m| > l or l < 0. Formulas: DLMF 14.10.3 (recurrence in degree), Wikipedia "Associated Legendre
    !! polynomials" (P_l^l and P_l^{l-1} identities). Recurrence: (l-m)P_l^m = (2l-1)x P_{l-1}^m - (l+m-1)P_{l-2}^m.
    !! @param x argument (typically cos(theta)), should be in [-1,1]
    !! @param l degree (>= 0)
    !! @param m_order order (0 <= m_order <= l)
    recursive function associated_legendre(x, l, m_order) result(result_P)
 
        integer, intent(in)  :: l, m_order
        real(wp), intent(in) :: x
        real(wp)             :: result_p
        real(wp)             :: one_minus_x2
 
        ! Out-of-domain: P_l^m = 0 for |m| > l or l < 0 (standard convention)
 
        if (l < 0 .or. m_order < 0 .or. m_order > l) then
            result_p = 0._wp
            return
        end if
 
        if (m_order <= 0 .and. l <= 0) then
            result_p = 1._wp
        else if (l == 1 .and. m_order <= 0) then
            result_p = x
        else if (l == 1 .and. m_order == 1) then
            one_minus_x2 = max(0._wp, 1._wp - x**2)
            result_p = -sqrt(one_minus_x2)
        else if (m_order == l) then
            ! P_l^l(x) = (-1)^l (2l-1)!! (1-x^2)^(l/2). Use real exponent for odd l
            one_minus_x2 = max(0._wp, 1._wp - x**2)
            result_p = (-1)**l*real(double_factorial(2*l - 1), wp)*one_minus_x2**(0.5_wp*real(l, wp))
        else if (m_order == l - 1) then
            result_p = x*(2*l - 1)*associated_legendre(x, l - 1, l - 1)
        else
            result_p = ((2*l - 1)*x*associated_legendre(x, l - 1, m_order) - (l + m_order - 1)*associated_legendre(x, l - 2, &
                        & m_order))/(l - m_order)
        end if
 
    end function associated_legendre
 
    !> Calculate the double factorial of an integer
    elemental function double_factorial(n_in) result(R_result)
 
        integer, intent(in)      :: n_in
        integer, parameter       :: int64_kind = selected_int_kind(18)  !< 18 bytes for 64-bit integer
        integer(kind=int64_kind) :: r_result
        integer                  :: i
 
        r_result = product((/(i, i=n_in, 1, -2)/))
 
    end function double_factorial
 
    !> Calculate the factorial of an integer
    elemental function factorial(n_in) result(R_result)
 
        integer, intent(in)      :: n_in
        integer, parameter       :: int64_kind = selected_int_kind(18)  !< 18 bytes for 64-bit integer
        integer(kind=int64_kind) :: r_result
        integer                  :: i
 
        r_result = product((/(i, i=n_in, 1, -1)/))
 
    end function factorial
 
    !> Calculate a smooth cut-on function that is zero for x values smaller than zero and goes to one, for generating smooth initial
    !! conditions
    function f_cut_on(x, eps) result(fx)
 
        real(wp), intent(in) :: x, eps
        real(wp)             :: fx
 
        fx = 1 - f_gx(x/eps)/(f_gx(x/eps) + f_gx(1 - x/eps))
 
    end function f_cut_on
 
    !> Calculate a smooth cut-off function that is one for x values smaller than zero and goes to zero, for generating smooth
    !! initial conditions
    function f_cut_off(x, eps) result(fx)
 
        real(wp), intent(in) :: x, eps
        real(wp)             :: fx
 
        fx = f_gx(x/eps)/(f_gx(x/eps) + f_gx(1 - x/eps))
 
    end function f_cut_off
 
    !> Helper function for f_cut_on and f_cut_off
    function f_gx(x) result(gx)
 
        real(wp), intent(in) :: x
        real(wp)             :: gx
 
        if (x > 0) then
            gx = exp(-1._wp/x)
        else
            gx = 0._wp
        end if
 
    end function f_gx
 
    !> Downsample conservative variable fields by a factor of 3 in each direction using volume averaging.
    subroutine s_downsample_data(q_cons_vf, q_cons_temp, m_ds, n_ds, p_ds, m_glb_ds, n_glb_ds, p_glb_ds)
 
        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf, q_cons_temp
 
        ! Down sampling variables
        integer                :: i, j, k, l
        integer                :: ix, iy, iz, x_id, y_id, z_id
        integer, intent(inout) :: m_ds, n_ds, p_ds, m_glb_ds, n_glb_ds, p_glb_ds
 
        m_ds = int((m + 1)/3) - 1
        n_ds = int((n + 1)/3) - 1
        p_ds = int((p + 1)/3) - 1
 
        m_glb_ds = int((m_glb + 1)/3) - 1
        n_glb_ds = int((n_glb + 1)/3) - 1
        p_glb_ds = int((p_glb + 1)/3) - 1
 
        do l = -1, p_ds + 1
            do k = -1, n_ds + 1
                do j = -1, m_ds + 1
                    x_id = 3*j + 1
                    y_id = 3*k + 1
                    z_id = 3*l + 1
                    do i = 1, sys_size
                        q_cons_temp(i)%sf(j, k, l) = 0
 
                        do iz = -1, 1
                            do iy = -1, 1
                                do ix = -1, 1
                                    q_cons_temp(i)%sf(j, k, l) = q_cons_temp(i)%sf(j, k, &
                                                & l) + (1._wp/27._wp)*q_cons_vf(i)%sf(x_id + ix, y_id + iy, z_id + iz)
                                end do
                            end do
                        end do
                    end do
                end do
            end do
        end do
 
    end subroutine s_downsample_data
 
    !> Upsample conservative variable fields from a coarsened grid back to the original resolution using interpolation.
    subroutine s_upsample_data(q_cons_vf, q_cons_temp)
 
        type(scalar_field), intent(inout), dimension(sys_size) :: q_cons_vf, q_cons_temp
        integer                                                :: i, j, k, l
        integer                                                :: ix, iy, iz
        integer                                                :: x_id, y_id, z_id
        real(wp), dimension(4)                                 :: temp
 
        do l = 0, p
            do k = 0, n
                do j = 0, m
                    do i = 1, sys_size
                        ix = int(j/3._wp)
                        iy = int(k/3._wp)
                        iz = int(l/3._wp)
 
                        x_id = j - int(3*ix) - 1
                        y_id = k - int(3*iy) - 1
                        z_id = l - int(3*iz) - 1
 
                        temp(1) = (2._wp/3._wp)*q_cons_temp(i)%sf(ix, iy, iz) + (1._wp/3._wp)*q_cons_temp(i)%sf(ix + x_id, iy, iz)
                        temp(2) = (2._wp/3._wp)*q_cons_temp(i)%sf(ix, iy + y_id, iz) + (1._wp/3._wp)*q_cons_temp(i)%sf(ix + x_id, &
                             & iy + y_id, iz)
                        temp(3) = (2._wp/3._wp)*temp(1) + (1._wp/3._wp)*temp(2)
 
                        temp(1) = (2._wp/3._wp)*q_cons_temp(i)%sf(ix, iy, iz + z_id) + (1._wp/3._wp)*q_cons_temp(i)%sf(ix + x_id, &
                             & iy, iz + z_id)
                        temp(2) = (2._wp/3._wp)*q_cons_temp(i)%sf(ix, iy + y_id, &
                             & iz + z_id) + (1._wp/3._wp)*q_cons_temp(i)%sf(ix + x_id, iy + y_id, iz + z_id)
                        temp(4) = (2._wp/3._wp)*temp(1) + (1._wp/3._wp)*temp(2)
 
                        q_cons_vf(i)%sf(j, k, l) = (2._wp/3._wp)*temp(3) + (1._wp/3._wp)*temp(4)
                    end do
                end do
            end do
        end do
 
    end subroutine s_upsample_data
 
end module m_helper