Going through second pass of marking subroutines

prathi-wind · prathi-wind · commit 0cdc1e1b114d · 2025-05-28T15:13:00.000-04:00
diff --git a/src/common/m_checker_common.fpp b/src/common/m_checker_common.fpp
@@ -170,7 +170,7 @@ contains
 
     !> Checks constraints on the Immersed Boundaries parameters.
         !! Called by s_check_inputs_common for pre-processing and simulation
-    subroutine s_check_inputs_ibm
+    impure subroutine s_check_inputs_ibm
         @:PROHIBIT(ib .and. n <= 0, "Immersed Boundaries do not work in 1D")
         @:PROHIBIT(ib .and. (num_ibs <= 0 .or. num_ibs > num_patches_max), "num_ibs must be between 1 and num_patches_max")
         @:PROHIBIT((.not. ib) .and. num_ibs > 0, "num_ibs is set, but ib is not enabled")
diff --git a/src/post_process/m_data_input.f90 b/src/post_process/m_data_input.f90
@@ -35,7 +35,7 @@ module m_data_input
 
         !> Subroutine for reading data files
         !!  @param t_step Current time-step to input
-        subroutine s_read_abstract_data_files(t_step)
+        impure subroutine s_read_abstract_data_files(t_step)
 
             implicit none
 
diff --git a/src/pre_process/m_compute_levelset.fpp b/src/pre_process/m_compute_levelset.fpp
@@ -23,16 +23,6 @@ module m_compute_levelset
  s_cuboid_levelset, &
  s_sphere_levelset
 
-    real(wp) :: x_centroid, y_centroid, z_centroid
-    real(wp) :: length_x, length_y, length_z
-    real(wp) :: radius
-
-    type(bounds_info) :: x_boundary, y_boundary, z_boundary  !<
-        !! These variables combine the centroid and length parameters associated with
-        !! a particular patch to yield the locations of the patch boundaries in the
-        !! x-, y- and z-coordinate directions. They are used as a means to concisely
-        !! perform the actions necessary to lay out a particular patch on the grid.
-
 contains
 
     pure subroutine s_circle_levelset(levelset, levelset_norm, ib_patch_id)
@@ -165,6 +155,8 @@ contains
         real(wp) :: x_centroid, y_centroid, z_centroid, lz, z_max, z_min, x_act, y_act, theta
         real(wp), dimension(3) :: dist_vec
 
+        real(wp) :: length_z
+
         integer :: i, j, k, l !< Loop index variables
 
         x_centroid = patch_ib(ib_patch_id)%x_centroid
@@ -256,7 +248,7 @@ contains
     end subroutine s_3D_airfoil_levelset
 
     !>  Initialize IBM module
-    subroutine s_rectangle_levelset(levelset, levelset_norm, ib_patch_id)
+    pure subroutine s_rectangle_levelset(levelset, levelset_norm, ib_patch_id)
 
         type(levelset_field), intent(INOUT) :: levelset
         type(levelset_norm_field), intent(INOUT) :: levelset_norm
@@ -266,6 +258,9 @@ contains
         real(wp) :: x, y, min_dist
         real(wp) :: side_dists(4)
 
+        real(wp) :: x_centroid, y_centroid
+        real(wp) :: length_x, length_y
+
         integer :: i, j, k !< Loop index variables
         integer :: idx !< Shortest path direction indicator
 
@@ -350,7 +345,7 @@ contains
 
     end subroutine s_rectangle_levelset
 
-    subroutine s_cuboid_levelset(levelset, levelset_norm, ib_patch_id)
+    pure subroutine s_cuboid_levelset(levelset, levelset_norm, ib_patch_id)
 
         type(levelset_field), intent(INOUT) :: levelset
         type(levelset_norm_field), intent(INOUT) :: levelset_norm
@@ -359,6 +354,9 @@ contains
         real(wp) :: Right, Left, Bottom, Top, Front, Back
         real(wp) :: x, y, z, min_dist
         real(wp) :: side_dists(6)
