diff --git a/CMakeLists.txt b/CMakeLists.txt
index 22268e2502..de84d7915d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -147,11 +147,11 @@ if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
     endif()
 elseif (CMAKE_Fortran_COMPILER_ID STREQUAL "Cray")
     add_compile_options(
-        "SHELL:-h nomessage=296:878:1391:1069"
+        "SHELL:-M 296,878,1391,1069,5025"
         "SHELL:-h static" "SHELL:-h keepfiles"
         "SHELL:-h acc_model=auto_async_none"
         "SHELL: -h acc_model=no_fast_addr"
-        "SHELL: -h list=adm" "-DCRAY_ACC_SIMPLIFY" "-DCRAY_ACC_WAR"
+        "SHELL: -h list=adm"
     )
 
     add_link_options("SHELL:-hkeepfiles")
diff --git a/src/common/include/macros.fpp b/src/common/include/macros.fpp
index 27a1bf3385..bcf7eaf78e 100644
--- a/src/common/include/macros.fpp
+++ b/src/common/include/macros.fpp
@@ -13,71 +13,31 @@
 #:def ALLOCATE(*args)
     @:LOG({'@:ALLOCATE(${re.sub(' +', ' ', ', '.join(args))}$)'})
     allocate (${', '.join(args)}$)
-#ifndef CRAY_ACC_WAR
-!$acc enter data create(${', '.join(args)}$)
-#endif
+    !$acc enter data create(${', '.join(args)}$)
 #:enddef ALLOCATE
 
 #:def DEALLOCATE(*args)
     @:LOG({'@:DEALLOCATE(${re.sub(' +', ' ', ', '.join(args))}$)'})
     deallocate (${', '.join(args)}$)
-#ifndef CRAY_ACC_WAR
-!$acc exit data delete(${', '.join(args)}$)
-#endif
+    !$acc exit data delete(${', '.join(args)}$)
 #:enddef DEALLOCATE
 
 #:def ALLOCATE_GLOBAL(*args)
     @:LOG({'@:ALLOCATE_GLOBAL(${re.sub(' +', ' ', ', '.join(args))}$)'})
-#ifdef CRAY_ACC_WAR
-    allocate (${', '.join(('p_' + arg.strip() for arg in args))}$)
-    #:for arg in args
-        ${re.sub('\\(.*\\)','',arg)}$ => ${ 'p_' + re.sub('\\(.*\\)','',arg.strip()) }$
-    #:endfor
-    !$acc enter data create(${', '.join(('p_' + re.sub('\\(.*\\)','',arg.strip()) for arg in args))}$) &
-    !$acc& attach(${', '.join(map(lambda x: re.sub('\\(.*\\)','',x), args))}$)
-#else
+
     allocate (${', '.join(args)}$)
-    !$acc enter data create(${', '.join(args)}$)
-#endif
 
 #:enddef ALLOCATE_GLOBAL
 
 #:def DEALLOCATE_GLOBAL(*args)
     @:LOG({'@:DEALLOCATE_GLOBAL(${re.sub(' +', ' ', ', '.join(args))}$)'})
-#ifdef CRAY_ACC_WAR
-    !$acc exit data delete(${', '.join(('p_' + arg.strip() for arg in args))}$) &
-    !$acc& detach(${', '.join(args)}$)
-    #:for arg in args
-        nullify (${arg}$)
-    #:endfor
-    deallocate (${', '.join(('p_' + arg.strip() for arg in args))}$)
-#else
+
     deallocate (${', '.join(args)}$)
-    !$acc exit data delete(${', '.join(args)}$)
-#endif
 
 #:enddef DEALLOCATE_GLOBAL
 
-#:def CRAY_DECLARE_GLOBAL(intype, dim, *args)
-#ifdef CRAY_ACC_WAR
-    ${intype}$, ${dim}$, allocatable, target :: ${', '.join(('p_' + arg.strip() for arg in args))}$
-    ${intype}$, ${dim}$, pointer :: ${', '.join(args)}$
-#else
-    ${intype}$, ${dim}$, allocatable :: ${', '.join(args)}$
-#endif
-#:enddef CRAY_DECLARE_GLOBAL
-
-#:def CRAY_DECLARE_GLOBAL_SCALAR(intype, *args)
-#ifdef CRAY_ACC_WAR
-    ${intype}$, target :: ${', '.join(('p_' + arg.strip() for arg in args))}$
-    ${intype}$, pointer :: ${', '.join(args)}$
-#else
-    ${intype}$::${', '.join(args)}$
-#endif
-#:enddef CRAY_DECLARE_GLOBAL_SCALAR
-
 #:def ACC_SETUP_VFs(*args)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
     block
         integer :: macros_setup_vfs_i
 
@@ -100,7 +60,7 @@
 #:enddef
 
 #:def ACC_SETUP_SFs(*args)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
     block
 
         @:LOG({'@:ACC_SETUP_SFs(${', '.join(args)}$)'})
@@ -116,7 +76,7 @@
 #:enddef
 
 #:def ACC_SETUP_source_spatials(*args)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
     block
 
         @:LOG({'@:ACC_SETUP_source_spatials(${', '.join(args)}$)'})
diff --git a/src/common/m_phase_change.fpp b/src/common/m_phase_change.fpp
index 9ae303ff3f..dc03e8a5d7 100644
--- a/src/common/m_phase_change.fpp
+++ b/src/common/m_phase_change.fpp
@@ -34,21 +34,6 @@ module m_phase_change
               s_infinite_relaxation_k, &
               s_finalize_relaxation_solver_module
 
-    !> @name Abstract interface for creating function pointers
-    !> @{
-    abstract interface
-
-        !> @name Abstract subroutine for the infinite relaxation solver
-        !> @{
-        subroutine s_abstract_relaxation_solver(q_cons_vf)
-            import :: scalar_field, sys_size
-            type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
-        end subroutine
-        !> @}
-
-    end interface
-    !> @}
-
     !> @name Parameters for the first order transition phase change
     !> @{
     integer, parameter :: max_iter = 1e8        !< max # of iterations
@@ -66,10 +51,18 @@ module m_phase_change
 
     !$acc declare create(max_iter,pCr,TCr,mixM,lp,vp,A,B,C,D)
 
-    procedure(s_abstract_relaxation_solver), pointer :: s_relaxation_solver => null()
-
 contains
 
+    !> This subroutine should dispatch to the correct relaxation solver based
+        !!      some parameter. It replaces the procedure pointer, which CCE
+        !!      is breaking on.
+    subroutine s_relaxation_solver(q_cons_vf)
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
+        ! This is empty because in current master the procedure pointer
+        ! was never assigned
+        @:ASSERT(.false., "s_relaxation_solver called but it currently does nothing")
+    end subroutine s_relaxation_solver
+
     !>  The purpose of this subroutine is to initialize the phase change module
         !!      by setting the parameters needed for phase change and
         !!      selecting the phase change module that will be used
@@ -298,8 +291,9 @@ contains
         !!  @param rhoe mixture energy
         !!  @param TS equilibrium temperature at the interface
     subroutine s_infinite_pt_relaxation_k(j, k, l, MFL, pS, p_infpT, rM, q_cons_vf, rhoe, TS)
-#ifdef CRAY_ACC_WAR
-        !DIR$ INLINEALWAYS s_compute_speed_of_sound
+
+#ifdef _CRAYFTN
+        !DIR$ INLINEALWAYS s_infinite_pt_relaxation_k
 #else
         !$acc routine seq
 #endif
@@ -406,7 +400,7 @@ contains
         !!  @param TS equilibrium temperature at the interface
     subroutine s_infinite_ptg_relaxation_k(j, k, l, pS, p_infpT, rhoe, q_cons_vf, TS)
 
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_infinite_ptg_relaxation_k
 #else
         !$acc routine seq
@@ -530,7 +524,8 @@ contains
         !!  @param k generic loop iterator for y direction
         !!  @param l generic loop iterator for z direction
     subroutine s_correct_partial_densities(MCT, q_cons_vf, rM, j, k, l)
-#ifdef CRAY_ACC_WAR
+
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_correct_partial_densities
 #else
         !$acc routine seq
@@ -593,7 +588,7 @@ contains
         !!  @param TJac Transpose of the Jacobian Matrix
     subroutine s_compute_jacobian_matrix(InvJac, j, Jac, k, l, mCPD, mCVGP, mCVGP2, pS, q_cons_vf, TJac)
 
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_compute_jacobian_matrix
 #else
         !$acc routine seq
@@ -700,7 +695,7 @@ contains
         !!  @param R2D (2D) residue array
     subroutine s_compute_pTg_residue(j, k, l, mCPD, mCVGP, mQD, q_cons_vf, pS, rhoe, R2D)
 
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_compute_pTg_residue
 #else
         !$acc routine seq
@@ -750,8 +745,9 @@ contains
         !!  @param TSat Saturation Temperature
         !!  @param TSIn equilibrium Temperature
     subroutine s_TSat(pSat, TSat, TSIn)
-#ifdef CRAY_ACC_WAR
-        !DIR$ INLINEALWAYS s_compute_speed_of_sound
+
+#ifdef _CRAYFTN
+        !DIR$ INLINEALWAYS s_TSat
 #else
         !$acc routine seq
 #endif
diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp
index 09462ce2ca..b919d6e0b1 100644
--- a/src/common/m_variables_conversion.fpp
+++ b/src/common/m_variables_conversion.fpp
@@ -49,57 +49,17 @@ module m_variables_conversion
 #endif
               s_finalize_variables_conversion_module
 
-    !> Abstract interface to two subroutines designed for the transfer/conversion
-    !! of the mixture/species variables to the mixture variables
-
-    abstract interface ! =======================================================
-
-        !> Structure of the s_convert_mixture_to_mixture_variables
-        !!      and s_convert_species_to_mixture_variables subroutines
-        !!  @param q_vf Conservative or primitive variables
-        !!  @param i First-coordinate cell index
-        !!  @param j First-coordinate cell index
-        !!  @param k First-coordinate cell index
-        !!  @param rho Density
-        !!  @param gamma Specific heat ratio function
-        !!  @param pi_inf Liquid stiffness function
-        !!  @param qv Fluid reference energy
-        subroutine s_convert_xxxxx_to_mixture_variables(q_vf, i, j, k, &
-                                                        rho, gamma, pi_inf, qv, Re_K, G_K, G)
-
-            ! Importing the derived type scalar_field from m_derived_types.f90
-            ! and global variable sys_size, from m_global_variables.f90, as
-            ! the abstract interface does not inherently have access to them
-            import :: scalar_field, sys_size, num_fluids
-
-            type(scalar_field), dimension(sys_size), intent(in) :: q_vf
-            integer, intent(in) :: i, j, k
-            real(kind(0d0)), intent(out), target :: rho, gamma, pi_inf, qv
-            real(kind(0d0)), optional, dimension(2), intent(out) :: Re_K
-            real(kind(0d0)), optional, intent(out) :: G_K
-            real(kind(0d0)), optional, dimension(num_fluids), intent(in) :: G
-
-        end subroutine s_convert_xxxxx_to_mixture_variables
-
-    end interface ! ============================================================
-
     !! In simulation, gammas, pi_infs, and qvs are already declared in m_global_variables
 #ifndef MFC_SIMULATION
     real(kind(0d0)), allocatable, public, dimension(:) :: gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps
     !$acc declare create(gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps)
 #endif
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), Gs)
-    @:CRAY_DECLARE_GLOBAL(integer,         dimension(:), bubrs)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :), Res)
-    !$acc declare link(bubrs, Gs, Res)
-#else
     real(kind(0d0)), allocatable, dimension(:) :: Gs
     integer, allocatable, dimension(:) :: bubrs
     real(kind(0d0)), allocatable, dimension(:, :) :: Res
     !$acc declare create(bubrs, Gs, Res)
-#endif
+
     integer :: is1b, is2b, is3b, is1e, is2e, is3e
     !$acc declare create(is1b, is2b, is3b, is1e, is2e, is3e)
 
@@ -108,13 +68,44 @@ module m_variables_conversion
     real(kind(0d0)), allocatable, dimension(:, :, :), public :: pi_inf_sf !< Scalar liquid stiffness function
     real(kind(0d0)), allocatable, dimension(:, :, :), public :: qv_sf !< Scalar liquid energy reference function
 
-    procedure(s_convert_xxxxx_to_mixture_variables), &
-        pointer :: s_convert_to_mixture_variables => null() !<
-    !! Pointer referencing the subroutine s_convert_mixture_to_mixture_variables
-    !! or s_convert_species_to_mixture_variables, based on model equations choice
-
 contains
 
+    !> Dispatch to the s_convert_mixture_to_mixture_variables
+        !!      and s_convert_species_to_mixture_variables subroutines.
+        !!      Replaces a procedure pointer.
+        !!  @param q_vf Conservative or primitive variables
+        !!  @param i First-coordinate cell index
+        !!  @param j First-coordinate cell index
+        !!  @param k First-coordinate cell index
+        !!  @param rho Density
+        !!  @param gamma Specific heat ratio function
+        !!  @param pi_inf Liquid stiffness function
+        !!  @param qv Fluid reference energy
+    subroutine s_convert_to_mixture_variables(q_vf, i, j, k, &
+                                              rho, gamma, pi_inf, qv, Re_K, G_K, G)
+
+        type(scalar_field), dimension(sys_size), intent(in) :: q_vf
+        integer, intent(in) :: i, j, k
+        real(kind(0d0)), intent(out), target :: rho, gamma, pi_inf, qv
+        real(kind(0d0)), optional, dimension(2), intent(out) :: Re_K
+        real(kind(0d0)), optional, intent(out) :: G_K
+        real(kind(0d0)), optional, dimension(num_fluids), intent(in) :: G
+
+        if (model_eqns == 1) then        ! Gamma/pi_inf model
+            call s_convert_mixture_to_mixture_variables(q_vf, i, j, k, &
+                                                        rho, gamma, pi_inf, qv, Re_K, G_K, G)
+
+        else if (bubbles) then
+            call s_convert_species_to_mixture_variables_bubbles(q_vf, i, j, k, &
+                                                                rho, gamma, pi_inf, qv, Re_K, G_K, G)
+        else
+            ! Volume fraction model
+            call s_convert_species_to_mixture_variables(q_vf, i, j, k, &
+                                                        rho, gamma, pi_inf, qv, Re_K, G_K, G)
+        end if
+
+    end subroutine s_convert_to_mixture_variables
+
     !>  This procedure conditionally calculates the appropriate pressure
         !! @param energy Energy
         !! @param alf Void Fraction
