From 1f4d0453a4d5c225be028a20ba8cfe44f3238336 Mon Sep 17 00:00:00 2001
From: huebleruwm <37674341+huebleruwm@users.noreply.github.com>
Date: Tue, 18 Jul 2023 20:27:07 -0500
Subject: [PATCH] Openacc tweaks and cleanup 2 (#1099)

* Making all end parallel directives specify end loop

* Replacing last acc declare in a procedure with acc enter/exit data commands

* Removing pure declarations, turns out they dont really improve performance and openmp isn't allowed within them.

* Adding explicit directives to penta_lu and tridiag_lu, as opposed to the previously used kernels directive. Also splitting up a loop to improve GPU performance.

* Splitting up a couple loops for performance reasons.

* Missed one somehow
---
 src/CLUBB_core/adg1_adg2_3d_luhar_pdf.F90     |   2 +-
 src/CLUBB_core/advance_clubb_core_module.F90  |   4 +-
 src/CLUBB_core/advance_helper_module.F90      |   2 +-
 .../advance_windm_edsclrm_module.F90          | 108 +++++++++---------
 src/CLUBB_core/advance_wp2_wp3_module.F90     |  57 +++++----
 src/CLUBB_core/advance_xm_wpxp_module.F90     |  12 +-
 src/CLUBB_core/advance_xp2_xpyp_module.F90    |   8 +-
 src/CLUBB_core/advance_xp3_module.F90         |   4 +-
 src/CLUBB_core/calc_roots.F90                 |   6 +-
 src/CLUBB_core/clip_explicit.F90              |  14 +--
 src/CLUBB_core/clip_semi_implicit.F90         |   2 +-
 src/CLUBB_core/clubb_api_module.F90           |   2 +-
 src/CLUBB_core/corr_varnce_module.F90         |   2 +-
 src/CLUBB_core/diffusion.F90                  |   6 +-
 src/CLUBB_core/fill_holes.F90                 |  18 ++-
 src/CLUBB_core/grid_class.F90                 |  18 +--
 src/CLUBB_core/mean_adv.F90                   |   4 +-
 src/CLUBB_core/mono_flux_limiter.F90          |   2 +-
 src/CLUBB_core/numerical_check.F90            |   2 +-
 src/CLUBB_core/penta_lu_solver.F90            | 102 ++++++++++-------
 src/CLUBB_core/tridiag_lu_solver.F90          |  48 +++++---
 src/CLUBB_core/turbulent_adv_pdf.F90          |   8 +-
 22 files changed, 245 insertions(+), 186 deletions(-)

diff --git a/src/CLUBB_core/adg1_adg2_3d_luhar_pdf.F90 b/src/CLUBB_core/adg1_adg2_3d_luhar_pdf.F90
index 09334a800..acf9fe634 100644
--- a/src/CLUBB_core/adg1_adg2_3d_luhar_pdf.F90
+++ b/src/CLUBB_core/adg1_adg2_3d_luhar_pdf.F90
@@ -1360,7 +1360,7 @@ subroutine backsolve_Luhar_params( Sk_max, Skx,          & ! In
   end subroutine backsolve_Luhar_params
 
   !=============================================================================
-  pure function max_cubic_root( a_coef, b_coef, c_coef, d_coef ) &
+  function max_cubic_root( a_coef, b_coef, c_coef, d_coef ) &
   result( max_root )
 
     ! Description:
diff --git a/src/CLUBB_core/advance_clubb_core_module.F90 b/src/CLUBB_core/advance_clubb_core_module.F90
index 925448535..4ac548303 100644
--- a/src/CLUBB_core/advance_clubb_core_module.F90
+++ b/src/CLUBB_core/advance_clubb_core_module.F90
@@ -5157,7 +5157,7 @@ subroutine calc_trapezoid_zm( nz, ngrdcol, gr, variable_zm, variable_zt, &
                                  * ( gr%zm(i,k) - gr%zt(i,k) ) * gr%invrs_dzm(i,k)
         end do
       end do 
-      !$acc end parallel
+      !$acc end parallel loop
 
       return
     end subroutine calc_trapezoid_zm
@@ -5514,7 +5514,7 @@ subroutine set_Lscale_max( ngrdcol, l_implemented, host_dx, host_dy, &
     end subroutine set_Lscale_max
 
 !===============================================================================
-  pure subroutine calculate_thlp2_rad &
+  subroutine calculate_thlp2_rad &
                   ( nz, rcm_zm, thlprcp, radht_zm, & ! Intent(in)
                     clubb_params,                  & ! Intent(in)
                     thlp2_forcing )                  ! Intent(inout)
diff --git a/src/CLUBB_core/advance_helper_module.F90 b/src/CLUBB_core/advance_helper_module.F90
index 19602e375..5fb4c88b2 100644
--- a/src/CLUBB_core/advance_helper_module.F90
+++ b/src/CLUBB_core/advance_helper_module.F90
@@ -2023,7 +2023,7 @@ function vertical_avg( total_idx, rho_ds, field, dz )
   end function vertical_avg
 
   !=============================================================================
-  pure function vertical_integral( total_idx, rho_ds, &
+  function vertical_integral( total_idx, rho_ds, &
                                        field, dz )
 
     ! Description:
diff --git a/src/CLUBB_core/advance_windm_edsclrm_module.F90 b/src/CLUBB_core/advance_windm_edsclrm_module.F90
index 37c639c24..3e502d531 100644
--- a/src/CLUBB_core/advance_windm_edsclrm_module.F90
+++ b/src/CLUBB_core/advance_windm_edsclrm_module.F90
@@ -294,7 +294,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
     do i = 1, ngrdcol
       nu_zero(i) = zero
     end do
-    !$acc end parallel
+    !$acc end parallel loop
     
     !$acc parallel loop gang vector collapse(2) default(present)
     do k = 1, nz
@@ -302,7 +302,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
         Km_zm_p_nu10(i,k) = Km_zm(i,k) + nu_vert_res_dep%nu10(i)
       end do
     end do
-    !$acc end parallel
+    !$acc end parallel loop
 
     l_perturbed_wind = ( .not. l_predict_upwp_vpwp ) .and. l_linearize_pbl_winds
     
@@ -320,7 +320,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           end do
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
     end if
 
     if ( .not. l_predict_upwp_vpwp ) then
@@ -333,7 +333,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           Km_zt(i,k) = max( Km_zt(i,k), zero )
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       ! Calculate diffusion terms
       call diffusion_zt_lhs( nz, ngrdcol, gr, Km_zm, Km_zt, nu_vert_res_dep%nu10, & ! In
@@ -357,7 +357,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           ! Thermodynamic subdiagonal: [ x xm(k-1,<t+1>) ]
           lhs_diff(km1_tdiag,i,2) = zero
         end do
-        !$acc end parallel
+        !$acc end parallel loop
       else
         !$acc parallel loop gang vector default(present)
         do i = 1, ngrdcol
@@ -374,7 +374,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           ! Thermodynamic subdiagonal: [ x xm(k-1,<t+1>) ]
           lhs_diff(km1_tdiag,i,2) = zero
         end do
-        !$acc end parallel
+        !$acc end parallel loop
       end if
 
       if ( l_lmm_stepping ) then
@@ -385,7 +385,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
             vm_old(i,k) = vm(i,k)
           end do
         end do
-        !$acc end parallel
+        !$acc end parallel loop
       end if ! l_lmm_stepping
 
       !----------------------------------------------------------------
@@ -419,14 +419,14 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           wind_speed(i,k) = max( sqrt( um(i,k)**2 + vm(i,k)**2 ), eps )
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       ! Compute u_star_sqd according to the definition of u_star.
       !$acc parallel loop gang vector default(present)
       do i = 1, ngrdcol
         u_star_sqd(i) = sqrt( upwp(i,1)**2 + vpwp(i,1)**2 )
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       ! Compute the explicit portion of the um equation.
       ! Build the right-hand side vector.
@@ -464,7 +464,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           upwp(i,k) = -one_half * xpwp(i,k)
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
       
       call calc_xpwp( nz, ngrdcol, gr, &
                       Km_zm_p_nu10, vm, &
@@ -476,7 +476,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           vpwp(i,k) = -one_half * xpwp(i,k)
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       ! A zero-flux boundary condition at the top of the model, d(xm)/dz = 0,
       ! means that x'w' at the top model level is 0,
@@ -486,7 +486,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
         upwp(i,nz) = zero
         vpwp(i,nz) = zero
       end do
-      !$acc end parallel
+      !$acc end parallel loop
       
       ! Compute the implicit portion of the um and vm equations.
       ! Build the left-hand side matrix.
@@ -522,7 +522,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           um(i,k) = solution(i,k,windm_edsclrm_um)
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       !----------------------------------------------------------------
       ! Update meridional (south-to-north) component of mean wind, vm
@@ -533,7 +533,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           vm(i,k) = solution(i,k,windm_edsclrm_vm)
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       if ( l_stats_samp ) then
 
@@ -571,7 +571,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
         um(i,1) = um(i,2)
         vm(i,1) = vm(i,2)
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       if ( l_lmm_stepping ) then
         !$acc parallel loop gang vector collapse(2) default(present)
@@ -581,7 +581,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
             vm(i,k) = one_half * ( vm_old(i,k) + vm(i,k) )
           end do
         end do
-        !$acc end parallel
+        !$acc end parallel loop
       endif ! l_lmm_stepping ) then
 
       if ( uv_sponge_damp_settings%l_sponge_damping ) then
@@ -645,7 +645,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           upwp(i,k) = upwp(i,k) - one_half * xpwp(i,k)
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
                                                       
       call calc_xpwp( nz, ngrdcol, gr, &
                       Km_zm_p_nu10, vm, &
@@ -657,7 +657,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           vpwp(i,k) = vpwp(i,k) - one_half * xpwp(i,k)
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       ! Adjust um and vm if nudging is turned on.
       if ( l_uv_nudge ) then
@@ -680,7 +680,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
             vm(i,k) = vm(i,k) - ( ( vm(i,k) - vm_ref(i,k) ) * (dt/ts_nudge) )
           end do
         end do
-        !$acc end parallel
+        !$acc end parallel loop
 
         if ( l_stats_samp ) then
           !$acc update host( um, vm )
@@ -800,14 +800,14 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           wind_speed_pert(i,k) = max( sqrt( (um_pert(i,k))**2 + (vm_pert(i,k))**2 ), eps )
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       ! Compute u_star_sqd according to the definition of u_star.
       !$acc parallel loop gang vector default(present)
       do i = 1, ngrdcol
         u_star_sqd_pert(i) = sqrt( upwp_pert(i,1)**2 + vpwp_pert(i,1)**2 )
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       ! Compute the explicit portion of the um equation.
       ! Build the right-hand side vector.
@@ -845,7 +845,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           upwp_pert(i,k) = -one_half * xpwp(i,k)
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
                                                   
       call calc_xpwp( nz, ngrdcol, gr, &
                       Km_zm_p_nu10, vm_pert, &
@@ -857,7 +857,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           vpwp_pert(i,k) = -one_half * xpwp(i,k)
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
       
       ! A zero-flux boundary condition at the top of the model, d(xm)/dz = 0,
       ! means that x'w' at the top model level is 0,
@@ -867,7 +867,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
         upwp_pert(i,nz) = zero
         vpwp_pert(i,nz) = zero
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       ! Compute the implicit portion of the um and vm equations.
       ! Build the left-hand side matrix.
@@ -903,7 +903,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           um_pert(i,k) = solution(i,k,windm_edsclrm_um)
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       !----------------------------------------------------------------
       ! Update meridional (south-to-north) component of mean wind, vm
@@ -914,7 +914,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           vm_pert(i,k) = solution(i,k,windm_edsclrm_vm)
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
   
       ! The values of um(1) and vm(1) are located below the model surface and
       ! do not affect the rest of the model.  The values of um(1) or vm(1) are
@@ -927,7 +927,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
         um_pert(i,1) = um_pert(i,2)
         vm_pert(i,1) = vm_pert(i,2)
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       ! Second part of momentum (implicit component)
 
@@ -943,7 +943,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           upwp_pert(i,k) = upwp_pert(i,k) - one_half * xpwp(i,k)
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
       
       call calc_xpwp( nz, ngrdcol, gr, &
                       Km_zm_p_nu10, vm_pert, &
@@ -955,7 +955,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           vpwp_pert(i,k) = vpwp_pert(i,k) - one_half * xpwp(i,k)
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
       
       if ( l_tke_aniso ) then
 
@@ -1034,7 +1034,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           Kmh_zt(i,k) = max( Kmh_zt(i,k), zero )
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       ! Calculate diffusion terms
       call diffusion_zt_lhs( nz, ngrdcol, gr, Kmh_zm, Kmh_zt, nu_zero,  & ! intent(in)
@@ -1059,7 +1059,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           ! Thermodynamic subdiagonal: [ x xm(k-1,<t+1>) ]
           lhs_diff(km1_tdiag,i,2) = zero
         end do
-        !$acc end parallel
+        !$acc end parallel loop
 
       else
 
@@ -1076,7 +1076,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           ! Thermodynamic subdiagonal: [ x xm(k-1,<t+1>) ]
           lhs_diff(km1_tdiag,i,2) = zero
         end do
-        !$acc end parallel
+        !$acc end parallel loop
 
       end if
 
@@ -1089,7 +1089,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
             end do
           end do
         end do
-        !$acc end parallel
+        !$acc end parallel loop
       endif ! l_lmm_stepping
 
       ! Eddy-scalar surface fluxes, x'w'|_sfc, are applied through an explicit
@@ -1130,7 +1130,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
             wpedsclrp(i,k,edsclr) = -one_half * xpwp(i,k)
           end do
         end do
-        !$acc end parallel
+        !$acc end parallel loop
       end do
 
       ! A zero-flux boundary condition at the top of the model, d(xm)/dz = 0,
@@ -1142,7 +1142,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           wpedsclrp(i,nz,j) = zero
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       ! Compute the implicit portion of the xm (eddy-scalar) equations.
       ! Build the left-hand side matrix.
@@ -1177,7 +1177,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           end do
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       ! The value of edsclrm(1) is located below the model surface and does not
       ! effect the rest of the model.  The value of edsclrm(1) is simply set to
@@ -1188,7 +1188,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
           edsclrm(i,1,edsclr) = edsclrm(i,2,edsclr)
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       if ( l_lmm_stepping ) then
         !$acc parallel loop gang vector collapse(3) default(present)
@@ -1199,7 +1199,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
             end do
           end do
         end do
-        !$acc end parallel
+        !$acc end parallel loop
       endif ! l_lmm_stepping
 
       ! Second part of momentum (implicit component)
@@ -1218,7 +1218,7 @@ subroutine advance_windm_edsclrm( nz, ngrdcol, gr, dt, &
             wpedsclrp(i,k,edsclr) = -one_half * xpwp(i,k)
           end do
         end do
-        !$acc end parallel
+        !$acc end parallel loop
       end do
 
       ! Note that the w'edsclr' terms are not clipped, since we don't compute
@@ -2153,7 +2153,7 @@ subroutine compute_uv_tndcy( nz, ngrdcol, solve_type, &
             xm_gf(i,k) = - fcor(i) * perp_wind_g(i,k)
           end do
         end do
-        !$acc end parallel
+        !$acc end parallel loop
         
         !$acc parallel loop gang vector collapse(2) default(present)
         do k = 1, nz
@@ -2161,7 +2161,7 @@ subroutine compute_uv_tndcy( nz, ngrdcol, solve_type, &
             xm_cf(i,k) = fcor(i) * perp_wind_m(i,k)
           end do
         end do
-        !$acc end parallel
+        !$acc end parallel loop
 
       case ( windm_edsclrm_vm )
 
@@ -2175,7 +2175,7 @@ subroutine compute_uv_tndcy( nz, ngrdcol, solve_type, &
             xm_gf(i,k) = fcor(i) * perp_wind_g(i,k)
           end do
         end do
-        !$acc end parallel
+        !$acc end parallel loop
 
         !$acc parallel loop gang vector collapse(2) default(present)
         do k = 1, nz
@@ -2183,7 +2183,7 @@ subroutine compute_uv_tndcy( nz, ngrdcol, solve_type, &
             xm_cf(i,k) = -fcor(i) * perp_wind_m(i,k)
           end do
         end do
-        !$acc end parallel
+        !$acc end parallel loop
 
       case default
 
@@ -2198,7 +2198,7 @@ subroutine compute_uv_tndcy( nz, ngrdcol, solve_type, &
             xm_cf(i,k) = 0._core_rknd
           end do
         end do
-        !$acc end parallel
+        !$acc end parallel loop
 
       end select
 
@@ -2208,7 +2208,7 @@ subroutine compute_uv_tndcy( nz, ngrdcol, solve_type, &
           xm_tndcy(i,k) = xm_gf(i,k) + xm_cf(i,k) + xm_forcing(i,k)
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       if ( l_stats_samp ) then
 
@@ -2237,7 +2237,7 @@ subroutine compute_uv_tndcy( nz, ngrdcol, solve_type, &
           xm_tndcy(i,k) = 0.0_core_rknd
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
     endif
 
@@ -2329,7 +2329,7 @@ subroutine windm_edsclrm_lhs( nz, ngrdcol, gr, dt, &
       lhs(2,i,1) = 1.0_core_rknd
       lhs(3,i,1) = 0.0_core_rknd
     end do
-    !$acc end parallel
+    !$acc end parallel loop
 
     ! Add terms to lhs
     !$acc parallel loop gang vector collapse(2) default(present)
@@ -2345,7 +2345,7 @@ subroutine windm_edsclrm_lhs( nz, ngrdcol, gr, dt, &
         lhs(3,i,k) = 0.5_core_rknd * lhs_diff(3,i,k)
       end do
     end do
-    !$acc end parallel 
+    !$acc end parallel loop
 
     ! LHS mean advection term.
     if ( .not. l_implemented ) then
@@ -2355,7 +2355,7 @@ subroutine windm_edsclrm_lhs( nz, ngrdcol, gr, dt, &
           lhs(1:3,i,k) = lhs(1:3,i,k) + lhs_ma_zt(:,i,k)
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
     endif
 
@@ -2367,7 +2367,7 @@ subroutine windm_edsclrm_lhs( nz, ngrdcol, gr, dt, &
         lhs(2,i,2) = lhs(2,i,2) + invrs_rho_ds_zt(i,2) * gr%invrs_dzt(i,2) &
                                   * rho_ds_zm(i,1) * ( u_star_sqd(i) / wind_speed(i,2) )
       end do
-      !$acc end parallel
+      !$acc end parallel loop
     end if ! l_imp_sfc_momentum_flux
 
     return
@@ -2488,7 +2488,7 @@ subroutine windm_edsclrm_rhs( nz, ngrdcol, gr, solve_type, dt, &
     do i = 1, ngrdcol
       rhs(i,1) = 0.0_core_rknd
     end do
-    !$acc end parallel
+    !$acc end parallel loop
 
     ! Non-boundary rhs calculation, this is a highly vectorized loop
     !$acc parallel loop gang vector collapse(2) default(present)
@@ -2502,7 +2502,7 @@ subroutine windm_edsclrm_rhs( nz, ngrdcol, gr, solve_type, dt, &
                    + invrs_dt * xm(i,k)                     ! RHS time tendency
       end do
     end do
-    !$acc end parallel
+    !$acc end parallel loop
 
     ! Upper boundary calculation
     !$acc parallel loop gang vector default(present)
@@ -2513,7 +2513,7 @@ subroutine windm_edsclrm_rhs( nz, ngrdcol, gr, solve_type, dt, &
                   + xm_tndcy(i,nz)                     & ! RHS forcings
                   + invrs_dt * xm(i,nz)                  ! RHS time tendency
     end do
-    !$acc end parallel
+    !$acc end parallel loop
 
     if ( l_stats_samp .and. ixm_ta > 0 ) then
 
@@ -2555,7 +2555,7 @@ subroutine windm_edsclrm_rhs( nz, ngrdcol, gr, solve_type, dt, &
                               * gr%invrs_dzt(i,2)  &
                               * rho_ds_zm(i,1) * xpwp_sfc(i)
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       if ( l_stats_samp .and. ixm_ta > 0 ) then
 
diff --git a/src/CLUBB_core/advance_wp2_wp3_module.F90 b/src/CLUBB_core/advance_wp2_wp3_module.F90
index a79562f8f..921e4309e 100644
--- a/src/CLUBB_core/advance_wp2_wp3_module.F90
+++ b/src/CLUBB_core/advance_wp2_wp3_module.F90
@@ -1985,7 +1985,6 @@ subroutine wp23_lhs( nz, ngrdcol, gr, dt, &
       do i = 1, ngrdcol
 
         k_wp3 = 2*k - 1
-        k_wp2 = 2*k
 
         ! ------ w'3 ------
 
@@ -2018,8 +2017,16 @@ subroutine wp23_lhs( nz, ngrdcol, gr, dt, &
 
         ! LHS mean advection (ma) and diffusion (diff) terms
         lhs(5,i,k_wp3) = lhs(5,i,k_wp3) + lhs_ma_zt(3,i,k) + lhs_diff_zt(3,i,k)
+      end do
+    end do
+    !$acc end parallel loop
 
+    !$acc parallel loop gang vector collapse(2) default(present)
+    do k = 2, nz-1, 1
+      do i = 1, ngrdcol
 
+        k_wp2 = 2*k
+        
         ! ------ w'2 ------
 
         ! LHS mean advection (ma) and diffusion (diff) terms
@@ -2473,7 +2480,6 @@ subroutine wp23_rhs( nz, ngrdcol, gr, dt, &
       do i = 1, ngrdcol
         
         k_wp3 = 2*k - 1
-        k_wp2 = 2*k
 
         ! ------ Combine terms for 3rd moment of vertical velocity, <w'^3> ------ !
 
@@ -2494,8 +2500,17 @@ subroutine wp23_rhs( nz, ngrdcol, gr, dt, &
         ! RHS "over implicit" pressure term 1 (pr1).
         rhs(i,k_wp3)  = rhs(i,k_wp3) + ( one - gamma_over_implicit_ts ) &
                                        * ( - lhs_pr1_wp3(i,k) * wp3(i,k) )
+      end do
+    end do
+    !$acc end parallel loop
 
 
+    !$acc parallel loop gang vector collapse(2) default(present)
+    do k = 2, nz-1
+      do i = 1, ngrdcol
+
+        k_wp2 = 2*k
+
         ! ------ Combine terms for 2nd moment of vertical velocity, <w'^2> ------ !
 
         ! RHS time tendency.
@@ -2906,7 +2921,7 @@ subroutine wp23_rhs( nz, ngrdcol, gr, dt, &
   end subroutine wp23_rhs
 
   !=============================================================================
-  pure subroutine wp2_term_ta_lhs( nz, ngrdcol, gr, &
+  subroutine wp2_term_ta_lhs( nz, ngrdcol, gr, &
                                    rho_ds_zt, invrs_rho_ds_zm, &
                                    lhs_ta_wp2 )
 
@@ -3027,7 +3042,7 @@ pure subroutine wp2_term_ta_lhs( nz, ngrdcol, gr, &
   end subroutine wp2_term_ta_lhs
 
   !=============================================================================
-  pure subroutine wp2_terms_ac_pr2_lhs( nz, ngrdcol, gr, C_uu_shr, wm_zt, &
+  subroutine wp2_terms_ac_pr2_lhs( nz, ngrdcol, gr, C_uu_shr, wm_zt, &
                                         lhs_ac_pr2_wp2 )
 
     ! Description:
@@ -3147,7 +3162,7 @@ pure subroutine wp2_terms_ac_pr2_lhs( nz, ngrdcol, gr, C_uu_shr, wm_zt, &
   end subroutine wp2_terms_ac_pr2_lhs
 
   !=============================================================================
-  pure subroutine wp2_term_dp1_lhs( nz, ngrdcol, &
+  subroutine wp2_term_dp1_lhs( nz, ngrdcol, &
                                     C1_Skw_fnc, invrs_tau1m, &
                                     lhs_dp1_wp2 )
 
@@ -3241,7 +3256,7 @@ pure subroutine wp2_term_dp1_lhs( nz, ngrdcol, &
   end subroutine wp2_term_dp1_lhs
 
   !=============================================================================
-  pure subroutine wp2_term_pr1_lhs( nz, ngrdcol, C4, invrs_tau_C4_zm, &
+  subroutine wp2_term_pr1_lhs( nz, ngrdcol, C4, invrs_tau_C4_zm, &
                                     lhs_pr1_wp2 )
 
     ! Description
@@ -3342,7 +3357,7 @@ pure subroutine wp2_term_pr1_lhs( nz, ngrdcol, C4, invrs_tau_C4_zm, &
   end subroutine wp2_term_pr1_lhs
 
   !=============================================================================
-  pure subroutine wp2_terms_bp_pr2_rhs( nz, ngrdcol, C_uu_buoy, &
+  subroutine wp2_terms_bp_pr2_rhs( nz, ngrdcol, C_uu_buoy, &
                                         thv_ds_zm, wpthvp, &
                                         rhs_bp_pr2_wp2 )
 
@@ -3435,7 +3450,7 @@ pure subroutine wp2_terms_bp_pr2_rhs( nz, ngrdcol, C_uu_buoy, &
   end subroutine wp2_terms_bp_pr2_rhs
 
   !=============================================================================
-  pure subroutine wp2_term_dp1_rhs( nz, ngrdcol, C1_Skw_fnc, &
+  subroutine wp2_term_dp1_rhs( nz, ngrdcol, C1_Skw_fnc, &
                                     invrs_tau1m, threshold, up2, vp2, &
                                     l_damp_wp2_using_em, &
                                     rhs_dp1_wp2 )
@@ -3547,7 +3562,7 @@ pure subroutine wp2_term_dp1_rhs( nz, ngrdcol, C1_Skw_fnc, &
   end subroutine wp2_term_dp1_rhs
 
   !=============================================================================
-  pure subroutine wp2_term_pr3_rhs( nz, ngrdcol, gr, C_uu_shr, C_uu_buoy, &
+  subroutine wp2_term_pr3_rhs( nz, ngrdcol, gr, C_uu_shr, C_uu_buoy, &
                                     thv_ds_zm, wpthvp, upwp, &
                                     um, vpwp, vm, &
                                     rhs_pr3_wp2 )
@@ -3685,7 +3700,7 @@ pure subroutine wp2_term_pr3_rhs( nz, ngrdcol, gr, C_uu_shr, C_uu_buoy, &
   end subroutine wp2_term_pr3_rhs
 
   !=============================================================================
-  pure subroutine wp2_term_pr1_rhs( nz, ngrdcol, C4, &
+  subroutine wp2_term_pr1_rhs( nz, ngrdcol, C4, &
                                     up2, vp2, invrs_tau_C4_zm, &
                                     rhs_pr1_wp2 )
 
@@ -3778,7 +3793,7 @@ pure subroutine wp2_term_pr1_rhs( nz, ngrdcol, C4, &
   end subroutine wp2_term_pr1_rhs
 
   !=============================================================================
-  pure subroutine wp2_term_pr_dfsn_rhs( nz, ngrdcol, gr, C_wp2_pr_dfsn, &
+  subroutine wp2_term_pr_dfsn_rhs( nz, ngrdcol, gr, C_wp2_pr_dfsn, &
                                         rho_ds_zt, invrs_rho_ds_zm, &
                                         wpup2, wpvp2, wp3, &
                                         rhs_pr_dfsn_wp2 )
@@ -3887,7 +3902,7 @@ pure subroutine wp2_term_pr_dfsn_rhs( nz, ngrdcol, gr, C_wp2_pr_dfsn, &
   end subroutine wp2_term_pr_dfsn_rhs
 
   !=============================================================================
-  pure subroutine wp3_term_ta_new_pdf_lhs( nz, ngrdcol, gr, coef_wp4_implicit, &
+  subroutine wp3_term_ta_new_pdf_lhs( nz, ngrdcol, gr, coef_wp4_implicit, &
                                            wp2, rho_ds_zm, invrs_rho_ds_zt, &
                                            lhs_ta_wp3 )
 
@@ -4047,7 +4062,7 @@ pure subroutine wp3_term_ta_new_pdf_lhs( nz, ngrdcol, gr, coef_wp4_implicit, &
   end subroutine wp3_term_ta_new_pdf_lhs
 
   !=============================================================================
-  pure subroutine wp3_term_ta_ADG1_lhs( nz, ngrdcol, gr, &
+  subroutine wp3_term_ta_ADG1_lhs( nz, ngrdcol, gr, &
                                         wp2, a1, a1_zt, a3, a3_zt, &
                                         wp3_on_wp2, rho_ds_zm, &
                                         rho_ds_zt, invrs_rho_ds_zt, &
@@ -4384,7 +4399,7 @@ pure subroutine wp3_term_ta_ADG1_lhs( nz, ngrdcol, gr, &
   end subroutine wp3_term_ta_ADG1_lhs
 
   !=============================================================================
-  pure subroutine wp3_term_tp_lhs( nz, ngrdcol, gr, coef_wp3_tp, & 
+  subroutine wp3_term_tp_lhs( nz, ngrdcol, gr, coef_wp3_tp, & 
                                    wp2, rho_ds_zm, invrs_rho_ds_zt, &
                                    lhs_tp_wp3 )
 
@@ -4536,7 +4551,7 @@ pure subroutine wp3_term_tp_lhs( nz, ngrdcol, gr, coef_wp3_tp, &
   end subroutine wp3_term_tp_lhs
 
   !=============================================================================
-  pure subroutine wp3_terms_ac_pr2_lhs( nz, ngrdcol, gr, C11_Skw_fnc, wm_zm, &
+  subroutine wp3_terms_ac_pr2_lhs( nz, ngrdcol, gr, C11_Skw_fnc, wm_zm, &
                                         lhs_ac_pr2_wp3 )
 
     ! Description:
@@ -4655,7 +4670,7 @@ pure subroutine wp3_terms_ac_pr2_lhs( nz, ngrdcol, gr, C11_Skw_fnc, wm_zm, &
   end subroutine wp3_terms_ac_pr2_lhs
 
   !=============================================================================
-  pure subroutine wp3_term_pr1_lhs( nz, ngrdcol, C8, C8b, &
+  subroutine wp3_term_pr1_lhs( nz, ngrdcol, C8, C8b, &
                                     invrs_tau_wp3_zt, Skw_zt, &
                                     l_damp_wp3_Skw_squared, &
                                     lhs_pr1_wp3 )
@@ -4787,7 +4802,7 @@ pure subroutine wp3_term_pr1_lhs( nz, ngrdcol, C8, C8b, &
   end subroutine wp3_term_pr1_lhs
 
   !=============================================================================
-  pure subroutine wp3_term_ta_explicit_rhs( nz, ngrdcol, gr, &
+  subroutine wp3_term_ta_explicit_rhs( nz, ngrdcol, gr, &
                                             wp4, rho_ds_zm, invrs_rho_ds_zt, &
                                             rhs_ta_wp3 )
 
@@ -4898,7 +4913,7 @@ pure subroutine wp3_term_ta_explicit_rhs( nz, ngrdcol, gr, &
   end subroutine wp3_term_ta_explicit_rhs
 
   !=============================================================================
-  pure subroutine wp3_terms_bp1_pr2_rhs( nz, ngrdcol, C11_Skw_fnc, &
+  subroutine wp3_terms_bp1_pr2_rhs( nz, ngrdcol, C11_Skw_fnc, &
                                          thv_ds_zt, wp2thvp, &
                                          rhs_bp1_pr2_wp3 )
 
@@ -4986,7 +5001,7 @@ pure subroutine wp3_terms_bp1_pr2_rhs( nz, ngrdcol, C11_Skw_fnc, &
   end subroutine wp3_terms_bp1_pr2_rhs
 
   !=============================================================================
-  pure subroutine wp3_term_pr_turb_rhs( nz, ngrdcol, gr, C_wp3_pr_turb, Kh_zt, wpthvp, &
+  subroutine wp3_term_pr_turb_rhs( nz, ngrdcol, gr, C_wp3_pr_turb, Kh_zt, wpthvp, &
                                         dum_dz, dvm_dz, &
                                         upwp, vpwp, &
                                         thv_ds_zt, &
@@ -5106,7 +5121,7 @@ pure subroutine wp3_term_pr_turb_rhs( nz, ngrdcol, gr, C_wp3_pr_turb, Kh_zt, wpt
   end subroutine wp3_term_pr_turb_rhs
 
   !=============================================================================
-  pure subroutine wp3_term_pr_dfsn_rhs( nz, ngrdcol, gr, C_wp3_pr_dfsn, &
+  subroutine wp3_term_pr_dfsn_rhs( nz, ngrdcol, gr, C_wp3_pr_dfsn, &
                                         rho_ds_zm, invrs_rho_ds_zt, &
                                         wp2up2, wp2vp2, wp4, &
                                         up2, vp2, wp2, &
@@ -5220,7 +5235,7 @@ pure subroutine wp3_term_pr_dfsn_rhs( nz, ngrdcol, gr, C_wp3_pr_dfsn, &
   end subroutine wp3_term_pr_dfsn_rhs
 
   !=============================================================================
-  pure subroutine wp3_term_pr1_rhs( nz, ngrdcol, gr, C8, C8b, &
+  subroutine wp3_term_pr1_rhs( nz, ngrdcol, gr, C8, C8b, &
                                     invrs_tau_wp3_zt, Skw_zt, wp3, &
                                     l_damp_wp3_Skw_squared, &
                                     rhs_pr1_wp3 )
diff --git a/src/CLUBB_core/advance_xm_wpxp_module.F90 b/src/CLUBB_core/advance_xm_wpxp_module.F90
index 78648bbad..ca7d671c7 100644
--- a/src/CLUBB_core/advance_xm_wpxp_module.F90
+++ b/src/CLUBB_core/advance_xm_wpxp_module.F90
@@ -5052,7 +5052,7 @@ subroutine xm_wpxp_clipping_and_stats( &
   end subroutine xm_wpxp_clipping_and_stats
 
   !=============================================================================
-  pure subroutine xm_term_ta_lhs( nz, ngrdcol, gr, &
+  subroutine xm_term_ta_lhs( nz, ngrdcol, gr, &
                                   rho_ds_zm, invrs_rho_ds_zt, &
                                   lhs_ta_xm )
 
@@ -5164,7 +5164,7 @@ pure subroutine xm_term_ta_lhs( nz, ngrdcol, gr, &
   end subroutine xm_term_ta_lhs
 
   !=============================================================================
-  pure subroutine wpxp_term_tp_lhs( nz, ngrdcol, gr, wp2, & 
+  subroutine wpxp_term_tp_lhs( nz, ngrdcol, gr, wp2, & 
                                     lhs_tp )
 
     ! Description:
@@ -5277,7 +5277,7 @@ pure subroutine wpxp_term_tp_lhs( nz, ngrdcol, gr, wp2, &
   end subroutine wpxp_term_tp_lhs
     
   !=============================================================================
-  pure subroutine wpxp_terms_ac_pr2_lhs( nz, ngrdcol, C7_Skw_fnc, &
+  subroutine wpxp_terms_ac_pr2_lhs( nz, ngrdcol, C7_Skw_fnc, &
                                          wm_zt, invrs_dzm, &
                                          lhs_ac_pr2  ) 
 
@@ -5392,7 +5392,7 @@ pure subroutine wpxp_terms_ac_pr2_lhs( nz, ngrdcol, C7_Skw_fnc, &
   end subroutine wpxp_terms_ac_pr2_lhs
 
   !=============================================================================
-  pure subroutine wpxp_term_pr1_lhs( nz, ngrdcol, C6rt_Skw_fnc, C6thl_Skw_fnc, C7_Skw_fnc, &
+  subroutine wpxp_term_pr1_lhs( nz, ngrdcol, C6rt_Skw_fnc, C6thl_Skw_fnc, C7_Skw_fnc, &
                                      invrs_tau_C6_zm, l_scalar_calc, &
                                      lhs_pr1_wprtp, lhs_pr1_wpthlp, &
                                      lhs_pr1_wpsclrp )
@@ -5522,7 +5522,7 @@ pure subroutine wpxp_term_pr1_lhs( nz, ngrdcol, C6rt_Skw_fnc, C6thl_Skw_fnc, C7_
   end subroutine wpxp_term_pr1_lhs
 
   !=============================================================================
-  pure subroutine wpxp_terms_bp_pr3_rhs( nz, ngrdcol, C7_Skw_fnc, thv_ds_zm, xpthvp, &
+  subroutine wpxp_terms_bp_pr3_rhs( nz, ngrdcol, C7_Skw_fnc, thv_ds_zm, xpthvp, &
                                          rhs_bp_pr3 )
 
     ! Description:
@@ -5902,7 +5902,7 @@ end subroutine damp_coefficient
   !-----------------------------------------------------------------------
 
   !=====================================================================================
-  pure subroutine diagnose_upxp( nz, ngrdcol, gr, ypwp, xm, wpxp, ym, &
+  subroutine diagnose_upxp( nz, ngrdcol, gr, ypwp, xm, wpxp, ym, &
                                  C6x_Skw_fnc, tau_C6_zm, C7_Skw_fnc, &
                                  ypxp )
     ! Description:
diff --git a/src/CLUBB_core/advance_xp2_xpyp_module.F90 b/src/CLUBB_core/advance_xp2_xpyp_module.F90
index 3c91399a2..b6a09be50 100644
--- a/src/CLUBB_core/advance_xp2_xpyp_module.F90
+++ b/src/CLUBB_core/advance_xp2_xpyp_module.F90
@@ -5193,7 +5193,7 @@ subroutine calc_xp2_xpyp_ta_terms( nz, ngrdcol, gr, wprtp, wprtp2, wpthlp, wpthl
   end subroutine calc_xp2_xpyp_ta_terms
 
   !=============================================================================
-  pure function term_tp( xamp1, xam, xbmp1, xbm,  & 
+  function term_tp( xamp1, xam, xbmp1, xbm,  & 
                          wpxbp, wpxap, invrs_dzm ) & 
   result( rhs )
   !$acc routine seq
@@ -5255,7 +5255,7 @@ pure function term_tp( xamp1, xam, xbmp1, xbm,  &
   end function term_tp
 
   !=============================================================================
-  pure function term_dp1_lhs( Cn, invrs_tau_zm )  & 
+  function term_dp1_lhs( Cn, invrs_tau_zm )  & 
   result( lhs )
   !$acc routine seq
 
@@ -5327,7 +5327,7 @@ pure function term_dp1_lhs( Cn, invrs_tau_zm )  &
   end function term_dp1_lhs
 
   !=============================================================================
-  pure function term_dp1_rhs( Cn, invrs_tau_zm, threshold ) &
+  function term_dp1_rhs( Cn, invrs_tau_zm, threshold ) &
   result( rhs )
   !$acc routine seq
 
@@ -5390,7 +5390,7 @@ pure function term_dp1_rhs( Cn, invrs_tau_zm, threshold ) &
   end function term_dp1_rhs
 
   !=============================================================================
-  pure function term_pr1( C4, C14, xbp2, wp2, invrs_tau_C4_zm, invrs_tau_C14_zm ) & 
+  function term_pr1( C4, C14, xbp2, wp2, invrs_tau_C4_zm, invrs_tau_C14_zm ) & 
   result( rhs )
   !$acc routine seq
 
diff --git a/src/CLUBB_core/advance_xp3_module.F90 b/src/CLUBB_core/advance_xp3_module.F90
index 1ab734e79..a239c3391 100644
--- a/src/CLUBB_core/advance_xp3_module.F90
+++ b/src/CLUBB_core/advance_xp3_module.F90
@@ -464,7 +464,7 @@ subroutine advance_xp3_simplified( nz, ngrdcol, gr, solve_type, dt, & ! Intent(i
   end subroutine advance_xp3_simplified
 
   !=============================================================================
-  pure function term_tp_rhs( xp2_zt, wpxp, wpxpm1, &
+  function term_tp_rhs( xp2_zt, wpxp, wpxpm1, &
                              rho_ds_zm, rho_ds_zmm1, &
                              invrs_rho_ds_zt, &
                              invrs_dzt ) &
@@ -541,7 +541,7 @@ pure function term_tp_rhs( xp2_zt, wpxp, wpxpm1, &
   end function term_tp_rhs
 
   !=============================================================================
-  pure function term_ac_rhs( xm_zm, xm_zmm1, wpxp2, &
+  function term_ac_rhs( xm_zm, xm_zmm1, wpxp2, &
                              invrs_dzt ) &
   result( term_ac )
 
diff --git a/src/CLUBB_core/calc_roots.F90 b/src/CLUBB_core/calc_roots.F90
index 9c2f663b4..15eff0c5f 100644
--- a/src/CLUBB_core/calc_roots.F90
+++ b/src/CLUBB_core/calc_roots.F90
@@ -14,7 +14,7 @@ module calc_roots
   contains
 
   !=============================================================================
-  pure function cubic_solve( nz, a_coef, b_coef, c_coef, d_coef ) &
+  function cubic_solve( nz, a_coef, b_coef, c_coef, d_coef ) &
   result( roots )
 
     ! Description:
@@ -172,7 +172,7 @@ pure function cubic_solve( nz, a_coef, b_coef, c_coef, d_coef ) &
   end function cubic_solve
 
   !=============================================================================
-  pure function quadratic_solve( nz, a_coef, b_coef, c_coef ) &
+  function quadratic_solve( nz, a_coef, b_coef, c_coef ) &
   result( roots )
 
     ! Description:
@@ -261,7 +261,7 @@ pure function quadratic_solve( nz, a_coef, b_coef, c_coef ) &
   end function quadratic_solve
 
   !=============================================================================
-  pure function cube_root( x )
+  function cube_root( x )
 
     ! Description:
     ! Calculates the cube root of x.
diff --git a/src/CLUBB_core/clip_explicit.F90 b/src/CLUBB_core/clip_explicit.F90
index 40825cad1..22e3842a1 100644
--- a/src/CLUBB_core/clip_explicit.F90
+++ b/src/CLUBB_core/clip_explicit.F90
@@ -970,7 +970,7 @@ subroutine clip_variance( nz, ngrdcol, gr, solve_type, dt, threshold, &
         end if
       end do
     end do
-    !$acc end parallel
+    !$acc end parallel loop
 
     if ( l_stats_samp ) then
       !$acc update host( xp2 )
@@ -1205,7 +1205,7 @@ subroutine clip_skewness_core( nz, ngrdcol, gr, sfc_elevation, &
         wp2_zt_cubed(i,k) = wp2_zt(i,k)**3
       end do
     end do
-    !$acc end parallel
+    !$acc end parallel loop
 
     if ( l_use_wp3_lim_with_smth_Heaviside ) then 
 
@@ -1217,7 +1217,7 @@ subroutine clip_skewness_core( nz, ngrdcol, gr, sfc_elevation, &
           zagl_thresh(i,k) = zagl_thresh(i,k)  - 1.0_core_rknd 
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
       H_zagl(:,:) = smooth_heaviside_peskin(nz, ngrdcol, zagl_thresh(:,:), 0.6_core_rknd) 
 
@@ -1230,7 +1230,7 @@ subroutine clip_skewness_core( nz, ngrdcol, gr, sfc_elevation, &
                                      * 0.0021_core_rknd *Skw_max_mag**2 )
         end do
       end do
-     !$acc end parallel
+     !$acc end parallel loop
 
     else ! default method 
 
@@ -1248,7 +1248,7 @@ subroutine clip_skewness_core( nz, ngrdcol, gr, sfc_elevation, &
           endif
         end do
       end do
-      !$acc end parallel
+      !$acc end parallel loop
 
     end if
   
@@ -1263,7 +1263,7 @@ subroutine clip_skewness_core( nz, ngrdcol, gr, sfc_elevation, &
         end if
       end do
     end do
-    !$acc end parallel
+    !$acc end parallel loop
 
     ! Clipping abs(wp3) to 100. This keeps wp3 from growing too large in some 
     ! deep convective cases, which helps prevent these cases from blowing up.
@@ -1275,7 +1275,7 @@ subroutine clip_skewness_core( nz, ngrdcol, gr, sfc_elevation, &
         end if
       end do
     end do
-    !$acc end parallel
+    !$acc end parallel loop
 
     !$acc exit data delete( wp2_zt_cubed, wp3_lim_sqd, zagl_thresh, H_zagl )
 
diff --git a/src/CLUBB_core/clip_semi_implicit.F90 b/src/CLUBB_core/clip_semi_implicit.F90
index ae1f53cad..65d7947df 100644
--- a/src/CLUBB_core/clip_semi_implicit.F90
+++ b/src/CLUBB_core/clip_semi_implicit.F90
@@ -366,7 +366,7 @@ function clip_semi_imp_lhs( dt, f_unclipped,  &
   end function clip_semi_imp_lhs
 
   !=============================================================================
-  pure function compute_clip_lhs( dt_clip, B_fnc ) & 
+  function compute_clip_lhs( dt_clip, B_fnc ) & 
   result( lhs_contribution )
 
     ! Description:
diff --git a/src/CLUBB_core/clubb_api_module.F90 b/src/CLUBB_core/clubb_api_module.F90
index 4a54552dc..afd3038a2 100644
--- a/src/CLUBB_core/clubb_api_module.F90
+++ b/src/CLUBB_core/clubb_api_module.F90
@@ -3877,7 +3877,7 @@ end function calculate_spurious_source_api
   !================================================================================================
   ! calculate_thlp2_rad - Computes the contribution of radiative cooling to thlp2
   !================================================================================================
-  pure subroutine calculate_thlp2_rad_api &
+  subroutine calculate_thlp2_rad_api &
                   ( nz, rcm_zm, thlprcp, radht_zm, & ! Intent(in)
                     clubb_params,                  & ! Intent(in)
                     thlp2_forcing )                  ! Intent(inout)
diff --git a/src/CLUBB_core/corr_varnce_module.F90 b/src/CLUBB_core/corr_varnce_module.F90
index cb1269c96..ea254b639 100644
--- a/src/CLUBB_core/corr_varnce_module.F90
+++ b/src/CLUBB_core/corr_varnce_module.F90
@@ -224,7 +224,7 @@ subroutine init_default_corr_arrays(  )
   end subroutine init_default_corr_arrays
 
   !-----------------------------------------------------------------------------
-  pure function def_corr_idx( iiPDF_x ) result(ii_def_corr)
+  function def_corr_idx( iiPDF_x ) result(ii_def_corr)
 
     ! Description:
     !   Map from a iiPDF index to the corresponding index in the default 
diff --git a/src/CLUBB_core/diffusion.F90 b/src/CLUBB_core/diffusion.F90
index f0ce296e6..3ab8793f1 100644
--- a/src/CLUBB_core/diffusion.F90
+++ b/src/CLUBB_core/diffusion.F90
@@ -31,7 +31,7 @@ module diffusion
   contains
 
   !=============================================================================
-  pure subroutine diffusion_zt_lhs( nz, ngrdcol, gr, K_zm, K_zt, nu,  & ! In
+  subroutine diffusion_zt_lhs( nz, ngrdcol, gr, K_zm, K_zt, nu,  & ! In
                                     invrs_rho_ds_zt, rho_ds_zm,       & ! In
                                     lhs )                               ! Out
 
@@ -522,7 +522,7 @@ pure subroutine diffusion_zt_lhs( nz, ngrdcol, gr, K_zm, K_zt, nu,  & ! In
   end subroutine diffusion_zt_lhs
 
   !=============================================================================
-  pure function diffusion_cloud_frac_zt_lhs &
+  function diffusion_cloud_frac_zt_lhs &
                 ( gr, K_zm, K_zmm1, cloud_frac_zt, cloud_frac_ztm1, &
                   cloud_frac_ztp1, cloud_frac_zm, &
                   cloud_frac_zmm1, &
@@ -693,7 +693,7 @@ pure function diffusion_cloud_frac_zt_lhs &
   end function diffusion_cloud_frac_zt_lhs
 
   !=============================================================================
-  pure subroutine diffusion_zm_lhs( nz, ngrdcol, gr, K_zt, K_zm, nu,   & ! In
+  subroutine diffusion_zm_lhs( nz, ngrdcol, gr, K_zt, K_zm, nu,   & ! In
                                     invrs_rho_ds_zm, rho_ds_zt, & ! In 
                                     lhs )                   ! Out
 
diff --git a/src/CLUBB_core/fill_holes.F90 b/src/CLUBB_core/fill_holes.F90
index 81b49d48a..96f6496ae 100644
--- a/src/CLUBB_core/fill_holes.F90
+++ b/src/CLUBB_core/fill_holes.F90
@@ -120,10 +120,8 @@ subroutine fill_holes_vertical( nz, ngrdcol, num_draw_pts, threshold, upper_hf_l
 
     ! --------------------- Begin Code --------------------- 
 
-    !$acc declare copyin( dz, rho_ds ) &
-    !$acc           copy( field ) &
-    !$acc         create( invrs_denom_integral, field_clipped, denom_integral_global, rho_ds_dz, &
-    !$acc                 numer_integral_global, field_avg_global, mass_fraction_global )
+    !$acc enter data create( invrs_denom_integral, field_clipped, denom_integral_global, rho_ds_dz, &
+    !$acc                    numer_integral_global, field_avg_global, mass_fraction_global )
 
     l_field_below_threshold = .false.
 
@@ -140,6 +138,8 @@ subroutine fill_holes_vertical( nz, ngrdcol, num_draw_pts, threshold, upper_hf_l
 
     ! If all field values are above the specified threshold, no hole filling is required
     if ( .not. l_field_below_threshold ) then
+      !$acc exit data delete( invrs_denom_integral, field_clipped, denom_integral_global, rho_ds_dz, &
+      !$acc                   numer_integral_global, field_avg_global, mass_fraction_global )
       return
     end if
 
@@ -236,6 +236,8 @@ subroutine fill_holes_vertical( nz, ngrdcol, num_draw_pts, threshold, upper_hf_l
 
     ! If all field values are above the threshold, no further hole filling is required
     if ( .not. l_field_below_threshold ) then
+      !$acc exit data delete( invrs_denom_integral, field_clipped, denom_integral_global, rho_ds_dz, &
+      !$acc                   numer_integral_global, field_avg_global, mass_fraction_global )
       return
     end if
 
@@ -322,6 +324,9 @@ subroutine fill_holes_vertical( nz, ngrdcol, num_draw_pts, threshold, upper_hf_l
 
     end do
     !$acc end parallel loop
+    
+    !$acc exit data delete( invrs_denom_integral, field_clipped, denom_integral_global, rho_ds_dz, &
+    !$acc                   numer_integral_global, field_avg_global, mass_fraction_global )
 
     return
 
@@ -822,12 +827,17 @@ subroutine fill_holes_driver( gr, nz, dt, hydromet_dim,        & ! Intent(in)
 
          if ( hydromet_name(1:1) == "r" .and. l_hole_fill ) then
 
+            !$acc data copyin( gr, gr%dzt, rho_ds_zt ) &
+            !$acc        copy( hydromet(:,i) )
+
             ! Apply the hole filling algorithm
             ! upper_hf_level = nz since we are filling the zt levels
             call fill_holes_vertical( gr%nz, 1, num_hf_draw_points, zero_threshold, gr%nz, & ! In
                                       gr%dzt, rho_ds_zt,                                   & ! In
                                       hydromet(:,i) )                                      ! InOut
 
+            !$acc end data
+
          endif ! Variable is a mixing ratio and l_hole_fill is true
 
       endif ! hydromet(:,i) < 0
diff --git a/src/CLUBB_core/grid_class.F90 b/src/CLUBB_core/grid_class.F90
index a94de605c..7bdb7f905 100644
--- a/src/CLUBB_core/grid_class.F90
+++ b/src/CLUBB_core/grid_class.F90
@@ -1447,7 +1447,7 @@ function redirect_interpolated_azt_2D( nz, ngrdcol, gr, azm )
   end function redirect_interpolated_azt_2D
   
   !=============================================================================
-  pure subroutine linear_interpolated_azm_2D( nz, ngrdcol, gr, azt, &
+  subroutine linear_interpolated_azm_2D( nz, ngrdcol, gr, azt, &
                                               linear_interpolated_azm )
 
     ! Description:
@@ -1708,7 +1708,7 @@ function cubic_interpolated_azm_2D( nz, ngrdcol, gr, azt )
   end function cubic_interpolated_azm_2D
 
   !=============================================================================
-  pure subroutine calc_zt2zm_weights( nz, ngrdcol, &
+  subroutine calc_zt2zm_weights( nz, ngrdcol, &
                                       gr ) 
 
     ! Description:
@@ -1918,7 +1918,7 @@ pure subroutine calc_zt2zm_weights( nz, ngrdcol, &
   end subroutine calc_zt2zm_weights
   
   !=============================================================================
-  pure subroutine linear_interpolated_azt_2D( nz, ngrdcol, gr, azm, &
+  subroutine linear_interpolated_azt_2D( nz, ngrdcol, gr, azm, &
                                               linear_interpolated_azt )
 
     ! Description:
@@ -2073,7 +2073,7 @@ function cubic_interpolated_azt_2D( nz, ngrdcol, gr, azm )
   end function cubic_interpolated_azt_2D
 
   !=============================================================================
-  pure subroutine calc_zm2zt_weights( nz, ngrdcol, &
+  subroutine calc_zm2zt_weights( nz, ngrdcol, &
                                       gr )
 
     ! Description:
@@ -2286,7 +2286,7 @@ end subroutine calc_zm2zt_weights
   
   !=============================================================================
   ! Wrapped in interface ddzm
-  pure function gradzm_2D( nz, ngrdcol, gr, azm )
+  function gradzm_2D( nz, ngrdcol, gr, azm )
 
     ! Description:
     !  2D version of gradzm
@@ -2339,7 +2339,7 @@ end function gradzm_2D
   
   !=============================================================================
   ! Wrapped in interface ddzm
-  pure function gradzm_1D( gr, azm )
+  function gradzm_1D( gr, azm )
 
     ! Description:
     !  2D version of gradzm
@@ -2380,7 +2380,7 @@ end function gradzm_1D
   
   !=============================================================================
   ! Wrapped in interface ddzt
-  pure function gradzt_2D( nz, ngrdcol, gr, azt )
+  function gradzt_2D( nz, ngrdcol, gr, azt )
 
     ! Description:
     !  2D version of gradzt
@@ -2433,7 +2433,7 @@ end function gradzt_2D
   
   !=============================================================================
   ! Wrapped in interface ddzt
-  pure function gradzt_1D( gr, azt )
+  function gradzt_1D( gr, azt )
 
     ! Description:
     !  2D version of gradzt
@@ -2474,7 +2474,7 @@ pure function gradzt_1D( gr, azt )
   end function gradzt_1D
 
   !=============================================================================
-  pure function flip( x, xdim )
+  function flip( x, xdim )
 
     ! Description:
     !   Flips a single dimension array (i.e. a vector), so the first element
diff --git a/src/CLUBB_core/mean_adv.F90 b/src/CLUBB_core/mean_adv.F90
index 69a5e2f20..8535cab6a 100644
--- a/src/CLUBB_core/mean_adv.F90
+++ b/src/CLUBB_core/mean_adv.F90
@@ -30,7 +30,7 @@ module mean_adv
   contains
 
   !=============================================================================
-  pure subroutine term_ma_zt_lhs( nz, ngrdcol, wm_zt, weights_zt2zm, & ! Intent(in)
+  subroutine term_ma_zt_lhs( nz, ngrdcol, wm_zt, weights_zt2zm, & ! Intent(in)
                                   invrs_dzt, invrs_dzm,     & ! Intent(in)
                                   l_upwind_xm_ma,           & ! Intent(in)
                                   lhs_ma )                    ! Intent(out)
@@ -354,7 +354,7 @@ pure subroutine term_ma_zt_lhs( nz, ngrdcol, wm_zt, weights_zt2zm, & ! Intent(in
   end subroutine term_ma_zt_lhs
 
   !=============================================================================
-  pure subroutine term_ma_zm_lhs( nz, ngrdcol, wm_zm, &
+  subroutine term_ma_zm_lhs( nz, ngrdcol, wm_zm, &
                                   invrs_dzm, weights_zm2zt, & 
                                   lhs_ma )
 
diff --git a/src/CLUBB_core/mono_flux_limiter.F90 b/src/CLUBB_core/mono_flux_limiter.F90
index 4b95fce9a..dcf64118e 100644
--- a/src/CLUBB_core/mono_flux_limiter.F90
+++ b/src/CLUBB_core/mono_flux_limiter.F90
@@ -1342,7 +1342,7 @@ subroutine mfl_xm_solve( nz, ngrdcol, solve_type, tridiag_solve_method, &
     do i = 1, ngrdcol
       xm(i,1) = xm(i,2)
     end do
-    !$acc end parallel
+    !$acc end parallel loop
 
     return
   end subroutine mfl_xm_solve
diff --git a/src/CLUBB_core/numerical_check.F90 b/src/CLUBB_core/numerical_check.F90
index d749a4ca1..a1569eb95 100644
--- a/src/CLUBB_core/numerical_check.F90
+++ b/src/CLUBB_core/numerical_check.F90
@@ -1035,7 +1035,7 @@ end subroutine check_nan_sclr
 !-------------------------------------------------------------------------
 
 !-----------------------------------------------------------------------
-  pure function calculate_spurious_source( integral_after, integral_before, &
+  function calculate_spurious_source( integral_after, integral_before, &
                                            flux_top, flux_sfc, & 
                                            integral_forcing, dt ) &
   result( spurious_source )
diff --git a/src/CLUBB_core/penta_lu_solver.F90 b/src/CLUBB_core/penta_lu_solver.F90
index f7f503d8d..2bb2aad8b 100644
--- a/src/CLUBB_core/penta_lu_solver.F90
+++ b/src/CLUBB_core/penta_lu_solver.F90
@@ -125,18 +125,18 @@ subroutine penta_lu_solve_single_rhs_multiple_lhs( ndim, ngrdcol, lhs, rhs, &
       upper_1,          & ! First U band 
       upper_2,          & ! Second U band 
       lower_diag_invrs, & ! Inverse of the diagonal of L
-      lower_1             ! First L band
+      lower_1,          & ! First L band
+      lower_2             ! Second L band
 
     integer :: i, k, j    ! Loop variables
 
     ! ----------------------- Begin Code -----------------------
        
-    !$acc data create( upper_1, upper_2, lower_1, lower_diag_invrs ) &
+    !$acc data create( upper_1, upper_2, lower_1, lower_2, lower_diag_invrs ) &
     !$acc      copyin( rhs, lhs ) &
     !$acc      copyout( soln )
-    
-    !$acc kernels
 
+    !$acc parallel loop gang vector default(present)
     do i = 1, ngrdcol
       lower_diag_invrs(i,1) = 1.0_core_rknd / lhs(0,i,1)
       upper_1(i,1)          = lower_diag_invrs(i,1) * lhs(-1,i,1) 
@@ -147,38 +147,45 @@ subroutine penta_lu_solve_single_rhs_multiple_lhs( ndim, ngrdcol, lhs, rhs, &
       upper_1(i,2)          = lower_diag_invrs(i,2) * ( lhs(-1,i,2) - lower_1(i,2) * upper_2(i,1) )
       upper_2(i,2)          = lower_diag_invrs(i,2) * lhs(-2,i,2)
     end do
+    !$acc end parallel loop
 
+    !$acc parallel loop gang vector default(present)
+    do i = 1, ngrdcol
+      do k = 3, ndim-2
+        lower_2(i,k) = lhs(2,i,k)
+        lower_1(i,k) = lhs(1,i,k) - lower_2(i,k) * upper_1(i,k-2)
 
-    do k = 3, ndim-2
-      do i = 1, ngrdcol
-        lower_1(i,k)          = lhs(1,i,k) - lhs(2,i,k) * upper_1(i,k-2)
-
-        lower_diag_invrs(i,k) = 1.0_core_rknd / ( lhs(0,i,k) - lhs(2,i,k) * upper_2(i,k-2) &
+        lower_diag_invrs(i,k) = 1.0_core_rknd / ( lhs(0,i,k) - lower_2(i,k) * upper_2(i,k-2) &
                                                              - lower_1(i,k) * upper_1(i,k-1) )
 
-        upper_1(i,k)          = lower_diag_invrs(i,k) * ( lhs(-1,i,k) - lower_1(i,k) * upper_2(i,k-1) )
-        upper_2(i,k)          = lower_diag_invrs(i,k) * lhs(-2,i,k)
+        upper_1(i,k) = lower_diag_invrs(i,k) * ( lhs(-1,i,k) - lower_1(i,k) * upper_2(i,k-1) )
+        upper_2(i,k) = lower_diag_invrs(i,k) * lhs(-2,i,k)
       end do
     end do
+    !$acc end parallel loop
 
+    !$acc parallel loop gang vector default(present)
     do i = 1, ngrdcol
-      lower_1(i,ndim-1) = lhs(1,i,ndim-1) - lhs(2,i,ndim-1) * upper_1(i,ndim-3)
+      lower_2(i,ndim-1) = lhs(2,i,ndim-1)
+      lower_1(i,ndim-1) = lhs(1,i,ndim-1) - lower_2(i,ndim-1) * upper_1(i,ndim-3)
 
       lower_diag_invrs(i,ndim-1) = 1.0_core_rknd  &
-                                   / ( lhs(0,i,ndim-1) - lhs(2,i,ndim-1) * upper_2(i,ndim-3) &
+                                   / ( lhs(0,i,ndim-1) - lower_2(i,ndim-1) * upper_2(i,ndim-3) &
                                                        - lower_1(i,ndim-1) * upper_1(i,ndim-2) )
 
       upper_1(i,ndim-1)  = lower_diag_invrs(i,ndim-1) * ( lhs(-1,i,ndim-1) - lower_1(i,ndim-1) &
-                                                      * upper_2(i,ndim-2) )
+                                                                             * upper_2(i,ndim-2) )
 
-      lower_1(i,ndim) = lhs(1,i,ndim) - lhs(2,i,ndim) * upper_1(i,ndim-2)
+      lower_2(i,ndim) = lhs(2,i,ndim)
+      lower_1(i,ndim) = lhs(1,i,ndim) - lower_2(i,ndim) * upper_1(i,ndim-2)
 
       lower_diag_invrs(i,ndim) = 1.0_core_rknd  &
-                                 / ( lhs(0,i,ndim-1) - lhs(2,i,ndim) * upper_2(i,ndim-2) &
+                                 / ( lhs(0,i,ndim-1) - lower_2(i,ndim) * upper_2(i,ndim-2) &
                                                      - lower_1(i,ndim) * upper_1(i,ndim-1) )
     end do
-
+    !$acc end parallel loop
     
+    !$acc parallel loop gang vector default(present)
     do i = 1, ngrdcol 
 
       soln(i,1)   = lower_diag_invrs(i,1) * rhs(i,1) 
@@ -186,10 +193,14 @@ subroutine penta_lu_solve_single_rhs_multiple_lhs( ndim, ngrdcol, lhs, rhs, &
       soln(i,2)   = lower_diag_invrs(i,2) * ( rhs(i,2) - lower_1(i,2) * soln(i,1) )
 
       do k = 3, ndim
-        soln(i,k) = lower_diag_invrs(i,k) * ( rhs(i,k) - lhs(2,i,k) * soln(i,k-2) &
+        soln(i,k) = lower_diag_invrs(i,k) * ( rhs(i,k) - lower_2(i,k) * soln(i,k-2) &
                                                        - lower_1(i,k) * soln(i,k-1) )
       end do
+    end do
+    !$acc end parallel loop
 
+    !$acc parallel loop gang vector default(present)
+    do i = 1, ngrdcol 
       soln(i,ndim-1) = soln(i,ndim-1) - upper_1(i,ndim-1) * soln(i,ndim)
 
       do k = ndim-2, 1, -1
@@ -197,8 +208,7 @@ subroutine penta_lu_solve_single_rhs_multiple_lhs( ndim, ngrdcol, lhs, rhs, &
       end do
 
     end do
-
-    !$acc end kernels
+    !$acc end parallel loop
 
     !$acc end data
 
@@ -237,18 +247,18 @@ subroutine penta_lu_solve_multiple_rhs_lhs( ndim, nrhs, ngrdcol, lhs, rhs, &
       upper_1,          & ! First U band 
       upper_2,          & ! Second U band 
       lower_diag_invrs, & ! Inverse of the diagonal of L
-      lower_1             ! First L band
+      lower_1,          & ! First L band
+      lower_2             ! Second L band
 
     integer :: i, k, j    ! Loop variables
 
     ! ----------------------- Begin Code -----------------------
        
-    !$acc data create( upper_1, upper_2, lower_1, lower_diag_invrs ) &
+    !$acc data create( upper_1, upper_2, lower_1, lower_2, lower_diag_invrs ) &
     !$acc      copyin( rhs, lhs ) &
     !$acc      copyout( soln )
 
-    !$acc kernels
-
+    !$acc parallel loop gang vector default(present)
     do i = 1, ngrdcol
       lower_diag_invrs(i,1) = 1.0_core_rknd / lhs(0,i,1)
       upper_1(i,1)          = lower_diag_invrs(i,1) * lhs(-1,i,1) 
@@ -259,38 +269,45 @@ subroutine penta_lu_solve_multiple_rhs_lhs( ndim, nrhs, ngrdcol, lhs, rhs, &
       upper_1(i,2)          = lower_diag_invrs(i,2) * ( lhs(-1,i,2) - lower_1(i,2) * upper_2(i,1) )
       upper_2(i,2)          = lower_diag_invrs(i,2) * lhs(-2,i,2)
     end do
+    !$acc end parallel loop
 
+    !$acc parallel loop gang vector default(present)
+    do i = 1, ngrdcol
+      do k = 3, ndim-2
+        lower_2(i,k) = lhs(2,i,k)
+        lower_1(i,k) = lhs(1,i,k) - lower_2(i,k) * upper_1(i,k-2)
 
-    do k = 3, ndim-2
-      do i = 1, ngrdcol
-        lower_1(i,k)          = lhs(1,i,k) - lhs(2,i,k) * upper_1(i,k-2)
-
-        lower_diag_invrs(i,k) = 1.0_core_rknd / ( lhs(0,i,k) - lhs(2,i,k) * upper_2(i,k-2) &
+        lower_diag_invrs(i,k) = 1.0_core_rknd / ( lhs(0,i,k) - lower_2(i,k) * upper_2(i,k-2) &
                                                              - lower_1(i,k) * upper_1(i,k-1) )
 
-        upper_1(i,k)          = lower_diag_invrs(i,k) * ( lhs(-1,i,k) - lower_1(i,k) * upper_2(i,k-1) )
-        upper_2(i,k)          = lower_diag_invrs(i,k) * lhs(-2,i,k)
+        upper_1(i,k) = lower_diag_invrs(i,k) * ( lhs(-1,i,k) - lower_1(i,k) * upper_2(i,k-1) )
+        upper_2(i,k) = lower_diag_invrs(i,k) * lhs(-2,i,k)
       end do
     end do
+    !$acc end parallel loop
 
+    !$acc parallel loop gang vector default(present)
     do i = 1, ngrdcol
-      lower_1(i,ndim-1) = lhs(1,i,ndim-1) - lhs(2,i,ndim-1) * upper_1(i,ndim-3)
+      lower_2(i,ndim-1) = lhs(2,i,ndim-1)
+      lower_1(i,ndim-1) = lhs(1,i,ndim-1) - lower_2(i,ndim-1) * upper_1(i,ndim-3)
 
       lower_diag_invrs(i,ndim-1) = 1.0_core_rknd  &
-                                   / ( lhs(0,i,ndim-1) - lhs(2,i,ndim-1) * upper_2(i,ndim-3) &
+                                   / ( lhs(0,i,ndim-1) - lower_2(i,ndim-1) * upper_2(i,ndim-3) &
                                                        - lower_1(i,ndim-1) * upper_1(i,ndim-2) )
 
       upper_1(i,ndim-1) = lower_diag_invrs(i,ndim-1) * ( lhs(-1,i,ndim-1) - lower_1(i,ndim-1) &
-                                                     * upper_2(i,ndim-2) )
+                                                                            * upper_2(i,ndim-2) )
 
-      lower_1(i,ndim) = lhs(1,i,ndim) - lhs(2,i,ndim) * upper_1(i,ndim-2)
+      lower_2(i,ndim) = lhs(2,i,ndim)
+      lower_1(i,ndim) = lhs(1,i,ndim) - lower_2(i,ndim) * upper_1(i,ndim-2)
 
       lower_diag_invrs(i,ndim) = 1.0_core_rknd  &
-                                 / ( lhs(0,i,ndim-1) - lhs(2,i,ndim) * upper_2(i,ndim-2) &
+                                 / ( lhs(0,i,ndim-1) - lower_2(i,ndim) * upper_2(i,ndim-2) &
                                                      - lower_1(  i,ndim) * upper_1(i,ndim-1) )
     end do
+    !$acc end parallel loop
 
-    
+    !$acc parallel loop gang vector collapse(2) default(present)
     do j = 1, nrhs
       do i = 1, ngrdcol 
 
@@ -299,10 +316,16 @@ subroutine penta_lu_solve_multiple_rhs_lhs( ndim, nrhs, ngrdcol, lhs, rhs, &
         soln(i,2,j)   = lower_diag_invrs(i,2) * ( rhs(i,2,j) - lower_1(i,2) * soln(i,1,j) )
 
         do k = 3, ndim
-          soln(i,k,j) = lower_diag_invrs(i,k) * ( rhs(i,k,j) - lhs(2,i,k) * soln(i,k-2,j) &
+          soln(i,k,j) = lower_diag_invrs(i,k) * ( rhs(i,k,j) - lower_2(i,k) * soln(i,k-2,j) &
                                                              - lower_1(i,k) * soln(i,k-1,j) )
         end do
+      end do
+    end do
+    !$acc end parallel loop
 
+    !$acc parallel loop gang vector collapse(2) default(present)
+    do j = 1, nrhs
+      do i = 1, ngrdcol 
         soln(i,ndim-1,j) = soln(i,ndim-1,j) - upper_1(i,ndim-1) * soln(i,ndim,j)
 
         do k = ndim-2, 1, -1
@@ -311,8 +334,7 @@ subroutine penta_lu_solve_multiple_rhs_lhs( ndim, nrhs, ngrdcol, lhs, rhs, &
 
       end do
     end do
-
-    !$acc end kernels
+    !$acc end parallel loop
 
     !$acc end data
 
diff --git a/src/CLUBB_core/tridiag_lu_solver.F90 b/src/CLUBB_core/tridiag_lu_solver.F90
index 3aefdfaf3..6bcd094b2 100644
--- a/src/CLUBB_core/tridiag_lu_solver.F90
+++ b/src/CLUBB_core/tridiag_lu_solver.F90
@@ -123,8 +123,8 @@ subroutine tridiag_lu_solve_single_rhs_lhs( ndim, lhs, rhs, &
     upper(1)            = lower_diag_invrs(1) * lhs(-1,1) 
 
     do k = 2, ndim-1
-        lower_diag_invrs(k) = 1.0_core_rknd / ( lhs(0,k) - lhs(1,k) * upper(k-1)  )
-        upper(k)            = lower_diag_invrs(k) * lhs(-1,k) 
+      lower_diag_invrs(k) = 1.0_core_rknd / ( lhs(0,k) - lhs(1,k) * upper(k-1)  )
+      upper(k)            = lower_diag_invrs(k) * lhs(-1,k) 
     end do
 
     lower_diag_invrs(ndim) = 1.0_core_rknd / ( lhs(0,ndim) - lhs(1,ndim) * upper(ndim-1)  )
@@ -139,7 +139,6 @@ subroutine tridiag_lu_solve_single_rhs_lhs( ndim, lhs, rhs, &
       soln(k) = soln(k) - upper(k) * soln(k+1)
     end do
 
-
     !$acc end kernels
 
     !$acc end data
@@ -185,26 +184,29 @@ subroutine tridiag_lu_solve_single_rhs_multiple_lhs( ndim, ngrdcol, lhs, rhs, &
     !$acc      copyin( rhs, lhs ) &
     !$acc      copyout( soln )
     
-    !$acc kernels
-    
+    !$acc parallel loop gang vector default(present)
     do i = 1, ngrdcol
       lower_diag_invrs(i,1) = 1.0_core_rknd / lhs(0,i,1)
       upper(i,1)            = lower_diag_invrs(i,1) * lhs(-1,i,1) 
     end do
+    !$acc end parallel loop
 
-
+    !$acc parallel loop gang vector default(present)
     do k = 2, ndim-1
       do i = 1, ngrdcol
         lower_diag_invrs(i,k) = 1.0_core_rknd / ( lhs(0,i,k) - lhs(1,i,k) * upper(i,k-1)  )
         upper(i,k)            = lower_diag_invrs(i,k) * lhs(-1,i,k) 
       end do
     end do
+    !$acc end parallel loop
 
+    !$acc parallel loop gang vector default(present)
     do i = 1, ngrdcol
       lower_diag_invrs(i,ndim) = 1.0_core_rknd / ( lhs(0,i,ndim) - lhs(1,i,ndim) * upper(i,ndim-1)  )
     end do
+    !$acc end parallel loop
 
-    
+    !$acc parallel loop gang vector default(present)
     do i = 1, ngrdcol 
 
       soln(i,1)   = lower_diag_invrs(i,1) * rhs(i,1) 
@@ -212,14 +214,17 @@ subroutine tridiag_lu_solve_single_rhs_multiple_lhs( ndim, ngrdcol, lhs, rhs, &
       do k = 2, ndim
         soln(i,k) = lower_diag_invrs(i,k) * ( rhs(i,k) - lhs(1,i,k) * soln(i,k-1) )
       end do
+    end do
+    !$acc end parallel loop
 
+    !$acc parallel loop gang vector default(present)
+    do i = 1, ngrdcol 
       do k = ndim-1, 1, -1
         soln(i,k) = soln(i,k) - upper(i,k) * soln(i,k+1)
       end do
 
     end do
-
-    !$acc end kernels
+    !$acc end parallel loop
 
     !$acc end data
 
@@ -266,26 +271,29 @@ subroutine tridiag_lu_solve_multiple_rhs_lhs( ndim, nrhs, ngrdcol, lhs, rhs, &
     !$acc      copyin( rhs, lhs ) &
     !$acc      copyout( soln )
 
-    !$acc kernels
-
+    !$acc parallel loop gang vector default(present)
     do i = 1, ngrdcol
       lower_diag_invrs(i,1) = 1.0_core_rknd / lhs(0,i,1)
       upper(i,1)            = lower_diag_invrs(i,1) * lhs(-1,i,1) 
     end do
+    !$acc end parallel loop
 
-
-    do k = 2, ndim-1
-      do i = 1, ngrdcol
+    !$acc parallel loop gang vector default(present)
+    do i = 1, ngrdcol
+      do k = 2, ndim-1
         lower_diag_invrs(i,k) = 1.0_core_rknd / ( lhs(0,i,k) - lhs(1,i,k) * upper(i,k-1)  )
         upper(i,k)            = lower_diag_invrs(i,k) * lhs(-1,i,k) 
       end do
     end do
+    !$acc end parallel loop
 
+    !$acc parallel loop gang vector default(present)
     do i = 1, ngrdcol
       lower_diag_invrs(i,ndim) = 1.0_core_rknd / ( lhs(0,i,ndim) - lhs(1,i,ndim) * upper(i,ndim-1)  )
     end do
+    !$acc end parallel loop
 
-    
+    !$acc parallel loop gang vector collapse(2) default(present)
     do j = 1, nrhs
       do i = 1, ngrdcol 
 
@@ -294,15 +302,19 @@ subroutine tridiag_lu_solve_multiple_rhs_lhs( ndim, nrhs, ngrdcol, lhs, rhs, &
         do k = 2, ndim
           soln(i,k,j) = lower_diag_invrs(i,k) * ( rhs(i,k,j) - lhs(1,i,k) * soln(i,k-1,j) )
         end do
+      end do
+    end do
+    !$acc end parallel loop
 
+    !$acc parallel loop gang vector collapse(2) default(present)
+    do j = 1, nrhs
+      do i = 1, ngrdcol 
         do k = ndim-1, 1, -1
           soln(i,k,j) = soln(i,k,j) - upper(i,k) * soln(i,k+1,j)
         end do
-
       end do
     end do
-
-    !$acc end kernels
+    !$acc end parallel loop
 
     !$acc end data
 
diff --git a/src/CLUBB_core/turbulent_adv_pdf.F90 b/src/CLUBB_core/turbulent_adv_pdf.F90
index 46eb4241c..2d340b82d 100644
--- a/src/CLUBB_core/turbulent_adv_pdf.F90
+++ b/src/CLUBB_core/turbulent_adv_pdf.F90
@@ -30,7 +30,7 @@ module turbulent_adv_pdf
   contains
 
   !=============================================================================
-  pure subroutine xpyp_term_ta_pdf_lhs( nz, ngrdcol, gr, coef_wpxpyp_implicit, & ! In
+  subroutine xpyp_term_ta_pdf_lhs( nz, ngrdcol, gr, coef_wpxpyp_implicit, & ! In
                                         rho_ds_zt, rho_ds_zm,            & ! In
                                         invrs_rho_ds_zm,                 & ! In
                                         l_upwind_xpyp_turbulent_adv,     & ! In
@@ -478,7 +478,7 @@ pure subroutine xpyp_term_ta_pdf_lhs( nz, ngrdcol, gr, coef_wpxpyp_implicit, & !
   end subroutine xpyp_term_ta_pdf_lhs
 
   !=============================================================================================
-  pure subroutine xpyp_term_ta_pdf_lhs_godunov( nz, ngrdcol, gr, & ! Intent(in)
+  subroutine xpyp_term_ta_pdf_lhs_godunov( nz, ngrdcol, gr, & ! Intent(in)
                                                 coef_wpxpyp_implicit, & ! Intent(in)
                                                 invrs_rho_ds_zm, rho_ds_zm,  & ! Intent(in)
                                                 lhs_ta )
@@ -589,7 +589,7 @@ pure subroutine xpyp_term_ta_pdf_lhs_godunov( nz, ngrdcol, gr, & ! Intent(in)
   end subroutine xpyp_term_ta_pdf_lhs_godunov
 
   !=============================================================================
-  pure subroutine xpyp_term_ta_pdf_rhs( nz, ngrdcol, gr, term_wpxpyp_explicit,  & ! In
+  subroutine xpyp_term_ta_pdf_rhs( nz, ngrdcol, gr, term_wpxpyp_explicit,  & ! In
                                         rho_ds_zt, rho_ds_zm,                   & ! In
                                         invrs_rho_ds_zm,                        & ! In
                                         l_upwind_xpyp_turbulent_adv,            & ! In
@@ -960,7 +960,7 @@ pure subroutine xpyp_term_ta_pdf_rhs( nz, ngrdcol, gr, term_wpxpyp_explicit,  &
   end subroutine xpyp_term_ta_pdf_rhs
 
   !=============================================================================
-  pure subroutine xpyp_term_ta_pdf_rhs_godunov( nz, ngrdcol, gr, & ! Intent(in)
+  subroutine xpyp_term_ta_pdf_rhs_godunov( nz, ngrdcol, gr, & ! Intent(in)
                                                 term_wpxpyp_explicit_zm, & ! Intent(in)
                                                 invrs_rho_ds_zm, & ! Intent(in)
                                                 sgn_turbulent_vel, & ! Intent(in)