+        
+        real(wp) :: x_centroid, y_centroid, z_centroid
+        real(wp) :: length_x, length_y, length_z
 
         integer :: i, j, k !< Loop index variables
 
diff --git a/src/pre_process/m_grid.f90 b/src/pre_process/m_grid.f90
@@ -35,7 +35,7 @@ module m_grid
 
     abstract interface
 
-        subroutine s_generate_abstract_grid
+        impure subroutine s_generate_abstract_grid
 
         end subroutine s_generate_abstract_grid
 
diff --git a/src/pre_process/m_initial_condition.fpp b/src/pre_process/m_initial_condition.fpp
@@ -173,7 +173,7 @@ contains
         !!              on the grid using the primitive variables included with
         !!              the patch parameters. The subroutine is complete once the
         !!              primitive variables are converted to conservative ones.
-    subroutine s_generate_initial_condition
+    impure subroutine s_generate_initial_condition
 
         integer :: i  !< Generic loop operator
 
diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp
@@ -76,13 +76,13 @@ module m_start_up
 
     abstract interface
 
-        subroutine s_read_abstract_grid_data_files
+        impure subroutine s_read_abstract_grid_data_files
 
         end subroutine s_read_abstract_grid_data_files
 
         !! @param q_cons_vf Conservative variables
         !! @param ib_markers track if a cell is within the immersed boundary
-        subroutine s_read_abstract_ic_data_files(q_cons_vf, ib_markers)
+        impure subroutine s_read_abstract_ic_data_files(q_cons_vf, ib_markers)
 
             import :: scalar_field, integer_field, sys_size, pres_field
 
diff --git a/src/simulation/m_bubbles.fpp b/src/simulation/m_bubbles.fpp
@@ -314,11 +314,13 @@ contains
     !>  Subroutine that computes bubble wall properties for vapor bubbles
         !!  @param pb Internal bubble pressure
         !!  @param iR0 Current bubble size index
-    subroutine s_bwproperty(pb, iR0)
+    pure elemental subroutine s_bwproperty(pb, iR0, chi_vw, k_mw, rho_mw)
         !$acc routine seq
         real(wp), intent(in) :: pb
         integer, intent(in) :: iR0
-
+        real(wp), intent(out) :: chi_vw  !< Bubble wall properties (Ando 2010)
+        real(wp), intent(out) :: k_mw    !< Bubble wall properties (Ando 2010)
+        real(wp), intent(out) :: rho_mw  !< Bubble wall properties (Ando 2010)
         real(wp) :: x_vw
 
         ! mass fraction of vapor
@@ -342,21 +344,21 @@ contains
         !!  @param fbeta_c Mass transfer coefficient (EL)
         !!  @param fR_m Mixture gas constant (EL)
         !!  @param fgamma_m Mixture gamma (EL)
-    function f_vflux(fR, fV, fpb, fmass_v, iR0, fmass_n, fbeta_c, fR_m, fgamma_m)
+    pure elemental subroutine s_vflux(fR, fV, fpb, fmass_v, iR0, vflux, fmass_n, fbeta_c, fR_m, fgamma_m)
         !$acc routine seq
         real(wp), intent(in) :: fR
         real(wp), intent(in) :: fV
         real(wp), intent(in) :: fpb
         real(wp), intent(in) :: fmass_v
         integer, intent(in) :: iR0
+        real(wp), intent(out) :: vflux
         real(wp), intent(in), optional :: fmass_n, fbeta_c
         real(wp), intent(out), optional :: fR_m, fgamma_m
 
         real(wp) :: chi_bar
         real(wp) :: rho_mw_lag
         real(wp) :: grad_chi
         real(wp) :: conc_v
-        real(wp) :: f_vflux
 
         if (thermal == 3) then !transfer
             ! constant transfer model
@@ -373,21 +375,21 @@ contains
                 chi_bar = fmass_v/(fmass_v + fmass_n)
                 grad_chi = (chi_bar - conc_v)
                 rho_mw_lag = (fmass_n + fmass_v)/(4._wp/3._wp*pi*fR**3._wp)