@@ -128,7 +119,7 @@ contains
         !! @param mom Momentum
     subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, qv, rhoYks, pres, T, stress, mom, G)
 
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_compute_pressure
 #else
         !$acc routine seq
@@ -473,7 +464,7 @@ contains
                                                           gamma_K, pi_inf_K, qv_K, &
                                                           alpha_K, alpha_rho_K, Re_K, k, l, r, &
                                                           G_K, G)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_convert_species_to_mixture_variables_acc
 #else
         !$acc routine seq
@@ -555,7 +546,7 @@ contains
     subroutine s_convert_species_to_mixture_variables_bubbles_acc(rho_K, &
                                                                   gamma_K, pi_inf_K, qv_K, &
                                                                   alpha_K, alpha_rho_K, Re_K, k, l, r)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_convert_species_to_mixture_variables_bubbles_acc
 #else
         !$acc routine seq
@@ -748,18 +739,6 @@ contains
         end if
 #endif
 
-        if (model_eqns == 1) then        ! Gamma/pi_inf model
-            s_convert_to_mixture_variables => &
-                s_convert_mixture_to_mixture_variables
-
-        else if (bubbles) then
-            s_convert_to_mixture_variables => &
-                s_convert_species_to_mixture_variables_bubbles
-        else
-            ! Volume fraction model
-            s_convert_to_mixture_variables => &
-                s_convert_species_to_mixture_variables
-        end if
     end subroutine s_initialize_variables_conversion_module
 
     !Initialize mv at the quadrature nodes based on the initialized moments and sigma
@@ -1393,15 +1372,11 @@ contains
         end if
 #endif
 
-        ! Nullifying the procedure pointer to the subroutine transferring/
-        ! computing the mixture/species variables to the mixture variables
-        s_convert_to_mixture_variables => null()
-
     end subroutine s_finalize_variables_conversion_module
 
 #ifndef MFC_PRE_PROCESS
-    subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c)
-#ifdef CRAY_ACC_WAR
+    pure subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c)
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_compute_speed_of_sound
 #else
         !$acc routine seq
diff --git a/src/simulation/m_acoustic_src.fpp b/src/simulation/m_acoustic_src.fpp
index 02ee735091..6f78f2edcc 100644
--- a/src/simulation/m_acoustic_src.fpp
+++ b/src/simulation/m_acoustic_src.fpp
@@ -23,41 +23,6 @@ module m_acoustic_src
     implicit none
     private; public :: s_initialize_acoustic_src, s_precalculate_acoustic_spatial_sources, s_acoustic_src_calculations
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(integer, dimension(:), pulse, support)
-    !$acc declare link(pulse, support)
-
-    @:CRAY_DECLARE_GLOBAL(logical, dimension(:), dipole)
-    !$acc declare link(dipole)
-
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :), loc_acoustic)
-    !$acc declare link(loc_acoustic)
-
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), mag, length, height, wavelength, frequency, gauss_sigma_dist, gauss_sigma_time, npulse, dir, delay)
-    !$acc declare link(mag, length, height, wavelength, frequency, gauss_sigma_dist, gauss_sigma_time, npulse, dir, delay)
-
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), foc_length, aperture)
-    !$acc declare link(foc_length, aperture)
-
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), element_spacing_angle, element_polygon_ratio, rotate_angle)
-    !$acc declare link(element_spacing_angle, element_polygon_ratio, rotate_angle)
-
-    @:CRAY_DECLARE_GLOBAL(integer, dimension(:), num_elements, element_on)
-    !$acc declare link(num_elements, element_on)
-
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), mass_src, e_src)
-    !$acc declare link(mass_src, e_src)
-
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), mom_src)
-    !$acc declare link(mom_src)
-
-    @:CRAY_DECLARE_GLOBAL(integer, dimension(:), source_spatials_num_points)
-    !$acc declare link(source_spatials_num_points)
-
-    @:CRAY_DECLARE_GLOBAL(type(source_spatial_type), dimension(:), source_spatials)
-    !$acc declare link(source_spatials)
-
-#else
     integer, allocatable, dimension(:) :: pulse, support
     !$acc declare create(pulse, support)
 
@@ -92,8 +57,6 @@ module m_acoustic_src
     type(source_spatial_type), dimension(:), allocatable :: source_spatials !< Data of non-zero source grid points for each source
     !$acc declare create(source_spatials)
 
-#endif
-
 contains
 
     !> This subroutine initializes the acoustic source module
diff --git a/src/simulation/m_body_forces.fpp b/src/simulation/m_body_forces.fpp
index 490f0f45bb..719a4f56af 100644
--- a/src/simulation/m_body_forces.fpp
+++ b/src/simulation/m_body_forces.fpp
@@ -24,13 +24,8 @@ module m_body_forces
               s_initialize_body_forces_module, &
               s_finalize_body_forces_module
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), rhoM)
-    !$acc declare link(rhoM)
-#else
     real(kind(0d0)), allocatable, dimension(:, :, :) :: rhoM
     !$acc declare create(rhoM)
-#endif
 
 contains
 
diff --git a/src/simulation/m_bubbles.fpp b/src/simulation/m_bubbles.fpp
index 804b473b8f..1dc0fdcf79 100644
--- a/src/simulation/m_bubbles.fpp
+++ b/src/simulation/m_bubbles.fpp
@@ -24,24 +24,8 @@ module m_bubbles
     real(kind(0.d0)) :: chi_vw  !< Bubble wall properties (Ando 2010)
     real(kind(0.d0)) :: k_mw    !< Bubble wall properties (Ando 2010)
     real(kind(0.d0)) :: rho_mw  !< Bubble wall properties (Ando 2010)
-!$acc declare create(chi_vw, k_mw, rho_mw)
+    !$acc declare create(chi_vw, k_mw, rho_mw)
 
-#ifdef CRAY_ACC_WAR
-    !> @name Bubble dynamic source terms
-    !> @{
-
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), bub_adv_src)
-    !$acc declare link(bub_adv_src)
-
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), bub_r_src, bub_v_src, bub_p_src, bub_m_src)
-    !$acc declare link(bub_r_src, bub_v_src, bub_p_src, bub_m_src)
-
-    type(scalar_field) :: divu !< matrix for div(u)
-    !$acc declare create(divu)
-
-    @:CRAY_DECLARE_GLOBAL(integer, dimension(:), rs, vs, ms, ps)
-    !$acc declare link(rs, vs, ms, ps)
-#else
     real(kind(0d0)), allocatable, dimension(:, :, :) :: bub_adv_src
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: bub_r_src, bub_v_src, bub_p_src, bub_m_src
     !$acc declare create(bub_adv_src, bub_r_src, bub_v_src, bub_p_src, bub_m_src)
@@ -51,7 +35,6 @@ module m_bubbles
 
     integer, allocatable, dimension(:) :: rs, vs, ms, ps
     !$acc declare create(rs, vs, ms, ps)
-#endif
 
 contains
 
diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp
index 05ad8e9da1..07e01d41da 100644
--- a/src/simulation/m_cbc.fpp
+++ b/src/simulation/m_cbc.fpp
@@ -39,87 +39,55 @@ module m_cbc
     !! The cell-average primitive variables. They are obtained by reshaping (RS)
     !! q_prim_vf in the coordinate direction normal to the domain boundary along
     !! which the CBC is applied.
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), q_prim_rsx_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), q_prim_rsy_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), q_prim_rsz_vf)
-    !$acc declare link(q_prim_rsx_vf, q_prim_rsy_vf, q_prim_rsz_vf)
-#else
+
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: q_prim_rsx_vf
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: q_prim_rsy_vf
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: q_prim_rsz_vf
-#endif
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(type(scalar_field), dimension(:), F_rs_vf, F_src_rs_vf)
-    !$acc declare link(F_rs_vf, F_src_rs_vf)
-#else
     type(scalar_field), allocatable, dimension(:) :: F_rs_vf, F_src_rs_vf !<
-#endif
+
     !! Cell-average fluxes (src - source). These are directly determined from the
     !! cell-average primitive variables, q_prims_rs_vf, and not a Riemann solver.
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), F_rsx_vf, F_src_rsx_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), F_rsy_vf, F_src_rsy_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), F_rsz_vf, F_src_rsz_vf)
-    !$acc declare link(F_rsx_vf, F_src_rsx_vf, F_rsy_vf, F_src_rsy_vf, F_rsz_vf, F_src_rsz_vf)
-#else
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: F_rsx_vf, F_src_rsx_vf !<
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: F_rsy_vf, F_src_rsy_vf !<
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: F_rsz_vf, F_src_rsz_vf !<
-#endif
-
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), flux_rsx_vf, flux_src_rsx_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), flux_rsy_vf, flux_src_rsy_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), flux_rsz_vf, flux_src_rsz_vf)
-    !$acc declare link(flux_rsx_vf, flux_src_rsx_vf, flux_rsy_vf, flux_src_rsy_vf, flux_rsz_vf, flux_src_rsz_vf)
-#else
-    real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_rsx_vf, flux_src_rsx_vf !<
-    real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_rsy_vf, flux_src_rsy_vf
-    real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_rsz_vf, flux_src_rsz_vf
-#endif
+
+    !! There is a CCE bug that is causing some subset of these variables to interfere
+    !! with variables of the same name in m_riemann_solvers.fpp, and giving this versions
+    !! unique "_l" names works around the bug. Other private module allocatable arrays
+    !! in `acc declare create` clauses don't have this problem, so we still need to
+    !! isolate this bug.
+
+    real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_rsx_vf_l, flux_src_rsx_vf_l !<
+    real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_rsy_vf_l, flux_src_rsy_vf_l
+    real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_rsz_vf_l, flux_src_rsz_vf_l
 
     real(kind(0d0)) :: c           !< Cell averaged speed of sound
     real(kind(0d0)), dimension(2) :: Re          !< Cell averaged Reynolds numbers
     !$acc declare create(c, Re)
 
     real(kind(0d0)) :: dpres_ds !< Spatial derivatives in s-dir of pressure
-!$acc declare create(dpres_ds)
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), ds)
-    !$acc declare link(ds)
-#else
+    !$acc declare create(dpres_ds)
+
     real(kind(0d0)), allocatable, dimension(:) :: ds !< Cell-width distribution in the s-direction
-#endif
 
     ! CBC Coefficients =========================================================
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :), fd_coef_x)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :), fd_coef_y)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :), fd_coef_z)
-    !$acc declare link(fd_coef_x, fd_coef_y, fd_coef_z)
-#else
+
     real(kind(0d0)), allocatable, dimension(:, :) :: fd_coef_x !< Finite diff. coefficients x-dir
     real(kind(0d0)), allocatable, dimension(:, :) :: fd_coef_y !< Finite diff. coefficients y-dir
     real(kind(0d0)), allocatable, dimension(:, :) :: fd_coef_z !< Finite diff. coefficients z-dir
-#endif
+
     !! The first dimension identifies the location of a coefficient in the FD
     !! formula, while the last dimension denotes the location of the CBC.
 
     ! Bug with NVHPC when using nullified pointers in a declare create
     !    real(kind(0d0)), pointer, dimension(:, :) :: fd_coef => null()
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), pi_coef_x)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), pi_coef_y)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), pi_coef_z)
-    !$acc declare link(pi_coef_x, pi_coef_y, pi_coef_z)
-#else
+
     real(kind(0d0)), allocatable, dimension(:, :, :) :: pi_coef_x !< Polynomial interpolant coefficients in x-dir
     real(kind(0d0)), allocatable, dimension(:, :, :) :: pi_coef_y !< Polynomial interpolant coefficients in y-dir
     real(kind(0d0)), allocatable, dimension(:, :, :) :: pi_coef_z !< Polynomial interpolant coefficients in z-dir
-#endif
+
     !! The first dimension of the array identifies the polynomial, the
     !! second dimension identifies the position of its coefficients and the last
     !! dimension denotes the location of the CBC.
@@ -132,14 +100,12 @@ module m_cbc
     integer :: dj
     integer :: bcxb, bcxe, bcyb, bcye, bczb, bcze
     integer :: cbc_dir, cbc_loc
-!$acc declare create(dj, bcxb, bcxe, bcyb, bcye, bczb, bcze, cbc_dir, cbc_loc)
+    !$acc declare create(dj, bcxb, bcxe, bcyb, bcye, bczb, bcze, cbc_dir, cbc_loc)
 
-#ifndef CRAY_ACC_WAR
-!$acc declare create(q_prim_rsx_vf, q_prim_rsy_vf, q_prim_rsz_vf,  F_rsx_vf, F_src_rsx_vf,flux_rsx_vf, flux_src_rsx_vf, &
-!$acc                 F_rsy_vf, F_src_rsy_vf,flux_rsy_vf, flux_src_rsy_vf, F_rsz_vf, F_src_rsz_vf,flux_rsz_vf, flux_src_rsz_vf,Re, &
-!$acc                 ds,fd_coef_x,fd_coef_y,fd_coef_z,      &
-!$acc                 pi_coef_x,pi_coef_y,pi_coef_z)
-#endif
+    !$acc declare create(q_prim_rsx_vf, q_prim_rsy_vf, q_prim_rsz_vf,  F_rsx_vf, F_src_rsx_vf,flux_rsx_vf_l, flux_src_rsx_vf_l, &
+    !$acc                 F_rsy_vf, F_src_rsy_vf,flux_rsy_vf_l, flux_src_rsy_vf_l, F_rsz_vf, F_src_rsz_vf,flux_rsz_vf_l, flux_src_rsz_vf_l, &
+    !$acc                 ds,fd_coef_x,fd_coef_y,fd_coef_z,      &
+    !$acc                 pi_coef_x,pi_coef_y,pi_coef_z)
 
 contains
 
@@ -188,11 +154,11 @@ contains
 
         end if
 
