diff --git a/sources/evolution.F90 b/sources/evolution.F90
index 1bc0460..a7245fa 100644
--- a/sources/evolution.F90
+++ b/sources/evolution.F90
@@ -1135,12 +1135,6 @@ module evolution
 
     end if ! toplev > 1
 
-#ifdef DEBUG
-! check variables for NaNs
-!
-    call check_variables()
-#endif /* DEBUG */
-
 ! error entry point
 !
     100 continue
@@ -4319,73 +4313,6 @@ module evolution
 !-------------------------------------------------------------------------------
 !
   end subroutine update_errors_max
-#ifdef DEBUG
-!
-!===============================================================================
-!
-! subroutine CHECK_VARIABLES:
-! --------------------------
-!
-!   Subroutine iterates over all data blocks and converts the conservative
-!   variables to their primitive representation.
-!
-!
-!===============================================================================
-!
-  subroutine check_variables()
-
-    use coordinates    , only : nn => bcells
-    use equations      , only : nv, pvars, cvars
-    use ieee_arithmetic, only : ieee_is_nan
-
-    use blocks        , only : block_meta
-    use blocks        , only : block_data, list_data
-
-    implicit none
-
-    integer :: i, j, k, p
-
-    type(block_meta), pointer :: pmeta
-    type(block_data), pointer :: pdata
-
-!-------------------------------------------------------------------------------
-!
-#if NDIMS == 2
-    k = 1
-#endif /* NDIMS == 2 */
-
-    pdata => list_data
-
-    do while (associated(pdata))
-
-      pmeta => pdata%meta
-
-#if NDIMS == 3
-      do k = 1, nn
-#endif /* NDIMS == 3 */
-        do j = 1, nn
-          do i = 1, nn
-            do p = 1, nv
-              if (ieee_is_nan(pdata%u(p,i,j,k))) then
-                print *, 'U NaN:', cvars(p), pdata%meta%id, i, j, k
-              end if
-              if (ieee_is_nan(pdata%q(p,i,j,k))) then
-                print *, 'Q NaN:', pvars(p), pdata%meta%id, i, j, k
-              end if
-            end do ! p = 1, nv
-          end do ! i = 1, nn
-        end do ! j = 1, nn
-#if NDIMS == 3
-      end do ! k = 1, nn
-#endif /* NDIMS == 3 */
-
-      pdata => pdata%next
-    end do
-
-!-------------------------------------------------------------------------------
-!
-  end subroutine check_variables
-#endif /* DEBUG */
 
 !===============================================================================
 !
diff --git a/sources/forcing.F90 b/sources/forcing.F90
index e19b09b..e93e7c3 100644
--- a/sources/forcing.F90
+++ b/sources/forcing.F90
@@ -1212,33 +1212,28 @@ module forcing
 !
   subroutine inject_eddy(xp, ap)
 
-    use blocks, only : block_data, list_data
+    use blocks, only : block_data, data_blocks, get_dblocks
 
     implicit none
 
     real(kind=8), dimension(3), intent(in) :: xp, ap
 
+    integer :: m, n
+
     type(block_data), pointer :: pdata
 
 !-------------------------------------------------------------------------------
 !
-! assign pdata with the first block on the data block list
-!
-    pdata => list_data
+    n = get_dblocks()
 
-! iterate over all data blocks
-!
-    do while (associated(pdata))
+!$omp parallel do default(shared) private(pdata)
+    do m = 1, n
+      pdata => data_blocks(m)%ptr
 
-! inject eddy into the current block
-!
       call inject_eddy_block(pdata, xp(:), ap(:))
 
-! assign pdata to the next block
-!
-      pdata => pdata%next
-
-    end do ! over data blocks
+    end do
+!$omp end parallel do
 
 !-------------------------------------------------------------------------------
 !
@@ -1695,45 +1690,50 @@ module forcing
   subroutine get_vcoefs()
 
     use blocks  , only : block_data, data_blocks, get_dblocks
+    use helpers , only : print_message
 #ifdef MPI
     use mpitools, only : reduce_sum
 #endif /* MPI */
 
     implicit none
 
-    integer :: m, n, nt
+    integer :: m, n, status
 
     type(block_data), pointer :: pdata
 
     complex(kind=8), dimension(:,:,:), allocatable :: vc
 
-!$  integer :: omp_get_num_threads, omp_get_thread_num
+    character(len=*), parameter :: loc = 'FORCING:get_vcoefs()'
 
 !-------------------------------------------------------------------------------
 !
-    nt = 0
-!$omp parallel
-!$  nt = omp_get_num_threads() - 1
-!$omp end parallel
-    allocate(vc(nmodes, NDIMS, 0:nt))
+    n = get_dblocks()
+
+    allocate(vc(nmodes, NDIMS, n), stat=status)
+    if (status /= 0) then
+      call print_message(loc, &
+                      "Could not allocate memory for the Fourier coefficients!")
+      go to 100
+    end if
 
     vc(:,:,:) = cmplx(0.0d+00, 0.0d+00, kind=8)
 
-    n = get_dblocks()
-!$omp parallel default(shared) private(pdata,nt)
-!$  nt = omp_get_thread_num()
-!$omp do
+!$omp parallel do default(shared) private(pdata)
     do m = 1, n
       pdata => data_blocks(m)%ptr
 
-      call get_vcoefs_block(pdata, vc(:,:,nt))
+      call get_vcoefs_block(pdata, vc(:,:,m))
     end do
-!$omp end do
-!$omp end parallel
+!$omp end parallel do
 
     vcoefs = sum(vc, 3)
 
-    deallocate(vc)
+    deallocate(vc, stat=status)
+    if (status /= 0) &
+      call print_message(loc, &
+                        "Could not release memory of the Fourier coefficients!")
+
+    100 continue
 
 #ifdef MPI
     call reduce_sum(vcoefs)