-                f_vflux = 0._wp
+                vflux = 0._wp
                 if (lag_params%massTransfer_model) then
-                    f_vflux = -fbeta_c*rho_mw_lag*grad_chi/(1._wp - conc_v)/fR
+                    vflux = -fbeta_c*rho_mw_lag*grad_chi/(1._wp - conc_v)/fR
                 end if
             else
                 chi_bar = fmass_v/(fmass_v + mass_n0(iR0))
                 grad_chi = -Re_trans_c(iR0)*(chi_bar - chi_vw)
-                f_vflux = rho_mw*grad_chi/Pe_c/(1._wp - chi_vw)/fR
+                vflux = rho_mw*grad_chi/Pe_c/(1._wp - chi_vw)/fR
             end if
         else
             ! polytropic
-            f_vflux = pv*fV/(R_v*Tw)
+            vflux = pv*fV/(R_v*Tw)
         end if
 
-    end function f_vflux
+    end subroutine s_vflux
 
     !>  Function that computes the time derivative of
         !!  the internal bubble pressure
@@ -457,7 +459,7 @@ contains
         !!  @param fbeta_t Heat transfer coefficient (EL)
         !!  @param fCson Speed of sound (EL)
         !!  @param adap_dt_stop Fail-safe exit if max iteration count reached