-        @:ALLOCATE_GLOBAL(flux_rsx_vf(-1:buff_size, &
+        @:ALLOCATE_GLOBAL(flux_rsx_vf_l(-1:buff_size, &
             is2%beg:is2%end, &
             is3%beg:is3%end, 1:adv_idx%end))
 
-        @:ALLOCATE_GLOBAL(flux_src_rsx_vf(-1:buff_size, &
+        @:ALLOCATE_GLOBAL(flux_src_rsx_vf_l(-1:buff_size, &
             is2%beg:is2%end, &
             is3%beg:is3%end, adv_idx%beg:adv_idx%end))
 
@@ -231,11 +197,11 @@ contains
 
             end if
 
-            @:ALLOCATE_GLOBAL(flux_rsy_vf(-1:buff_size, &
+            @:ALLOCATE_GLOBAL(flux_rsy_vf_l(-1:buff_size, &
                 is2%beg:is2%end, &
                 is3%beg:is3%end, 1:adv_idx%end))
 
-            @:ALLOCATE_GLOBAL(flux_src_rsy_vf(-1:buff_size, &
+            @:ALLOCATE_GLOBAL(flux_src_rsy_vf_l(-1:buff_size, &
                 is2%beg:is2%end, &
                 is3%beg:is3%end, adv_idx%beg:adv_idx%end))
 
@@ -276,11 +242,11 @@ contains
 
             end if
 
-            @:ALLOCATE_GLOBAL(flux_rsz_vf(-1:buff_size, &
+            @:ALLOCATE_GLOBAL(flux_rsz_vf_l(-1:buff_size, &
                 is2%beg:is2%end, &
                 is3%beg:is3%end, 1:adv_idx%end))
 
-            @:ALLOCATE_GLOBAL(flux_src_rsz_vf(-1:buff_size, &
+            @:ALLOCATE_GLOBAL(flux_src_rsz_vf_l(-1:buff_size, &
                 is2%beg:is2%end, &
                 is3%beg:is3%end, adv_idx%beg:adv_idx%end))
 
@@ -692,10 +658,10 @@ contains
                     do i = 1, advxe
                         do r = is3%beg, is3%end
                             do k = is2%beg, is2%end
-                                flux_rs${XYZ}$_vf(0, k, r, i) = F_rs${XYZ}$_vf(0, k, r, i) &
-                                                                + pi_coef_${XYZ}$ (0, 0, cbc_loc)* &
-                                                                (F_rs${XYZ}$_vf(1, k, r, i) - &
-                                                                 F_rs${XYZ}$_vf(0, k, r, i))
+                                flux_rs${XYZ}$_vf_l(0, k, r, i) = F_rs${XYZ}$_vf(0, k, r, i) &
+                                                                  + pi_coef_${XYZ}$ (0, 0, cbc_loc)* &
+                                                                  (F_rs${XYZ}$_vf(1, k, r, i) - &
+                                                                   F_rs${XYZ}$_vf(0, k, r, i))
                             end do
                         end do
                     end do
@@ -704,10 +670,10 @@ contains
                     do i = advxb, advxe
                         do r = is3%beg, is3%end
                             do k = is2%beg, is2%end
-                                flux_src_rs${XYZ}$_vf(0, k, r, i) = F_src_rs${XYZ}$_vf(0, k, r, i) + &
-                                                                    (F_src_rs${XYZ}$_vf(1, k, r, i) - &
-                                                                     F_src_rs${XYZ}$_vf(0, k, r, i)) &
-                                                                    *pi_coef_${XYZ}$ (0, 0, cbc_loc)
+                                flux_src_rs${XYZ}$_vf_l(0, k, r, i) = F_src_rs${XYZ}$_vf(0, k, r, i) + &
+                                                                      (F_src_rs${XYZ}$_vf(1, k, r, i) - &
+                                                                       F_src_rs${XYZ}$_vf(0, k, r, i)) &
+                                                                      *pi_coef_${XYZ}$ (0, 0, cbc_loc)
                             end do
                         end do
                     end do
@@ -725,16 +691,16 @@ contains
                         do j = 0, 1
                             do r = is3%beg, is3%end
                                 do k = is2%beg, is2%end
-                                    flux_rs${XYZ}$_vf(j, k, r, i) = F_rs${XYZ}$_vf(j, k, r, i) &
-                                                                    + pi_coef_${XYZ}$ (j, 0, cbc_loc)* &
-                                                                    (F_rs${XYZ}$_vf(3, k, r, i) - &
-                                                                     F_rs${XYZ}$_vf(2, k, r, i)) &
-                                                                    + pi_coef_${XYZ}$ (j, 1, cbc_loc)* &
-                                                                    (F_rs${XYZ}$_vf(2, k, r, i) - &
-                                                                     F_rs${XYZ}$_vf(1, k, r, i)) &
-                                                                    + pi_coef_${XYZ}$ (j, 2, cbc_loc)* &
-                                                                    (F_rs${XYZ}$_vf(1, k, r, i) - &
-                                                                     F_rs${XYZ}$_vf(0, k, r, i))
+                                    flux_rs${XYZ}$_vf_l(j, k, r, i) = F_rs${XYZ}$_vf(j, k, r, i) &
+                                                                      + pi_coef_${XYZ}$ (j, 0, cbc_loc)* &
+                                                                      (F_rs${XYZ}$_vf(3, k, r, i) - &
+                                                                       F_rs${XYZ}$_vf(2, k, r, i)) &
+                                                                      + pi_coef_${XYZ}$ (j, 1, cbc_loc)* &
+                                                                      (F_rs${XYZ}$_vf(2, k, r, i) - &
+                                                                       F_rs${XYZ}$_vf(1, k, r, i)) &
+                                                                      + pi_coef_${XYZ}$ (j, 2, cbc_loc)* &
+                                                                      (F_rs${XYZ}$_vf(1, k, r, i) - &
+                                                                       F_rs${XYZ}$_vf(0, k, r, i))
                                 end do
                             end do
                         end do
@@ -745,16 +711,16 @@ contains
                         do j = 0, 1
                             do r = is3%beg, is3%end
                                 do k = is2%beg, is2%end
-                                    flux_src_rs${XYZ}$_vf(j, k, r, i) = F_src_rs${XYZ}$_vf(j, k, r, i) + &
-                                                                        (F_src_rs${XYZ}$_vf(3, k, r, i) - &
-                                                                         F_src_rs${XYZ}$_vf(2, k, r, i)) &
-                                                                        *pi_coef_${XYZ}$ (j, 0, cbc_loc) + &
-                                                                        (F_src_rs${XYZ}$_vf(2, k, r, i) - &
-                                                                         F_src_rs${XYZ}$_vf(1, k, r, i)) &
-                                                                        *pi_coef_${XYZ}$ (j, 1, cbc_loc) + &
-                                                                        (F_src_rs${XYZ}$_vf(1, k, r, i) - &
-                                                                         F_src_rs${XYZ}$_vf(0, k, r, i)) &
-                                                                        *pi_coef_${XYZ}$ (j, 2, cbc_loc)
+                                    flux_src_rs${XYZ}$_vf_l(j, k, r, i) = F_src_rs${XYZ}$_vf(j, k, r, i) + &
+                                                                          (F_src_rs${XYZ}$_vf(3, k, r, i) - &
+                                                                           F_src_rs${XYZ}$_vf(2, k, r, i)) &
+                                                                          *pi_coef_${XYZ}$ (j, 0, cbc_loc) + &
+                                                                          (F_src_rs${XYZ}$_vf(2, k, r, i) - &
+                                                                           F_src_rs${XYZ}$_vf(1, k, r, i)) &
+                                                                          *pi_coef_${XYZ}$ (j, 1, cbc_loc) + &
+                                                                          (F_src_rs${XYZ}$_vf(1, k, r, i) - &
+                                                                           F_src_rs${XYZ}$_vf(0, k, r, i)) &
+                                                                          *pi_coef_${XYZ}$ (j, 2, cbc_loc)
                                 end do
                             end do
                         end do
@@ -937,42 +903,42 @@ contains
                         end if
                         ! ============================================================
 
-                        ! flux_rs_vf and flux_src_rs_vf at j = -1/2 ==================
+                        ! flux_rs_vf_l and flux_src_rs_vf_l at j = -1/2 ==================
                         !$acc loop seq
                         do i = 1, contxe
-                            flux_rs${XYZ}$_vf(-1, k, r, i) = flux_rs${XYZ}$_vf(0, k, r, i) &
-                                                             + ds(0)*dalpha_rho_dt(i)
+                            flux_rs${XYZ}$_vf_l(-1, k, r, i) = flux_rs${XYZ}$_vf_l(0, k, r, i) &
+                                                               + ds(0)*dalpha_rho_dt(i)
                         end do
 
                         !$acc loop seq
                         do i = momxb, momxe
-                            flux_rs${XYZ}$_vf(-1, k, r, i) = flux_rs${XYZ}$_vf(0, k, r, i) &
-                                                             + ds(0)*(vel(i - contxe)*drho_dt &
-                                                                      + rho*dvel_dt(i - contxe))
+                            flux_rs${XYZ}$_vf_l(-1, k, r, i) = flux_rs${XYZ}$_vf_l(0, k, r, i) &
+                                                               + ds(0)*(vel(i - contxe)*drho_dt &
+                                                                        + rho*dvel_dt(i - contxe))
                         end do
 
-                        flux_rs${XYZ}$_vf(-1, k, r, E_idx) = flux_rs${XYZ}$_vf(0, k, r, E_idx) &
-                                                             + ds(0)*(pres*dgamma_dt &
-                                                                      + gamma*dpres_dt &
-                                                                      + dpi_inf_dt &
-                                                                      + dqv_dt &
-                                                                      + rho*vel_dv_dt_sum &
-                                                                      + 5d-1*drho_dt*vel_K_sum)
+                        flux_rs${XYZ}$_vf_l(-1, k, r, E_idx) = flux_rs${XYZ}$_vf_l(0, k, r, E_idx) &
+                                                               + ds(0)*(pres*dgamma_dt &
+                                                                        + gamma*dpres_dt &
+                                                                        + dpi_inf_dt &
+                                                                        + dqv_dt &
+                                                                        + rho*vel_dv_dt_sum &
+                                                                        + 5d-1*drho_dt*vel_K_sum)
 
                         if (riemann_solver == 1) then
                             !$acc loop seq
                             do i = advxb, advxe
-                                flux_rs${XYZ}$_vf(-1, k, r, i) = 0d0
+                                flux_rs${XYZ}$_vf_l(-1, k, r, i) = 0d0
                             end do
 
                             !$acc loop seq
                             do i = advxb, advxe
-                                flux_src_rs${XYZ}$_vf(-1, k, r, i) = &
+                                flux_src_rs${XYZ}$_vf_l(-1, k, r, i) = &
                                     1d0/max(abs(vel(dir_idx(1))), sgm_eps) &
                                     *sign(1d0, vel(dir_idx(1))) &
-                                    *(flux_rs${XYZ}$_vf(0, k, r, i) &
+                                    *(flux_rs${XYZ}$_vf_l(0, k, r, i) &
                                       + vel(dir_idx(1)) &
-                                      *flux_src_rs${XYZ}$_vf(0, k, r, i) &
+                                      *flux_src_rs${XYZ}$_vf_l(0, k, r, i) &
                                       + ds(0)*dadv_dt(i - E_idx))
                             end do
 
@@ -980,17 +946,17 @@ contains
 
                             !$acc loop seq
                             do i = advxb, advxe
-                                flux_rs${XYZ}$_vf(-1, k, r, i) = flux_rs${XYZ}$_vf(0, k, r, i) + &
-                                                                 ds(0)*dadv_dt(i - E_idx)
+                                flux_rs${XYZ}$_vf_l(-1, k, r, i) = flux_rs${XYZ}$_vf_l(0, k, r, i) + &
+                                                                   ds(0)*dadv_dt(i - E_idx)
                             end do
 
                             !$acc loop seq
                             do i = advxb, advxe
-                                flux_src_rs${XYZ}$_vf(-1, k, r, i) = flux_src_rs${XYZ}$_vf(0, k, r, i)
+                                flux_src_rs${XYZ}$_vf_l(-1, k, r, i) = flux_src_rs${XYZ}$_vf_l(0, k, r, i)
                             end do
 
                         end if
-                        ! END: flux_rs_vf and flux_src_rs_vf at j = -1/2 =============
+                        ! END: flux_rs_vf_l and flux_src_rs_vf_l at j = -1/2 =============
 
                     end do
                 end do
@@ -1083,7 +1049,7 @@ contains
                 do r = is3%beg, is3%end
                     do k = is2%beg, is2%end
                         do j = -1, buff_size
-                            flux_rsx_vf(j, k, r, i) = &
+                            flux_rsx_vf_l(j, k, r, i) = &
                                 flux_vf(i)%sf(dj*((m - 1) - 2*j) + j, k, r)* &
                                 sign(1d0, -real(cbc_loc, kind(0d0)))
                         end do
@@ -1095,7 +1061,7 @@ contains
             do r = is3%beg, is3%end
                 do k = is2%beg, is2%end
                     do j = -1, buff_size
-                        flux_rsx_vf(j, k, r, momxb) = &
+                        flux_rsx_vf_l(j, k, r, momxb) = &
                             flux_vf(momxb)%sf(dj*((m - 1) - 2*j) + j, k, r)
                     end do
                 end do
@@ -1107,7 +1073,7 @@ contains
                     do r = is3%beg, is3%end
                         do k = is2%beg, is2%end
                             do j = -1, buff_size
-                                flux_src_rsx_vf(j, k, r, i) = &
+                                flux_src_rsx_vf_l(j, k, r, i) = &
                                     flux_src_vf(i)%sf(dj*((m - 1) - 2*j) + j, k, r)
                             end do
                         end do
@@ -1118,7 +1084,7 @@ contains
                 do r = is3%beg, is3%end
                     do k = is2%beg, is2%end
                         do j = -1, buff_size
-                            flux_src_rsx_vf(j, k, r, advxb) = &
+                            flux_src_rsx_vf_l(j, k, r, advxb) = &
                                 flux_src_vf(advxb)%sf(dj*((m - 1) - 2*j) + j, k, r)* &
                                 sign(1d0, -real(cbc_loc, kind(0d0)))
                         end do
@@ -1159,7 +1125,7 @@ contains
                 do r = is3%beg, is3%end
                     do k = is2%beg, is2%end
                         do j = -1, buff_size
-                            flux_rsy_vf(j, k, r, i) = &
+                            flux_rsy_vf_l(j, k, r, i) = &
                                 flux_vf(i)%sf(k, dj*((n - 1) - 2*j) + j, r)* &
                                 sign(1d0, -real(cbc_loc, kind(0d0)))
                         end do
@@ -1171,7 +1137,7 @@ contains
             do r = is3%beg, is3%end
                 do k = is2%beg, is2%end
                     do j = -1, buff_size
-                        flux_rsy_vf(j, k, r, momxb + 1) = &
+                        flux_rsy_vf_l(j, k, r, momxb + 1) = &
                             flux_vf(momxb + 1)%sf(k, dj*((n - 1) - 2*j) + j, r)
                     end do
                 end do
@@ -1183,7 +1149,7 @@ contains
                     do r = is3%beg, is3%end
                         do k = is2%beg, is2%end
                             do j = -1, buff_size
-                                flux_src_rsy_vf(j, k, r, i) = &
+                                flux_src_rsy_vf_l(j, k, r, i) = &
                                     flux_src_vf(i)%sf(k, dj*((n - 1) - 2*j) + j, r)
                             end do
                         end do
@@ -1194,7 +1160,7 @@ contains
                 do r = is3%beg, is3%end
                     do k = is2%beg, is2%end
                         do j = -1, buff_size
-                            flux_src_rsy_vf(j, k, r, advxb) = &
+                            flux_src_rsy_vf_l(j, k, r, advxb) = &
                                 flux_src_vf(advxb)%sf(k, dj*((n - 1) - 2*j) + j, r)* &
                                 sign(1d0, -real(cbc_loc, kind(0d0)))
                         end do
@@ -1235,7 +1201,7 @@ contains
                 do r = is3%beg, is3%end
                     do k = is2%beg, is2%end
                         do j = -1, buff_size
-                            flux_rsz_vf(j, k, r, i) = &
+                            flux_rsz_vf_l(j, k, r, i) = &
                                 flux_vf(i)%sf(r, k, dj*((p - 1) - 2*j) + j)* &
                                 sign(1d0, -real(cbc_loc, kind(0d0)))
                         end do
@@ -1247,7 +1213,7 @@ contains
             do r = is3%beg, is3%end
                 do k = is2%beg, is2%end
                     do j = -1, buff_size
-                        flux_rsz_vf(j, k, r, momxe) = &
+                        flux_rsz_vf_l(j, k, r, momxe) = &
                             flux_vf(momxe)%sf(r, k, dj*((p - 1) - 2*j) + j)
                     end do
                 end do
@@ -1259,7 +1225,7 @@ contains
                     do r = is3%beg, is3%end
                         do k = is2%beg, is2%end
                             do j = -1, buff_size
-                                flux_src_rsz_vf(j, k, r, i) = &
+                                flux_src_rsz_vf_l(j, k, r, i) = &
                                     flux_src_vf(i)%sf(r, k, dj*((p - 1) - 2*j) + j)
                             end do
                         end do
@@ -1270,7 +1236,7 @@ contains
                 do r = is3%beg, is3%end
                     do k = is2%beg, is2%end
                         do j = -1, buff_size
-                            flux_src_rsz_vf(j, k, r, advxb) = &
+                            flux_src_rsz_vf_l(j, k, r, advxb) = &
                                 flux_src_vf(advxb)%sf(r, k, dj*((p - 1) - 2*j) + j)* &
                                 sign(1d0, -real(cbc_loc, kind(0d0)))
                         end do
@@ -1321,7 +1287,7 @@ contains
                     do k = is2%beg, is2%end
                         do j = -1, buff_size
                             flux_vf(i)%sf(dj*((m - 1) - 2*j) + j, k, r) = &
-                                flux_rsx_vf(j, k, r, i)* &
+                                flux_rsx_vf_l(j, k, r, i)* &
                                 sign(1d0, -real(cbc_loc, kind(0d0)))
                         end do
                     end do
@@ -1332,7 +1298,7 @@ contains
                 do k = is2%beg, is2%end
                     do j = -1, buff_size
                         flux_vf(momxb)%sf(dj*((m - 1) - 2*j) + j, k, r) = &
-                            flux_rsx_vf(j, k, r, momxb)
+                            flux_rsx_vf_l(j, k, r, momxb)
                     end do
                 end do
             end do
@@ -1344,7 +1310,7 @@ contains
                         do k = is2%beg, is2%end
                             do j = -1, buff_size
                                 flux_src_vf(i)%sf(dj*((m - 1) - 2*j) + j, k, r) = &
-                                    flux_src_rsx_vf(j, k, r, i)
+                                    flux_src_rsx_vf_l(j, k, r, i)
                             end do
                         end do
                     end do
@@ -1355,7 +1321,7 @@ contains
                     do k = is2%beg, is2%end
                         do j = -1, buff_size
                             flux_src_vf(advxb)%sf(dj*((m - 1) - 2*j) + j, k, r) = &
-                                flux_src_rsx_vf(j, k, r, advxb)* &
+                                flux_src_rsx_vf_l(j, k, r, advxb)* &
                                 sign(1d0, -real(cbc_loc, kind(0d0)))
                         end do
                     end do
@@ -1372,7 +1338,7 @@ contains
                     do k = is2%beg, is2%end
                         do j = -1, buff_size
                             flux_vf(i)%sf(k, dj*((n - 1) - 2*j) + j, r) = &
-                                flux_rsy_vf(j, k, r, i)* &
+                                flux_rsy_vf_l(j, k, r, i)* &
                                 sign(1d0, -real(cbc_loc, kind(0d0)))
                         end do
                     end do
@@ -1384,7 +1350,7 @@ contains
                 do k = is2%beg, is2%end
                     do j = -1, buff_size
                         flux_vf(momxb + 1)%sf(k, dj*((n - 1) - 2*j) + j, r) = &
-                            flux_rsy_vf(j, k, r, momxb + 1)
+                            flux_rsy_vf_l(j, k, r, momxb + 1)
                     end do
                 end do
             end do
@@ -1396,7 +1362,7 @@ contains
                         do k = is2%beg, is2%end
                             do j = -1, buff_size
                                 flux_src_vf(i)%sf(k, dj*((n - 1) - 2*j) + j, r) = &
-                                    flux_src_rsy_vf(j, k, r, i)
+                                    flux_src_rsy_vf_l(j, k, r, i)
                             end do
                         end do
                     end do
@@ -1407,7 +1373,7 @@ contains
                     do k = is2%beg, is2%end
                         do j = -1, buff_size
                             flux_src_vf(advxb)%sf(k, dj*((n - 1) - 2*j) + j, r) = &
-                                flux_src_rsy_vf(j, k, r, advxb)* &
+                                flux_src_rsy_vf_l(j, k, r, advxb)* &
                                 sign(1d0, -real(cbc_loc, kind(0d0)))
                         end do
                     end do
@@ -1425,7 +1391,7 @@ contains
                     do k = is2%beg, is2%end
                         do j = -1, buff_size
                             flux_vf(i)%sf(r, k, dj*((p - 1) - 2*j) + j) = &
-                                flux_rsz_vf(j, k, r, i)* &
+                                flux_rsz_vf_l(j, k, r, i)* &
                                 sign(1d0, -real(cbc_loc, kind(0d0)))
                         end do
                     end do
@@ -1437,7 +1403,7 @@ contains
                 do k = is2%beg, is2%end
                     do j = -1, buff_size
                         flux_vf(momxe)%sf(r, k, dj*((p - 1) - 2*j) + j) = &
-                            flux_rsz_vf(j, k, r, momxe)
+                            flux_rsz_vf_l(j, k, r, momxe)
                     end do
                 end do
             end do
@@ -1449,7 +1415,7 @@ contains
                         do k = is2%beg, is2%end
                             do j = -1, buff_size
                                 flux_src_vf(i)%sf(r, k, dj*((p - 1) - 2*j) + j) = &
-                                    flux_src_rsz_vf(j, k, r, i)
+                                    flux_src_rsz_vf_l(j, k, r, i)
                             end do
                         end do
                     end do
@@ -1460,7 +1426,7 @@ contains
                     do k = is2%beg, is2%end
                         do j = -1, buff_size
                             flux_src_vf(advxb)%sf(r, k, dj*((p - 1) - 2*j) + j) = &
-                                flux_src_rsz_vf(j, k, r, advxb)* &
+                                flux_src_rsz_vf_l(j, k, r, advxb)* &
                                 sign(1d0, -real(cbc_loc, kind(0d0)))
                         end do
                     end do
@@ -1507,21 +1473,21 @@ contains
         if (weno_order > 1) then
             @:DEALLOCATE_GLOBAL(F_rsx_vf, F_src_rsx_vf)
         end if
-        @:DEALLOCATE_GLOBAL(flux_rsx_vf, flux_src_rsx_vf)
+        @:DEALLOCATE_GLOBAL(flux_rsx_vf_l, flux_src_rsx_vf_l)
 
         if (n > 0) then
             @:DEALLOCATE_GLOBAL(q_prim_rsy_vf)
             if (weno_order > 1) then
                 @:DEALLOCATE_GLOBAL(F_rsy_vf, F_src_rsy_vf)
             end if
-            @:DEALLOCATE_GLOBAL(flux_rsy_vf, flux_src_rsy_vf)
+            @:DEALLOCATE_GLOBAL(flux_rsy_vf_l, flux_src_rsy_vf_l)
         end if
         if (p > 0) then
             @:DEALLOCATE_GLOBAL(q_prim_rsz_vf)
             if (weno_order > 1) then
                 @:DEALLOCATE_GLOBAL(F_rsz_vf, F_src_rsz_vf)
             end if
-            @:DEALLOCATE_GLOBAL(flux_rsz_vf, flux_src_rsz_vf)
+            @:DEALLOCATE_GLOBAL(flux_rsz_vf_l, flux_src_rsz_vf_l)
         end if
 
         ! Deallocating the cell-width distribution in the s-direction
diff --git a/src/simulation/m_compute_cbc.fpp b/src/simulation/m_compute_cbc.fpp
index c4b369945e..6a18c1eb41 100644
--- a/src/simulation/m_compute_cbc.fpp
+++ b/src/simulation/m_compute_cbc.fpp
@@ -28,7 +28,7 @@ contains
         !!      the normal component of velocity is zero at all times,
         !!      while the transverse velocities may be nonzero.
     subroutine s_compute_slip_wall_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_compute_slip_wall_L
 #else
         !$acc routine seq
@@ -58,7 +58,7 @@ contains
         !!      buffer reduces the amplitude of any reflections caused by
         !!      outgoing waves.
     subroutine s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_compute_nonreflecting_subsonic_buffer_L
 #else
         !$acc routine seq
@@ -100,7 +100,7 @@ contains
         !!      CBC assumes an incoming flow and reduces the amplitude of
         !!      any reflections caused by outgoing waves.
     subroutine s_compute_nonreflecting_subsonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_compute_nonreflecting_subsonic_inflow_L
 #else
         !$acc routine seq
@@ -128,7 +128,7 @@ contains
         !!      subsonic CBC presumes an outgoing flow and reduces the
         !!      amplitude of any reflections caused by outgoing waves.
     subroutine s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_compute_nonreflecting_subsonic_outflow_L
 #else
         !$acc routine seq
@@ -170,7 +170,7 @@ contains
         !!      at the boundary is simply advected outward at the fluid
         !!      velocity.
     subroutine s_compute_force_free_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_compute_force_free_subsonic_outflow_L
 #else
         !$acc routine seq
@@ -208,7 +208,7 @@ contains
         !!      subsonic outflow maintains a fixed pressure at the CBC
         !!      boundary in absence of any transverse effects.
     subroutine s_compute_constant_pressure_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_compute_constant_pressure_subsonic_outflow_L
 #else
         !$acc routine seq
@@ -247,7 +247,7 @@ contains
         !!      transverse terms may generate a time dependence at the
         !!      inflow boundary.
     subroutine s_compute_supersonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_compute_supersonic_inflow_L
 #else
         !$acc routine seq
@@ -272,7 +272,7 @@ contains
         !!      flow evolution at the boundary is determined completely
         !!      by the interior data.
     subroutine s_compute_supersonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_compute_supersonic_outflow_L
 #else
         !$acc routine seq
diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp
index dcba87d601..d5f1ed0f36 100644
--- a/src/simulation/m_data_output.fpp
+++ b/src/simulation/m_data_output.fpp
@@ -48,41 +48,11 @@ module m_data_output
               s_close_probe_files, &
               s_finalize_data_output_module
 
-    abstract interface ! ===================================================
-
-        !> Write data files
-        !! @param q_cons_vf Conservative variables
-        !! @param q_prim_vf Primitive variables
-        !! @param t_step Current time step
-        subroutine s_write_abstract_data_files(q_cons_vf, q_prim_vf, t_step)
-
-            import :: scalar_field, sys_size, pres_field
-
-            type(scalar_field), &
-                dimension(sys_size), &
-                intent(in) :: q_cons_vf
-
-            type(scalar_field), &
-                dimension(sys_size), &
-                intent(inout) :: q_prim_vf
-
-            integer, intent(in) :: t_step
-
-        end subroutine s_write_abstract_data_files
-    end interface ! ========================================================
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), icfl_sf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), vcfl_sf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), ccfl_sf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), Rc_sf)
-    !$acc declare link(icfl_sf, vcfl_sf, ccfl_sf, Rc_sf)
-#else
     real(kind(0d0)), allocatable, dimension(:, :, :) :: icfl_sf  !< ICFL stability criterion
     real(kind(0d0)), allocatable, dimension(:, :, :) :: vcfl_sf  !< VCFL stability criterion
     real(kind(0d0)), allocatable, dimension(:, :, :) :: ccfl_sf  !< CCFL stability criterion
     real(kind(0d0)), allocatable, dimension(:, :, :) :: Rc_sf  !< Rc stability criterion
     !$acc declare create(icfl_sf, vcfl_sf, ccfl_sf, Rc_sf)
-#endif
 
     real(kind(0d0)) :: icfl_max_loc, icfl_max_glb !< ICFL stability extrema on local and global grids
     real(kind(0d0)) :: vcfl_max_loc, vcfl_max_glb !< VCFL stability extrema on local and global grids
@@ -98,10 +68,31 @@ module m_data_output
     real(kind(0d0)) :: Rc_min !< Rc criterion maximum
     !> @}
 
-    procedure(s_write_abstract_data_files), pointer :: s_write_data_files => null()
-
 contains
 
+    !> Write data files. Dispatch subroutine that replaces procedure pointer.
+        !! @param q_cons_vf Conservative variables
+        !! @param q_prim_vf Primitive variables
+        !! @param t_step Current time step
+    subroutine s_write_data_files(q_cons_vf, q_prim_vf, t_step)
+
+        type(scalar_field), &
+            dimension(sys_size), &
+            intent(in) :: q_cons_vf
+
+        type(scalar_field), &
+            dimension(sys_size), &
+            intent(inout) :: q_prim_vf
+
+        integer, intent(in) :: t_step
+
+        if (.not. parallel_io) then
+            call s_write_serial_data_files(q_cons_vf, q_prim_vf, t_step)
+        else
+            call s_write_parallel_data_files(q_cons_vf, q_prim_vf, t_step)
+        end if
+    end subroutine s_write_data_files
+
     !>  The purpose of this subroutine is to open a new or pre-
         !!          existing run-time information file and append to it the
         !!      basic header information relevant to current simulation.
@@ -275,7 +266,7 @@ contains
 
         ! Determining local stability criteria extrema at current time-step
 
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !$acc update host(icfl_sf)
 
         if (viscous) then
@@ -1619,26 +1610,6 @@ contains
             Rc_min = 1d3
         end if
 
-        ! Associating the procedural pointer to the appropriate subroutine
-        ! that will be utilized in the conversion to the mixture variables
-
-        if (model_eqns == 1) then        ! Gamma/pi_inf model
-            s_convert_to_mixture_variables => &
-                s_convert_mixture_to_mixture_variables
-        elseif (bubbles) then           ! Volume fraction for bubbles
-            s_convert_to_mixture_variables => &
-                s_convert_species_to_mixture_variables_bubbles
-        else                            ! Volume fraction model
-            s_convert_to_mixture_variables => &
-                s_convert_species_to_mixture_variables
-        end if
-
-        if (parallel_io .neqv. .true.) then
-            s_write_data_files => s_write_serial_data_files
-        else
-            s_write_data_files => s_write_parallel_data_files
-        end if
-
     end subroutine s_initialize_data_output_module
 
     !> Module deallocation and/or disassociation procedures
@@ -1652,11 +1623,6 @@ contains
             @:DEALLOCATE_GLOBAL(vcfl_sf, Rc_sf)
         end if
 
-        ! Disassociating the pointer to the procedure that was utilized to
-        ! to convert mixture or species variables to the mixture variables
-        s_convert_to_mixture_variables => null()
-        s_write_data_files => null()
-
     end subroutine s_finalize_data_output_module
 
 end module m_data_output
diff --git a/src/simulation/m_fftw.fpp b/src/simulation/m_fftw.fpp
index 58fb51be77..3ba557a139 100644
--- a/src/simulation/m_fftw.fpp
+++ b/src/simulation/m_fftw.fpp
@@ -49,19 +49,12 @@ module m_fftw
     !! Filtered complex data in Fourier space
 
 #if defined(MFC_OpenACC)
-!$acc declare create(real_size, cmplx_size, x_size, batch_size, Nfq)
+    !$acc declare create(real_size, cmplx_size, x_size, batch_size, Nfq)
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:),  data_real_gpu)
-    @:CRAY_DECLARE_GLOBAL(complex(kind(0d0)), dimension(:), data_cmplx_gpu)
-    @:CRAY_DECLARE_GLOBAL(complex(kind(0d0)), dimension(:), data_fltr_cmplx_gpu)
-    !$acc declare link(data_real_gpu, data_cmplx_gpu, data_fltr_cmplx_gpu)
-#else
     real(kind(0d0)), allocatable, target :: data_real_gpu(:)
     complex(kind(0d0)), allocatable, target :: data_cmplx_gpu(:)
     complex(kind(0d0)), allocatable, target :: data_fltr_cmplx_gpu(:)
-    !$acc declare create(data_real_gpu, data_cmplx_gpu, data_fltr_cmplx_gpu)
-#endif
+!$acc declare create(data_real_gpu, data_cmplx_gpu, data_fltr_cmplx_gpu)
 
 #if defined(__PGI)
     integer :: fwd_plan_gpu, bwd_plan_gpu
diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp
index 63fc21c73f..519a57886c 100644
--- a/src/simulation/m_global_parameters.fpp
+++ b/src/simulation/m_global_parameters.fpp
@@ -61,40 +61,24 @@ module m_global_parameters
 
     !> @name Cell-boundary (CB) locations in the x-, y- and z-directions, respectively
     !> @{
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), x_cb, y_cb, z_cb)
-#else
     real(kind(0d0)), target, allocatable, dimension(:) :: x_cb, y_cb, z_cb
-#endif
     !> @}
 
     !> @name Cell-center (CC) locations in the x-, y- and z-directions, respectively
     !> @{
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), x_cc, y_cc, z_cc)
-#else
     real(kind(0d0)), target, allocatable, dimension(:) :: x_cc, y_cc, z_cc
-#endif
     !> @}
     !type(bounds_info) :: x_domain, y_domain, z_domain !<
     !! Locations of the domain bounds in the x-, y- and z-coordinate directions
     !> @name Cell-width distributions in the x-, y- and z-directions, respectively
     !> @{
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), dx, dy, dz)
-#else
     real(kind(0d0)), target, allocatable, dimension(:) :: dx, dy, dz
