diff --git a/Source/driver/Castro_util.H b/Source/driver/Castro_util.H
index 4202132493..e7b8751d12 100644
--- a/Source/driver/Castro_util.H
+++ b/Source/driver/Castro_util.H
@@ -157,10 +157,12 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE
 Real volume(const int& i, const int& j, const int& k,
             const GeometryData& geomdata)
 {
+    //
     // Given 3D cell-centered indices (i,j,k), return the volume of the zone.
-    // Note that Castro has no support for angular coordinates, so
-    // this function only provides Cartesian in 1D/2D/3D, Cylindrical (R-Z)
-    // in 2D, and Spherical in 1D.
+    // this function only provides Cartesian in 1D/2D/3D,
+    // Cylindrical (R-Z) in 2D,
+    // and Spherical in 1D and 2D (R-THETA).
+    //
 
     amrex::ignore_unused(i);
     amrex::ignore_unused(j);
@@ -210,8 +212,7 @@ Real volume(const int& i, const int& j, const int& k,
 
         vol = dx[0] * dx[1];
 
-    }
-    else {
+    } else if (coord == 1) {
 
         // Cylindrical
 
@@ -220,6 +221,20 @@ Real volume(const int& i, const int& j, const int& k,
 
         vol = M_PI * (r_l + r_r) * dx[0] * dx[1];
 
+    } else {
+
+        // Spherical
+
+        constexpr Real twoThirdsPi = 2.0_rt * M_PI / 3.0_rt;
+
+        Real r_l = geomdata.ProbLo()[0] + static_cast<Real>(i) * dx[0];
+        Real r_r = geomdata.ProbLo()[0] + static_cast<Real>(i+1) * dx[0];
+        Real theta_l = geomdata.ProbLo()[1] + static_cast<Real>(j) * dx[1];
+        Real theta_r = geomdata.ProbLo()[1] + static_cast<Real>(j+1) * dx[1];
+
+        vol = twoThirdsPi * (std::cos(theta_l) - std::cos(theta_r)) * dx[0] *
+            (r_r * r_r + r_r * r_l + r_l * r_l);
+
     }
 
 #else
@@ -290,8 +305,7 @@ Real area(const int& i, const int& j, const int& k,
             a = dx[0];
         }
 
-    }
-    else  {
+    } else  if (coord == 1) {
 
         // Cylindrical
 
@@ -304,6 +318,24 @@ Real area(const int& i, const int& j, const int& k,
             a = 2.0_rt * M_PI * r * dx[0];
         }
 
+    } else {
+
+        // Spherical
+
+        if (idir == 0) {
+            Real r = geomdata.ProbLo()[0] + static_cast<Real>(i) * dx[0];
+            Real theta_l = geomdata.ProbLo()[1] + static_cast<Real>(j) * dx[1];
+            Real theta_r = geomdata.ProbLo()[1] + static_cast<Real>(j+1) * dx[1];
+
+            a = 2.0_rt * M_PI * r * r * (std::cos(theta_l) - std::cos(theta_r));
+        }
+        else {
+            Real r = geomdata.ProbLo()[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0];
+            Real theta = geomdata.ProbLo()[1] + static_cast<Real>(j) * dx[1];
+
+            a = 2.0_rt * M_PI * std::sin(theta) * r * dx[0];
+        }
+
     }
 
 #else