diff --git a/src/common/m_checker_common.fpp b/src/common/m_checker_common.fpp
index 5ec438a7dd..1c3a238941 100644
--- a/src/common/m_checker_common.fpp
+++ b/src/common/m_checker_common.fpp
@@ -322,7 +322,7 @@ contains
         @:PROHIBIT(surface_tension .and. sigma < 0._wp, &
             "sigma must be greater than or equal to zero")
 
-        @:PROHIBIT(surface_tension .and. sigma == dflt_real, &
+        @:PROHIBIT(surface_tension .and. f_approx_equal(sigma, dflt_real), &
             "sigma must be set if surface_tension is enabled")
 
         @:PROHIBIT(.not. f_is_default(sigma) .and. .not. surface_tension, &
@@ -347,9 +347,12 @@ contains
         !! Called by s_check_inputs_common for all three stages
     impure subroutine s_check_inputs_moving_bc
         #:for X, VB2, VB3 in [('x', 'vb2', 'vb3'), ('y', 'vb3', 'vb1'), ('z', 'vb1', 'vb2')]
-            if (any((/bc_${X}$%vb1, bc_${X}$%vb2, bc_${X}$%vb3/) /= 0._wp)) then
+            if (.not. (f_approx_equal(bc_${X}$%vb1, 0._wp) .and. &
+                       f_approx_equal(bc_${X}$%vb2, 0._wp) .and. &
+                       f_approx_equal(bc_${X}$%vb3, 0._wp))) then
                 if (bc_${X}$%beg == BC_SLIP_WALL) then
-                    if (any((/bc_${X}$%${VB2}$, bc_${X}$%${VB3}$/) /= 0._wp)) then
+                    if (.not. (f_approx_equal(bc_${X}$%${VB2}$, 0._wp) .and. &
+                               f_approx_equal(bc_${X}$%${VB3}$, 0._wp))) then
                         call s_mpi_abort("bc_${X}$%beg must be -15 if "// &
                                          "bc_${X}$%${VB2}$ or bc_${X}$%${VB3}$ "// &
                                          "is set. Exiting.", CASE_FILE_ERROR_CODE)
@@ -362,9 +365,12 @@ contains
         #:endfor
 
         #:for X, VE2, VE3 in [('x', 've2', 've3'), ('y', 've3', 've1'), ('z', 've1', 've2')]
-            if (any((/bc_${X}$%ve1, bc_${X}$%ve2, bc_${X}$%ve3/) /= 0._wp)) then
+            if (.not. (f_approx_equal(bc_${X}$%ve1, 0._wp) .and. &
+                       f_approx_equal(bc_${X}$%ve2, 0._wp) .and. &
+                       f_approx_equal(bc_${X}$%ve3, 0._wp))) then
                 if (bc_${X}$%end == BC_SLIP_WALL) then
-                    if (any((/bc_${X}$%${VE2}$, bc_${X}$%${VE3}$/) /= 0._wp)) then
+                    if (.not. (f_approx_equal(bc_${X}$%${VE2}$, 0._wp) .and. &
+                               f_approx_equal(bc_${X}$%${VE3}$, 0._wp))) then
                         call s_mpi_abort("bc_${X}$%end must be -15 if "// &
                                          "bc_${X}$%${VE2}$ or bc_${X}$%${VE3}$ "// &
                                          "is set. Exiting.", CASE_FILE_ERROR_CODE)
diff --git a/src/common/m_eigen_solver.f90 b/src/common/m_eigen_solver.f90
index 292ba1f730..2bf38d04bd 100644
--- a/src/common/m_eigen_solver.f90
+++ b/src/common/m_eigen_solver.f90
@@ -10,6 +10,8 @@ module m_eigen_solver
 
     use m_precision_select
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     implicit none
 
     private; 
@@ -124,7 +126,7 @@ pure subroutine cbal(nm, nl, ar, ai, low, igh, scale)
 
             do 110 i = 1, l
                 if (i == j) go to 110
-                if (ar(j, i) /= 0.0_wp .or. ai(j, i) /= 0.0_wp) go to 120
+                if (.not. f_approx_equal(ar(j, i), 0.0_wp) .or. .not. f_approx_equal(ai(j, i), 0.0_wp)) go to 120
 110         end do
 
             ml = l
@@ -140,7 +142,7 @@ pure subroutine cbal(nm, nl, ar, ai, low, igh, scale)
 
             do 150 i = k, l
                 if (i == j) go to 150
-                if (ar(i, j) /= 0.0_wp .or. ai(i, j) /= 0.0_wp) go to 170
+                if (.not. f_approx_equal(ar(i, j), 0.0_wp) .or. .not. f_approx_equal(ai(i, j), 0.0_wp)) go to 170
 150         end do
 
             ml = k
@@ -164,7 +166,7 @@ pure subroutine cbal(nm, nl, ar, ai, low, igh, scale)
                 r = r + abs(ar(i, j)) + abs(ai(i, j))
 200         end do
 !     guard against zero c or r due to underflow
-            if (c == 0.0_wp .or. r == 0.0_wp) go to 270
+            if (f_approx_equal(c, 0.0_wp) .or. f_approx_equal(r, 0.0_wp)) go to 270
             g = r/radix
             f = 1.0_wp
             s = c + r
@@ -242,7 +244,7 @@ pure subroutine corth(nm, nl, low, igh, ar, ai, ortr, orti)
             do 90 i = ml, igh
                 scale = scale + abs(ar(i, ml - 1)) + abs(ai(i, ml - 1))
 90          end do
-            if (scale == 0._wp) go to 180
+            if (f_approx_equal(scale, 0._wp)) go to 180
             mp = ml + igh
 !     for i=igh step -1 until ml do
             do 100 ii = ml, igh
@@ -254,7 +256,7 @@ pure subroutine corth(nm, nl, low, igh, ar, ai, ortr, orti)
 
             g = sqrt(h)
             call pythag(ortr(ml), orti(ml), f)
-            if (f == 0._wp) go to 103
+            if (f_approx_equal(f, 0._wp)) go to 103
             h = h + f*g
             g = g/f
             ortr(ml) = (1.0_wp + g)*ortr(ml)
@@ -374,8 +376,8 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
 !     for i=igh-1 step -1 until low+1 do
 105     do 140 ii = 1, iend
             i = igh - ii
-            if (abs(ortr(i)) == 0._wp .and. abs(orti(i)) == 0._wp) go to 140
-            if (abs(hr(i, i - 1)) == 0._wp .and. abs(hi(i, i - 1)) == 0._wp) go to 140
+            if (f_approx_equal(abs(ortr(i)), 0._wp) .and. f_approx_equal(abs(orti(i)), 0._wp)) go to 140
+            if (f_approx_equal(abs(hr(i, i - 1)), 0._wp) .and. f_approx_equal(abs(hi(i, i - 1)), 0._wp)) go to 140
 !     norm below is negative of h formed in corth
             norm = hr(i, i - 1)*ortr(i) + hi(i, i - 1)*orti(i)
             ip1 = i + 1
@@ -410,7 +412,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
 
         do 170 i = l, igh
             ll = min0(i + 1, igh)
-            if (abs(hi(i, i - 1)) == 0._wp) go to 170
+            if (f_approx_equal(abs(hi(i, i - 1)), 0._wp)) go to 170
             call pythag(hr(i, i - 1), hi(i, i - 1), norm)
             yr = hr(i, i - 1)/norm
             yi = hi(i, i - 1)/norm
@@ -458,7 +460,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
             tst1 = abs(hr(l - 1, l - 1)) + abs(hi(l - 1, l - 1)) &
                    + abs(hr(l, l)) + abs(hi(l, l))
             tst2 = tst1 + abs(hr(l, l - 1))
-            if (tst2 == tst1) go to 300
+            if (f_approx_equal(tst2, tst1)) go to 300
 260     end do
 !     form shift
 300     if (l == en) go to 660
@@ -468,7 +470,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
         si = hi(en, en)
         xr = hr(enm1, en)*hr(en, enm1)
         xi = hi(enm1, en)*hr(en, enm1)
-        if (xr == 0.0_wp .and. xi == 0.0_wp) go to 340
+        if (f_approx_equal(xr, 0.0_wp) .and. f_approx_equal(xi, 0.0_wp)) go to 340
         yr = (hr(enm1, enm1) - sr)/2.0_wp
         yi = (hi(enm1, enm1) - si)/2.0_wp
         call csroot(yr**2 - yi**2 + xr, 2.0_wp*yr*yi + xi, zzr, zzi)
@@ -522,7 +524,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
 500     end do
 
         si = hi(en, en)
-        if (abs(si) == 0._wp) go to 540
+        if (f_approx_equal(abs(si), 0._wp)) go to 540
         call pythag(hr(en, en), si, norm)
         sr = hr(en, en)/norm
         si = si/norm
@@ -567,7 +569,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
 590         end do
 600     end do
 
-        if (abs(si) == 0._wp) go to 240
+        if (f_approx_equal(abs(si), 0._wp)) go to 240
 
         do 630 i = 1, en
             yr = hr(i, en)
@@ -602,7 +604,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
             end do
         end do
 
-        if (nl == 1 .or. norm == 0._wp) go to 1001
+        if (nl == 1 .or. f_approx_equal(norm, 0._wp)) go to 1001
 !     for en=nl step -1 until 2 do
         do 800 nn = 2, nl
             en = nl + 2 - nn
@@ -625,7 +627,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
 
                 yr = xr - wr(i)
                 yi = xi - wi(i)
-                if (yr /= 0.0_wp .or. yi /= 0.0_wp) go to 765
+                if (.not. f_approx_equal(yr, 0.0_wp) .or. .not. f_approx_equal(yi, 0.0_wp)) go to 765
                 tst1 = norm
                 yr = tst1
 760             yr = 0.01_wp*yr
@@ -635,7 +637,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
                 call cdiv(zzr, zzi, yr, yi, hr(i, en), hi(i, en))
 !     overflow control
                 tr = abs(hr(i, en)) + abs(hi(i, en))
-                if (tr == 0.0_wp) go to 780
+                if (f_approx_equal(tr, 0.0_wp)) go to 780
                 tst1 = tr
                 tst2 = tst1 + 1.0_wp/tst1
                 if (tst2 > tst1) go to 780
@@ -796,11 +798,11 @@ pure elemental subroutine pythag(a, b, c)
 
         real(wp) :: p, r, s, t, u
         p = max(abs(a), abs(b))
-        if (p == 0.0_wp) go to 20
+        if (f_approx_equal(p, 0.0_wp)) go to 20
         r = (min(abs(a), abs(b))/p)**2
 10      continue
         t = 4.0_wp + r
-        if (t == 4.0_wp) go to 20
+        if (f_approx_equal(t, 4.0_wp)) go to 20
         s = r/t
         u = 1.0_wp + 2.0_wp*s
         p = u*p
diff --git a/src/common/m_helper.fpp b/src/common/m_helper.fpp
index 4be5974c2e..835413cb80 100644
--- a/src/common/m_helper.fpp
+++ b/src/common/m_helper.fpp
@@ -163,7 +163,7 @@ contains
         phi_nv = (1._wp + sqrt(mu_n/mu_v)*(M_v/M_n)**(0.25_wp))**2 &
                  /(sqrt(8._wp)*sqrt(1._wp + M_n/M_v))
         ! internal bubble pressure
-        pb0 = pl0 + 2._wp*ss/(R0ref*R0)
+        pb0(:) = pl0 + 2._wp*ss/(R0ref*R0(:))
 
         ! mass fraction of vapor
         chi_vw0 = 1._wp/(1._wp + R_v/R_n*(pb0/pv - 1._wp))
@@ -179,10 +179,10 @@ contains
         rho_m0 = pv/(chi_vw0*R_v*temp)
 
         ! mass of gas/vapor computed using dimensional quantities
-        mass_n0 = 4._wp*(pb0 - pv)*pi/(3._wp*R_n*temp*rhol0)*R0**3
-        mass_v0 = 4._wp*pv*pi/(3._wp*R_v*temp*rhol0)*R0**3
+        mass_n0(:) = 4._wp*(pb0(:) - pv)*pi/(3._wp*R_n*temp*rhol0)*R0(:)**3
+        mass_v0(:) = 4._wp*pv*pi/(3._wp*R_v*temp*rhol0)*R0(:)**3
         ! Peclet numbers
-        Pe_T = rho_m0*cp_m0*uu*R0ref/k_m0
+        Pe_T(:) = rho_m0*cp_m0(:)*uu*R0ref/k_m0(:)
         Pe_c = uu*R0ref/D_m
 
         Tw = temp
@@ -191,8 +191,8 @@ contains
         !if(.not. qbmm) then
         R_n = rhol0*R_n*temp/pl0
         R_v = rhol0*R_v*temp/pl0
-        k_n = k_n/k_m0
-        k_v = k_v/k_m0
+        k_n(:) = k_n(:)/k_m0(:)
+        k_v(:) = k_v(:)/k_m0(:)
         pb0 = pb0/pl0
         pv = pv/pl0
         Tw = 1._wp
@@ -203,7 +203,7 @@ contains
         !end if
 
         ! natural frequencies
-        omegaN = sqrt(3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/(Web*R0))/R0
+        omegaN(:) = sqrt(3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/(Web*R0))/R0
         do ir = 1, Nb
             call s_transcoeff(omegaN(ir)*R0(ir), Pe_T(ir)*R0(ir), &
                               Re_trans_T(ir), Im_trans_T(ir))
diff --git a/src/common/m_helper_basic.f90 b/src/common/m_helper_basic.f90
index 4205279cd4..47db999d49 100644
--- a/src/common/m_helper_basic.f90
+++ b/src/common/m_helper_basic.f90
@@ -10,6 +10,7 @@ module m_helper_basic
 
     private; 
     public :: f_approx_equal, &
+              f_approx_in_array, &
               f_is_default, &
               f_all_default, &
               f_is_integer, &
@@ -20,7 +21,7 @@ module m_helper_basic
     !> This procedure checks if two floating point numbers of wp are within tolerance.
     !! @param a First number.
     !! @param b Second number.
-    !! @param tol_input Relative error (default = 1e-6_wp).
+    !! @param tol_input Relative error (default = 1e-10_wp).
     !! @return Result of the comparison.
     logical pure elemental function f_approx_equal(a, b, tol_input) result(res)
         !$acc routine seq
@@ -31,7 +32,7 @@ logical pure elemental function f_approx_equal(a, b, tol_input) result(res)
         if (present(tol_input)) then
             tol = tol_input
         else
-            tol = 1e-6_wp
+            tol = 1e-10_wp
         end if
 
         if (a == b) then
@@ -43,6 +44,35 @@ logical pure elemental function f_approx_equal(a, b, tol_input) result(res)
         end if
     end function f_approx_equal
 
+    !> This procedure checks if the point numbers of wp belongs to another array are within tolerance.
+    !! @param a First number.
+    !! @param b Array that contains several point numbers.
+    !! @param tol_input Relative error (default = 1e-10_wp).
+    !! @return Result of the comparison.
+    logical pure function f_approx_in_array(a, b, tol_input) result(res)
+        !$acc routine seq
+        real(wp), intent(in) :: a
+        real(wp), intent(in) :: b(:)
+        real(wp), optional, intent(in) :: tol_input
+        real(wp) :: tol
+        integer :: i
+
+        res = .false.
+
+        if (present(tol_input)) then
+            tol = tol_input
+        else
+            tol = 1e-10_wp
+        end if
+
+        do i = 1, size(b)
+            if (f_approx_equal(a, b(i), tol)) then
+                res = .true.
+                exit
+            end if
+        end do
+    end function f_approx_in_array
+
     !> Checks if a real(wp) variable is of default value.
     !! @param var Variable to check.
     logical pure elemental function f_is_default(var) result(res)
diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp
index 4a96787e15..b920151488 100644
--- a/src/common/m_mpi_common.fpp
+++ b/src/common/m_mpi_common.fpp
@@ -608,8 +608,6 @@ contains
 
         integer :: pack_offset, unpack_offset
 
-        real(wp), pointer :: p_send, p_recv
-
 #ifdef MFC_MPI
 
         call nvtxStartRange("RHS-COMM-PACKBUF")
diff --git a/src/common/m_phase_change.fpp b/src/common/m_phase_change.fpp
index b520d700ba..2ccfa63659 100644
--- a/src/common/m_phase_change.fpp
+++ b/src/common/m_phase_change.fpp
@@ -17,6 +17,8 @@ module m_phase_change
 
     use ieee_arithmetic
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     implicit none
 
     private; 
@@ -748,7 +750,7 @@ contains
         ! Generic loop iterators
         integer :: ns
 
-        if ((pSat == 0.0_wp) .and. (TSIn == 0.0_wp)) then
+        if ((f_approx_equal(pSat, 0.0_wp)) .and. (f_approx_equal(TSIn, 0.0_wp))) then
 
             ! assigning Saturation temperature
             TSat = 0.0_wp
diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp
index 382f8387ae..ce8eebeffc 100644
--- a/src/post_process/m_data_output.fpp
+++ b/src/post_process/m_data_output.fpp
@@ -183,28 +183,28 @@ contains
         if (format == 1 .and. n > 0) then
             if (p > 0) then
                 if (grid_geometry == 3) then
-                    lo_offset = (/offset_y%beg, offset_z%beg, offset_x%beg/)
-                    hi_offset = (/offset_y%end, offset_z%end, offset_x%end/)
+                    lo_offset(:) = (/offset_y%beg, offset_z%beg, offset_x%beg/)
+                    hi_offset(:) = (/offset_y%end, offset_z%end, offset_x%end/)
                 else
-                    lo_offset = (/offset_x%beg, offset_y%beg, offset_z%beg/)
-                    hi_offset = (/offset_x%end, offset_y%end, offset_z%end/)
+                    lo_offset(:) = (/offset_x%beg, offset_y%beg, offset_z%beg/)
+                    hi_offset(:) = (/offset_x%end, offset_y%end, offset_z%end/)
                 end if
 
                 if (grid_geometry == 3) then
-                    dims = (/n + offset_y%beg + offset_y%end + 2, &
-                             p + offset_z%beg + offset_z%end + 2, &
-                             m + offset_x%beg + offset_x%end + 2/)
+                    dims(:) = (/n + offset_y%beg + offset_y%end + 2, &
+                                p + offset_z%beg + offset_z%end + 2, &
+                                m + offset_x%beg + offset_x%end + 2/)
                 else
-                    dims = (/m + offset_x%beg + offset_x%end + 2, &
-                             n + offset_y%beg + offset_y%end + 2, &
-                             p + offset_z%beg + offset_z%end + 2/)
+                    dims(:) = (/m + offset_x%beg + offset_x%end + 2, &
+                                n + offset_y%beg + offset_y%end + 2, &
+                                p + offset_z%beg + offset_z%end + 2/)
                 end if
             else
-                lo_offset = (/offset_x%beg, offset_y%beg/)
-                hi_offset = (/offset_x%end, offset_y%end/)
+                lo_offset(:) = (/offset_x%beg, offset_y%beg/)
+                hi_offset(:) = (/offset_x%end, offset_y%end/)
 
-                dims = (/m + offset_x%beg + offset_x%end + 2, &
-                         n + offset_y%beg + offset_y%end + 2/)
+                dims(:) = (/m + offset_x%beg + offset_x%end + 2, &
+                            n + offset_y%beg + offset_y%end + 2/)
             end if
         end if
 
@@ -710,10 +710,10 @@ contains
 
             if (precision == 1) then
                 if (p > 0) then
-                    z_cb_s = real(z_cb, sp)
+                    z_cb_s(:) = real(z_cb(:), sp)
                 end if
-                x_cb_s = real(x_cb, sp)
-                y_cb_s = real(y_cb, sp)
+                x_cb_s(:) = real(x_cb(:), sp)
+                y_cb_s(:) = real(y_cb(:), sp)
             end if
 
             #:for PRECISION, SFX, DBT in [(1,'_s','DB_FLOAT'),(2,'',"DB_DOUBLE")]
@@ -805,7 +805,7 @@ contains
                 if (num_procs > 1) then
                     call s_mpi_defragment_1d_grid_variable()
                 else
-                    x_root_cb = x_cb
+                    x_root_cb(:) = x_cb(:)
                 end if
 
                 if (proc_rank == 0) then
@@ -872,11 +872,11 @@ contains
             if (n == 0) then
 
                 if (precision == 1 .and. wp == dp) then
-                    x_cc_s = real(x_cc, sp)
-                    q_sf_s = real(q_sf, sp)
+                    x_cc_s(:) = real(x_cc(:), sp)
+                    q_sf_s(:, :, :) = real(q_sf(:, :, :), sp)
                 elseif (precision == 1 .and. wp == sp) then
-                    x_cc_s = x_cc
-                    q_sf_s = q_sf
+                    x_cc_s(:) = x_cc(:)
+                    q_sf_s(:, :, :) = q_sf(:, :, :)
                 end if
 
                 ! Writing the curve object associated with the local process
@@ -897,16 +897,16 @@ contains
                     call s_mpi_defragment_1d_flow_variable(q_sf, q_root_sf)
 
                     if (precision == 1) then
-                        x_root_cc_s = real(x_root_cc, sp)
-                        q_root_sf_s = real(q_root_sf, sp)
+                        x_root_cc_s(:) = real(x_root_cc(:), sp)
+                        q_root_sf_s(:, :, :) = real(q_root_sf(:, :, :), sp)
                     end if
                 else
                     if (precision == 1) then
-                        x_root_cc_s = real(x_cc, sp)
-                        q_root_sf_s = real(q_sf, sp)
+                        x_root_cc_s(:) = real(x_cc(:), sp)
+                        q_root_sf_s(:, :, :) = real(q_sf(:, :, :), sp)
                     else
-                        x_root_cc = x_cc
-                        q_root_sf = q_sf
+                        x_root_cc(:) = x_cc(:)
+                        q_root_sf(:, :, :) = q_sf(:, :, :)
                     end if
                 end if
 
@@ -997,7 +997,7 @@ contains
                             end do
                         end if
                     end if
-                elseif (wp == dp) then
+                elseif (wp == sp) then
                     do i = -offset_x%beg, m + offset_x%end
                         do j = -offset_y%beg, n + offset_y%end
                             do k = -offset_z%beg, p + offset_z%end
@@ -1069,7 +1069,7 @@ contains
                 if (num_procs > 1) then
                     call s_mpi_defragment_1d_flow_variable(q_sf, q_root_sf)
                 else
-                    q_root_sf = q_sf
+                    q_root_sf(:, :, :) = q_sf(:, :, :)
                 end if
 
                 if (proc_rank == 0) then
@@ -1196,7 +1196,6 @@ contains
         type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf
         integer :: i, j, k, l, cent !< Generic loop iterators
         integer :: counter, root !< number of data points extracted to fit shape to SH perturbations
-        real(wp), parameter :: pi = 4._wp*tan(1._wp)
         real(wp), allocatable :: x_td(:), y_td(:), x_d1(:), y_d1(:), y_d(:), x_d(:)
         real(wp) :: axp, axm, ayp, aym, tgp, euc_d, thres, maxalph_loc, maxalph_glb
 
diff --git a/src/post_process/m_derived_variables.fpp b/src/post_process/m_derived_variables.fpp
index 0deffe61ba..6eaefe30a3 100644
--- a/src/post_process/m_derived_variables.fpp
+++ b/src/post_process/m_derived_variables.fpp
@@ -16,6 +16,8 @@ module m_derived_variables
 
     use m_mpi_proxy             !< Message passing interface (MPI) module proxy
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     use m_variables_conversion
 
     implicit none
@@ -286,7 +288,7 @@ contains
                     if (abs(top) < 1e-8_wp) top = 0._wp
                     if (abs(bottom) < 1e-8_wp) bottom = 0._wp
 
-                    if (top == bottom) then
+                    if (f_approx_equal(top, bottom)) then
                         slope = 1._wp
                         !       ELSEIF((top == 0._wp .AND. bottom /= 0._wp) &
                         !               .OR.            &
diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90
index 0b9435b0b8..b733f43dd9 100644
--- a/src/post_process/m_start_up.f90
+++ b/src/post_process/m_start_up.f90
@@ -262,7 +262,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
         if ((model_eqns == 2) .or. (model_eqns == 3) .or. (model_eqns == 4)) then
             do i = 1, num_fluids
                 if (alpha_rho_wrt(i) .or. (cons_vars_wrt .or. prim_vars_wrt)) then
-                    q_sf = q_cons_vf(i)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                    q_sf(:, :, :) = q_cons_vf(i)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
                     if (model_eqns /= 4) then
                         write (varname, '(A,I0)') 'alpha_rho', i
                     else
@@ -278,7 +278,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
 
         ! Adding the density to the formatted database file
         if ((rho_wrt .or. (model_eqns == 1 .and. (cons_vars_wrt .or. prim_vars_wrt))) .and. (.not. relativity)) then
-            q_sf = rho_sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+            q_sf(:, :, :) = rho_sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
             write (varname, '(A)') 'rho'
             call s_write_variable_to_formatted_database_file(varname, t_step)
 
@@ -286,7 +286,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
         end if
 
         if (relativity .and. (rho_wrt .or. prim_vars_wrt)) then
-            q_sf = q_prim_vf(1)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+            q_sf(:, :, :) = q_prim_vf(1)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
             write (varname, '(A)') 'rho'
             call s_write_variable_to_formatted_database_file(varname, t_step)
 
@@ -296,7 +296,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
         if (relativity .and. (rho_wrt .or. cons_vars_wrt)) then
             ! For relativistic flow, conservative and primitive densities are different
             ! Hard-coded single-component for now
-            q_sf = q_cons_vf(1)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+            q_sf(:, :, :) = q_cons_vf(1)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
             write (varname, '(A)') 'D'
             call s_write_variable_to_formatted_database_file(varname, t_step)
 
@@ -306,7 +306,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
         ! Adding the momentum to the formatted database file
         do i = 1, E_idx - mom_idx%beg
             if (mom_wrt(i) .or. cons_vars_wrt) then
-                q_sf = q_cons_vf(i + cont_idx%end)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                q_sf(:, :, :) = q_cons_vf(i + cont_idx%end)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
                 write (varname, '(A,I0)') 'mom', i
                 call s_write_variable_to_formatted_database_file(varname, t_step)
 
@@ -318,7 +318,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
         ! Adding the velocity to the formatted database file
         do i = 1, E_idx - mom_idx%beg
             if (vel_wrt(i) .or. prim_vars_wrt) then
-                q_sf = q_prim_vf(i + cont_idx%end)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                q_sf(:, :, :) = q_prim_vf(i + cont_idx%end)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
                 write (varname, '(A,I0)') 'vel', i
                 call s_write_variable_to_formatted_database_file(varname, t_step)
 
@@ -331,7 +331,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
         if (chemistry) then
             do i = 1, num_species
                 if (chem_wrt_Y(i) .or. prim_vars_wrt) then
-                    q_sf = q_prim_vf(chemxb + i - 1)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                    q_sf(:, :, :) = q_prim_vf(chemxb + i - 1)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
                     write (varname, '(A,A)') 'Y_', trim(species_names(i))
                     call s_write_variable_to_formatted_database_file(varname, t_step)
 
@@ -341,7 +341,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
             end do
 
             if (chem_wrt_T) then
-                q_sf = q_T_sf%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                q_sf(:, :, :) = q_T_sf%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
                 write (varname, '(A)') 'T'
                 call s_write_variable_to_formatted_database_file(varname, t_step)
 
@@ -364,7 +364,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
 
         ! Adding the energy to the formatted database file
         if (E_wrt .or. cons_vars_wrt) then
-            q_sf = q_cons_vf(E_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+            q_sf(:, :, :) = q_cons_vf(E_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
             write (varname, '(A)') 'E'
             call s_write_variable_to_formatted_database_file(varname, t_step)
 
@@ -375,7 +375,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
         ! Adding the magnetic field to the formatted database file
         if (mhd .and. prim_vars_wrt) then
             do i = B_idx%beg, B_idx%end
-                q_sf = q_prim_vf(i)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                q_sf(:, :, :) = q_prim_vf(i)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
 
                 ! 1D: output By, Bz
                 if (n == 0) then
@@ -404,7 +404,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
         if (elasticity) then
             do i = 1, stress_idx%end - stress_idx%beg + 1
                 if (prim_vars_wrt) then
-                    q_sf = q_prim_vf(i - 1 + stress_idx%beg)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                    q_sf(:, :, :) = q_prim_vf(i - 1 + stress_idx%beg)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
                     write (varname, '(A,I0)') 'tau', i
                     call s_write_variable_to_formatted_database_file(varname, t_step)
                 end if
@@ -415,7 +415,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
         if (hyperelasticity) then
             do i = 1, xiend - xibeg + 1
                 if (prim_vars_wrt) then
-                    q_sf = q_prim_vf(i - 1 + xibeg)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                    q_sf(:, :, :) = q_prim_vf(i - 1 + xibeg)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
                     write (varname, '(A,I0)') 'xi', i
                     call s_write_variable_to_formatted_database_file(varname, t_step)
                 end if
@@ -424,7 +424,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
         end if
 
         if (cont_damage) then
-            q_sf = q_cons_vf(damage_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+            q_sf(:, :, :) = q_cons_vf(damage_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
             write (varname, '(A)') 'damage_state'
             call s_write_variable_to_formatted_database_file(varname, t_step)
 
@@ -433,7 +433,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
 
         ! Adding the pressure to the formatted database file
         if (pres_wrt .or. prim_vars_wrt) then
-            q_sf = q_prim_vf(E_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+            q_sf(:, :, :) = q_prim_vf(E_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
             write (varname, '(A)') 'pres'
             call s_write_variable_to_formatted_database_file(varname, t_step)
 
@@ -448,7 +448,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
 
             do i = 1, num_fluids - 1
                 if (alpha_wrt(i) .or. (cons_vars_wrt .or. prim_vars_wrt)) then
-                    q_sf = q_cons_vf(i + E_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                    q_sf(:, :, :) = q_cons_vf(i + E_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
                     write (varname, '(A,I0)') 'alpha', i
                     call s_write_variable_to_formatted_database_file(varname, t_step)
 
@@ -460,7 +460,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
             if (alpha_wrt(num_fluids) &
                 .or. &
                 (cons_vars_wrt .or. prim_vars_wrt)) then
-                q_sf = q_cons_vf(adv_idx%end)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                q_sf(:, :, :) = q_cons_vf(adv_idx%end)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
                 write (varname, '(A,I0)') 'alpha', num_fluids
                 call s_write_variable_to_formatted_database_file(varname, t_step)
 
@@ -474,7 +474,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
         if (gamma_wrt &
             .or. &
             (model_eqns == 1 .and. (cons_vars_wrt .or. prim_vars_wrt))) then
-            q_sf = gamma_sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+            q_sf(:, :, :) = gamma_sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
             write (varname, '(A)') 'gamma'
             call s_write_variable_to_formatted_database_file(varname, t_step)
 
@@ -498,7 +498,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
         if (pi_inf_wrt &
             .or. &
             (model_eqns == 1 .and. (cons_vars_wrt .or. prim_vars_wrt))) then
-            q_sf = pi_inf_sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+            q_sf(:, :, :) = pi_inf_sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
             write (varname, '(A)') 'pi_inf'
             call s_write_variable_to_formatted_database_file(varname, t_step)
 
@@ -562,7 +562,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
         end do
 
         if (ib) then
-            q_sf = real(ib_markers%sf(-offset_x%beg:m + offset_x%end, -offset_y%beg:n + offset_y%end, -offset_z%beg:p + offset_z%end))
+            q_sf(:, :, :) = real(ib_markers%sf(-offset_x%beg:m + offset_x%end, -offset_y%beg:n + offset_y%end, -offset_z%beg:p + offset_z%end))
             varname = 'ib_markers'
             call s_write_variable_to_formatted_database_file(varname, t_step)
         end if
@@ -591,7 +591,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
 
         ! Adding the color function to formatted database file
         if (cf_wrt) then
-            q_sf = q_cons_vf(c_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+            q_sf(:, :, :) = q_cons_vf(c_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
             write (varname, '(A,I0)') 'color_function'
             call s_write_variable_to_formatted_database_file(varname, t_step)
             varname(:) = ' '
@@ -601,7 +601,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
         ! Adding the volume fraction(s) to the formatted database file
         if (bubbles_euler) then
             do i = adv_idx%beg, adv_idx%end
-                q_sf = q_cons_vf(i)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                q_sf(:, :, :) = q_cons_vf(i)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
                 write (varname, '(A,I0)') 'alpha', i - E_idx
                 call s_write_variable_to_formatted_database_file(varname, t_step)
                 varname(:) = ' '
@@ -612,7 +612,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
         if (bubbles_euler) then
             !nR
             do i = 1, nb
-                q_sf = q_cons_vf(bub_idx%rs(i))%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                q_sf(:, :, :) = q_cons_vf(bub_idx%rs(i))%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
                 write (varname, '(A,I3.3)') 'nR', i
                 call s_write_variable_to_formatted_database_file(varname, t_step)
                 varname(:) = ' '
@@ -620,7 +620,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
 
             !nRdot
             do i = 1, nb
-                q_sf = q_cons_vf(bub_idx%vs(i))%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                q_sf(:, :, :) = q_cons_vf(bub_idx%vs(i))%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
                 write (varname, '(A,I3.3)') 'nV', i
                 call s_write_variable_to_formatted_database_file(varname, t_step)
                 varname(:) = ' '
@@ -628,7 +628,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
             if ((polytropic .neqv. .true.) .and. (.not. qbmm)) then
                 !nP
                 do i = 1, nb
-                    q_sf = q_cons_vf(bub_idx%ps(i))%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                    q_sf(:, :, :) = q_cons_vf(bub_idx%ps(i))%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
                     write (varname, '(A,I3.3)') 'nP', i
                     call s_write_variable_to_formatted_database_file(varname, t_step)
                     varname(:) = ' '
@@ -636,7 +636,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
 
                 !nM
                 do i = 1, nb
-                    q_sf = q_cons_vf(bub_idx%ms(i))%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                    q_sf(:, :, :) = q_cons_vf(bub_idx%ms(i))%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
                     write (varname, '(A,I3.3)') 'nM', i
                     call s_write_variable_to_formatted_database_file(varname, t_step)
                     varname(:) = ' '
@@ -645,7 +645,7 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
 
             ! number density
             if (adv_n) then
-                q_sf = q_cons_vf(n_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
+                q_sf(:, :, :) = q_cons_vf(n_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end)
                 write (varname, '(A)') 'n'
                 call s_write_variable_to_formatted_database_file(varname, t_step)
                 varname(:) = ' '
@@ -655,10 +655,10 @@ impure subroutine s_save_data(t_step, varname, pres, c, H)
         ! Adding the lagrangian subgrid variables  to the formatted database file
         if (bubbles_lagrange) then
             !! Void fraction field
-            q_sf = 1._wp - q_cons_vf(beta_idx)%sf( &
-                   -offset_x%beg:m + offset_x%end, &
-                   -offset_y%beg:n + offset_y%end, &
-                   -offset_z%beg:p + offset_z%end)
+            q_sf(:, :, :) = 1._wp - q_cons_vf(beta_idx)%sf( &
+                            -offset_x%beg:m + offset_x%end, &
+                            -offset_y%beg:n + offset_y%end, &
+                            -offset_z%beg:p + offset_z%end)
             write (varname, '(A)') 'voidFraction'
             call s_write_variable_to_formatted_database_file(varname, t_step)
             varname(:) = ' '
diff --git a/src/pre_process/m_assign_variables.fpp b/src/pre_process/m_assign_variables.fpp
index 11eb62f9e6..e2398d4c84 100644
--- a/src/pre_process/m_assign_variables.fpp
+++ b/src/pre_process/m_assign_variables.fpp
@@ -244,7 +244,7 @@ contains
         rhoH = (1._wp - vfH)/ratio
         deno = 1._wp - (1._wp - q_prim_vf(alf_idx)%sf(j, k, l))/rhoH
 
-        if (deno == 0._wp) then
+        if (f_approx_equal(deno, 0._wp)) then
             velH = 0._wp
         else
             velH = (q_prim_vf(E_idx)%sf(j, k, l) - 1._wp)/(1._wp - q_prim_vf(alf_idx)%sf(j, k, l))/deno
@@ -502,7 +502,7 @@ contains
                 theta = atan2(y_cc(k), x_cc(j))
                 phi = atan2(sqrt(x_cc(j)**2 + y_cc(k)**2), z_cc(l))
                 !spherical coord, assuming Rmax=1
-                xi_sph = (rcoord**3 - R0ref**3 + 1_wp)**(1_wp/3_wp)
+                xi_sph = (rcoord**3 - R0ref**3 + 1_wp)**(1._wp/3._wp)
                 xi_cart(1) = xi_sph*sin(phi)*cos(theta)
                 xi_cart(2) = xi_sph*sin(phi)*sin(theta)
                 xi_cart(3) = xi_sph*cos(phi)
diff --git a/src/pre_process/m_check_patches.fpp b/src/pre_process/m_check_patches.fpp
index d1a8aade96..969b6cd471 100644
--- a/src/pre_process/m_check_patches.fpp
+++ b/src/pre_process/m_check_patches.fpp
@@ -318,7 +318,7 @@ contains
         @:PROHIBIT(f_is_default(patch_icpp(patch_id)%x_centroid), "Spherical harmonic patch "//trim(iStr)//": x_centroid must be set")
         @:PROHIBIT(f_is_default(patch_icpp(patch_id)%y_centroid), "Spherical harmonic patch "//trim(iStr)//": y_centroid must be set")
         @:PROHIBIT(f_is_default(patch_icpp(patch_id)%z_centroid), "Spherical harmonic patch "//trim(iStr)//": z_centroid must be set")
-        @:PROHIBIT(all(patch_icpp(patch_id)%epsilon /= (/1._wp, 2._wp, 3._wp, 4._wp, 5._wp/)), &
+        @:PROHIBIT(.not. f_approx_in_array(patch_icpp(patch_id)%epsilon, (/1._wp, 2._wp, 3._wp, 4._wp, 5._wp/)), &
             "Spherical harmonic patch "//trim(iStr)//": epsilon must be one of 1, 2, 3, 4, 5")
         @:PROHIBIT(patch_icpp(patch_id)%beta < 0._wp, &
             "Spherical harmonic patch "//trim(iStr)//": beta must be greater than or equal to zero")
@@ -514,11 +514,11 @@ contains
 
         @:PROHIBIT(f_is_default(patch_icpp(patch_id)%vel(1)), &
             "Patch "//trim(iStr)//": vel(1) must be set")
-        @:PROHIBIT(n == 0 .and. (.not. f_is_default(patch_icpp(patch_id)%vel(2))) .and. patch_icpp(patch_id)%vel(2) /= 0 .and. (.not. mhd), &
+        @:PROHIBIT(n == 0 .and. (.not. f_is_default(patch_icpp(patch_id)%vel(2))) .and. (.not. f_approx_equal(patch_icpp(patch_id)%vel(2) , 0._wp)) .and. (.not. mhd), &
             "Patch "//trim(iStr)//": vel(2) must not be set when n = 0")
         @:PROHIBIT(n > 0 .and. f_is_default(patch_icpp(patch_id)%vel(2)), &
             "Patch "//trim(iStr)//": vel(2) must be set when n > 0")
-        @:PROHIBIT(p == 0 .and. (.not. f_is_default(patch_icpp(patch_id)%vel(3))) .and. patch_icpp(patch_id)%vel(3) /= 0 .and. (.not. mhd), &
+        @:PROHIBIT(p == 0 .and. (.not. f_is_default(patch_icpp(patch_id)%vel(3))) .and. (.not. f_approx_equal(patch_icpp(patch_id)%vel(3) , 0._wp)) .and. (.not. mhd), &
             "Patch "//trim(iStr)//": vel(3) must not be set when p = 0")
         @:PROHIBIT(p > 0 .and. f_is_default(patch_icpp(patch_id)%vel(3)), &
             "Patch "//trim(iStr)//": vel(3) must be set when p > 0")
diff --git a/src/pre_process/m_checker.fpp b/src/pre_process/m_checker.fpp
index 46a7c55941..972f3bcb70 100644
--- a/src/pre_process/m_checker.fpp
+++ b/src/pre_process/m_checker.fpp
@@ -88,11 +88,11 @@ contains
             "n must be positive (2D or 3D) for cylindrical coordinates")
         @:PROHIBIT(cyl_coord .and. (f_is_default(y_domain%beg) .or. f_is_default(y_domain%end)), &
             "y_domain%beg and y_domain%end must be set for n = 0 (2D cylindrical coordinates)")
-        @:PROHIBIT(cyl_coord .and. (y_domain%beg /= 0._wp .or. y_domain%end <= 0._wp), &
+        @:PROHIBIT(cyl_coord .and. ((.not. f_approx_equal(y_domain%beg , 0._wp)) .or. y_domain%end <= 0._wp), &
             "y_domain%beg must be 0 and y_domain%end must be positive for cylindrical coordinates")
         @:PROHIBIT(cyl_coord .and. p == 0 .and. ((.not. f_is_default(z_domain%beg)) .or. (.not. f_is_default(z_domain%end))), &
             "z_domain%beg and z_domain%end are not supported for p = 0 (2D cylindrical coordinates)")
-        @:PROHIBIT(cyl_coord .and. p > 0 .and. (z_domain%beg /= 0._wp .or. z_domain%end /= 2._wp*pi), &
+        @:PROHIBIT(cyl_coord .and. p > 0 .and. (.not. (f_approx_equal(z_domain%beg, 0._wp) .and. f_approx_equal(z_domain%end, 2._wp*pi))), &
             "z_domain%beg must be 0 and z_domain%end must be 2*pi for 3D cylindrical coordinates")
 
         @:PROHIBIT(num_patches < 0)
diff --git a/src/pre_process/m_compute_levelset.fpp b/src/pre_process/m_compute_levelset.fpp
index 25958a36b3..17f66f8d68 100644
--- a/src/pre_process/m_compute_levelset.fpp
+++ b/src/pre_process/m_compute_levelset.fpp
@@ -14,6 +14,8 @@ module m_compute_levelset
 
     use m_mpi_proxy            !< Message passing interface (MPI) module proxy
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     implicit none
 
     private; public :: s_cylinder_levelset, s_circle_levelset, &
@@ -49,7 +51,7 @@ contains
                 dist_vec(3) = 0
                 dist = sqrt(sum(dist_vec**2))
                 levelset%sf(i, j, 0, ib_patch_id) = dist - radius
-                if (dist == 0) then
+                if (f_approx_equal(dist, 0._wp)) then
                     levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0
                 else
                     levelset_norm%sf(i, j, 0, ib_patch_id, :) = &
@@ -132,7 +134,7 @@ contains
                 end if
 
                 levelset%sf(i, j, 0, ib_patch_id) = dist
-                if (dist == 0) then
+                if (f_approx_equal(dist, 0._wp)) then
                     levelset_norm%sf(i, j, 0, ib_patch_id, :) = 0
                 else
                     levelset_norm%sf(i, j, 0, ib_patch_id, :) = &
@@ -226,14 +228,14 @@ contains
 
                     if (dist_side < dist_surf) then
                         levelset%sf(i, j, l, ib_patch_id) = dist_side
-                        if (dist_side == abs(z_cc(l) - z_min)) then
+                        if (f_approx_equal(dist_side, abs(z_cc(l) - z_min))) then
                             levelset_norm%sf(i, j, l, ib_patch_id, :) = (/0, 0, -1/)
                         else
                             levelset_norm%sf(i, j, l, ib_patch_id, :) = (/0, 0, 1/)
                         end if
                     else
                         levelset%sf(i, j, l, ib_patch_id) = dist_surf
-                        if (dist == 0) then
+                        if (f_approx_equal(dist, 0._wp)) then
                             levelset_norm%sf(i, j, l, ib_patch_id, :) = 0
                         else
                             levelset_norm%sf(i, j, l, ib_patch_id, :) = &
@@ -300,7 +302,7 @@ contains
 
                     if (idx == 1) then
                         levelset%sf(i, j, 0, ib_patch_id) = side_dists(1)
-                        if (side_dists(1) == 0) then
+                        if (f_approx_equal(side_dists(1), 0._wp)) then
                             levelset_norm%sf(i, j, 0, ib_patch_id, 1) = 0._wp
                         else
                             levelset_norm%sf(i, j, 0, ib_patch_id, 1) = side_dists(1)/ &
@@ -309,7 +311,7 @@ contains
 
                     else if (idx == 2) then
                         levelset%sf(i, j, 0, ib_patch_id) = side_dists(2)
-                        if (side_dists(2) == 0) then
+                        if (f_approx_equal(side_dists(2), 0._wp)) then
                             levelset_norm%sf(i, j, 0, ib_patch_id, 1) = 0._wp
                         else
                             levelset_norm%sf(i, j, 0, ib_patch_id, 1) = -side_dists(2)/ &
@@ -318,7 +320,7 @@ contains
 
                     else if (idx == 3) then
                         levelset%sf(i, j, 0, ib_patch_id) = side_dists(3)
-                        if (side_dists(3) == 0) then
+                        if (f_approx_equal(side_dists(3), 0._wp)) then
 
                             levelset_norm%sf(i, j, 0, ib_patch_id, 2) = 0._wp
                         else
@@ -328,7 +330,7 @@ contains
 
                     else if (idx == 4) then
                         levelset%sf(i, j, 0, ib_patch_id) = side_dists(4)
-                        if (side_dists(4) == 0) then
+                        if (f_approx_equal(side_dists(4), 0._wp)) then
 
                             levelset_norm%sf(i, j, 0, ib_patch_id, 2) = 0._wp
                         else
@@ -398,54 +400,54 @@ contains
 
                         min_dist = minval(abs(side_dists))
 
-                        if (min_dist == abs(side_dists(1))) then
+                        if (f_approx_equal(min_dist, abs(side_dists(1)))) then
                             levelset%sf(i, j, k, ib_patch_id) = side_dists(1)
-                            if (side_dists(1) == 0) then
+                            if (f_approx_equal(side_dists(1), 0._wp)) then
                                 levelset_norm%sf(i, j, k, ib_patch_id, 1) = 0._wp
                             else
                                 levelset_norm%sf(i, j, k, ib_patch_id, 1) = side_dists(1)/ &
                                                                             abs(side_dists(1))
                             end if
 
-                        else if (min_dist == abs(side_dists(2))) then
+                        else if (f_approx_equal(min_dist, abs(side_dists(2)))) then
                             levelset%sf(i, j, k, ib_patch_id) = side_dists(2)
-                            if (side_dists(2) == 0) then
+                            if (f_approx_equal(side_dists(2), 0._wp)) then
                                 levelset_norm%sf(i, j, k, ib_patch_id, 1) = 0._wp
                             else
                                 levelset_norm%sf(i, j, k, ib_patch_id, 1) = -side_dists(2)/ &
                                                                             abs(side_dists(2))
                             end if
 
-                        else if (min_dist == abs(side_dists(3))) then
+                        else if (f_approx_equal(min_dist, abs(side_dists(3)))) then
                             levelset%sf(i, j, k, ib_patch_id) = side_dists(3)
-                            if (side_dists(3) == 0) then
+                            if (f_approx_equal(side_dists(3), 0._wp)) then
                                 levelset_norm%sf(i, j, k, ib_patch_id, 2) = 0._wp
                             else
                                 levelset_norm%sf(i, j, k, ib_patch_id, 2) = side_dists(3)/ &
                                                                             abs(side_dists(3))
                             end if
 
-                        else if (min_dist == abs(side_dists(4))) then
+                        else if (f_approx_equal(min_dist, abs(side_dists(4)))) then
                             levelset%sf(i, j, k, ib_patch_id) = side_dists(4)
-                            if (side_dists(4) == 0) then
+                            if (f_approx_equal(side_dists(4), 0._wp)) then
                                 levelset_norm%sf(i, j, k, ib_patch_id, 2) = 0._wp
                             else
                                 levelset_norm%sf(i, j, k, ib_patch_id, 2) = -side_dists(4)/ &
                                                                             abs(side_dists(4))
                             end if
 
-                        else if (min_dist == abs(side_dists(5))) then
+                        else if (f_approx_equal(min_dist, abs(side_dists(5)))) then
                             levelset%sf(i, j, k, ib_patch_id) = side_dists(5)
-                            if (side_dists(5) == 0) then
+                            if (f_approx_equal(side_dists(5), 0._wp)) then
                                 levelset_norm%sf(i, j, k, ib_patch_id, 3) = 0._wp
                             else
                                 levelset_norm%sf(i, j, k, ib_patch_id, 3) = side_dists(5)/ &
                                                                             abs(side_dists(5))
                             end if
 
-                        else if (min_dist == abs(side_dists(6))) then
+                        else if (f_approx_equal(min_dist, abs(side_dists(6)))) then
                             levelset%sf(i, j, k, ib_patch_id) = side_dists(6)
-                            if (side_dists(6) == 0) then
+                            if (f_approx_equal(side_dists(6), 0._wp)) then
                                 levelset_norm%sf(i, j, k, ib_patch_id, 3) = 0._wp
                             else
                                 levelset_norm%sf(i, j, k, ib_patch_id, 3) = -side_dists(6)/ &
@@ -487,7 +489,7 @@ contains
                     dist_vec(3) = z_cc(k) - z_centroid
                     dist = sqrt(sum(dist_vec**2))
                     levelset%sf(i, j, k, ib_patch_id) = dist - radius
-                    if (dist == 0) then
+                    if (f_approx_equal(dist, 0._wp)) then
                         levelset_norm%sf(i, j, k, ib_patch_id, :) = (/1, 0, 0/)
                     else
                         levelset_norm%sf(i, j, k, ib_patch_id, :) = &
@@ -521,17 +523,17 @@ contains
         length_y = patch_ib(ib_patch_id)%length_y
         length_z = patch_ib(ib_patch_id)%length_z
 
-        if (length_x /= 0._wp) then
+        if (.not. f_approx_equal(length_x, 0._wp)) then
             boundary%beg = x_centroid - 0.5_wp*length_x
             boundary%end = x_centroid + 0.5_wp*length_x
             dist_sides_vec = (/1, 0, 0/)
             dist_surface_vec = (/0, 1, 1/)
-        else if (length_y /= 0._wp) then
+        else if (.not. f_approx_equal(length_y, 0._wp)) then
             boundary%beg = y_centroid - 0.5_wp*length_y
             boundary%end = y_centroid + 0.5_wp*length_y
             dist_sides_vec = (/0, 1, 0/)
             dist_surface_vec = (/1, 0, 1/)
-        else if (length_z /= 0._wp) then
+        else if (.not. f_approx_equal(length_z, 0._wp)) then
             boundary%beg = z_centroid - 0.5_wp*length_z
             boundary%end = z_centroid + 0.5_wp*length_z
             dist_sides_vec = (/0, 0, 1/)
@@ -552,7 +554,7 @@ contains
 
                     if (dist_side < abs(dist_surface)) then
                         levelset%sf(i, j, k, ib_patch_id) = -dist_side
-                        if (dist_side == abs(side_pos - boundary%beg)) then
+                        if (f_approx_equal(dist_side, abs(side_pos - boundary%beg))) then
                             levelset_norm%sf(i, j, k, ib_patch_id, :) = -dist_sides_vec
                         else
                             levelset_norm%sf(i, j, k, ib_patch_id, :) = dist_sides_vec
diff --git a/src/pre_process/m_grid.f90 b/src/pre_process/m_grid.f90
index 2c73bc4fe1..07277433d2 100644
--- a/src/pre_process/m_grid.f90
+++ b/src/pre_process/m_grid.f90
@@ -20,6 +20,8 @@ module m_grid
 
     use m_mpi_proxy             ! Message passing interface (MPI) module proxy
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
 #ifdef MFC_MPI
     use mpi                     ! Message passing interface (MPI) module
 #endif
@@ -83,7 +85,7 @@ impure subroutine s_generate_serial_grid
             end do
             x_cb = x_cb*length
 
-            x_cc = (x_cb(0:m) + x_cb(-1:m - 1))/2._wp
+            x_cc(0:m) = (x_cb(0:m) + x_cb(-1:m - 1))/2._wp
 
             dx = minval(x_cb(0:m) - x_cb(-1:m - 1))
             print *, 'Stretched grid: min/max x grid: ', minval(x_cc(:)), maxval(x_cc(:))
@@ -94,7 +96,7 @@ impure subroutine s_generate_serial_grid
         ! Grid Generation in the y-direction
         if (n == 0) return
 
-        if (grid_geometry == 2 .and. y_domain%beg == 0.0_wp) then
+        if (grid_geometry == 2 .and. f_approx_equal(y_domain%beg, 0.0_wp)) then
             !IF (grid_geometry == 2) THEN
 
             dy = (y_domain%end - y_domain%beg)/real(2*n + 1, wp)
@@ -137,7 +139,7 @@ impure subroutine s_generate_serial_grid
             end do
 
             y_cb = y_cb*length
-            y_cc = (y_cb(0:n) + y_cb(-1:n - 1))/2._wp
+            y_cc(0:m) = (y_cb(0:n) + y_cb(-1:n - 1))/2._wp
 
             dy = minval(y_cb(0:n) - y_cb(-1:n - 1))
 
@@ -174,7 +176,7 @@ impure subroutine s_generate_serial_grid
             end do
 
             z_cb = z_cb*length
-            z_cc = (z_cb(0:p) + z_cb(-1:p - 1))/2._wp
+            z_cc(0:m) = (z_cb(0:p) + z_cb(-1:p - 1))/2._wp
 
             dz = minval(z_cb(0:p) - z_cb(-1:p - 1))
 
@@ -241,7 +243,7 @@ impure subroutine s_generate_parallel_grid
         ! Grid generation in the y-direction
         if (n_glb > 0) then
 
-            if (grid_geometry == 2 .and. y_domain%beg == 0.0_wp) then
+            if (grid_geometry == 2 .and. f_approx_equal(y_domain%beg, 0.0_wp)) then
                 dy = (y_domain%end - y_domain%beg)/real(2*n_glb + 1, wp)
                 y_cb_glb(-1) = y_domain%beg
                 do i = 1, n_glb
diff --git a/src/pre_process/m_model.fpp b/src/pre_process/m_model.fpp
index 2936fa0a3c..e42c94a82f 100644
--- a/src/pre_process/m_model.fpp
+++ b/src/pre_process/m_model.fpp
@@ -327,7 +327,7 @@ contains
 
         character(kind=c_char, len=80), parameter :: header = "Model file written by MFC."
         integer(kind=c_int32_t) :: nTriangles
-        real(kind=c_float) :: normal(3), v(3)
+        real(wp) :: normal(3), v(3)
         integer(kind=c_int16_t) :: attribute
 
         open (newunit=iunit, file=filepath, action='WRITE', &
diff --git a/src/pre_process/m_patches.fpp b/src/pre_process/m_patches.fpp
index ec9b2baf05..91cb219e2d 100644
--- a/src/pre_process/m_patches.fpp
+++ b/src/pre_process/m_patches.fpp
@@ -48,7 +48,6 @@ module m_patches
     !! is to act as a pseudo volume fraction to indicate the contribution of each
     !! patch toward the composition of a cell's fluid state.
 
-    real(wp) :: r_cyl, theta_cyl, x_cart, y_cart, z_cart
     real(wp) :: cart_x, cart_y, cart_z
     real(wp) :: sph_phi !<
     !! Variables to be used to hold cell locations in Cartesian coordinates if
@@ -551,7 +550,7 @@ contains
                         do while (airfoil_grid_u(k)%x < x_act)
                             k = k + 1
                         end do
-                        if (airfoil_grid_u(k)%x == x_act) then
+                        if (f_approx_equal(airfoil_grid_u(k)%x, x_act)) then
                             if (y_act <= airfoil_grid_u(k)%y) then
                                 !!IB
                                 !call s_assign_patch_primitive_variables(patch_id, i, j, 0, &
@@ -572,7 +571,7 @@ contains
                         do while (airfoil_grid_l(k)%x < x_act)
                             k = k + 1
                         end do
-                        if (airfoil_grid_l(k)%x == x_act) then
+                        if (f_approx_equal(airfoil_grid_l(k)%x, x_act)) then
                             if (y_act >= airfoil_grid_l(k)%y) then
                                 !!IB
                                 !call s_assign_patch_primitive_variables(patch_id, i, j, 0, &
@@ -720,7 +719,7 @@ contains
                                 do while (airfoil_grid_u(k)%x < x_act)
                                     k = k + 1
                                 end do
-                                if (airfoil_grid_u(k)%x == x_act) then
+                                if (f_approx_equal(airfoil_grid_u(k)%x, x_act)) then
                                     if (y_act <= airfoil_grid_u(k)%y) then
                                         !!IB
                                         !call s_assign_patch_primitive_variables(patch_id, i, j, 0, &
@@ -741,7 +740,7 @@ contains
                                 do while (airfoil_grid_l(k)%x < x_act)
                                     k = k + 1
                                 end do
-                                if (airfoil_grid_l(k)%x == x_act) then
+                                if (f_approx_equal(airfoil_grid_l(k)%x, x_act)) then
                                     if (y_act >= airfoil_grid_l(k)%y) then
                                         !!IB
                                         !call s_assign_patch_primitive_variables(patch_id, i, j, 0, &
@@ -1570,8 +1569,6 @@ contains
 
         integer :: i, j, k !< generic loop iterators
 
-        complex(wp), parameter :: cmplx_i = (0._wp, 1._wp)
-
         ! Transferring the patch's centroid and radius information
         x_centroid = patch_icpp(patch_id)%x_centroid
         y_centroid = patch_icpp(patch_id)%y_centroid
diff --git a/src/pre_process/m_perturbation.fpp b/src/pre_process/m_perturbation.fpp
index 209c752287..15e12dfe6e 100644
--- a/src/pre_process/m_perturbation.fpp
+++ b/src/pre_process/m_perturbation.fpp
@@ -16,6 +16,8 @@ module m_perturbation
 
     use m_boundary_common   ! Boundary conditions module
 
+    use m_helper_basic           !< Functions to compare floating point numbers
+
     use ieee_arithmetic
 
     implicit none
@@ -72,7 +74,7 @@ contains
 
                     ! Perturb partial density fields to match perturbed volume fraction fields
                     !    IF ((perturb_alpha >= 25e-2_wp) .AND. (perturb_alpha <= 75e-2_wp)) THEN
-                    if ((perturb_alpha /= 0._wp) .and. (perturb_alpha /= 1._wp)) then
+                    if ((.not. f_approx_equal(perturb_alpha, 0._wp)) .and. (.not. f_approx_equal(perturb_alpha, 1._wp))) then
 
                         ! Derive new partial densities
                         do l = 1, num_fluids
@@ -350,7 +352,7 @@ contains
         call cg(mixlayer_nvar*n - n_bc_skip, mixlayer_nvar*n - n_bc_skip, hr, hi, wr, wi, zr, zi, fv1, fv2, fv3, ierr)
 
         ! Generate instability wave
-        call s_generate_wave(wr, wi, zr, zi, rho_mean, mach, alpha, beta, wave, shift)
+        call s_generate_wave(wi, zr, zi, rho_mean, mach, alpha, beta, wave, shift)
 
     end subroutine s_solve_linear_system
 
@@ -492,8 +494,8 @@ contains
     !>  This subroutine generates an instability wave using the most unstable
         !!              eigenvalue and corresponding eigenvector among the
         !!              given set of eigenvalues and eigenvectors.
-    pure subroutine s_generate_wave(wr, wi, zr, zi, rho_mean, mach, alpha, beta, wave, shift)
-        real(wp), dimension(0:mixlayer_nvar*n - n_bc_skip - 1), intent(in) :: wr, wi !< eigenvalues
+    pure subroutine s_generate_wave(wi, zr, zi, rho_mean, mach, alpha, beta, wave, shift)
+        real(wp), dimension(0:mixlayer_nvar*n - n_bc_skip - 1), intent(in) :: wi !< eigenvalues
         real(wp), dimension(0:mixlayer_nvar*n - n_bc_skip - 1, 0:mixlayer_nvar*n - n_bc_skip - 1), intent(in) :: zr, zi !< eigenvectors
         real(wp), intent(in) :: rho_mean
         real(wp), dimension(mixlayer_nvar, 0:m, 0:n, 0:p), intent(inout) :: wave
@@ -597,7 +599,7 @@ contains
         do i = 0, m
             do j = 0, n
                 do k = 0, p
-                    if (beta == 0) then
+                    if (f_approx_equal(beta, 0._wp)) then
                         ang = alpha*(x_cc(i)*xratio)
                     else
                         ang = alpha*(x_cc(i)*xratio) + beta*(z_cc(k)*xratio) + shift
diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp
index e33491576f..4cd5bf0349 100644
--- a/src/pre_process/m_start_up.fpp
+++ b/src/pre_process/m_start_up.fpp
@@ -775,8 +775,8 @@ contains
         end if
         !Initialize pb based on surface tension for qbmm (polytropic)
         if (qbmm .and. polytropic .and. (.not. f_is_default(Web))) then
-            pb0 = pref + 2._wp*fluid_pp(1)%ss/(R0*R0ref)
-            pb0 = pb0/pref
+            pb0(:) = pref + 2._wp*fluid_pp(1)%ss/(R0(:)*R0ref)
+            pb0(:) = pb0(:)/pref
             pref = 1._wp
         end if
         call s_initialize_mpi_common_module()
diff --git a/src/simulation/m_bubbles.fpp b/src/simulation/m_bubbles.fpp
index 4ff8a54df0..187d78c2c0 100644
--- a/src/simulation/m_bubbles.fpp
+++ b/src/simulation/m_bubbles.fpp
@@ -15,6 +15,8 @@ module m_bubbles
 
     use m_variables_conversion !< State variables type conversion procedures
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     implicit none
 
     real(wp) :: chi_vw  !< Bubble wall properties (Ando 2010)
@@ -569,7 +571,7 @@ contains
             end do
 
             ! Exit the loop if the final time reached dt
-            if (t_new == 0.5_wp*dt .or. iter_count >= adap_dt_max_iters) exit
+            if (f_approx_equal(t_new, 0.5_wp*dt) .or. iter_count >= adap_dt_max_iters) exit
 
         end do
 
diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
index 008840caa1..854e3f63a0 100644
--- a/src/simulation/m_bubbles_EL.fpp
+++ b/src/simulation/m_bubbles_EL.fpp
@@ -21,6 +21,8 @@ module m_bubbles_EL
 
     use m_boundary_common
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     use m_sim_helpers
 
     use m_helper
@@ -348,7 +350,7 @@ contains
 
         ! Initial particle pressure
         gas_p(bub_id, 1) = pliq + 2._wp*(1._wp/Web)/bub_R0(bub_id)
-        if ((1._wp/Web) /= 0._wp) then
+        if (.not. f_approx_equal((1._wp/Web), 0._wp)) then
             pcrit = pv - 4._wp*(1._wp/Web)/(3._wp*sqrt(3._wp*gas_p(bub_id, 1)*bub_R0(bub_id)**3._wp/(2._wp*(1._wp/Web))))
             pref = gas_p(bub_id, 1)
         else
diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp
index 1645ffef25..665fc6b34b 100644
--- a/src/simulation/m_cbc.fpp
+++ b/src/simulation/m_cbc.fpp
@@ -49,8 +49,6 @@ module m_cbc
     real(wp), allocatable, dimension(:, :, :, :) :: q_prim_rsy_vf
     real(wp), allocatable, dimension(:, :, :, :) :: q_prim_rsz_vf
 
-    type(scalar_field), allocatable, dimension(:) :: F_rs_vf, F_src_rs_vf !<
-
     !! Cell-average fluxes (src - source). These are directly determined from the
     !! cell-average primitive variables, q_prims_rs_vf, and not a Riemann solver.
 
@@ -68,10 +66,6 @@ module m_cbc
     real(wp), allocatable, dimension(:, :, :, :) :: flux_rsy_vf_l, flux_src_rsy_vf_l
     real(wp), allocatable, dimension(:, :, :, :) :: flux_rsz_vf_l, flux_src_rsz_vf_l
 
-    real(wp) :: c           !< Cell averaged speed of sound
-    real(wp), dimension(2) :: Re          !< Cell averaged Reynolds numbers
-    !$acc declare create(c, Re)
-
     real(wp) :: dpres_ds !< Spatial derivatives in s-dir of pressure
     !$acc declare create(dpres_ds)
 
diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp
index 66ee9b9921..3879bd8850 100644
--- a/src/simulation/m_checker.fpp
+++ b/src/simulation/m_checker.fpp
@@ -14,6 +14,8 @@ module m_checker
 
     use m_helper
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     implicit none
 
     private; public :: s_check_inputs
@@ -73,9 +75,9 @@ contains
             "weno_order != 1, but weno_eps is not set. A typical value of weno_eps is 1e-6")
 
         @:PROHIBIT(weno_eps <= 0._wp, "weno_eps must be positive. A typical value of weno_eps is 1e-6")
-        @:PROHIBIT(wenoz .and. weno_order == real(7, wp) .and. f_is_default(real(wenoz_q, wp)), &
+        @:PROHIBIT(wenoz .and. weno_order == 7 .and. f_is_default(real(wenoz_q, wp)), &
             "wenoz is used at 7th order, but wenoz_q is not set. It should be either 2, 3, or 4")
-        @:PROHIBIT(wenoz .and. weno_order == real(7, wp) .and. .not. (f_approx_equal(real(wenoz_q, wp), real(2, wp)) .or. &
+        @:PROHIBIT(wenoz .and. weno_order == 7 .and. .not. (f_approx_equal(real(wenoz_q, wp), real(2, wp)) .or. &
             f_approx_equal(real(wenoz_q, wp), real(3, wp)) .or. f_approx_equal(real(wenoz_q, wp), real(4, wp))), &
             "wenoz_q must be either 2, 3, or 4")
         @:PROHIBIT(teno .and. f_is_default(teno_CT), "teno is used, but teno_CT is not set. A typical value of teno_CT is 1e-6")
diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp
index 98a8a73046..a0cc35fbbd 100644
--- a/src/simulation/m_data_output.fpp
+++ b/src/simulation/m_data_output.fpp
@@ -25,6 +25,8 @@ module m_data_output
 
     use m_helper
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     use m_sim_helpers
 
     use m_delay_file_access
@@ -360,7 +362,7 @@ contains
                     t_step, dt, t_step*dt, icfl_max_glb
             end if
 
-            if (icfl_max_glb /= icfl_max_glb) then
+            if (.not. f_approx_equal(icfl_max_glb, icfl_max_glb)) then
                 call s_mpi_abort('ICFL is NaN. Exiting.')
             elseif (icfl_max_glb > 1._wp) then
                 print *, 'icfl', icfl_max_glb
@@ -368,7 +370,7 @@ contains
             end if
 
             if (viscous) then
-                if (vcfl_max_glb /= vcfl_max_glb) then
+                if (.not. f_approx_equal(vcfl_max_glb, vcfl_max_glb)) then
                     call s_mpi_abort('VCFL is NaN. Exiting.')
                 elseif (vcfl_max_glb > 1._wp) then
                     print *, 'vcfl', vcfl_max_glb
diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp
index 0314aad864..806b45c2da 100644
--- a/src/simulation/m_ibm.fpp
+++ b/src/simulation/m_ibm.fpp
@@ -18,6 +18,8 @@ module m_ibm
 
     use m_helper
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     use m_constants
 
     implicit none
@@ -145,7 +147,7 @@ contains
         real(wp) :: qv_K
         real(wp), dimension(num_fluids) :: Gs
 
-        real(wp) :: pres_IP, coeff
+        real(wp) :: pres_IP
         real(wp), dimension(3) :: vel_IP, vel_norm_IP
         real(wp) :: c_IP
         real(wp), dimension(num_fluids) :: alpha_rho_IP, alpha_IP
@@ -164,7 +166,7 @@ contains
         type(ghost_point) :: gp
         type(ghost_point) :: innerp
 
-        !$acc parallel loop gang vector private(physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, vel_g, vel_norm_IP, r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, G_K, Gs, gp, innerp, norm, buf, j, k, l, q, coeff)
+        !$acc parallel loop gang vector private(physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, vel_g, vel_norm_IP, r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, G_K, Gs, gp, innerp, norm, buf, j, k, l, q)
         do i = 1, num_gps
 
             gp = ghost_points(i)
@@ -408,7 +410,7 @@ contains
                     bound = p
                 end if
 
-                if (norm(dim) == 0) then
+                if (f_approx_equal(norm(dim), 0._wp)) then
                     ghost_points(q)%ip_grid(dim) = ghost_points(q)%loc(dim)
                 else
                     if (norm(dim) > 0) then
diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp
index ca4a4e3e8e..93d864c5e8 100644
--- a/src/simulation/m_mpi_proxy.fpp
+++ b/src/simulation/m_mpi_proxy.fpp
@@ -253,8 +253,6 @@ contains
 
         integer :: pack_offset, unpack_offset
 
-        integer, pointer :: p_i_send, p_i_recv
-
 #ifdef MFC_MPI
 
         call nvtxStartRange("IB-MARKER-COMM-PACKBUF")
diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp
index 93bcd55113..6a7a96cd08 100644
--- a/src/simulation/m_rhs.fpp
+++ b/src/simulation/m_rhs.fpp
@@ -600,7 +600,6 @@ contains
         real(wp), intent(inout) :: time_avg
         integer, intent(in) :: stage
 
-        real(wp), dimension(0:m, 0:n, 0:p) :: nbub
         real(wp) :: t_start, t_finish
         integer :: i, j, k, l, id !< Generic loop iterators
 
diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp
index 36d831a439..7468fa605e 100644
--- a/src/simulation/m_riemann_solvers.fpp
+++ b/src/simulation/m_riemann_solvers.fpp
@@ -39,6 +39,8 @@ module m_riemann_solvers
 
     use m_surface_tension      !< To get the capilary fluxes
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     use m_chemistry
 
     use m_thermochem, only: &
@@ -2187,7 +2189,7 @@ contains
                                                                                                           rho_R*R3V2Rbar/R3Rbar)
                                     end if
 
-                                    if ((ptilde_L /= ptilde_L) .or. (ptilde_R /= ptilde_R)) then
+                                    if ((.not. f_approx_equal(ptilde_L, ptilde_L)) .or. (.not. f_approx_equal(ptilde_R, ptilde_R))) then
                                     end if
 
                                     rho_avg = 5e-1_wp*(rho_L + rho_R)
diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp
index 330cb75712..22665f71c8 100644
--- a/src/simulation/m_start_up.fpp
+++ b/src/simulation/m_start_up.fpp
@@ -103,8 +103,6 @@ module m_start_up
  s_save_performance_metrics
 
 
-    type(scalar_field), allocatable, dimension(:) :: grad_x_vf, grad_y_vf, grad_z_vf, norm_vf
-
     real(wp) :: dt_init
 
 contains
@@ -1251,7 +1249,7 @@ contains
         end if
         !Initialize pb based on surface tension for qbmm (polytropic)
         if (qbmm .and. polytropic .and. (.not. f_is_default(Web))) then
-            pb0 = pref + 2._wp*fluid_pp(1)%ss/(R0*R0ref)
+            pb0(:) = pref + 2._wp*fluid_pp(1)%ss/(R0(:)*R0ref)
             pb0 = pb0/pref
             pref = 1._wp
         end if
diff --git a/src/simulation/m_weno.fpp b/src/simulation/m_weno.fpp
index 331479a5cb..b4dc5ad0c9 100644
--- a/src/simulation/m_weno.fpp
+++ b/src/simulation/m_weno.fpp
@@ -95,9 +95,6 @@ module m_weno
     !
     !> @}
 
-    real(wp) :: test
-    !$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, &