-#endif
     !> @}
 
     real(kind(0d0)) :: dt !< Size of the time-step
 
-#ifdef CRAY_ACC_WAR
-    !$acc declare link(x_cb, y_cb, z_cb, x_cc, y_cc, z_cc, dx, dy, dz)
-    !$acc declare create(m, n, p, dt)
-#else
     !$acc declare create(x_cb, y_cb, z_cb, x_cc, y_cc, z_cc, dx, dy, dz, dt, m, n, p)
-#endif
+
     !> @name Starting time-step iteration, stopping time-step iteration and the number
     !! of time-step iterations between successive solution backups, respectively
     !> @{
@@ -263,18 +247,10 @@ module m_global_parameters
     !! numbers, will be non-negligible.
     !> @{
     integer, dimension(2) :: Re_size
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(integer, dimension(:, :), Re_idx)
-#else
     integer, allocatable, dimension(:, :) :: Re_idx
-#endif
     !> @}
-#ifdef CRAY_ACC_WAR
-    !$acc declare create(Re_size)
-    !$acc declare link(Re_idx)
-#else
+
     !$acc declare create(Re_size, Re_idx)
-#endif
 
     ! The WENO average (WA) flag regulates whether the calculation of any cell-
     ! average spatial derivatives is carried out in each cell by utilizing the
@@ -375,17 +351,12 @@ module m_global_parameters
     real(kind(0d0)) :: Ca       !< Cavitation number
     real(kind(0d0)) :: Web      !< Weber number
     real(kind(0d0)) :: Re_inv   !< Inverse Reynolds number
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), weight)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), R0)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), V0)
-    !$acc declare link(weight, R0, V0)
-#else
+
     real(kind(0d0)), dimension(:), allocatable :: weight !< Simpson quadrature weights
     real(kind(0d0)), dimension(:), allocatable :: R0     !< Bubble sizes
     real(kind(0d0)), dimension(:), allocatable :: V0     !< Bubble velocities
     !$acc declare create(weight, R0, V0)