-    subroutine s_advance_step(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
+    pure subroutine s_advance_step(fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
                               fntait, fBtait, f_bub_adv_src, f_divu, &
                               bub_id, fmass_v, fmass_n, fbeta_c, &
                               fbeta_t, fCson, adap_dt_stop)
@@ -676,7 +678,7 @@ contains
         !!  @param myV_tmp Bubble radial velocity at each stage
         !!  @param myPb_tmp Internal bubble pressure at each stage (EL)
         !!  @param myMv_tmp Mass of vapor in the bubble at each stage (EL)
-    subroutine s_advance_substep(err, fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
+    pure subroutine s_advance_substep(err, fRho, fP, fR, fV, fR0, fpb, fpbdot, alf, &
                                  fntait, fBtait, f_bub_adv_src, f_divu, &
                                  bub_id, fmass_v, fmass_n, fbeta_c, &
                                  fbeta_t, fCson, h, &
@@ -705,8 +707,8 @@ contains
         if (bubbles_lagrange) then
             myPb_tmp(1) = fpb
             myMv_tmp(1) = fmass_v
-            mydMvdt_tmp(1) = f_advance_EL(myR_tmp(1), myV_tmp(1), myPb_tmp(1), myMv_tmp(1), bub_id, &
-                                          fmass_n, fbeta_c, fbeta_t, mydPbdt_tmp(1))
+            call s_advance_EL(myR_tmp(1), myV_tmp(1), myPb_tmp(1), myMv_tmp(1), bub_id, &
+                                          fmass_n, fbeta_c, fbeta_t, mydPbdt_tmp(1), mydMvdt_tmp(1))
         end if
         myA_tmp(1) = f_rddot(fRho, fP, myR_tmp(1), myV_tmp(1), fR0, &
                              myPb_tmp(1), mydPbdt_tmp(1), alf, fntait, fBtait, &
@@ -719,8 +721,8 @@ contains
         if (bubbles_lagrange) then
             myPb_tmp(2) = myPb_tmp(1) + h*mydPbdt_tmp(1)
             myMv_tmp(2) = myMv_tmp(1) + h*mydMvdt_tmp(1)
-            mydMvdt_tmp(2) = f_advance_EL(myR_tmp(2), myV_tmp(2), myPb_tmp(2), myMv_tmp(2), &
-                                          bub_id, fmass_n, fbeta_c, fbeta_t, mydPbdt_tmp(2))
+            call s_advance_EL(myR_tmp(2), myV_tmp(2), myPb_tmp(2), myMv_tmp(2), &
+                                          bub_id, fmass_n, fbeta_c, fbeta_t, mydPbdt_tmp(2), mydMvdt_tmp(2))
         end if
         myA_tmp(2) = f_rddot(fRho, fP, myR_tmp(2), myV_tmp(2), fR0, &
                              myPb_tmp(2), mydPbdt_tmp(2), alf, fntait, fBtait, &
@@ -733,8 +735,8 @@ contains
         if (bubbles_lagrange) then
             myPb_tmp(3) = myPb_tmp(1) + (h/4._wp)*(mydPbdt_tmp(1) + mydPbdt_tmp(2))
             myMv_tmp(3) = myMv_tmp(1) + (h/4._wp)*(mydMvdt_tmp(1) + mydMvdt_tmp(2))
-            mydMvdt_tmp(3) = f_advance_EL(myR_tmp(3), myV_tmp(3), myPb_tmp(3), myMv_tmp(3), &
-                                          bub_id, fmass_n, fbeta_c, fbeta_t, mydPbdt_tmp(3))
+            call s_advance_EL(myR_tmp(3), myV_tmp(3), myPb_tmp(3), myMv_tmp(3), &
+                                          bub_id, fmass_n, fbeta_c, fbeta_t, mydPbdt_tmp(3), mydMvdt_tmp(3))
         end if
         myA_tmp(3) = f_rddot(fRho, fP, myR_tmp(3), myV_tmp(3), fR0, &
                              myPb_tmp(3), mydPbdt_tmp(3), alf, fntait, fBtait, &
@@ -747,8 +749,8 @@ contains
         if (bubbles_lagrange) then
             myPb_tmp(4) = myPb_tmp(1) + (h/6._wp)*(mydPbdt_tmp(1) + mydPbdt_tmp(2) + 4._wp*mydPbdt_tmp(3))
             myMv_tmp(4) = myMv_tmp(1) + (h/6._wp)*(mydMvdt_tmp(1) + mydMvdt_tmp(2) + 4._wp*mydMvdt_tmp(3))
-            mydMvdt_tmp(4) = f_advance_EL(myR_tmp(4), myV_tmp(4), myPb_tmp(4), myMv_tmp(4), &
-                                          bub_id, fmass_n, fbeta_c, fbeta_t, mydPbdt_tmp(4))
+            call s_advance_EL(myR_tmp(4), myV_tmp(4), myPb_tmp(4), myMv_tmp(4), &
+                                          bub_id, fmass_n, fbeta_c, fbeta_t, mydPbdt_tmp(4), mydMvdt_tmp(4))
         end if
         myA_tmp(4) = f_rddot(fRho, fP, myR_tmp(4), myV_tmp(4), fR0, &
                              myPb_tmp(4), mydPbdt_tmp(4), alf, fntait, fBtait, &
@@ -780,20 +782,20 @@ contains
         !!  @param fMv_tmp Mass of vapor in the bubble
         !!  @param fdPbdt_tmp Rate of change of the internal bubble pressure
         !!  @param fdMvdt_tmp Rate of change of the mass of vapor in the bubble
-    function f_advance_EL(fR_tmp, fV_tmp, fPb_tmp, fMv_tmp, bub_id, &
-                          fmass_n, fbeta_c, fbeta_t, fdPbdt_tmp)
+    pure elemental subroutine s_advance_EL(fR_tmp, fV_tmp, fPb_tmp, fMv_tmp, bub_id, &
+                          fmass_n, fbeta_c, fbeta_t, fdPbdt_tmp, advance_EL)
         !$acc routine seq
         real(wp), intent(IN) :: fR_tmp, fV_tmp, fPb_tmp, fMv_tmp
         real(wp), intent(IN) :: fmass_n, fbeta_c, fbeta_t
         integer, intent(IN) :: bub_id
         real(wp), intent(INOUT) :: fdPbdt_tmp
+        real(wp), intent(out) :: advance_EL
+        real(wp) :: fVapFlux, myR_m, mygamma_m
 
-        real(wp) :: fVapFlux, f_advance_EL, myR_m, mygamma_m
-
-        fVapFlux = f_vflux(fR_tmp, fV_tmp, fPb_tmp, fMv_tmp, bub_id, fmass_n, fbeta_c, myR_m, mygamma_m)
+        call s_vflux(fR_tmp, fV_tmp, fPb_tmp, fMv_tmp, bub_id, fVapFlux, fmass_n, fbeta_c, myR_m, mygamma_m)
         fdPbdt_tmp = f_bpres_dot(fVapFlux, fR_tmp, fV_tmp, fPb_tmp, fMv_tmp, bub_id, fbeta_t, myR_m, mygamma_m)
-        f_advance_EL = 4._wp*pi*fR_tmp**2._wp*fVapFlux
+        advance_EL = 4._wp*pi*fR_tmp**2._wp*fVapFlux
 
-    end function f_advance_EL
+    end subroutine s_advance_EL
 
 end module m_bubbles
diff --git a/src/simulation/m_bubbles_EE.fpp b/src/simulation/m_bubbles_EE.fpp
@@ -92,10 +92,11 @@ contains
 
     end subroutine s_comp_alpha_from_n
 
-    subroutine s_compute_bubbles_EE_rhs(idir, q_prim_vf)
+    pure subroutine s_compute_bubbles_EE_rhs(idir, q_prim_vf, divu)
 
         integer, intent(in) :: idir
         type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
+        type(scalar_field), intent(inout) :: divu !< matrix for div(u)
 
         integer :: j, k, l
 
@@ -152,7 +153,7 @@ contains
         !!      that are needed for the bubble modeling
         !!  @param q_prim_vf Primitive variables
         !!  @param q_cons_vf Conservative variables
-    subroutine s_compute_bubble_EE_source(q_cons_vf, q_prim_vf, t_step, rhs_vf)
+    impure subroutine s_compute_bubble_EE_source(q_cons_vf, q_prim_vf, t_step, rhs_vf)
         type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
         type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
         integer, intent(in) :: t_step
@@ -271,8 +272,8 @@ contains
                         if (.not. polytropic) then
                             pb = q_prim_vf(ps(q))%sf(j, k, l)
                             mv = q_prim_vf(ms(q))%sf(j, k, l)
-                            call s_bwproperty(pb, q)
-                            vflux = f_vflux(myR, myV, pb, mv, q)
+                            call s_bwproperty(pb, q, chi_vw, k_mw, rho_mw)
+                            call s_vflux(myR, myV, pb, mv, q, vflux)
                             pbdot = f_bpres_dot(vflux, myR, myV, pb, mv, q)
 
                             bub_p_src(j, k, l, q) = nbub*pbdot
diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp
@@ -561,7 +561,7 @@ contains
             myR0 = bub_R0(k)
 
             ! Vapor and heat fluxes
-            myVapFlux = f_vflux(myR, myV, myPb, myMass_v, k, myMass_n, myBeta_c, myR_m, mygamma_m)
+            call s_vflux(myR, myV, myPb, myMass_v, k, myVapFlux, myMass_n, myBeta_c, myR_m, mygamma_m)
             myPbdot = f_bpres_dot(myVapFlux, myR, myV, myPb, myMass_v, k, myBeta_t, myR_m, mygamma_m)
             myMvdot = 4._wp*pi*myR**2._wp*myVapFlux
 
@@ -1017,7 +1017,7 @@ contains
     !>  This subroutine updates the Lagrange variables using the tvd RK time steppers.
         !!      The time derivative of the bubble variables must be stored at every stage to avoid precision errors.
         !! @param stage Current tvd RK stage
-    subroutine s_update_lagrange_tdv_rk(stage)
+    impure subroutine s_update_lagrange_tdv_rk(stage)
 
         integer, intent(in) :: stage
 
diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp
@@ -838,7 +838,7 @@ contains
             ! RHS additions for sub-grid bubbles_euler
             if (bubbles_euler) then
                 call nvtxStartRange("RHS-BUBBLES-COMPUTE")
-                call s_compute_bubbles_EE_rhs(id, q_prim_qp%vf)
+                call s_compute_bubbles_EE_rhs(id, q_prim_qp%vf, divu)
                 call nvtxEndRange
             end if
 
