Convert sphe to C++ (#2579)

Source/radiation/RAD_1D.F90

@@ -11,21 +11,6 @@ module rad_module
 
 contains
 
-subroutine sphe(r, s, n, &
-                DIMS(reg), dx) bind(C, name="sphe")
-  use amrex_fort_module, only : rt => amrex_real
-  implicit none
-  integer :: DIMDEC(reg)
-  real(rt)         :: r(reg_l1:reg_h1)
-  real(rt)         :: s(1)
-  integer :: n
-  real(rt)         :: dx(1)
-  integer :: i
-  do i = reg_l1, reg_h1
-     r(i) = r(i)**2
-  enddo
-end subroutine sphe
-
 subroutine rfface(fine, &
                   DIMS(fbox), &
                   crse, &

Source/radiation/RAD_2D.F90

@@ -11,38 +11,6 @@ module rad_module
 
 contains
 
-subroutine sphe(r, s, n, &
-                DIMS(reg), dx) bind(C, name="sphe")
-
-  use amrex_fort_module, only : rt => amrex_real
-  integer :: DIMDEC(reg)
-  real(rt)         :: r(reg_l1:reg_h1)
-  real(rt)         :: s(reg_l2:reg_h2)
-  integer :: n
-  real(rt)         :: dx(2)
-  real(rt)         :: h1, h2, d1, d2
-  integer :: i, j
-  if (n == 0) then
-     do i = reg_l1, reg_h1
-        r(i) = r(i)**2
-     enddo
-     h2 = 0.5e0_rt * dx(2)
-     d2 = 1.e0_rt / dx(2)
-     do j = reg_l2, reg_h2
-        s(j) = d2 * (cos(s(j) - h2) - cos(s(j) + h2))
-     enddo
-  else
-     h1 = 0.5e0_rt * dx(1)
-     d1 = 1.e0_rt / (3.e0_rt * dx(1))
-     do i = reg_l1, reg_h1
-        r(i) = d1 * ((r(i) + h1)**3 - (r(i) - h1)**3)
-     enddo
-     do j = reg_l2, reg_h2
-        s(j) = sin(s(j))
-     enddo
-  endif
-end subroutine sphe
-
 subroutine rfface(fine, &
                   DIMS(fbox), &
                   crse, &

Source/radiation/RAD_3D.F90

@@ -14,38 +14,6 @@ module rad_module
 ! The following routines implement metric terms in 2D and are included
 ! in the 3D source only to enable the code to link.
 
-subroutine sphe(r, s, n, &
-                DIMS(reg), dx) bind(C, name="sphe")
-
-  use amrex_fort_module, only : rt => amrex_real
-  integer :: DIMDEC(reg)
-  real(rt)         :: r(reg_l1:reg_h1)
-  real(rt)         :: s(reg_l2:reg_h2)
-  integer :: n
-  real(rt)         :: dx(2)
-  real(rt)         :: h1, h2, d1, d2
-  integer :: i, j
-  if (n == 0) then
-     do i = reg_l1, reg_h1
-        r(i) = r(i)**2
-     enddo
-     h2 = 0.5e0_rt * dx(2)
-     d2 = 1.e0_rt / dx(2)
-     do j = reg_l2, reg_h2
-        s(j) = d2 * (cos(s(j) - h2) - cos(s(j) + h2))
-     enddo
-  else
-     h1 = 0.5e0_rt * dx(1)
-     d1 = 1.e0_rt / (3.e0_rt * dx(1))
-     do i = reg_l1, reg_h1
-        r(i) = d1 * ((r(i) + h1)**3 - (r(i) - h1)**3)
-     enddo
-     do j = reg_l2, reg_h2
-        s(j) = sin(s(j))
-     enddo
-  endif
-end subroutine sphe
-
 subroutine rfface(fine, &
                   DIMS(fbox), &
                   crse, &

Source/radiation/RAD_F.H

@@ -126,9 +126,6 @@ extern "C"
 #ifdef __cplusplus
 extern "C" {
 #endif
-  void sphe(amrex::Real* r, amrex::Real* s, const int& idim, 
-            ARLIM_P(edgeboxlo), ARLIM_P(edgeboxhi), const amrex::Real* dx);
-
   void fkpn(const int* lo, const int* hi,
             BL_FORT_FAB_ARG_3D(fkp),       
             amrex::Real con, amrex::Real em, amrex::Real en,

Source/radiation/RadSolve.cpp

@@ -877,9 +877,18 @@ void RadSolve::levelDterm(int level, MultiFab& Dterm, MultiFab& Er, int igroup)
               parent->Geom(level).GetCellLoc(rc, reg, 0);
               parent->Geom(level).GetCellLoc(s, reg, 0);
               const Box &dbox = Dterm_face[0][fi].box();
-              sphe(re.dataPtr(), s.dataPtr(), 0,
-                   ARLIM(dbox.loVect()), ARLIM(dbox.hiVect()), dx.data());
-
+
+              for (int i = dbox.loVect()[0]; i <= dbox.hiVect()[0]; ++i) {
+                  re[i] *= re[i];
+              }
+#if AMREX_SPACEDIM >= 2
+              Real h2 = 0.5e0_rt * dx[1];
+              Real d2 = 1.e0_rt / dx[1];
+              for (int j = dbox.loVect()[1]; j <= dbox.hiVect()[1]; ++j) {
+                  s[j] = d2 * (std::cos(s[j] - h2) - std::cos(s[j] + h2));
+              }
+#endif
+
               ca_correct_dterm(D_DECL(BL_TO_FORTRAN(Dterm_face[0][fi]),
                                       BL_TO_FORTRAN(Dterm_face[1][fi]),
                                       BL_TO_FORTRAN(Dterm_face[2][fi])),
@@ -1160,8 +1169,31 @@ void RadSolve::getEdgeMetric(int idim, const Geometry& geom, const Box& edgebox,
         geom.GetEdgeLoc(s, reg, I);
       }
       const Real *dx = geom.CellSize();
-      sphe(r.dataPtr(), s.dataPtr(), idim,
-           ARLIM(edgebox.loVect()), ARLIM(edgebox.hiVect()), dx);
+
+      if (idim == 0) {
+          for (int i = edgebox.loVect()[0]; i <= edgebox.hiVect()[0]; ++i) {
+              r[i] *= r[i];
+          }
+#if AMREX_SPACEDIM >= 2
+          Real h2 = 0.5e0_rt * dx[1];
+          Real d2 = 1.e0_rt / dx[1];
+          for (int j = edgebox.loVect()[1]; j <= edgebox.hiVect()[1]; ++j) {
+              s[j] = d2 * (std::cos(s[j] - h2) - std::cos(s[j] + h2));
+          }
+#endif
+      }
+      else {
+          Real h1 = 0.5e0_rt * dx[0];
+          Real d1 = 1.e0_rt / (3.e0_rt * dx[0]);
+          for (int i = edgebox.loVect()[0]; i <= edgebox.hiVect()[0]; ++i) {
+              r[i] = d1 * (std::pow(r[i] + h1, 3) - std::pow(r[i] - h1, 3));
+          }
+#if AMREX_SPACEDIM >= 2
+          for (int j = edgebox.loVect()[1]; j <= edgebox.hiVect()[1]; ++j) {
+              s[j] = std::sin(s[j]);
+          }
+#endif
+      }
     }
 }