-#endif
+
     logical :: bubbles      !< Bubbles on/off
     logical :: polytropic   !< Polytropic  switch
     logical :: polydisperse !< Polydisperse bubbles
@@ -394,13 +365,10 @@ module m_global_parameters
 
     integer :: bubble_model !< Gilmore or Keller--Miksis bubble model
     integer :: thermal      !< Thermal behavior. 1 = adiabatic, 2 = isotherm, 3 = transfer
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), ptil)
-    !$acc declare link(ptil)
-#else
+
     real(kind(0d0)), allocatable, dimension(:, :, :) :: ptil  !< Pressure modification
     !$acc declare create(ptil)
-#endif
+
     real(kind(0d0)) :: poly_sigma  !< log normal sigma for polydisperse PDF
 
     logical :: qbmm      !< Quadrature moment method
@@ -415,17 +383,12 @@ module m_global_parameters
         !$acc declare create(nb)
     #:endif
 
-!$acc declare create(R0ref, Ca, Web, Re_inv, bubbles, polytropic, polydisperse, qbmm, nmomsp, nmomtot, R0_type, bubble_model, thermal, poly_sigma, adv_n, adap_dt, pi_fac)
+    !$acc declare create(R0ref, Ca, Web, Re_inv, bubbles, polytropic, polydisperse, qbmm, nmomsp, nmomtot, R0_type, bubble_model, thermal, poly_sigma, adv_n, adap_dt, pi_fac)
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(type(scalar_field), dimension(:), mom_sp)
-    @:CRAY_DECLARE_GLOBAL(type(scalar_field), dimension(:, :, :), mom_3d)
-    !$acc declare link(mom_sp, mom_3d)
-#else
     type(scalar_field), allocatable, dimension(:) :: mom_sp
     type(scalar_field), allocatable, dimension(:, :, :) :: mom_3d
     !$acc declare create(mom_sp, mom_3d)
-#endif
+
     !> @}
 
     type(chemistry_parameters) :: chem_params
@@ -434,16 +397,12 @@ module m_global_parameters
     !> @name Physical bubble parameters (see Ando 2010, Preston 2007)
     !> @{
     real(kind(0d0)) :: R_n, R_v, phi_vn, phi_nv, Pe_c, Tw, pv, M_n, M_v
-!$acc declare create(R_n, R_v, phi_vn, phi_nv, Pe_c, Tw, pv, M_n, M_v)
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), k_n, k_v, pb0, mass_n0, mass_v0, Pe_T)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c, omegaN)
-    !$acc declare link( k_n, k_v, pb0, mass_n0, mass_v0, Pe_T, Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c, omegaN)
-#else
+    !$acc declare create(R_n, R_v, phi_vn, phi_nv, Pe_c, Tw, pv, M_n, M_v)
+
     real(kind(0d0)), dimension(:), allocatable :: k_n, k_v, pb0, mass_n0, mass_v0, Pe_T
     real(kind(0d0)), dimension(:), allocatable :: Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c, omegaN
     !$acc declare create( k_n, k_v, pb0, mass_n0, mass_v0, Pe_T, Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c, omegaN)
-#endif
+
     real(kind(0d0)) :: mul0, ss, gamma_v, mu_v
     real(kind(0d0)) :: gamma_m, gamma_n, mu_n
     real(kind(0d0)) :: gam
@@ -474,34 +433,22 @@ module m_global_parameters
     integer :: strxb, strxe
     integer :: chemxb, chemxe
 
-!$acc declare create(momxb, momxe, advxb, advxe, contxb, contxe, intxb, intxe, bubxb, bubxe, strxb, strxe,  chemxb, chemxe)
+    !$acc declare create(momxb, momxe, advxb, advxe, contxb, contxe, intxb, intxe, bubxb, bubxe, strxb, strxe,  chemxb, chemxe)
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps)
-    !$acc declare link(gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps)
-#else
     real(kind(0d0)), allocatable, dimension(:) :: gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps
     !$acc declare create(gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps)
-#endif
 
     real(kind(0d0)) :: mytime       !< Current simulation time
     real(kind(0d0)) :: finaltime    !< Final simulation time
 
     logical :: weno_flat, riemann_flat, rdma_mpi
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(type(pres_field), dimension(:), pb_ts)
-
-    @:CRAY_DECLARE_GLOBAL(type(pres_field), dimension(:), mv_ts)
-
-    !$acc declare link(pb_ts, mv_ts)
-#else
     type(pres_field), allocatable, dimension(:) :: pb_ts
 
     type(pres_field), allocatable, dimension(:) :: mv_ts
 
     !$acc declare create(pb_ts, mv_ts)
-#endif
+
     ! ======================================================================
 
 contains
diff --git a/src/simulation/m_hypoelastic.fpp b/src/simulation/m_hypoelastic.fpp
index e3bb7ec08f..94a57b148e 100644
--- a/src/simulation/m_hypoelastic.fpp
+++ b/src/simulation/m_hypoelastic.fpp
@@ -22,18 +22,6 @@ module m_hypoelastic
     private; public :: s_initialize_hypoelastic_module, &
  s_compute_hypoelastic_rhs
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), Gs)
-    !$acc declare link(Gs)
-
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), du_dx, du_dy, du_dz)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), dv_dx, dv_dy, dv_dz)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), dw_dx, dw_dy, dw_dz)
-    !$acc declare link(du_dx,du_dy,du_dz,dv_dx,dv_dy,dv_dz,dw_dx,dw_dy,dw_dz)
-
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), rho_K_field, G_K_field)
-    !$acc declare link(rho_K_field, G_K_field)
-#else
     real(kind(0d0)), allocatable, dimension(:) :: Gs
     !$acc declare create(Gs)
 
@@ -45,8 +33,6 @@ module m_hypoelastic
     real(kind(0d0)), allocatable, dimension(:, :, :) :: rho_K_field, G_K_field
     !$acc declare create(rho_K_field, G_K_field)
 
-#endif
-
 contains
 
     subroutine s_initialize_hypoelastic_module
diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp
index cf64495355..a6dce0fb6d 100644
--- a/src/simulation/m_ibm.fpp
+++ b/src/simulation/m_ibm.fpp
@@ -37,16 +37,7 @@ module m_ibm
  s_finalize_ibm_module
 
     type(integer_field), public :: ib_markers
-!$acc declare create(ib_markers)
-
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), levelset)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :, :), levelset_norm)
-    @:CRAY_DECLARE_GLOBAL(type(ghost_point), dimension(:), ghost_points)
-    @:CRAY_DECLARE_GLOBAL(type(ghost_point), dimension(:), inner_points)
-
-    !$acc declare link(levelset, levelset_norm, ghost_points, inner_points)
-#else
+    !$acc declare create(ib_markers)
 
     !! Marker for solid cells. 0 if liquid, the patch id of its IB if solid
     real(kind(0d0)), dimension(:, :, :, :), allocatable :: levelset
@@ -58,7 +49,6 @@ module m_ibm
     !! Matrix of normal vector to IB
 
     !$acc declare create(levelset, levelset_norm, ghost_points, inner_points)
-#endif
 
     integer :: gp_layers !< Number of ghost point layers
     integer :: num_gps !< Number of ghost points
diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp
index bc0552b1d0..96860ef5e2 100644
--- a/src/simulation/m_mpi_proxy.fpp
+++ b/src/simulation/m_mpi_proxy.fpp
@@ -32,17 +32,6 @@ module m_mpi_proxy
 
     implicit none
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), q_cons_buff_send)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), q_cons_buff_recv)
-    @:CRAY_DECLARE_GLOBAL(integer, dimension(:), ib_buff_send)
-    @:CRAY_DECLARE_GLOBAL(integer, dimension(:), ib_buff_recv)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), c_divs_buff_send)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), c_divs_buff_recv)
-    !$acc declare link(q_cons_buff_recv, q_cons_buff_send)
-    !$acc declare link(ib_buff_send, ib_buff_recv)
-    !$acc declare link(c_divs_buff_send, c_divs_buff_recv)
-#else
     real(kind(0d0)), private, allocatable, dimension(:), target :: q_cons_buff_send !<
     !! This variable is utilized to pack and send the buffer of the cell-average
     !! conservative variables, for a single computational domain boundary at the
@@ -76,7 +65,7 @@ module m_mpi_proxy
     !$acc declare create(q_cons_buff_send, q_cons_buff_recv)
     !$acc declare create( ib_buff_send, ib_buff_recv)
     !$acc declare create(c_divs_buff_send, c_divs_buff_recv)
-#endif
+
     !> @name Generic flags used to identify and report MPI errors
     !> @{
     integer, private :: err_code, ierr, v_size
diff --git a/src/simulation/m_qbmm.fpp b/src/simulation/m_qbmm.fpp
index c478495a56..62af480e00 100644
--- a/src/simulation/m_qbmm.fpp
+++ b/src/simulation/m_qbmm.fpp
@@ -28,13 +28,9 @@ module m_qbmm
 
     private; public :: s_initialize_qbmm_module, s_mom_inv, s_coeff, s_compute_qbmm_rhs
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :, :), momrhs)
-    !$acc declare link(momrhs)
-#else
     real(kind(0d0)), allocatable, dimension(:, :, :, :, :) :: momrhs
     !$acc declare create(momrhs)
-#endif
+
     #:if MFC_CASE_OPTIMIZATION
         integer, parameter :: nterms = ${nterms}$
     #:else
@@ -43,17 +39,11 @@ module m_qbmm
     #:endif
 
     type(int_bounds_info) :: is1_qbmm, is2_qbmm, is3_qbmm
-!$acc declare create(is1_qbmm, is2_qbmm, is3_qbmm)
+    !$acc declare create(is1_qbmm, is2_qbmm, is3_qbmm)
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(integer, dimension(:), bubrs)
-    @:CRAY_DECLARE_GLOBAL(integer, dimension(:, :), bubmoms)
-    !$acc declare link(bubrs, bubmoms)
-#else
     integer, allocatable, dimension(:) :: bubrs
     integer, allocatable, dimension(:, :) :: bubmoms
     !$acc declare create(bubrs, bubmoms)
-#endif
 
 contains
 
@@ -694,7 +684,7 @@ contains
 !Coefficient array for non-polytropic model (pb and mv values are accounted in wght_pb and wght_mv)
 
     subroutine s_coeff_nonpoly(pres, rho, c, coeffs)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_coeff_nonpoly
 #else
         !$acc routine seq
@@ -767,7 +757,7 @@ contains
 
 !Coefficient array for polytropic model (pb for each R0 bin accounted for in wght_pb)
     subroutine s_coeff(pres, rho, c, coeffs)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_coeff
 #else
         !$acc routine seq
@@ -1046,7 +1036,7 @@ contains
     end subroutine s_mom_inv
 
     subroutine s_chyqmom(momin, wght, abscX, abscY)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_chyqmom
 #else
         !$acc routine seq
@@ -1113,7 +1103,7 @@ contains
     end subroutine s_chyqmom
 
     subroutine s_hyqmom(frho, fup, fmom)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_hyqmom
 #else
         !$acc routine seq
diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp
index 0846498375..a624425942 100644
--- a/src/simulation/m_rhs.fpp
+++ b/src/simulation/m_rhs.fpp
@@ -84,39 +84,23 @@ module m_rhs
     !! of the divergence theorem on the integral-average cell-boundary values
     !! of the primitive variables, located in qK_prim_n, where K = L or R.
     !> @{
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(type(vector_field), dimension(:), dq_prim_dx_qp, dq_prim_dy_qp, dq_prim_dz_qp)
-    !$acc declare link(dq_prim_dx_qp, dq_prim_dy_qp, dq_prim_dz_qp)
-#else
     type(vector_field), allocatable, dimension(:) :: dq_prim_dx_qp, dq_prim_dy_qp, dq_prim_dz_qp
     !$acc declare create(dq_prim_dx_qp, dq_prim_dy_qp, dq_prim_dz_qp)
-#endif
+    !> @}
 
     !> @name The left and right WENO-reconstructed cell-boundary values of the cell-
     !! average first-order spatial derivatives of the primitive variables. The
     !! cell-average of the first-order spatial derivatives may be found in the
     !! variables dq_prim_ds_qp, where s = x, y or z.
     !> @{
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(type(vector_field), dimension(:), dqL_prim_dx_n, dqL_prim_dy_n, dqL_prim_dz_n)
-    @:CRAY_DECLARE_GLOBAL(type(vector_field), dimension(:), dqR_prim_dx_n, dqR_prim_dy_n, dqR_prim_dz_n)
-    !$acc declare link(dqL_prim_dx_n, dqL_prim_dy_n, dqL_prim_dz_n)
-    !$acc declare link(dqR_prim_dx_n, dqR_prim_dy_n, dqR_prim_dz_n)
-#else
     type(vector_field), allocatable, dimension(:) :: dqL_prim_dx_n, dqL_prim_dy_n, dqL_prim_dz_n
     type(vector_field), allocatable, dimension(:) :: dqR_prim_dx_n, dqR_prim_dy_n, dqR_prim_dz_n
     !$acc declare create(dqL_prim_dx_n, dqL_prim_dy_n, dqL_prim_dz_n)
     !$acc declare create(dqR_prim_dx_n, dqR_prim_dy_n, dqR_prim_dz_n)
-#endif
     !> @}
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(type(scalar_field), dimension(:), tau_Re_vf)
-    !$acc declare link(tau_Re_vf)
-#else
     type(scalar_field), allocatable, dimension(:) :: tau_Re_vf
     !$acc declare create(tau_Re_vf)
-#endif
 
     type(vector_field) :: gm_alpha_qp  !<
     !! The gradient magnitude of the volume fractions at cell-interior Gaussian
@@ -128,41 +112,23 @@ module m_rhs
     !> @name The left and right WENO-reconstructed cell-boundary values of the cell-
     !! average gradient magnitude of volume fractions, located in gm_alpha_qp.
     !> @{
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(type(vector_field), dimension(:), gm_alphaL_n)
-    @:CRAY_DECLARE_GLOBAL(type(vector_field), dimension(:), gm_alphaR_n)
-    !$acc declare link(gm_alphaL_n, gm_alphaR_n)
-#else
     type(vector_field), allocatable, dimension(:) :: gm_alphaL_n
     type(vector_field), allocatable, dimension(:) :: gm_alphaR_n
     !$acc declare create(gm_alphaL_n, gm_alphaR_n)
-#endif
     !> @}
 
     !> @name The cell-boundary values of the fluxes (src - source, gsrc - geometrical
     !! source). These are computed by applying the chosen Riemann problem solver
     !! .on the left and right cell-boundary values of the primitive variables
     !> @{
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(type(vector_field), dimension(:), flux_n)
-    @:CRAY_DECLARE_GLOBAL(type(vector_field), dimension(:), flux_src_n)
-    @:CRAY_DECLARE_GLOBAL(type(vector_field), dimension(:), flux_gsrc_n)
-    !$acc declare link(flux_n, flux_src_n, flux_gsrc_n)
-#else
     type(vector_field), allocatable, dimension(:) :: flux_n
     type(vector_field), allocatable, dimension(:) :: flux_src_n
     type(vector_field), allocatable, dimension(:) :: flux_gsrc_n
     !$acc declare create(flux_n, flux_src_n, flux_gsrc_n)
-#endif
     !> @}
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(type(vector_field), dimension(:), qL_prim, qR_prim)
-    !$acc declare link(qL_prim, qR_prim)
-#else
     type(vector_field), allocatable, dimension(:) :: qL_prim, qR_prim
     !$acc declare create(qL_prim, qR_prim)
-#endif
 
     type(int_bounds_info) :: iv !< Vector field indical bounds
     !$acc declare create(iv)
@@ -179,49 +145,23 @@ module m_rhs
     !> @{
     type(scalar_field) :: alf_sum
     !> @}
-!$acc declare create(alf_sum)
-
-#ifdef CRAY_ACC_WAR
+    !$acc declare create(alf_sum)
 
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), blkmod1, blkmod2, alpha1, alpha2, Kterm)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, qR_rsx_vf, qR_rsy_vf, qR_rsz_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), dqL_rsx_vf, dqL_rsy_vf, dqL_rsz_vf, dqR_rsx_vf, dqR_rsy_vf, dqR_rsz_vf)
-    !$acc declare link(blkmod1, blkmod2, alpha1, alpha2, Kterm)
-    !$acc declare link(qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, qR_rsx_vf, qR_rsy_vf, qR_rsz_vf)
-    !$acc declare link(dqL_rsx_vf, dqL_rsy_vf, dqL_rsz_vf, dqR_rsx_vf, dqR_rsy_vf, dqR_rsz_vf)
-
-#else
     real(kind(0d0)), allocatable, dimension(:, :, :) :: blkmod1, blkmod2, alpha1, alpha2, Kterm
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, qR_rsx_vf, qR_rsy_vf, qR_rsz_vf
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: dqL_rsx_vf, dqL_rsy_vf, dqL_rsz_vf, dqR_rsx_vf, dqR_rsy_vf, dqR_rsz_vf
     !$acc declare create(blkmod1, blkmod2, alpha1, alpha2, Kterm)
     !$acc declare create(qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, qR_rsx_vf, qR_rsy_vf, qR_rsz_vf)
     !$acc declare create(dqL_rsx_vf, dqL_rsy_vf, dqL_rsz_vf, dqR_rsx_vf, dqR_rsy_vf, dqR_rsz_vf)
-#endif
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), gamma_min, pres_inf)
-    !$acc declare link(gamma_min, pres_inf)
-#else
     real(kind(0d0)), allocatable, dimension(:) :: gamma_min, pres_inf
     !$acc declare create(gamma_min, pres_inf)
-#endif
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :), Res)
-    !$acc declare link(Res)
-#else
     real(kind(0d0)), allocatable, dimension(:, :) :: Res
     !$acc declare create(Res)
-#endif
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), nbub)
-    !$acc declare link(nbub)
-#else
     real(kind(0d0)), allocatable, dimension(:, :, :) :: nbub !< Bubble number density
     !$acc declare create(nbub)
-#endif
 
 contains
 
@@ -649,27 +589,6 @@ contains
             !$acc update device(Res, Re_idx, Re_size)
         end if
 
-        ! Associating procedural pointer to the subroutine that will be
-        ! utilized to calculate the solution of a given Riemann problem
-        if (riemann_solver == 1) then
-            s_riemann_solver => s_hll_riemann_solver
-        elseif (riemann_solver == 2) then
-            s_riemann_solver => s_hllc_riemann_solver
-        end if
-
-        ! Associating the procedural pointer to the appropriate subroutine
-        ! that will be utilized in the conversion to the mixture variables
-        if (model_eqns == 1) then        ! Gamma/pi_inf model
-            s_convert_to_mixture_variables => &
-                s_convert_mixture_to_mixture_variables
-        else if (bubbles) then          ! Volume fraction for bubbles
-            s_convert_to_mixture_variables => &
-                s_convert_species_to_mixture_variables_bubbles
-        else                            ! Volume fraction model
-            s_convert_to_mixture_variables => &
-                s_convert_species_to_mixture_variables
-        end if
-
         !$acc parallel loop collapse(4) gang vector default(present)
         do id = 1, num_dims
             do i = 1, sys_size
@@ -2163,9 +2082,7 @@ contains
 
         end if
 
-#ifndef _CRAYFTN
-!$acc update device(is1, is2, is3, iv)
-#endif
+        !$acc update device(is1, is2, is3, iv)
 
         if (recon_dir == 1) then
             !$acc parallel loop collapse(4) default(present)
@@ -2361,9 +2278,6 @@ contains
             @:DEALLOCATE_GLOBAL(tau_re_vf)
         end if
 
-        s_riemann_solver => null()
-        s_convert_to_mixture_variables => null()
-
     end subroutine s_finalize_rhs_module
 
 end module m_rhs
diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp
index 99da073bba..58313ff597 100644
--- a/src/simulation/m_riemann_solvers.fpp
+++ b/src/simulation/m_riemann_solvers.fpp
@@ -54,197 +54,44 @@ module m_riemann_solvers
  s_hllc_riemann_solver, &
  s_finalize_riemann_solvers_module
 
-    abstract interface ! =======================================================
-
-        !> Abstract interface to the subroutines that are utilized to compute the
-        !! Riemann problem solution. For additional information please reference:
-        !!                        1) s_hll_riemann_solver
-        !!                        2) s_hllc_riemann_solver
-        !!                        3) s_exact_riemann_solver
-        !!  @param qL_prim_vf The  left WENO-reconstructed cell-boundary values of the
-        !!      cell-average primitive variables
-        !!  @param qR_prim_vf The right WENO-reconstructed cell-boundary values of the
-        !!      cell-average primitive variables
-        !!  @param dqL_prim_dx_vf The  left WENO-reconstructed cell-boundary values of the
-        !!      first-order x-dir spatial derivatives
-        !!  @param dqL_prim_dy_vf The  left WENO-reconstructed cell-boundary values of the
-        !!      first-order y-dir spatial derivatives
-        !!  @param dqL_prim_dz_vf The  left WENO-reconstructed cell-boundary values of the
-        !!      first-order z-dir spatial derivatives
-        !!  @param dqR_prim_dx_vf The right WENO-reconstructed cell-boundary values of the
-        !!      first-order x-dir spatial derivatives
-        !!  @param dqR_prim_dy_vf The right WENO-reconstructed cell-boundary values of the
-        !!      first-order y-dir spatial derivatives
-        !!  @param dqR_prim_dz_vf The right WENO-reconstructed cell-boundary values of the
-        !!      first-order z-dir spatial derivatives
-        !!  @param gm_alphaL_vf  Left averaged gradient magnitude
-        !!  @param gm_alphaR_vf Right averaged gradient magnitude
-        !!  @param flux_vf Intra-cell fluxes
-        !!  @param flux_src_vf Intra-cell fluxes sources
-        !!  @param flux_gsrc_vf Intra-cell geometric fluxes sources
-        !!  @param norm_dir Dir. splitting direction
-        !!  @param ix Index bounds in the x-dir
-        !!  @param iy Index bounds in the y-dir
-        !!  @param iz Index bounds in the z-dir
-        !!  @param q_prim_vf Cell-averaged primitive variables
-        subroutine s_abstract_riemann_solver(qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, dqL_prim_dx_vf, &
-                                             dqL_prim_dy_vf, &
-                                             dqL_prim_dz_vf, &
-                                             qL_prim_vf, &
-                                             qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf, dqR_prim_dx_vf, &
-                                             dqR_prim_dy_vf, &
-                                             dqR_prim_dz_vf, &
-                                             qR_prim_vf, &
-                                             q_prim_vf, &
-                                             flux_vf, flux_src_vf, &
-                                             flux_gsrc_vf, &
-                                             norm_dir, ix, iy, iz)
-
-            import :: scalar_field, int_bounds_info, sys_size, startx, starty, startz
-
-            real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
-            type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
-
-            type(scalar_field), allocatable, dimension(:), intent(inout) :: qL_prim_vf, qR_prim_vf
-
-            type(scalar_field), &
-                allocatable, dimension(:), &
-                intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, &
-                                 dqL_prim_dy_vf, dqR_prim_dy_vf, &
-                                 dqL_prim_dz_vf, dqR_prim_dz_vf
-
-            type(scalar_field), &
-                dimension(sys_size), &
-                intent(inout) :: flux_vf, flux_src_vf, flux_gsrc_vf
-
-            integer, intent(in) :: norm_dir
-
-            type(int_bounds_info), intent(in) :: ix, iy, iz
-
-        end subroutine s_abstract_riemann_solver
-
-        !> The abstract interface to the subroutines that are utilized to compute
-        !! the viscous source fluxes for either Cartesian or cylindrical geometries.
-        !! For more information please refer to:
-        !!      1) s_compute_cartesian_viscous_source_flux
-        !!      2) s_compute_cylindrical_viscous_source_flux
-        subroutine s_compute_abstract_viscous_source_flux(velL_vf, &
-                                                          dvelL_dx_vf, &
-                                                          dvelL_dy_vf, &
-                                                          dvelL_dz_vf, &
-                                                          velR_vf, &
-                                                          dvelR_dx_vf, &
-                                                          dvelR_dy_vf, &
-                                                          dvelR_dz_vf, &
-                                                          flux_src_vf, &
-                                                          norm_dir, &
-                                                          ix, iy, iz)
-
-            import :: scalar_field, int_bounds_info, num_dims, sys_size
-
-            type(scalar_field), &
-                dimension(num_dims), &
-                intent(in) :: velL_vf, velR_vf, &
-                              dvelL_dx_vf, dvelR_dx_vf, &
-                              dvelL_dy_vf, dvelR_dy_vf, &
-                              dvelL_dz_vf, dvelR_dz_vf
-
-            type(scalar_field), &
-                dimension(sys_size), &
-                intent(inout) :: flux_src_vf
-
-            integer, intent(in) :: norm_dir
-
-            type(int_bounds_info), intent(in) :: ix, iy, iz
-
-        end subroutine s_compute_abstract_viscous_source_flux
-
-    end interface ! ============================================================
-
     !> The cell-boundary values of the fluxes (src - source) that are computed
     !! through the chosen Riemann problem solver, and the direct evaluation of
     !! source terms, by using the left and right states given in qK_prim_rs_vf,
     !! dqK_prim_ds_vf where ds = dx, dy or dz.
     !> @{
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), flux_rsx_vf, flux_src_rsx_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), flux_rsy_vf, flux_src_rsy_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), flux_rsz_vf, flux_src_rsz_vf)
-    !$acc declare link( flux_rsx_vf, flux_src_rsx_vf, flux_rsy_vf,  &
-    !$acc   flux_src_rsy_vf, flux_rsz_vf, flux_src_rsz_vf )
-#else
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_rsx_vf, flux_src_rsx_vf
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_rsy_vf, flux_src_rsy_vf
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_rsz_vf, flux_src_rsz_vf
     !$acc declare create( flux_rsx_vf, flux_src_rsx_vf, flux_rsy_vf,  &
     !$acc   flux_src_rsy_vf, flux_rsz_vf, flux_src_rsz_vf )
-#endif
     !> @}
 
     !> The cell-boundary values of the geometrical source flux that are computed
     !! through the chosen Riemann problem solver by using the left and right
     !! states given in qK_prim_rs_vf. Currently 2D axisymmetric for inviscid only.
     !> @{
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), flux_gsrc_rsx_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), flux_gsrc_rsy_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), flux_gsrc_rsz_vf)
-    !$acc declare link( flux_gsrc_rsx_vf, flux_gsrc_rsy_vf, flux_gsrc_rsz_vf )
-#else
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_gsrc_rsx_vf !<
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_gsrc_rsy_vf !<
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: flux_gsrc_rsz_vf !<
     !$acc declare create( flux_gsrc_rsx_vf, flux_gsrc_rsy_vf, flux_gsrc_rsz_vf )
-#endif
     !> @}
 
     ! The cell-boundary values of the velocity. vel_src_rs_vf is determined as
     ! part of Riemann problem solution and is used to evaluate the source flux.
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), vel_src_rsx_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), vel_src_rsy_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), vel_src_rsz_vf)
-    !$acc declare link(vel_src_rsx_vf, vel_src_rsy_vf, vel_src_rsz_vf)
-#else
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: vel_src_rsx_vf
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: vel_src_rsy_vf
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: vel_src_rsz_vf
     !$acc declare create(vel_src_rsx_vf, vel_src_rsy_vf, vel_src_rsz_vf)
-#endif
-
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), mom_sp_rsx_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), mom_sp_rsy_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), mom_sp_rsz_vf)
-    !$acc declare link(mom_sp_rsx_vf, mom_sp_rsy_vf, mom_sp_rsz_vf)
-#else
+
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: mom_sp_rsx_vf
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: mom_sp_rsy_vf
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: mom_sp_rsz_vf
     !$acc declare create(mom_sp_rsx_vf, mom_sp_rsy_vf, mom_sp_rsz_vf)
-#endif
-
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), Re_avg_rsx_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), Re_avg_rsy_vf)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), Re_avg_rsz_vf)
-    !$acc declare link(Re_avg_rsx_vf, Re_avg_rsy_vf, Re_avg_rsz_vf)
-#else
+
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: Re_avg_rsx_vf
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: Re_avg_rsy_vf
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: Re_avg_rsz_vf
-    !$acc declare link(Re_avg_rsx_vf, Re_avg_rsy_vf, Re_avg_rsz_vf)
-#endif
-
-    procedure(s_abstract_riemann_solver), &
-        pointer :: s_riemann_solver => null() !<
-    !! Pointer to the procedure that is utilized to calculate either the HLL,
-    !! HLLC or exact intercell fluxes, based on the choice of Riemann solver
-
-    procedure(s_compute_abstract_viscous_source_flux), &
-        pointer :: s_compute_viscous_source_flux => null() !<
-    !! Pointer to the subroutine that is utilized to compute the viscous source
-    !! flux for either Cartesian or cylindrical geometries.
+    !$acc declare create(Re_avg_rsx_vf, Re_avg_rsy_vf, Re_avg_rsz_vf)
 
     !> @name Indical bounds in the s1-, s2- and s3-directions
     !> @{
@@ -252,26 +99,168 @@ module m_riemann_solvers
     type(int_bounds_info) :: isx, isy, isz
     !> @}
 
-!$acc declare create(is1, is2, is3, isx, isy, isz)
+    !$acc declare create(is1, is2, is3, isx, isy, isz)
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:),  Gs)
-    !$acc declare link(Gs)
-#else
     real(kind(0d0)), allocatable, dimension(:) :: Gs
     !$acc declare create(Gs)
-#endif
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :), Res)
-    !$acc declare link(Res)
-#else
     real(kind(0d0)), allocatable, dimension(:, :) :: Res
     !$acc declare create(Res)
-#endif
 
 contains
 
+    !> Dispatch to the subroutines that are utilized to compute the
+        !! Riemann problem solution. For additional information please reference:
+        !!                        1) s_hll_riemann_solver
+        !!                        2) s_hllc_riemann_solver
+        !!                        3) s_exact_riemann_solver
+        !!  @param qL_prim_vf The  left WENO-reconstructed cell-boundary values of the
+        !!      cell-average primitive variables
+        !!  @param qR_prim_vf The right WENO-reconstructed cell-boundary values of the
+        !!      cell-average primitive variables
+        !!  @param dqL_prim_dx_vf The  left WENO-reconstructed cell-boundary values of the
+        !!      first-order x-dir spatial derivatives
+        !!  @param dqL_prim_dy_vf The  left WENO-reconstructed cell-boundary values of the
+        !!      first-order y-dir spatial derivatives
+        !!  @param dqL_prim_dz_vf The  left WENO-reconstructed cell-boundary values of the
+        !!      first-order z-dir spatial derivatives
+        !!  @param dqR_prim_dx_vf The right WENO-reconstructed cell-boundary values of the
+        !!      first-order x-dir spatial derivatives
+        !!  @param dqR_prim_dy_vf The right WENO-reconstructed cell-boundary values of the
+        !!      first-order y-dir spatial derivatives
+        !!  @param dqR_prim_dz_vf The right WENO-reconstructed cell-boundary values of the
+        !!      first-order z-dir spatial derivatives
+        !!  @param gm_alphaL_vf  Left averaged gradient magnitude
+        !!  @param gm_alphaR_vf Right averaged gradient magnitude
+        !!  @param flux_vf Intra-cell fluxes
+        !!  @param flux_src_vf Intra-cell fluxes sources
+        !!  @param flux_gsrc_vf Intra-cell geometric fluxes sources
+        !!  @param norm_dir Dir. splitting direction
+        !!  @param ix Index bounds in the x-dir
+        !!  @param iy Index bounds in the y-dir
+        !!  @param iz Index bounds in the z-dir
+        !!  @param q_prim_vf Cell-averaged primitive variables
+    subroutine s_riemann_solver(qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, dqL_prim_dx_vf, &
+                                dqL_prim_dy_vf, &
+                                dqL_prim_dz_vf, &
+                                qL_prim_vf, &
+                                qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf, dqR_prim_dx_vf, &
+                                dqR_prim_dy_vf, &
+                                dqR_prim_dz_vf, &
+                                qR_prim_vf, &
+                                q_prim_vf, &
+                                flux_vf, flux_src_vf, &
+                                flux_gsrc_vf, &
+                                norm_dir, ix, iy, iz)
+
+        real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(INOUT) :: qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf
+        type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
+
+        type(scalar_field), allocatable, dimension(:), intent(INOUT) :: qL_prim_vf, qR_prim_vf
+
+        type(scalar_field), &
+            allocatable, dimension(:), &
+            intent(INOUT) :: dqL_prim_dx_vf, dqR_prim_dx_vf, &
+                             dqL_prim_dy_vf, dqR_prim_dy_vf, &
+                             dqL_prim_dz_vf, dqR_prim_dz_vf
+
+        type(scalar_field), &
+            dimension(sys_size), &
+            intent(INOUT) :: flux_vf, flux_src_vf, flux_gsrc_vf
+
+        integer, intent(IN) :: norm_dir
+
+        type(int_bounds_info), intent(IN) :: ix, iy, iz
+
+        if (riemann_solver == 1) then
+            call s_hll_riemann_solver(qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, dqL_prim_dx_vf, &
+                                      dqL_prim_dy_vf, &
+                                      dqL_prim_dz_vf, &
+                                      qL_prim_vf, &
+                                      qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf, dqR_prim_dx_vf, &
+                                      dqR_prim_dy_vf, &
+                                      dqR_prim_dz_vf, &
+                                      qR_prim_vf, &
+                                      q_prim_vf, &
+                                      flux_vf, flux_src_vf, &
+                                      flux_gsrc_vf, &
+                                      norm_dir, ix, iy, iz)
+        elseif (riemann_solver == 2) then
+            call s_hllc_riemann_solver(qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, dqL_prim_dx_vf, &
+                                       dqL_prim_dy_vf, &
+                                       dqL_prim_dz_vf, &
+                                       qL_prim_vf, &
+                                       qR_prim_rsx_vf, qR_prim_rsy_vf, qR_prim_rsz_vf, dqR_prim_dx_vf, &
+                                       dqR_prim_dy_vf, &
+                                       dqR_prim_dz_vf, &
+                                       qR_prim_vf, &
+                                       q_prim_vf, &
+                                       flux_vf, flux_src_vf, &
+                                       flux_gsrc_vf, &
+                                       norm_dir, ix, iy, iz)
+        end if
+
+    end subroutine s_riemann_solver
+
+    !> Dispatch to the subroutines that are utilized to compute
+        !! the viscous source fluxes for either Cartesian or cylindrical geometries.
+        !! For more information please refer to:
+        !!      1) s_compute_cartesian_viscous_source_flux
+        !!      2) s_compute_cylindrical_viscous_source_flux
+    subroutine s_compute_viscous_source_flux(velL_vf, & ! -------------
+                                             dvelL_dx_vf, &
+                                             dvelL_dy_vf, &
+                                             dvelL_dz_vf, &
+                                             velR_vf, &
+                                             dvelR_dx_vf, &
+                                             dvelR_dy_vf, &
+                                             dvelR_dz_vf, &
+                                             flux_src_vf, &
+                                             norm_dir, &
+                                             ix, iy, iz)
+
+        type(scalar_field), &
+            dimension(num_dims), &
+            intent(IN) :: velL_vf, velR_vf, &
+                          dvelL_dx_vf, dvelR_dx_vf, &
+                          dvelL_dy_vf, dvelR_dy_vf, &
+                          dvelL_dz_vf, dvelR_dz_vf
+
+        type(scalar_field), &
+            dimension(sys_size), &
+            intent(INOUT) :: flux_src_vf
+
+        integer, intent(IN) :: norm_dir
+
+        type(int_bounds_info), intent(IN) :: ix, iy, iz
+
+        if (grid_geometry == 3) then
+            call s_compute_cylindrical_viscous_source_flux(velL_vf, & ! -------------
+                                                           dvelL_dx_vf, &
+                                                           dvelL_dy_vf, &
+                                                           dvelL_dz_vf, &
+                                                           velR_vf, &
+                                                           dvelR_dx_vf, &
+                                                           dvelR_dy_vf, &
+                                                           dvelR_dz_vf, &
+                                                           flux_src_vf, &
+                                                           norm_dir, &
+                                                           ix, iy, iz)
+        else
+            call s_compute_cartesian_viscous_source_flux(velL_vf, & ! -------------
+                                                         dvelL_dx_vf, &
+                                                         dvelL_dy_vf, &
+                                                         dvelL_dz_vf, &
+                                                         velR_vf, &
+                                                         dvelR_dx_vf, &
+                                                         dvelR_dy_vf, &
+                                                         dvelR_dz_vf, &
+                                                         flux_src_vf, &
+                                                         norm_dir, &
+                                                         ix, iy, iz)
+        end if
+    end subroutine s_compute_viscous_source_flux
+
     subroutine s_hll_riemann_solver(qL_prim_rsx_vf, qL_prim_rsy_vf, qL_prim_rsz_vf, dqL_prim_dx_vf, & ! -------
                                     dqL_prim_dy_vf, &
                                     dqL_prim_dz_vf, &
@@ -2543,22 +2532,6 @@ contains
 
         !$acc enter data copyin(is1, is2, is3, isx, isy, isz)
 
-        ! Associating procedural pointer to the subroutine that will be
-        ! utilized to calculate the solution of a given Riemann problem
-        if (riemann_solver == 1) then
-            s_riemann_solver => s_hll_riemann_solver
-        elseif (riemann_solver == 2) then
-            s_riemann_solver => s_hllc_riemann_solver
-        end if
-
-        ! Associating procedural pointer to the subroutine that will be
-        ! utilized to compute the viscous source flux
-        if (grid_geometry == 3) then
-            s_compute_viscous_source_flux => s_compute_cylindrical_viscous_source_flux
-        else
-            s_compute_viscous_source_flux => s_compute_cartesian_viscous_source_flux
-        end if
-
         is1%beg = -1; is2%beg = 0; is3%beg = 0
         is1%end = m; is2%end = n; is3%end = p
 
@@ -4359,18 +4332,6 @@ contains
     !> Module deallocation and/or disassociation procedures
     subroutine s_finalize_riemann_solvers_module
 
-        ! Disassociating procedural pointer to the subroutine which was
-        ! utilized to calculate the solution of a given Riemann problem
-        s_riemann_solver => null()
-
-        ! Disassociating procedural pointer to the subroutine which was
-        ! utilized to calculate the viscous source flux
-        s_compute_viscous_source_flux => null()
-
-        ! Disassociating the pointer to the procedure that was utilized to
-        ! to convert mixture or species variables to the mixture variables
-        ! s_convert_to_mixture_variables => null()
-
         if (viscous) then
             @:DEALLOCATE_GLOBAL(Re_avg_rsx_vf)
         end if
diff --git a/src/simulation/m_sim_helpers.f90 b/src/simulation/m_sim_helpers.f90
index a4af912815..f6b1a3e381 100644
--- a/src/simulation/m_sim_helpers.f90
+++ b/src/simulation/m_sim_helpers.f90
@@ -31,7 +31,7 @@ module m_sim_helpers
         !! @param k y index
         !! @param l z index
     subroutine s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l)
-#ifdef CRAY_ACC_WAR
+#ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_compute_enthalpy
 #else
         !$acc routine seq
diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp
index 633f0e8ce0..6e8fbe2bfe 100644
--- a/src/simulation/m_start_up.fpp
+++ b/src/simulation/m_start_up.fpp
@@ -94,28 +94,28 @@ module m_start_up
  s_perform_time_step, s_save_data, &
  s_save_performance_metrics
 
-    abstract interface ! ===================================================
 
-        !! @param q_cons_vf  Conservative variables
-        subroutine s_read_abstract_data_files(q_cons_vf)
-
-            import :: scalar_field, sys_size, pres_field
-
-            type(scalar_field), &
-                dimension(sys_size), &
-                intent(inout) :: q_cons_vf
+    type(scalar_field), allocatable, dimension(:) :: grad_x_vf, grad_y_vf, grad_z_vf, norm_vf
 
-        end subroutine s_read_abstract_data_files
+    real(kind(0d0)) :: dt_init
 
-    end interface ! ========================================================
+contains
 
-    type(scalar_field), allocatable, dimension(:) :: grad_x_vf, grad_y_vf, grad_z_vf, norm_vf
+   !> Read data files. Dispatch subroutine that replaces procedure pointer.
+        !! @param q_cons_vf Conservative variables
+    subroutine s_read_data_files(q_cons_vf)
 
-    procedure(s_read_abstract_data_files), pointer :: s_read_data_files => null()
+        type(scalar_field), &
+            dimension(sys_size), &
+            intent(inout) :: q_cons_vf
 
-    real(kind(0d0)) :: dt_init
+        if (.not. parallel_io) then
+            call s_read_serial_data_files(q_cons_vf)
+        else
+            call s_read_parallel_data_files(q_cons_vf)
+        end if
 
-contains
+    end subroutine s_read_data_files
 
     !>  The purpose of this procedure is to first verify that an
         !!      input file has been made available by the user. Provided
@@ -1350,15 +1350,6 @@ contains
         call acc_present_dump()
 #endif
 
-        ! Associate pointers for serial or parallel I/O
-        if (parallel_io .neqv. .true.) then
-            s_read_data_files => s_read_serial_data_files
-            s_write_data_files => s_write_serial_data_files
-        else
-            s_read_data_files => s_read_parallel_data_files
-            s_write_data_files => s_write_parallel_data_files
-        end if
-
         ! Reading in the user provided initial condition and grid data
         call s_read_data_files(q_cons_ts(1)%vf)
 
@@ -1492,9 +1483,6 @@ contains
     end subroutine s_initialize_gpu_vars
 
     subroutine s_finalize_modules
-        ! Disassociate pointers for serial and parallel I/O
-        s_read_data_files => null()
-        s_write_data_files => null()
 
         call s_finalize_time_steppers_module()
         call s_finalize_derived_variables_module()
diff --git a/src/simulation/m_surface_tension.fpp b/src/simulation/m_surface_tension.fpp
index 587539b6b2..926d90f833 100644
--- a/src/simulation/m_surface_tension.fpp
+++ b/src/simulation/m_surface_tension.fpp
@@ -28,27 +28,17 @@ module m_surface_tension
  s_get_capilary, &
  s_finalize_surface_tension_module
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(type(scalar_field), dimension(:), c_divs)
-    !$acc declare link(c_divs)
-#else
     !> @name color function gradient components and magnitude
     !> @{
     type(scalar_field), allocatable, dimension(:) :: c_divs
     !> @)
     !$acc declare create(c_divs)
-#endif
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:,:,:,:), gL_x, gL_y, gL_z, gR_x, gR_y, gR_z)
-    !$acc declare link(gL_x, gL_y, gL_z, gR_x, gR_y, gR_z)
-#else
     !> @name cell boundary reconstructed gradient components and magnitude
     !> @{
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: gL_x, gR_x, gL_y, gR_y, gL_z, gR_z
     !> @}
     !$acc declare create(gL_x, gR_x, gL_y, gR_y, gL_z, gR_z)
-#endif
 
     type(int_bounds_info) :: is1, is2, is3, iv
     !$acc declare create(is1, is2, is3, iv)
diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp
index e9b2603708..28fa61e9d5 100644
--- a/src/simulation/m_time_steppers.fpp
+++ b/src/simulation/m_time_steppers.fpp
@@ -44,30 +44,6 @@ module m_time_steppers
 
     implicit none
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(type(vector_field), dimension(:), q_cons_ts)
-    !! Cell-average conservative variables at each time-stage (TS)
-
-    @:CRAY_DECLARE_GLOBAL(type(scalar_field), dimension(:), q_prim_vf)
-    !! Cell-average primitive variables at the current time-stage
-
-    @:CRAY_DECLARE_GLOBAL(type(scalar_field), dimension(:), rhs_vf)
-    !! Cell-average RHS variables at the current time-stage
-
-    @:CRAY_DECLARE_GLOBAL(type(vector_field), dimension(:), q_prim_ts)
-    !! Cell-average primitive variables at consecutive TIMESTEPS
-
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :, :), rhs_pb)
-
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :, :), rhs_mv)
-
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension( :, :, :), max_dt)
-
-    integer, private :: num_ts !<
-    !! Number of time stages in the time-stepping scheme
-
-    !$acc declare link(q_cons_ts,q_prim_vf,rhs_vf,q_prim_ts, rhs_mv, rhs_pb, max_dt)
-#else
     type(vector_field), allocatable, dimension(:) :: q_cons_ts !<
     !! Cell-average conservative variables at each time-stage (TS)
 
@@ -90,7 +66,6 @@ module m_time_steppers
     !! Number of time stages in the time-stepping scheme
 
     !$acc declare create(q_cons_ts,q_prim_vf,rhs_vf,q_prim_ts, rhs_mv, rhs_pb, max_dt)
-#endif
 
 contains
 
diff --git a/src/simulation/m_viscous.fpp b/src/simulation/m_viscous.fpp
index 3e2c3bf70d..55ff395cf1 100644
--- a/src/simulation/m_viscous.fpp
+++ b/src/simulation/m_viscous.fpp
@@ -26,15 +26,10 @@ module m_viscous
 
     type(int_bounds_info) :: iv
     type(int_bounds_info) :: is1_viscous, is2_viscous, is3_viscous
-!$acc declare create(is1_viscous, is2_viscous, is3_viscous, iv)
+    !$acc declare create(is1_viscous, is2_viscous, is3_viscous, iv)
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :), Res_viscous)
-    !$acc declare link(Res_viscous)
-#else
     real(kind(0d0)), allocatable, dimension(:, :) :: Res_viscous
-    !$acc declare create(Re_viscous)
-#endif
+    !$acc declare create(Res_viscous)
 
 contains
 
diff --git a/src/simulation/m_weno.fpp b/src/simulation/m_weno.fpp
index f0b5490b41..eef47b873c 100644
--- a/src/simulation/m_weno.fpp
+++ b/src/simulation/m_weno.fpp
@@ -44,12 +44,7 @@ module m_weno
     !! stencils (WS) that are annexed to each position of a given scalar field.
     !> @{
 
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :, :), v_rs_ws_x, v_rs_ws_y, v_rs_ws_z)
-    !$acc declare link(v_rs_ws_x, v_rs_ws_y, v_rs_ws_z)
-#else
     real(kind(0d0)), allocatable, dimension(:, :, :, :) :: v_rs_ws_x, v_rs_ws_y, v_rs_ws_z
-#endif
     !> @}
 
     ! WENO Coefficients ========================================================
@@ -60,17 +55,6 @@ module m_weno
     !! second dimension identifies the position of its coefficients and the last
     !! dimension denotes the cell-location in the relevant coordinate direction.
     !> @{
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), poly_coef_cbL_x)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), poly_coef_cbL_y)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), poly_coef_cbL_z)
-
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), poly_coef_cbR_x)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), poly_coef_cbR_y)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), poly_coef_cbR_z)
-    !$acc declare link(poly_coef_cbL_x, poly_coef_cbL_y, poly_coef_cbL_z)
-    !$acc declare link(poly_coef_cbR_x, poly_coef_cbR_y, poly_coef_cbR_z)
-#else
     real(kind(0d0)), target, allocatable, dimension(:, :, :) :: poly_coef_cbL_x
     real(kind(0d0)), target, allocatable, dimension(:, :, :) :: poly_coef_cbL_y
     real(kind(0d0)), target, allocatable, dimension(:, :, :) :: poly_coef_cbL_z
@@ -78,7 +62,6 @@ module m_weno
     real(kind(0d0)), target, allocatable, dimension(:, :, :) :: poly_coef_cbR_x
     real(kind(0d0)), target, allocatable, dimension(:, :, :) :: poly_coef_cbR_y
     real(kind(0d0)), target, allocatable, dimension(:, :, :) :: poly_coef_cbR_z
-#endif
 
     !    real(kind(0d0)), pointer, dimension(:, :, :) :: poly_coef_L => null()
     !    real(kind(0d0)), pointer, dimension(:, :, :) :: poly_coef_R => null()
@@ -89,16 +72,6 @@ module m_weno
     !! that the first dimension of the array identifies the weight, while the
     !! last denotes the cell-location in the relevant coordinate direction.
     !> @{
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :), d_cbL_y)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :), d_cbL_x)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :), d_cbL_z)
-
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :), d_cbR_x)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :), d_cbR_y)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :), d_cbR_z)
-    !$acc declare link(d_cbL_x, d_cbL_y, d_cbL_z, d_cbR_x, d_cbR_y, d_cbR_z)
-#else
     real(kind(0d0)), target, allocatable, dimension(:, :) :: d_cbL_x
     real(kind(0d0)), target, allocatable, dimension(:, :) :: d_cbL_y
     real(kind(0d0)), target, allocatable, dimension(:, :) :: d_cbL_z
@@ -106,7 +79,6 @@ module m_weno
     real(kind(0d0)), target, allocatable, dimension(:, :) :: d_cbR_x
     real(kind(0d0)), target, allocatable, dimension(:, :) :: d_cbR_y
     real(kind(0d0)), target, allocatable, dimension(:, :) :: d_cbR_z
-#endif
 !    real(kind(0d0)), pointer, dimension(:, :) :: d_L => null()
 !    real(kind(0d0)), pointer, dimension(:, :) :: d_R => null()
     !> @}
@@ -116,16 +88,9 @@ module m_weno
     !! second identifies the position of its coefficients and the last denotes
     !! the cell-location in the relevant coordinate direction.
     !> @{
-#ifdef CRAY_ACC_WAR
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), beta_coef_x)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), beta_coef_y)
-    @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), beta_coef_z)
-    !$acc declare link(beta_coef_x, beta_coef_y, beta_coef_z)
-#else
     real(kind(0d0)), target, allocatable, dimension(:, :, :) :: beta_coef_x
     real(kind(0d0)), target, allocatable, dimension(:, :, :) :: beta_coef_y
     real(kind(0d0)), target, allocatable, dimension(:, :, :) :: beta_coef_z
-#endif
 !    real(kind(0d0)), pointer, dimension(:, :, :) :: beta_coef => null()
     !> @}
 
@@ -143,15 +108,13 @@ module m_weno
     !> @}
 
     real(kind(0d0)) :: test
-!$acc declare create(test)
-
-#ifndef CRAY_ACC_WAR
-!$acc declare create( &
-!$acc                v_rs_ws_x, v_rs_ws_y, v_rs_ws_z, &
-!$acc                poly_coef_cbL_x,poly_coef_cbL_y,poly_coef_cbL_z, &
-!$acc                poly_coef_cbR_x,poly_coef_cbR_y,poly_coef_cbR_z,d_cbL_x,       &
-!$acc                d_cbL_y,d_cbL_z,d_cbR_x,d_cbR_y,d_cbR_z,beta_coef_x,beta_coef_y,beta_coef_z)
-#endif
+    !$acc declare create(test)
+
+    !$acc declare create( &
+    !$acc                v_rs_ws_x, v_rs_ws_y, v_rs_ws_z, &
+    !$acc                poly_coef_cbL_x,poly_coef_cbL_y,poly_coef_cbL_z, &
+    !$acc                poly_coef_cbR_x,poly_coef_cbR_y,poly_coef_cbR_z,d_cbL_x,       &
+    !$acc                d_cbL_y,d_cbL_z,d_cbR_x,d_cbR_y,d_cbR_z,beta_coef_x,beta_coef_y,beta_coef_z)
 
 contains