diff --git a/adflow/__init__.py b/adflow/__init__.py
index 172f6ad82..33baa2569 100644
--- a/adflow/__init__.py
+++ b/adflow/__init__.py
@@ -6,3 +6,6 @@
 from .pyADflow_C import ADFLOW_C
 from .oversetCheck import OversetCheck
 from .checkZipper import checkZipper
+
+
+from .actuatorRegion import UniformActuatorRegion, BSplineActuatorRegion
diff --git a/adflow/actuatorRegion.py b/adflow/actuatorRegion.py
new file mode 100644
index 000000000..10b588797
--- /dev/null
+++ b/adflow/actuatorRegion.py
@@ -0,0 +1,331 @@
+from abc import ABC, abstractmethod
+from typing import Tuple
+from types import ModuleType
+import os
+
+from mpi4py import MPI
+
+import numpy as np
+import numpy.typing as npt
+import scipy.interpolate as interpolate
+
+
+from baseclasses import AeroProblem
+
+
+class AbstractActuatorRegion(ABC):
+    def __init__(self, centerPoint: npt.NDArray, thrustVector: npt.NDArray, thrust: float, heat: float):
+        if centerPoint.shape != (3,):
+            raise ValueError('"centerPoint" must have shape "(3,)"')
+        if thrustVector.shape != (3,):
+            raise ValueError('"thrustVector" must have shape "(3,)"')
+        if np.linalg.norm(thrustVector) == 0:
+            raise ValueError('"trustVector" can not have a length of "0"')
+        if thrust < 0:
+            raise ValueError('"thrust" must not be smaller than 0.')
+        if heat < 0:
+            raise ValueError('"heat" must not be smaller than 0.')
+
+        self.centerPoint = centerPoint
+        self.thrustVector = thrustVector / np.linalg.norm(thrustVector)
+        self.iRegion = -1
+        self.nLocalCells = -1
+        self.familyName = ""
+
+        self._heatBcName = "Heat"
+        self._thrustBcName = "Thrust"
+
+        self._boundaryConditions = {
+            self._thrustBcName: thrust,
+            self._heatBcName: heat,
+        }
+
+    def setBCValue(self, bcName, bcValue):
+        if bcName not in self._boundaryConditions.keys():
+            raise ValueError(
+                f'"{bcName}" is not a valid Boundary Condition. Valid Conditions are: {list(self._boundaryConditions.keys())}'
+            )
+
+        self._boundaryConditions[bcName] = bcValue
+
+    def _computeScalingConstant(self, cellValues, totalValue, comm):
+        computedTotalvalue = comm.allreduce(np.sum(cellValues), op=MPI.SUM)
+
+        if np.isclose(computedTotalvalue, 0):
+            return 0.0
+
+        return totalValue / computedTotalvalue
+
+    @abstractmethod
+    def tagActiveCells(
+        self, distance2axis: npt.NDArray, distance2plane: npt.NDArray, tangent: npt.NDArray
+    ) -> npt.NDArray:
+        flags = np.zeros_like(distance2axis)
+        return flags
+
+    @abstractmethod
+    def computeCellForceVector(
+        self,
+        distance2axis: npt.NDArray,
+        distance2plane: npt.NDArray,
+        tangent: npt.NDArray,
+        cellVolume: npt.NDArray,
+        totalVolume: float,
+        comm: MPI.Intracomm,
+    ) -> npt.NDArray:
+        force = np.zeros((3, len(distance2axis)))
+        return force
+
+    @abstractmethod
+    def computeCellHeatVector(
+        self,
+        distance2axis: npt.NDArray,
+        distance2plane: npt.NDArray,
+        tangent: npt.NDArray,
+        cellVolume: npt.NDArray,
+        totalVolume: float,
+        comm: MPI.Intracomm,
+    ) -> npt.NDArray:
+        heat = np.zeros_like(distance2axis)
+        return heat
+
+
+class ActuatorRegionHandler:
+    def __init__(self, adflow_fortran: ModuleType, comm):
+        self._comm = comm
+        self._adflow_fortran = adflow_fortran
+
+        self._actuatorRegions = list()
+
+    def addActuatorRegion(
+        self,
+        actuatorRegion: AbstractActuatorRegion,
+        familyName: str,
+        familyID: int,
+        relaxStart: float,
+        relaxEnd: float,
+    ):
+        # compute the distance of each cell to the AR plane and axis
+        ncells = self._adflow_fortran.adjointvars.ncellslocal[0]
+
+        distance2plane, distance2axis, tangent = self._adflow_fortran.actuatorregion.computeinitialspatialmetrics(
+            actuatorRegion.centerPoint, actuatorRegion.thrustVector, ncells
+        )
+        tangent = tangent.T
+
+        # tag the active cells
+        flag = actuatorRegion.tagActiveCells(distance2axis, distance2plane, tangent)
+        iRegion, nLocalCells = self._adflow_fortran.actuatorregion.addactuatorregion(
+            flag,
+            familyName,
+            familyID,
+            relaxStart,
+            relaxEnd,
+        )
+        actuatorRegion.iRegion = iRegion
+        actuatorRegion.nLocalCells = nLocalCells
+        actuatorRegion.familyName = familyName
+
+        # book keep the new region
+        self._actuatorRegions.append(actuatorRegion)
+
+    def updateActuatorRegionsBC(self, aeroproblem: AeroProblem):
+        for actuatorRegion in self._actuatorRegions:
+
+            # set BC values
+            for (bcName, familyName), bcValue in aeroproblem.bcVarData.items():
+                if familyName != actuatorRegion.familyName:
+                    continue
+                actuatorRegion.setBCValue(bcName, bcValue)
+
+            # compute local metrics
+            distance2plane, distance2axis, tangent, volume = self._adflow_fortran.actuatorregion.computespatialmetrics(
+                actuatorRegion.iRegion,
+                actuatorRegion.nLocalCells,
+                actuatorRegion.centerPoint,
+                actuatorRegion.thrustVector,
+            )
+            tangent = tangent.T
+
+            totalVolume = self._comm.allreduce(np.sum(volume), op=MPI.SUM)
+
+            # compute the source terms
+            force = actuatorRegion.computeCellForceVector(
+                distance2axis, distance2plane, tangent, volume, totalVolume, self._comm
+            )
+            heat = actuatorRegion.computeCellHeatVector(
+                distance2axis, distance2plane, tangent, volume, totalVolume, self._comm
+            )
+
+            # skip this proc if no cells are active
+            if actuatorRegion.nLocalCells == 0:
+                continue
+
+            # apply the source terms in Fortran
+            self._adflow_fortran.actuatorregion.populatebcvalues(
+                actuatorRegion.iRegion,
+                actuatorRegion.nLocalCells,
+                force.T,
+                heat.T,
+            )
+
+    def writeActuatorRegions(self, fileName: str, outputDir: str):
+        # Join to get the actual filename root
+        fileName = os.path.join(outputDir, fileName)
+
+        # Ensure extension is .plt even if the user didn't specify
+        fileName, ext = os.path.splitext(fileName)
+        fileName += ".plt"
+
+        # just call the underlying fortran routine
+        self._adflow_fortran.actuatorregion.writeactuatorregions(fileName)
+
+
+class CircularActuatorRegion(AbstractActuatorRegion):
+    def __init__(
+        self,
+        centerPoint: npt.NDArray,
+        thrustVector: npt.NDArray,
+        innerDiameter: float,
+        outerDiameter: float,
+        regionDepth: float,
+        thrust: float = 0,
+        heat: float = 0,
+    ):
+        super().__init__(centerPoint, thrustVector, thrust, heat)
+
+        if regionDepth <= 0:
+            raise ValueError('"depth" must be greater than 0.')
+        if innerDiameter < 0:
+            raise ValueError('"innerDiameter" must not be smaller than 0.')
+        if outerDiameter <= 0:
+            raise ValueError('"outerDiameter" must be greater than 0.')
+        if outerDiameter <= innerDiameter:
+            raise ValueError('"outerDiameter" must be greater than "innerDiameter".')
+
+        self._innerDiameter = innerDiameter
+        self._outerDiameter = outerDiameter
+        self._regionDepth = regionDepth
+
+    def tagActiveCells(self, distance2axis, distance2plane, tangent):
+        flags = np.zeros_like(distance2axis)
+
+        indices = np.logical_and(
+            np.logical_and(
+                distance2axis >= self._innerDiameter / 2,
+                distance2axis <= self._outerDiameter / 2,
+            ),
+            np.logical_and(distance2plane >= 0, distance2plane <= self._regionDepth),
+        )
+
+        flags[indices] = 1
+
+        return flags
+
+
+class UniformActuatorRegion(CircularActuatorRegion):
+    def computeCellForceVector(self, distance2axis, distance2plane, tangent, cellVolume, totalVolume, comm):
+        thrustBC = self._boundaryConditions[self._thrustBcName]
+
+        cellThrusts = thrustBC * cellVolume / totalVolume
+        scaling_constant = self._computeScalingConstant(cellThrusts, thrustBC, comm)
+        force = np.outer(scaling_constant * cellThrusts, -self.thrustVector)
+
+        return force
+
+    def computeCellHeatVector(self, distance2axis, distance2plane, tangent, cellVolume, totalVolume, comm):
+        heatBC = self._boundaryConditions[self._heatBcName]
+
+        cellHeats = heatBC * cellVolume / totalVolume
+        scaling_constant = self._computeScalingConstant(cellHeats, heatBC, comm)
+
+        return cellHeats * scaling_constant
+
+
+class BSplineActuatorRegion(CircularActuatorRegion):
+    def __init__(
+        self,
+        centerPoint: npt.NDArray,
+        thrustVector: npt.NDArray,
+        innerDiameter: float,
+        outerDiameter: float,
+        regionDepth: float,
+        thrustDistribution: Tuple[npt.NDArray, npt.NDArray, float] = (
+            np.array([0.0, 0.0, 1.0, 1.0]),
+            np.array([0.0, 0.0, 0.0, 0.0]),
+            1,
+        ),
+        tangentDistribution: Tuple[npt.NDArray, npt.NDArray, float] = (
+            np.array([0.0, 0.0, 1.0, 1.0]),
+            np.array([0.0, 0.0, 0.0, 0.0]),
+            1,
+        ),
+        radialDistribution: Tuple[npt.NDArray, npt.NDArray, float] = (
+            np.array([0.0, 0.0, 1.0, 1.0]),
+            np.array([0.0, 0.0, 0.0, 0.0]),
+            1,
+        ),
+        heatDistribution: Tuple[npt.NDArray, npt.NDArray, float] = (
+            np.array([0.0, 0.0, 1.0, 1.0]),
+            np.array([0.0, 0.0, 0.0, 0.0]),
+            1,
+        ),
+        thrust: float = 0,
+        heat: float = 0,
+    ):
+        super().__init__(centerPoint, thrustVector, innerDiameter, outerDiameter, regionDepth, thrust, heat)
+
+        self._thrustSpline = interpolate.BSpline(*thrustDistribution)
+        self._tangentSpline = interpolate.BSpline(*tangentDistribution)
+        self._radialSpline = interpolate.BSpline(*radialDistribution)
+        self._heatSpline = interpolate.BSpline(*heatDistribution)
+
+    def _formSpline(self, distribution: Tuple[npt.NDArray, npt.NDArray, float], scalingConstant: float):
+        t = distribution[0]
+        c = distribution[1]
+        k = distribution[2]
+
+        return interpolate.BSpline(t, c * scalingConstant, k)
+
+    def _computeNormalizedRadius(self, distance2axis):
+        normalizedRadius = distance2axis - self._innerDiameter / 2
+        normalizedRadius /= self._outerDiameter / 2 - self._innerDiameter / 2
+        return normalizedRadius
+
+    def computeCellForceVector(self, distance2axis, distance2plane, tangent, cellVolume, totalVolume, comm):
+        thrustBC = self._boundaryConditions[self._thrustBcName]
+
+        normalizedRadius = self._computeNormalizedRadius(distance2axis)
+        thrustFactor = thrustBC * cellVolume / totalVolume
+
+        # add axial contribution (thrust)
+        cellThrusts = self._thrustSpline(normalizedRadius) * thrustFactor
+
+        # compute a scaling constant such that the summ of the thrust equals the prescribed thrust
+        scaling_constant = self._computeScalingConstant(cellThrusts, thrustBC, comm)
+        force = np.outer(scaling_constant * cellThrusts, -self.thrustVector)
+
+        # add tangential contribution (swirl)
+        force += np.multiply(
+            tangent.T, np.array([scaling_constant * self._tangentSpline(normalizedRadius) * thrustFactor])
+        ).T
+
+        # add radial contribution
+        radius = np.cross(self.thrustVector, tangent)
+        force += np.multiply(
+            radius.T, np.array([scaling_constant * self._radialSpline(normalizedRadius) * 2 * thrustFactor])
+        ).T
+
+        return force
+
+    def computeCellHeatVector(self, distance2axis, distance2plane, tangent, cellVolume, totalVolume, comm):
+        heatBC = self._boundaryConditions[self._heatBcName]
+
+        normalizedRadius = self._computeNormalizedRadius(distance2axis)
+
+        # compute a scaling constant such that the summ of the heat equals the prescribed heat
+
+        cellHeats = self._heatSpline(normalizedRadius) * heatBC * cellVolume / totalVolume
+        scaling_constant = self._computeScalingConstant(cellHeats, heatBC, comm)
+
+        return cellHeats * scaling_constant
diff --git a/adflow/pyADflow.py b/adflow/pyADflow.py
index 9f11d44ce..861bdaf11 100644
--- a/adflow/pyADflow.py
+++ b/adflow/pyADflow.py
@@ -33,6 +33,7 @@
 from . import MExt
 import hashlib
 from collections import OrderedDict
+from .actuatorRegion import ActuatorRegionHandler, AbstractActuatorRegion
 
 
 class ADFLOWWarning(object):
@@ -202,6 +203,8 @@ def __init__(self, comm=None, options=None, debug=False, dtype="d"):
         # had we read in a param file
         self.adflow.iteration.deforming_grid = True
 
+        self._actuatorRegionHandler = ActuatorRegionHandler(self.adflow, self.comm)
+
         # In order to properly initialize we need to have mach number
         # and a few other things set. Just create a dummy aeroproblem,
         # use it, and then it will be deleted.
@@ -832,115 +835,15 @@ def addIntegrationSurface(self, fileName, familyName, isInflow=True, coordXfer=N
 
     def addActuatorRegion(
         self,
-        fileName,
-        axis1,
-        axis2,
+        actuatorRegion: AbstractActuatorRegion,
         familyName,
-        thrust=0.0,
-        torque=0.0,
-        heat=0.0,
-        relaxStart=None,
-        relaxEnd=None,
-        coordXfer=None,
+        relaxStart=-1,
+        relaxEnd=-1,
     ):
-        """
-        Add an actuator disk zone defined by the supplied closed surface in the
-        plot3d file "fileName". This surface defines the physical extent of the
-        region over which to apply the source terms. The plot3d file may be
-        multi-block but all the surface normals must point outside and no
-        additional surfaces can be inside. Internally, we find all of the CFD
-        volume cells that have their center inside this closed surface and
-        apply the source terms over these cells. This surface is only used with
-        the original mesh coordinates to mark the internal CFD cells, and we
-        keep using these cells even if the geometric design and the mesh
-        coordinates change. For example, the marked cells can lie outside of
-        the original closed surface after a large design change, but we will
-        still use the cells that were inside the surface with the baseline
-        design.
-
-        axis1 and axis2 define the vector that we use to determine the
-        direction of the thrust addition. Internally, we compute a
-        vector by $axisVec = axis2-axis1$ and then we normalize this
-        vector. When adding the thrust terms, the direction of the thrust
-        is obtained by multiplying the thrust magnitude by this vector.
-
-        Optionally, the source terms in the actuator zone can be
-        gradually ramped up as the solution converges. This continuation
-        approach can be more robust but the users should be careful with
-        the side effects of partially converged solutions. This behavior
-        can be controlled by relaxStart and relaxEnd parameters. By
-        default, the full magnitude of the source terms are added.
-        relaxStart controls when the continuation starts ramping up.
-        The value represents the relative convergence in a log10 basis.
-        So relaxStart = 2 means the source terms will be inactive until
-        the residual is decreased by a factor of 100 compared to free
-        stream conditions. The source terms are ramped to the full
-        magnitude at relaxEnd. E.g., a relaxEnd value of 4 would
-        result in the full source terms after a relative reduction of
-        1e4 in the total residuals. If relaxStart is not provided, but
-        relaxEnd is provided, then the relaxStart is assumed to be 0.
-        If both are not provided, we do not do ramping and just activate
-        the full source terms from the beginning. When this continuation
-        is used, we internally ramp up the magnitude of the source terms
-        monotonically to prevent oscillations; i.e., decrease in the total
-        residuals increase the source term magnitudes, but an increase
-        in the residuals do not reduce the source terms back down.
-
-        Parameters
-        ----------
-
-        fileName : str
-           Surface Plot 3D file (multiblock ascii) defining the closed
-           region over which the integration is to be applied.
-
-        axis1 : numpy array, length 3
-           The physical location of the start of the axis
-
-        axis2 : numpy array, length 4
-           The physical location of the end of the axis
-
-        familyName : str
-           The name to be associated with the functions defined on this
-           region.
-
-        thrust : scalar, float
-           The total amount of axial force to apply to this region, in the direction
-           of axis1 -> axis2
-
-        torque : scalar, float
-           The total amount of torque to apply to the region, about the
-           specified axis.
-
-        heat : scalar, float
-            The total amount of head added in the actuator zone with source terms
-
-        relaxStart : scalar, float
-            The start of the relaxation in terms of
-            orders of magnitude of relative convergence
-
-        relaxEnd : scalar, float
-            The end of the relaxation in terms of
-            orders of magnitude of relative convergence
-
-        coordXfer : function
-            A callback function that performs a coordinate transformation
-            between the original reference frame and any other reference
-            frame relevant to the current CFD case. This allows user to apply
-            arbitrary modifications to the loaded plot3d surface. The call
-            signature is documented in DVGeometry's :meth:`addPointset <pygeo:pygeo.parameterization.DVGeo.DVGeometry.addPointSet>` method.
-
-        """
         # ActuatorDiskRegions cannot be used in timeSpectralMode
         if self.getOption("equationMode").lower() == "time spectral":
             raise Error("ActuatorRegions cannot be used in Time Spectral Mode.")
 
-        # Load in the user supplied surface file.
-        pts, conn = self._readPlot3DSurfFile(fileName, convertToTris=False, coordXfer=coordXfer)
-
-        # We should do a topology check to ensure that the surface the
-        # user supplied is actually closed.
-        # self._closedTopoCheck(pts, conn)
-
         # Check if the family name given is already used for something
         # else.
         if familyName.lower() in self.families:
@@ -955,25 +858,16 @@ def addActuatorRegion(
         for fam in self.families:
             if len(self.families[fam]) > 0:
                 maxInd = max(maxInd, numpy.max(self.families[fam]))
-        famID = maxInd + 1
-        self.families[familyName.lower()] = [famID]
-
-        if relaxStart is None and relaxEnd is None:
-            # No relaxation at all
-            relaxStart = -1.0
-            relaxEnd = -1.0
+        familyID = maxInd + 1
+        self.families[familyName.lower()] = [familyID]
 
-        if relaxStart is None and relaxEnd is not None:
-            # Start at 0 orders if start is not given
+        if relaxStart < 0 and relaxEnd > 0:
             relaxStart = 0.0
 
-        if relaxEnd is None and relaxStart is not None:
-            raise Error("relaxEnd must be given is relaxStart is specified")
+        if relaxEnd < 0 and relaxStart >= 0:
+            raise ValueError("relaxEnd must be given if relaxStart is specified")
 
-        # Now continue to fortran were we setup the actual actuator region.
-        self.adflow.actuatorregion.addactuatorregion(
-            pts.T, conn.T, axis1, axis2, familyName, famID, thrust, torque, heat, relaxStart, relaxEnd
-        )
+        self._actuatorRegionHandler.addActuatorRegion(actuatorRegion, familyName, familyID, relaxStart, relaxEnd)
 
     def writeActuatorRegions(self, fileName, outputDir=None):
         """
@@ -995,15 +889,7 @@ def writeActuatorRegions(self, fileName, outputDir=None):
         if outputDir is None:
             outputDir = self.getOption("outputDirectory")
 
-        # Join to get the actual filename root
-        fileName = os.path.join(outputDir, fileName)
-
-        # Ensure extension is .plt even if the user didn't specify
-        fileName, ext = os.path.splitext(fileName)
-        fileName += ".plt"
-
-        # just call the underlying fortran routine
-        self.adflow.actuatorregion.writeactuatorregions(fileName)
+        self._actuatorRegionHandler.writeActuatorRegions(fileName, outputDir)
 
     def addUserFunction(self, funcName, functions, callBack):
         """Add a new function to ADflow by combining existing functions in a
@@ -3696,6 +3582,8 @@ def _setAeroProblemData(self, aeroProblem, firstCall=False):
         if not empty:
             self.adflow.bcdata.setbcdata(nameArray, dataArray, groupArray, 1)
 
+        self._actuatorRegionHandler.updateActuatorRegionsBC(AP)
+
         if not firstCall:
             self.adflow.initializeflow.updatebcdataalllevels()
 
diff --git a/adflow/pyPerformanceClass.py b/adflow/pyPerformanceClass.py
index 5bdeb8bf1..e0f365bbb 100644
--- a/adflow/pyPerformanceClass.py
+++ b/adflow/pyPerformanceClass.py
@@ -29,7 +29,6 @@
 
 
 class PERFORMANCE(object):
-
     """
     Basic Performance class to handle dynamic handling qualities constraints
     """
diff --git a/adflow/pyWeightAndBalance.py b/adflow/pyWeightAndBalance.py
index 600df7135..4e81cb567 100644
--- a/adflow/pyWeightAndBalance.py
+++ b/adflow/pyWeightAndBalance.py
@@ -48,7 +48,6 @@
 
 
 class WEIGHTANDBALANCE(Base):
-
     """
     Implementation of basic preliminary level CG and inertia
     estimation methods.
@@ -405,12 +404,7 @@ def calculateWingInertias(self, acg, xcg):
                         # print 'centroids',Xs,Ys,Zs3
 
                         Ix = (
-                            Ix
-                            + I1x
-                            - W * (Ys_dot**2)
-                            - W * (Zs3**2)
-                            + W * (Ys_dot + Ysoff) ** 2
-                            + W * (Zs3 + Zs1) ** 2
+                            Ix + I1x - W * (Ys_dot**2) - W * (Zs3**2) + W * (Ys_dot + Ysoff) ** 2 + W * (Zs3 + Zs1) ** 2
                         )
                         # print 'w',W,(Ys_dot),W*(Ys_dot**2)
                         # Ix = Ix-W*(Ys_dot**2)#-W*(Zs3**2)+W*(Ys_dot+Ysoff)**2 + W*(Zs3+Zs1)**2
diff --git a/adflow/pyWingCG.py b/adflow/pyWingCG.py
index 9f8f0b751..feaf39232 100644
--- a/adflow/pyWingCG.py
+++ b/adflow/pyWingCG.py
@@ -367,9 +367,7 @@ def calculateWingInertias(acg):
                     # print 'zs3',Zs3
                     Zs4 = acg[i][j].z_Centroid
 
-                    Ix = (
-                        Ix + I1x - W * (Ys_dot**2) - W * (Zs3**2) + W * (Ys_dot + Ysoff) ** 2 + W * (Zs3 + Zs1) ** 2
-                    )
+                    Ix = Ix + I1x - W * (Ys_dot**2) - W * (Zs3**2) + W * (Ys_dot + Ysoff) ** 2 + W * (Zs3 + Zs1) ** 2
                     Iy = Iy + I1y - W * (Xs**2) - W * (Zs3**2) + W * (Xs + Xs4) ** 2 + W * (Zs3 + Zs1) ** 2
                     Iz = Iz + I1z - W * (Xs**2 + Ys_dot**2) + W * (Xs + Xs4) ** 2 + W * (Ys_dot + Ysoff) ** 2
                     # Ixz = ...
diff --git a/src/adjoint/Makefile_tapenade b/src/adjoint/Makefile_tapenade
index d0077a2ea..678ba2d62 100644
--- a/src/adjoint/Makefile_tapenade
+++ b/src/adjoint/Makefile_tapenade
@@ -104,9 +104,6 @@ bcData%BCDataSubsonicOutflow(bcVarArray, Pref) > \
 bcData%BCDataSupersonicInflow(bcVarArray, muRef, rhoRef, Pref, uRef, wInf, pInfCorr) > \
 	                     (bcData%rho, bcData%velx, bcData%vely, bcData%velz, bcData%ps, bcData%turbInlet, muRef, rhoRef, Pref, uRef, wInf, pInfCorr) \
 \
-actuatorRegion%computeActuatorRegionVolume(vol, actuatorRegions%volLocal) > \
-										 (vol, actuatorRegions%volLocal)\
-\
 adjointExtra%xhalo_block(x) > \
                         (x) \
 \
@@ -183,8 +180,8 @@ sa%saViscous(w, vol, si, sj, sk, rlv, scratch) > \
 sa%saResScale(scratch, dw) > \
              (dw) \
 \
-residuals%sourceTerms_block(w, pref, uref, plocal, dw, vol, actuatorRegions%force, actuatorRegions%heat, actuatorRegions%volume) > \
-                           (w, pref, uref, plocal, dw, vol, actuatorRegions%force, actuatorRegions%heat, actuatorRegions%volume) \
+residuals%sourceTerms_block(w, pref, uref, plocal, dw, vol, actuatorRegions%force, actuatorRegions%heat) > \
+                           (w, pref, uref, plocal, dw, vol, actuatorRegions%force, actuatorRegions%heat) \
 \
 residuals%initres_block(dw, fw, flowDoms%w, flowDoms%vol, dscalar, dvector) > \
                        (dw, fw, flowDoms%w, flowDoms%vol, dscalar, dvector) \
diff --git a/src/adjoint/adjointDebug.F90 b/src/adjoint/adjointDebug.F90
index 24119b4d8..ec0f25070 100644
--- a/src/adjoint/adjointDebug.F90
+++ b/src/adjoint/adjointDebug.F90
@@ -559,7 +559,6 @@ subroutine printADSeeds(nn, level, sps)
         ! And the reverse seeds in the actuator zones
         do i = 1, nActuatorRegions
             write (*, *) 'actuatorRegionsd(i)%Force ', actuatorRegionsd(i)%force
-            write (*, *) 'actuatorRegionsd(i)%Torque ', actuatorRegionsd(i)%torque
         end do
 
     end subroutine printADSeeds
diff --git a/src/adjoint/adjointUtils.F90 b/src/adjoint/adjointUtils.F90
index c497ddec8..8fb30fa4e 100644
--- a/src/adjoint/adjointUtils.F90
+++ b/src/adjoint/adjointUtils.F90
@@ -1075,9 +1075,7 @@ subroutine zeroADSeeds(nn, level, sps)
         ! And the reverse seeds in the actuator zones
         do i = 1, nActuatorRegions
             actuatorRegionsd(i)%force = zero
-            actuatorRegionsd(i)%torque = zero
             actuatorRegionsd(i)%heat = zero
-            actuatorRegionsd(i)%volume = zero
         end do
 
     end subroutine zeroADSeeds
diff --git a/src/adjoint/masterRoutines.F90 b/src/adjoint/masterRoutines.F90
index 050bfbbd2..f0f450fe4 100644
--- a/src/adjoint/masterRoutines.F90
+++ b/src/adjoint/masterRoutines.F90
@@ -39,7 +39,6 @@ subroutine master(useSpatial, famLists, funcValues, forces, &
         use oversetCommUtilities, only: updateOversetConnectivity
         use actuatorRegionData, only: nActuatorRegions, actuatorRegions
         use wallDistanceData, only: xSurfVec, xSurf
-        use actuatorRegion, only: computeActuatorRegionVolume
 
         implicit none
 
@@ -90,13 +89,6 @@ subroutine master(useSpatial, famLists, funcValues, forces, &
             end do
         end if
 
-        if (useSpatial) then
-            ! Zero out the local volume pointers for the actuator zone
-            do iRegion = 1, nActuatorRegions
-                actuatorRegions(iRegion)%volLocal = zero
-            end do
-        end if
-
         do sps = 1, nTimeIntervalsSpectral
             do nn = 1, nDom
                 call setPointers(nn, 1, sps)
@@ -108,11 +100,6 @@ subroutine master(useSpatial, famLists, funcValues, forces, &
 
                     call volume_block
 
-                    ! Compute the volume of each actuator region
-                    do iRegion = 1, nActuatorRegions
-                        call computeActuatorRegionVolume(nn, iRegion)
-                    end do
-
                     call metric_block
                     call boundaryNormals
 
@@ -143,16 +130,6 @@ subroutine master(useSpatial, famLists, funcValues, forces, &
             end do
         end do
 
-        ! Sum the local actuator zone volumes into the global actuator
-        ! zone volumes.
-        if (useSpatial) then
-            do iRegion = 1, nActuatorRegions
-                call mpi_allreduce(actuatorRegions(iRegion)%volLocal, actuatorRegions(iRegion)%volume, 1, &
-                                   adflow_real, MPI_SUM, adflow_comm_world, ierr)
-                call ECHK(ierr, __FILE__, __LINE__)
-            end do
-        end if
-
         ! Exchange values
         call whalo2(currentLevel, 1_intType, nw, .True., .True., .True.)
 
@@ -299,7 +276,6 @@ subroutine master_d(wdot, xdot, forcesDot, dwDot, famLists, funcValues, funcValu
         use inputOverset, only: oversetUpdateMode
         use oversetCommUtilities, only: updateOversetConnectivity_d
         use actuatorRegionData, only: nActuatorRegions, actuatorRegionsd
-        use actuatorregion_d, only: computeactuatorregionvolume_d
 #include <petsc/finclude/petsc.h>
         use petsc
         implicit none
@@ -409,11 +385,6 @@ subroutine master_d(wdot, xdot, forcesDot, dwDot, famLists, funcValues, funcValu
             call setBCDataFineGrid_d(.true.)
         end if
 
-        ! Zero out the local volume seeds for the actuator zones
-        do iRegion = 1, nActuatorRegions
-            actuatorRegionsd(iRegion)%volLocal = zero
-        end do
-
         do sps = 1, nTimeIntervalsSpectral
             do nn = 1, nDom
 
@@ -433,11 +404,6 @@ subroutine master_d(wdot, xdot, forcesDot, dwDot, famLists, funcValues, funcValu
                 call volume_block_d()
                 call metric_block_d()
 
-                ! Loop over the actuator regions to compute the local
-                ! volume seeds
-                do iRegion = 1, nActuatorRegions
-                    call computeActuatorRegionVolume_d(nn, iRegion)
-                end do
                 call boundaryNormals_d()
 
                 time = timeunsteadyrestart
@@ -481,14 +447,6 @@ subroutine master_d(wdot, xdot, forcesDot, dwDot, famLists, funcValues, funcValu
             end do
         end do
 
-        ! Loop over the actuator regions again to sum the local volume
-        ! seeds into the global volume seeds
-        do iRegion = 1, nActuatorRegions
-            call mpi_allreduce(actuatorRegionsd(iRegion)%volLocal, actuatorRegionsd(iRegion)%volume, 1, &
-                               adflow_real, MPI_SUM, adflow_comm_world, ierr)
-            call ECHK(ierr, __FILE__, __LINE__)
-        end do
-
         ! Just exchange the derivative values.
         call whalo2_d(1, 1, nw, .True., .True., .True.)
 
@@ -662,7 +620,6 @@ subroutine master_b(wbar, xbar, extraBar, forcesBar, dwBar, nState, famLists, &
         use oversetCommUtilities, only: updateOversetConnectivity_b
         use BCRoutines, only: applyAllBC_block
         use actuatorRegionData, only: nActuatorRegions, actuatorRegionsd
-        use actuatorregion_b, only: computeactuatorregionvolume_b
         use monitor, only: timeUnsteadyRestart
         use section, only: sections, nSections ! used in time-declaration
 
@@ -803,16 +760,6 @@ subroutine master_b(wbar, xbar, extraBar, forcesBar, dwBar, nState, famLists, &
             end do domainLoop1
         end do spsLoop1
 
-        ! All reduce the global AZ volume seeds and store them in the
-        ! local AZ volume seeds.
-        ! This is the inverse of the all reduce that is done in forward
-        ! mode to sum the local seeds into the global seeds.
-        do iRegion = 1, nActuatorRegions
-            call mpi_allreduce(actuatorRegionsd(iRegion)%volume, actuatorRegionsd(iRegion)%volLocal, 1, &
-                               adflow_real, mpi_sum, adflow_comm_world, ierr)
-            call ECHK(ierr, __FILE__, __LINE__)
-        end do
-
         ! Need to re-apply the BCs. The reason is that BC halos behind
         ! interpolated cells need to be recomputed with their new
         ! interpolated values from actual compute cells. Only needed for
@@ -889,9 +836,6 @@ subroutine master_b(wbar, xbar, extraBar, forcesBar, dwBar, nState, famLists, &
                 call gridvelocitiesfinelevel_block_b(useoldcoor, time, sps, nn)
 
                 call boundaryNormals_b
-                do iRegion = 1, nActuatorRegions
-                    call computeActuatorRegionVolume_b(nn, iRegion)
-                end do
                 call metric_block_b
                 call volume_block_b
 
@@ -917,11 +861,6 @@ subroutine master_b(wbar, xbar, extraBar, forcesBar, dwBar, nState, famLists, &
             end if
         end do spsLoop2
 
-        ! Zero out the local volume seeds of the actuator zone
-        do iRegion = 1, nActuatorRegions
-            actuatorRegionsd(iRegion)%volLocal = zero
-        end do
-
         if (present(bcDataNames)) then
             allocate (bcDataValuesdLocal(size(bcDataValuesd)))
             bcDataValuesdLocal = zero
diff --git a/src/adjoint/outputForward/actuatorRegion_d.f90 b/src/adjoint/outputForward/actuatorRegion_d.f90
index 0fd8c8992..608b96aba 100644
--- a/src/adjoint/outputForward/actuatorRegion_d.f90
+++ b/src/adjoint/outputForward/actuatorRegion_d.f90
@@ -8,52 +8,48 @@ module actuatorregion_d
   implicit none
 
 contains
-!  differentiation of computeactuatorregionvolume in forward (tangent) mode (with options i4 dr8 r8):
-!   variations   of useful results: actuatorregions.vollocal
-!   with respect to varying inputs: *vol actuatorregions.vollocal
-!   rw status of diff variables: *vol:in actuatorregions.vollocal:in-out
-!   plus diff mem management of: vol:in
-  subroutine computeactuatorregionvolume_d(nn, iregion)
-    use blockpointers, only : ndom, vol, vold
+  subroutine computecellspatialmetrics(i, j, k, centerpoint, &
+&   thrustvector, distance2plane, distance2axis, tangent)
+    use constants
+    use blockpointers, only : x
     implicit none
 ! inputs
-    integer(kind=inttype), intent(in) :: nn, iregion
+    real(kind=realtype), dimension(3), intent(in) :: centerpoint, &
+&   thrustvector
+    integer(kind=inttype), intent(in) :: i, j, k
+! outputs
+    real(kind=realtype), intent(out) :: distance2axis, distance2plane
+    real(kind=realtype), dimension(3), intent(out) :: tangent
 ! working
-    integer(kind=inttype) :: iii
-    integer(kind=inttype) :: i, j, k
-! loop over the region for this block
-    do iii=actuatorregions(iregion)%blkptr(nn-1)+1,actuatorregions(&
-&       iregion)%blkptr(nn)
-      i = actuatorregions(iregion)%cellids(1, iii)
-      j = actuatorregions(iregion)%cellids(2, iii)
-      k = actuatorregions(iregion)%cellids(3, iii)
-! sum the volume of each cell within the region on this proc
-      actuatorregionsd(iregion)%vollocal = actuatorregionsd(iregion)%&
-&       vollocal + vold(i, j, k)
-      actuatorregions(iregion)%vollocal = actuatorregions(iregion)%&
-&       vollocal + vol(i, j, k)
-    end do
-  end subroutine computeactuatorregionvolume_d
-
-  subroutine computeactuatorregionvolume(nn, iregion)
-    use blockpointers, only : ndom, vol
-    implicit none
-! inputs
-    integer(kind=inttype), intent(in) :: nn, iregion
-! working
-    integer(kind=inttype) :: iii
-    integer(kind=inttype) :: i, j, k
-! loop over the region for this block
-    do iii=actuatorregions(iregion)%blkptr(nn-1)+1,actuatorregions(&
-&       iregion)%blkptr(nn)
-      i = actuatorregions(iregion)%cellids(1, iii)
-      j = actuatorregions(iregion)%cellids(2, iii)
-      k = actuatorregions(iregion)%cellids(3, iii)
-! sum the volume of each cell within the region on this proc
-      actuatorregions(iregion)%vollocal = actuatorregions(iregion)%&
-&       vollocal + vol(i, j, k)
-    end do
-  end subroutine computeactuatorregionvolume
+    real(kind=realtype) :: thrustvectornorm, dotproduct, tangentnorm
+    real(kind=realtype), dimension(3) :: xcen, distance2center, &
+&   rawtangent
+    intrinsic sqrt
+    real(kind=realtype) :: arg1
+! compute the cell center
+    xcen = eighth*(x(i-1, j-1, k-1, :)+x(i, j-1, k-1, :)+x(i-1, j, k-1, &
+&     :)+x(i, j, k-1, :)+x(i-1, j-1, k, :)+x(i, j-1, k, :)+x(i-1, j, k, &
+&     :)+x(i, j, k, :))
+    arg1 = thrustvector(1)**2 + thrustvector(2)**2 + thrustvector(3)**2
+    thrustvectornorm = sqrt(arg1)
+    distance2center = xcen - centerpoint
+! compute distance to plane
+    dotproduct = thrustvector(1)*distance2center(1) + thrustvector(2)*&
+&     distance2center(2) + thrustvector(3)*distance2center(3)
+    distance2plane = dotproduct/thrustvectornorm
+! compute distance to axis
+    rawtangent(1) = distance2center(2)*thrustvector(3) - distance2center&
+&     (3)*thrustvector(2)
+    rawtangent(2) = distance2center(3)*thrustvector(1) - distance2center&
+&     (1)*thrustvector(3)
+    rawtangent(3) = distance2center(1)*thrustvector(2) - distance2center&
+&     (2)*thrustvector(1)
+    arg1 = rawtangent(1)**2 + rawtangent(2)**2 + rawtangent(3)**2
+    tangentnorm = sqrt(arg1)
+    distance2axis = tangentnorm/thrustvectornorm
+! compute tangential vector
+    tangent = rawtangent/tangentnorm
+  end subroutine computecellspatialmetrics
 ! ----------------------------------------------------------------------
 !                                                                      |
 !                    no tapenade routine below this line               |
diff --git a/src/adjoint/outputForward/residuals_d.f90 b/src/adjoint/outputForward/residuals_d.f90
index 8e12ed2e6..09e7c3771 100644
--- a/src/adjoint/outputForward/residuals_d.f90
+++ b/src/adjoint/outputForward/residuals_d.f90
@@ -340,12 +340,13 @@ end subroutine residual_block
 
 !  differentiation of sourceterms_block in forward (tangent) mode (with options i4 dr8 r8):
 !   variations   of useful results: *dw plocal
-!   with respect to varying inputs: uref pref *w *dw *vol actuatorregions.force
-!                actuatorregions.heat actuatorregions.volume plocal
+!   with respect to varying inputs: uref pref *w *dw *(actuatorregions.force)
+!                *(actuatorregions.heat) plocal
 !   rw status of diff variables: uref:in pref:in *w:in *dw:in-out
-!                *vol:in actuatorregions.force:in actuatorregions.heat:in
-!                actuatorregions.volume:in plocal:in-out
-!   plus diff mem management of: w:in dw:in vol:in
+!                *(actuatorregions.force):in *(actuatorregions.heat):in
+!                plocal:in-out
+!   plus diff mem management of: w:in dw:in actuatorregions.force:in
+!                actuatorregions.heat:in
   subroutine sourceterms_block_d(nn, res, iregion, plocal, plocald)
 ! apply the source terms for the given block. assume that the
 ! block pointers are already set.
@@ -363,13 +364,11 @@ subroutine sourceterms_block_d(nn, res, iregion, plocal, plocald)
     real(kind=realtype), intent(inout) :: plocald
 ! working
     integer(kind=inttype) :: i, j, k, ii, istart, iend
-    real(kind=realtype) :: ftmp(3), vx, vy, vz, f_fact(3), q_fact, qtmp&
-&   , redim, factor, ostart, oend
-    real(kind=realtype) :: ftmpd(3), vxd, vyd, vzd, f_factd(3), q_factd&
-&   , qtmpd, redimd
+    real(kind=realtype) :: vx, vy, vz, fx, fy, fz, q, redim, factor, &
+&   ostart, oend
+    real(kind=realtype) :: vxd, vyd, vzd, fxd, fyd, fzd, qd, redimd
     real(kind=realtype) :: temp
     real(kind=realtype) :: temp0
-    real(kind=realtype) :: temp1
     redimd = uref*prefd + pref*urefd
     redim = pref*uref
 ! compute the relaxation factor based on the ordersconverged
@@ -385,20 +384,6 @@ subroutine sourceterms_block_d(nn, res, iregion, plocal, plocald)
       oend = actuatorregions(iregion)%relaxend
       factor = (ordersconverged-ostart)/(oend-ostart)
     end if
-! compute the constant force factor
-    temp = actuatorregions(iregion)%volume*pref
-    f_factd = factor*(actuatorregionsd(iregion)%force-actuatorregions(&
-&     iregion)%force*(pref*actuatorregionsd(iregion)%volume+&
-&     actuatorregions(iregion)%volume*prefd)/temp)/temp
-    f_fact = factor*(actuatorregions(iregion)%force/temp)
-! heat factor. this is heat added per unit volume per unit time
-    temp = lref*lref*actuatorregions(iregion)%volume
-    temp0 = temp*pref*uref
-    temp1 = actuatorregions(iregion)%heat/temp0
-    q_factd = factor*(actuatorregionsd(iregion)%heat-temp1*(pref*uref*&
-&     lref**2*actuatorregionsd(iregion)%volume+temp*(uref*prefd+pref*&
-&     urefd)))/temp0
-    q_fact = factor*temp1
 ! loop over the ranges for this block
     istart = actuatorregions(iregion)%blkptr(nn-1) + 1
     iend = actuatorregions(iregion)%blkptr(nn)
@@ -408,8 +393,18 @@ subroutine sourceterms_block_d(nn, res, iregion, plocal, plocald)
       j = actuatorregions(iregion)%cellids(2, ii)
       k = actuatorregions(iregion)%cellids(3, ii)
 ! this actually gets the force
-      ftmpd = f_fact*vold(i, j, k) + vol(i, j, k)*f_factd
-      ftmp = vol(i, j, k)*f_fact
+      temp = actuatorregions(iregion)%force(1, ii)/pref
+      fxd = factor*(actuatorregionsd(iregion)%force(1, ii)-temp*prefd)/&
+&       pref
+      fx = factor*temp
+      temp = actuatorregions(iregion)%force(2, ii)/pref
+      fyd = factor*(actuatorregionsd(iregion)%force(2, ii)-temp*prefd)/&
+&       pref
+      fy = factor*temp
+      temp = actuatorregions(iregion)%force(3, ii)/pref
+      fzd = factor*(actuatorregionsd(iregion)%force(3, ii)-temp*prefd)/&
+&       pref
+      fz = factor*temp
       vxd = wd(i, j, k, ivx)
       vx = w(i, j, k, ivx)
       vyd = wd(i, j, k, ivy)
@@ -417,24 +412,30 @@ subroutine sourceterms_block_d(nn, res, iregion, plocal, plocald)
       vzd = wd(i, j, k, ivz)
       vz = w(i, j, k, ivz)
 ! this gets the heat addition rate
-      qtmpd = q_fact*vold(i, j, k) + vol(i, j, k)*q_factd
-      qtmp = vol(i, j, k)*q_fact
+      temp = lref*lref*pref*uref
+      temp0 = actuatorregions(iregion)%heat(ii)/temp
+      qd = factor*(actuatorregionsd(iregion)%heat(ii)-temp0*lref**2*(&
+&       uref*prefd+pref*urefd))/temp
+      q = factor*temp0
       if (res) then
 ! momentum residuals
-        dwd(i, j, k, imx:imz) = dwd(i, j, k, imx:imz) - ftmpd
-        dw(i, j, k, imx:imz) = dw(i, j, k, imx:imz) - ftmp
+        dwd(i, j, k, imx) = dwd(i, j, k, imx) - fxd
+        dw(i, j, k, imx) = dw(i, j, k, imx) - fx
+        dwd(i, j, k, imy) = dwd(i, j, k, imy) - fyd
+        dw(i, j, k, imy) = dw(i, j, k, imy) - fy
+        dwd(i, j, k, imz) = dwd(i, j, k, imz) - fzd
+        dw(i, j, k, imz) = dw(i, j, k, imz) - fz
 ! energy residuals
-        dwd(i, j, k, irhoe) = dwd(i, j, k, irhoe) - vx*ftmpd(1) - ftmp(1&
-&         )*vxd - vy*ftmpd(2) - ftmp(2)*vyd - vz*ftmpd(3) - ftmp(3)*vzd &
-&         - qtmpd
-        dw(i, j, k, irhoe) = dw(i, j, k, irhoe) - ftmp(1)*vx - ftmp(2)*&
-&         vy - ftmp(3)*vz - qtmp
+        dwd(i, j, k, irhoe) = dwd(i, j, k, irhoe) - vx*fxd - fx*vxd - vy&
+&         *fyd - fy*vyd - vz*fzd - fz*vzd - qd
+        dw(i, j, k, irhoe) = dw(i, j, k, irhoe) - fx*vx - fy*vy - fz*vz &
+&         - q
       else
 ! add in the local power contribution:
-        temp1 = vx*ftmp(1) + vy*ftmp(2) + vz*ftmp(3)
-        plocald = plocald + redim*(ftmp(1)*vxd+vx*ftmpd(1)+ftmp(2)*vyd+&
-&         vy*ftmpd(2)+ftmp(3)*vzd+vz*ftmpd(3)) + temp1*redimd
-        plocal = plocal + temp1*redim
+        temp0 = vx*fx + vy*fy + vz*fz
+        plocald = plocald + redim*(fx*vxd+vx*fxd+fy*vyd+vy*fyd+fz*vzd+vz&
+&         *fzd) + temp0*redimd
+        plocal = plocal + temp0*redim
       end if
     end do
   end subroutine sourceterms_block_d
@@ -455,8 +456,8 @@ subroutine sourceterms_block(nn, res, iregion, plocal)
     real(kind=realtype), intent(inout) :: plocal
 ! working
     integer(kind=inttype) :: i, j, k, ii, istart, iend
-    real(kind=realtype) :: ftmp(3), vx, vy, vz, f_fact(3), q_fact, qtmp&
-&   , redim, factor, ostart, oend
+    real(kind=realtype) :: vx, vy, vz, fx, fy, fz, q, redim, factor, &
+&   ostart, oend
     redim = pref*uref
 ! compute the relaxation factor based on the ordersconverged
 ! how far we are into the ramp:
@@ -471,12 +472,6 @@ subroutine sourceterms_block(nn, res, iregion, plocal)
       oend = actuatorregions(iregion)%relaxend
       factor = (ordersconverged-ostart)/(oend-ostart)
     end if
-! compute the constant force factor
-    f_fact = factor*actuatorregions(iregion)%force/actuatorregions(&
-&     iregion)%volume/pref
-! heat factor. this is heat added per unit volume per unit time
-    q_fact = factor*actuatorregions(iregion)%heat/actuatorregions(&
-&     iregion)%volume/(pref*uref*lref*lref)
 ! loop over the ranges for this block
     istart = actuatorregions(iregion)%blkptr(nn-1) + 1
     iend = actuatorregions(iregion)%blkptr(nn)
@@ -487,21 +482,25 @@ subroutine sourceterms_block(nn, res, iregion, plocal)
       j = actuatorregions(iregion)%cellids(2, ii)
       k = actuatorregions(iregion)%cellids(3, ii)
 ! this actually gets the force
-      ftmp = vol(i, j, k)*f_fact
+      fx = factor*actuatorregions(iregion)%force(1, ii)/pref
+      fy = factor*actuatorregions(iregion)%force(2, ii)/pref
+      fz = factor*actuatorregions(iregion)%force(3, ii)/pref
       vx = w(i, j, k, ivx)
       vy = w(i, j, k, ivy)
       vz = w(i, j, k, ivz)
 ! this gets the heat addition rate
-      qtmp = vol(i, j, k)*q_fact
+      q = factor*actuatorregions(iregion)%heat(ii)/(pref*uref*lref*lref)
       if (res) then
 ! momentum residuals
-        dw(i, j, k, imx:imz) = dw(i, j, k, imx:imz) - ftmp
+        dw(i, j, k, imx) = dw(i, j, k, imx) - fx
+        dw(i, j, k, imy) = dw(i, j, k, imy) - fy
+        dw(i, j, k, imz) = dw(i, j, k, imz) - fz
 ! energy residuals
-        dw(i, j, k, irhoe) = dw(i, j, k, irhoe) - ftmp(1)*vx - ftmp(2)*&
-&         vy - ftmp(3)*vz - qtmp
+        dw(i, j, k, irhoe) = dw(i, j, k, irhoe) - fx*vx - fy*vy - fz*vz &
+&         - q
       else
 ! add in the local power contribution:
-        plocal = plocal + (vx*ftmp(1)+vy*ftmp(2)+vz*ftmp(3))*redim
+        plocal = plocal + (vx*fx+vy*fy+vz*fz)*redim
       end if
     end do
   end subroutine sourceterms_block
diff --git a/src/adjoint/outputReverse/actuatorRegion_b.f90 b/src/adjoint/outputReverse/actuatorRegion_b.f90
index 069350fae..c74214599 100644
--- a/src/adjoint/outputReverse/actuatorRegion_b.f90
+++ b/src/adjoint/outputReverse/actuatorRegion_b.f90
@@ -8,47 +8,47 @@ module actuatorregion_b
   implicit none
 
 contains
-!  differentiation of computeactuatorregionvolume in reverse (adjoint) mode (with options noisize i4 dr8 r8):
-!   gradient     of useful results: *vol actuatorregions.vollocal
-!   with respect to varying inputs: *vol actuatorregions.vollocal
-!   rw status of diff variables: *vol:incr actuatorregions.vollocal:in-out
-!   plus diff mem management of: vol:in
-  subroutine computeactuatorregionvolume_b(nn, iregion)
-    use blockpointers, only : ndom, vol, vold
+  subroutine computecellspatialmetrics(i, j, k, centerpoint, &
+&   thrustvector, distance2plane, distance2axis, tangent)
+    use constants
+    use blockpointers, only : x
     implicit none
 ! inputs
-    integer(kind=inttype), intent(in) :: nn, iregion
+    real(kind=realtype), dimension(3), intent(in) :: centerpoint, &
+&   thrustvector
+    integer(kind=inttype), intent(in) :: i, j, k
+! outputs
+    real(kind=realtype), intent(out) :: distance2axis, distance2plane
+    real(kind=realtype), dimension(3), intent(out) :: tangent
 ! working
-    integer(kind=inttype) :: iii
-    integer(kind=inttype) :: i, j, k
-    do iii=actuatorregions(iregion)%blkptr(nn),actuatorregions(iregion)%&
-&       blkptr(nn-1)+1,-1
-      i = actuatorregions(iregion)%cellids(1, iii)
-      j = actuatorregions(iregion)%cellids(2, iii)
-      k = actuatorregions(iregion)%cellids(3, iii)
-      vold(i, j, k) = vold(i, j, k) + actuatorregionsd(iregion)%vollocal
-    end do
-  end subroutine computeactuatorregionvolume_b
-
-  subroutine computeactuatorregionvolume(nn, iregion)
-    use blockpointers, only : ndom, vol
-    implicit none
-! inputs
-    integer(kind=inttype), intent(in) :: nn, iregion
-! working
-    integer(kind=inttype) :: iii
-    integer(kind=inttype) :: i, j, k
-! loop over the region for this block
-    do iii=actuatorregions(iregion)%blkptr(nn-1)+1,actuatorregions(&
-&       iregion)%blkptr(nn)
-      i = actuatorregions(iregion)%cellids(1, iii)
-      j = actuatorregions(iregion)%cellids(2, iii)
-      k = actuatorregions(iregion)%cellids(3, iii)
-! sum the volume of each cell within the region on this proc
-      actuatorregions(iregion)%vollocal = actuatorregions(iregion)%&
-&       vollocal + vol(i, j, k)
-    end do
-  end subroutine computeactuatorregionvolume
+    real(kind=realtype) :: thrustvectornorm, dotproduct, tangentnorm
+    real(kind=realtype), dimension(3) :: xcen, distance2center, &
+&   rawtangent
+    intrinsic sqrt
+! compute the cell center
+    xcen = eighth*(x(i-1, j-1, k-1, :)+x(i, j-1, k-1, :)+x(i-1, j, k-1, &
+&     :)+x(i, j, k-1, :)+x(i-1, j-1, k, :)+x(i, j-1, k, :)+x(i-1, j, k, &
+&     :)+x(i, j, k, :))
+    thrustvectornorm = sqrt(thrustvector(1)**2 + thrustvector(2)**2 + &
+&     thrustvector(3)**2)
+    distance2center = xcen - centerpoint
+! compute distance to plane
+    dotproduct = thrustvector(1)*distance2center(1) + thrustvector(2)*&
+&     distance2center(2) + thrustvector(3)*distance2center(3)
+    distance2plane = dotproduct/thrustvectornorm
+! compute distance to axis
+    rawtangent(1) = distance2center(2)*thrustvector(3) - distance2center&
+&     (3)*thrustvector(2)
+    rawtangent(2) = distance2center(3)*thrustvector(1) - distance2center&
+&     (1)*thrustvector(3)
+    rawtangent(3) = distance2center(1)*thrustvector(2) - distance2center&
+&     (2)*thrustvector(1)
+    tangentnorm = sqrt(rawtangent(1)**2 + rawtangent(2)**2 + rawtangent(&
+&     3)**2)
+    distance2axis = tangentnorm/thrustvectornorm
+! compute tangential vector
+    tangent = rawtangent/tangentnorm
+  end subroutine computecellspatialmetrics
 ! ----------------------------------------------------------------------
 !                                                                      |
 !                    no tapenade routine below this line               |
diff --git a/src/adjoint/outputReverse/residuals_b.f90 b/src/adjoint/outputReverse/residuals_b.f90
index b16d723d8..883ed69a7 100644
--- a/src/adjoint/outputReverse/residuals_b.f90
+++ b/src/adjoint/outputReverse/residuals_b.f90
@@ -335,14 +335,15 @@ subroutine residual_block()
   end subroutine residual_block
 
 !  differentiation of sourceterms_block in reverse (adjoint) mode (with options noisize i4 dr8 r8):
-!   gradient     of useful results: uref pref *w *dw *vol actuatorregions.force
-!                actuatorregions.heat actuatorregions.volume plocal
-!   with respect to varying inputs: uref pref *w *dw *vol actuatorregions.force
-!                actuatorregions.heat actuatorregions.volume plocal
+!   gradient     of useful results: uref pref *w *dw *(actuatorregions.force)
+!                *(actuatorregions.heat) plocal
+!   with respect to varying inputs: uref pref *w *dw *(actuatorregions.force)
+!                *(actuatorregions.heat) plocal
 !   rw status of diff variables: uref:incr pref:incr *w:incr *dw:in-out
-!                *vol:incr actuatorregions.force:incr actuatorregions.heat:incr
-!                actuatorregions.volume:incr plocal:in-out
-!   plus diff mem management of: w:in dw:in vol:in
+!                *(actuatorregions.force):incr *(actuatorregions.heat):incr
+!                plocal:in-out
+!   plus diff mem management of: w:in dw:in actuatorregions.force:in
+!                actuatorregions.heat:in
   subroutine sourceterms_block_b(nn, res, iregion, plocal, plocald)
 ! apply the source terms for the given block. assume that the
 ! block pointers are already set.
@@ -360,16 +361,13 @@ subroutine sourceterms_block_b(nn, res, iregion, plocal, plocald)
     real(kind=realtype), intent(inout) :: plocald
 ! working
     integer(kind=inttype) :: i, j, k, ii, istart, iend
-    real(kind=realtype) :: ftmp(3), vx, vy, vz, f_fact(3), q_fact, qtmp&
-&   , redim, factor, ostart, oend
-    real(kind=realtype) :: ftmpd(3), vxd, vyd, vzd, f_factd(3), q_factd&
-&   , qtmpd, redimd
-    real(kind=realtype), dimension(3) :: tempd
+    real(kind=realtype) :: vx, vy, vz, fx, fy, fz, q, redim, factor, &
+&   ostart, oend
+    real(kind=realtype) :: vxd, vyd, vzd, fxd, fyd, fzd, qd, redimd
     real(kind=realtype) :: temp
+    real(kind=realtype) :: tempd
     real(kind=realtype) :: tempd0
-    real(kind=realtype) :: temp0
     real(kind=realtype) :: tempd1
-    real(kind=realtype) :: tempd2
     integer :: branch
     redim = pref*uref
 ! compute the relaxation factor based on the ordersconverged
@@ -388,18 +386,10 @@ subroutine sourceterms_block_b(nn, res, iregion, plocal, plocald)
       oend = actuatorregions(iregion)%relaxend
       factor = (ordersconverged-ostart)/(oend-ostart)
     end if
-! compute the constant force factor
-    f_fact = factor*actuatorregions(iregion)%force/actuatorregions(&
-&     iregion)%volume/pref
-! heat factor. this is heat added per unit volume per unit time
-    q_fact = factor*actuatorregions(iregion)%heat/actuatorregions(&
-&     iregion)%volume/(pref*uref*lref*lref)
 ! loop over the ranges for this block
     istart = actuatorregions(iregion)%blkptr(nn-1) + 1
     iend = actuatorregions(iregion)%blkptr(nn)
-    q_factd = 0.0_8
     redimd = 0.0_8
-    f_factd = 0.0_8
 !$bwd-of ii-loop 
     do ii=istart,iend
 ! extract the cell id.
@@ -407,56 +397,55 @@ subroutine sourceterms_block_b(nn, res, iregion, plocal, plocald)
       j = actuatorregions(iregion)%cellids(2, ii)
       k = actuatorregions(iregion)%cellids(3, ii)
 ! this actually gets the force
-      ftmp = vol(i, j, k)*f_fact
+      fx = factor*actuatorregions(iregion)%force(1, ii)/pref
+      fy = factor*actuatorregions(iregion)%force(2, ii)/pref
+      fz = factor*actuatorregions(iregion)%force(3, ii)/pref
       vx = w(i, j, k, ivx)
       vy = w(i, j, k, ivy)
       vz = w(i, j, k, ivz)
 ! this gets the heat addition rate
       if (res) then
-        ftmpd = 0.0_8
-        ftmpd(1) = ftmpd(1) - vx*dwd(i, j, k, irhoe)
-        vxd = -(ftmp(1)*dwd(i, j, k, irhoe))
-        ftmpd(2) = ftmpd(2) - vy*dwd(i, j, k, irhoe)
-        vyd = -(ftmp(2)*dwd(i, j, k, irhoe))
-        ftmpd(3) = ftmpd(3) - vz*dwd(i, j, k, irhoe)
-        vzd = -(ftmp(3)*dwd(i, j, k, irhoe))
-        qtmpd = -dwd(i, j, k, irhoe)
-        ftmpd = ftmpd - dwd(i, j, k, imx:imz)
+        fxd = -(vx*dwd(i, j, k, irhoe)) - dwd(i, j, k, imx)
+        vxd = -(fx*dwd(i, j, k, irhoe))
+        fyd = -(vy*dwd(i, j, k, irhoe)) - dwd(i, j, k, imy)
+        vyd = -(fy*dwd(i, j, k, irhoe))
+        fzd = -(vz*dwd(i, j, k, irhoe)) - dwd(i, j, k, imz)
+        vzd = -(fz*dwd(i, j, k, irhoe))
+        qd = -dwd(i, j, k, irhoe)
       else
-        ftmpd = 0.0_8
         tempd1 = redim*plocald
-        redimd = redimd + (vx*ftmp(1)+vy*ftmp(2)+vz*ftmp(3))*plocald
-        vxd = ftmp(1)*tempd1
-        ftmpd(1) = ftmpd(1) + vx*tempd1
-        vyd = ftmp(2)*tempd1
-        ftmpd(2) = ftmpd(2) + vy*tempd1
-        vzd = ftmp(3)*tempd1
-        ftmpd(3) = ftmpd(3) + vz*tempd1
-        qtmpd = 0.0_8
+        redimd = redimd + (vx*fx+vy*fy+vz*fz)*plocald
+        vxd = fx*tempd1
+        fxd = vx*tempd1
+        vyd = fy*tempd1
+        fyd = vy*tempd1
+        vzd = fz*tempd1
+        fzd = vz*tempd1
+        qd = 0.0_8
       end if
-      vold(i, j, k) = vold(i, j, k) + q_fact*qtmpd + sum(f_fact*ftmpd)
-      q_factd = q_factd + vol(i, j, k)*qtmpd
+      tempd = factor*fzd/pref
+      temp = lref*lref*pref*uref
+      tempd0 = factor*qd/temp
+      actuatorregionsd(iregion)%heat(ii) = actuatorregionsd(iregion)%&
+&       heat(ii) + tempd0
+      tempd1 = -(lref**2*actuatorregions(iregion)%heat(ii)*tempd0/temp)
+      prefd = prefd + uref*tempd1 - actuatorregions(iregion)%force(3, ii&
+&       )*tempd/pref
+      urefd = urefd + pref*tempd1
       wd(i, j, k, ivz) = wd(i, j, k, ivz) + vzd
       wd(i, j, k, ivy) = wd(i, j, k, ivy) + vyd
       wd(i, j, k, ivx) = wd(i, j, k, ivx) + vxd
-      f_factd = f_factd + vol(i, j, k)*ftmpd
+      actuatorregionsd(iregion)%force(3, ii) = actuatorregionsd(iregion)&
+&       %force(3, ii) + tempd
+      tempd = factor*fyd/pref
+      actuatorregionsd(iregion)%force(2, ii) = actuatorregionsd(iregion)&
+&       %force(2, ii) + tempd
+      prefd = prefd - actuatorregions(iregion)%force(2, ii)*tempd/pref
+      tempd = factor*fxd/pref
+      actuatorregionsd(iregion)%force(1, ii) = actuatorregionsd(iregion)&
+&       %force(1, ii) + tempd
+      prefd = prefd - actuatorregions(iregion)%force(1, ii)*tempd/pref
     end do
-    tempd = factor*f_factd/(actuatorregions(iregion)%volume*pref)
-    tempd0 = -(sum(actuatorregions(iregion)%force*tempd)/(&
-&     actuatorregions(iregion)%volume*pref))
-    temp = lref*lref*actuatorregions(iregion)%volume
-    temp0 = temp*pref*uref
-    tempd1 = factor*q_factd/temp0
-    actuatorregionsd(iregion)%heat = actuatorregionsd(iregion)%heat + &
-&     tempd1
-    tempd2 = -(actuatorregions(iregion)%heat*tempd1/temp0)
-    actuatorregionsd(iregion)%volume = actuatorregionsd(iregion)%volume &
-&     + lref**2*pref*uref*tempd2 + pref*tempd0
-    prefd = prefd + uref*temp*tempd2 + actuatorregions(iregion)%volume*&
-&     tempd0
-    urefd = urefd + pref*temp*tempd2
-    actuatorregionsd(iregion)%force = actuatorregionsd(iregion)%force + &
-&     tempd
     call popcontrol1b(branch)
     prefd = prefd + uref*redimd
     urefd = urefd + pref*redimd
@@ -478,8 +467,8 @@ subroutine sourceterms_block(nn, res, iregion, plocal)
     real(kind=realtype), intent(inout) :: plocal
 ! working
     integer(kind=inttype) :: i, j, k, ii, istart, iend
-    real(kind=realtype) :: ftmp(3), vx, vy, vz, f_fact(3), q_fact, qtmp&
-&   , redim, factor, ostart, oend
+    real(kind=realtype) :: vx, vy, vz, fx, fy, fz, q, redim, factor, &
+&   ostart, oend
     redim = pref*uref
 ! compute the relaxation factor based on the ordersconverged
 ! how far we are into the ramp:
@@ -494,12 +483,6 @@ subroutine sourceterms_block(nn, res, iregion, plocal)
       oend = actuatorregions(iregion)%relaxend
       factor = (ordersconverged-ostart)/(oend-ostart)
     end if
-! compute the constant force factor
-    f_fact = factor*actuatorregions(iregion)%force/actuatorregions(&
-&     iregion)%volume/pref
-! heat factor. this is heat added per unit volume per unit time
-    q_fact = factor*actuatorregions(iregion)%heat/actuatorregions(&
-&     iregion)%volume/(pref*uref*lref*lref)
 ! loop over the ranges for this block
     istart = actuatorregions(iregion)%blkptr(nn-1) + 1
     iend = actuatorregions(iregion)%blkptr(nn)
@@ -510,21 +493,25 @@ subroutine sourceterms_block(nn, res, iregion, plocal)
       j = actuatorregions(iregion)%cellids(2, ii)
       k = actuatorregions(iregion)%cellids(3, ii)
 ! this actually gets the force
-      ftmp = vol(i, j, k)*f_fact
+      fx = factor*actuatorregions(iregion)%force(1, ii)/pref
+      fy = factor*actuatorregions(iregion)%force(2, ii)/pref
+      fz = factor*actuatorregions(iregion)%force(3, ii)/pref
       vx = w(i, j, k, ivx)
       vy = w(i, j, k, ivy)
       vz = w(i, j, k, ivz)
 ! this gets the heat addition rate
-      qtmp = vol(i, j, k)*q_fact
+      q = factor*actuatorregions(iregion)%heat(ii)/(pref*uref*lref*lref)
       if (res) then
 ! momentum residuals
-        dw(i, j, k, imx:imz) = dw(i, j, k, imx:imz) - ftmp
+        dw(i, j, k, imx) = dw(i, j, k, imx) - fx
+        dw(i, j, k, imy) = dw(i, j, k, imy) - fy
+        dw(i, j, k, imz) = dw(i, j, k, imz) - fz
 ! energy residuals
-        dw(i, j, k, irhoe) = dw(i, j, k, irhoe) - ftmp(1)*vx - ftmp(2)*&
-&         vy - ftmp(3)*vz - qtmp
+        dw(i, j, k, irhoe) = dw(i, j, k, irhoe) - fx*vx - fy*vy - fz*vz &
+&         - q
       else
 ! add in the local power contribution:
-        plocal = plocal + (vx*ftmp(1)+vy*ftmp(2)+vz*ftmp(3))*redim
+        plocal = plocal + (vx*fx+vy*fy+vz*fz)*redim
       end if
     end do
   end subroutine sourceterms_block
diff --git a/src/adjoint/outputReverseFast/residuals_fast_b.f90 b/src/adjoint/outputReverseFast/residuals_fast_b.f90
index edf78f8d1..9b8be42e8 100644
--- a/src/adjoint/outputReverseFast/residuals_fast_b.f90
+++ b/src/adjoint/outputReverseFast/residuals_fast_b.f90
@@ -355,8 +355,8 @@ subroutine sourceterms_block_fast_b(nn, res, iregion, plocal)
     real(kind=realtype), intent(inout) :: plocal
 ! working
     integer(kind=inttype) :: i, j, k, ii, istart, iend
-    real(kind=realtype) :: ftmp(3), vx, vy, vz, f_fact(3), q_fact, qtmp&
-&   , redim, factor, ostart, oend
+    real(kind=realtype) :: vx, vy, vz, fx, fy, fz, q, redim, factor, &
+&   ostart, oend
     real(kind=realtype) :: vxd, vyd, vzd
     integer :: branch
 ! compute the relaxation factor based on the ordersconverged
@@ -378,10 +378,6 @@ subroutine sourceterms_block_fast_b(nn, res, iregion, plocal)
       oend = actuatorregions(iregion)%relaxend
       factor = (ordersconverged-ostart)/(oend-ostart)
     end if
-! compute the constant force factor
-    f_fact = factor*actuatorregions(iregion)%force/actuatorregions(&
-&     iregion)%volume/pref
-! heat factor. this is heat added per unit volume per unit time
 ! loop over the ranges for this block
     istart = actuatorregions(iregion)%blkptr(nn-1) + 1
     iend = actuatorregions(iregion)%blkptr(nn)
@@ -392,12 +388,14 @@ subroutine sourceterms_block_fast_b(nn, res, iregion, plocal)
       j = actuatorregions(iregion)%cellids(2, ii)
       k = actuatorregions(iregion)%cellids(3, ii)
 ! this actually gets the force
-      ftmp = vol(i, j, k)*f_fact
+      fx = factor*actuatorregions(iregion)%force(1, ii)/pref
+      fy = factor*actuatorregions(iregion)%force(2, ii)/pref
+      fz = factor*actuatorregions(iregion)%force(3, ii)/pref
 ! this gets the heat addition rate
       if (res) then
-        vxd = -(ftmp(1)*dwd(i, j, k, irhoe))
-        vyd = -(ftmp(2)*dwd(i, j, k, irhoe))
-        vzd = -(ftmp(3)*dwd(i, j, k, irhoe))
+        vxd = -(fx*dwd(i, j, k, irhoe))
+        vyd = -(fy*dwd(i, j, k, irhoe))
+        vzd = -(fz*dwd(i, j, k, irhoe))
       else
         vxd = 0.0_8
         vyd = 0.0_8
@@ -427,8 +425,8 @@ subroutine sourceterms_block(nn, res, iregion, plocal)
     real(kind=realtype), intent(inout) :: plocal
 ! working
     integer(kind=inttype) :: i, j, k, ii, istart, iend
-    real(kind=realtype) :: ftmp(3), vx, vy, vz, f_fact(3), q_fact, qtmp&
-&   , redim, factor, ostart, oend
+    real(kind=realtype) :: vx, vy, vz, fx, fy, fz, q, redim, factor, &
+&   ostart, oend
     redim = pref*uref
 ! compute the relaxation factor based on the ordersconverged
 ! how far we are into the ramp:
@@ -443,12 +441,6 @@ subroutine sourceterms_block(nn, res, iregion, plocal)
       oend = actuatorregions(iregion)%relaxend
       factor = (ordersconverged-ostart)/(oend-ostart)
     end if
-! compute the constant force factor
-    f_fact = factor*actuatorregions(iregion)%force/actuatorregions(&
-&     iregion)%volume/pref
-! heat factor. this is heat added per unit volume per unit time
-    q_fact = factor*actuatorregions(iregion)%heat/actuatorregions(&
-&     iregion)%volume/(pref*uref*lref*lref)
 ! loop over the ranges for this block
     istart = actuatorregions(iregion)%blkptr(nn-1) + 1
     iend = actuatorregions(iregion)%blkptr(nn)
@@ -459,21 +451,25 @@ subroutine sourceterms_block(nn, res, iregion, plocal)
       j = actuatorregions(iregion)%cellids(2, ii)
       k = actuatorregions(iregion)%cellids(3, ii)
 ! this actually gets the force
-      ftmp = vol(i, j, k)*f_fact
+      fx = factor*actuatorregions(iregion)%force(1, ii)/pref
+      fy = factor*actuatorregions(iregion)%force(2, ii)/pref
+      fz = factor*actuatorregions(iregion)%force(3, ii)/pref
       vx = w(i, j, k, ivx)
       vy = w(i, j, k, ivy)
       vz = w(i, j, k, ivz)
 ! this gets the heat addition rate
-      qtmp = vol(i, j, k)*q_fact
+      q = factor*actuatorregions(iregion)%heat(ii)/(pref*uref*lref*lref)
       if (res) then
 ! momentum residuals
-        dw(i, j, k, imx:imz) = dw(i, j, k, imx:imz) - ftmp
+        dw(i, j, k, imx) = dw(i, j, k, imx) - fx
+        dw(i, j, k, imy) = dw(i, j, k, imy) - fy
+        dw(i, j, k, imz) = dw(i, j, k, imz) - fz
 ! energy residuals
-        dw(i, j, k, irhoe) = dw(i, j, k, irhoe) - ftmp(1)*vx - ftmp(2)*&
-&         vy - ftmp(3)*vz - qtmp
+        dw(i, j, k, irhoe) = dw(i, j, k, irhoe) - fx*vx - fy*vy - fz*vz &
+&         - q
       else
 ! add in the local power contribution:
-        plocal = plocal + (vx*ftmp(1)+vy*ftmp(2)+vz*ftmp(3))*redim
+        plocal = plocal + (vx*fx+vy*fy+vz*fz)*redim
       end if
     end do
   end subroutine sourceterms_block
diff --git a/src/bcdata/BCData.F90 b/src/bcdata/BCData.F90
index 4ebb6dd75..14837f94e 100644
--- a/src/bcdata/BCData.F90
+++ b/src/bcdata/BCData.F90
@@ -1413,7 +1413,6 @@ subroutine setBCData(bcDataNamesIn, bcDataIn, famLists, sps, &
         use sorting, only: famInList
         use utils, only: setPointers, terminate, char2str
         use communication, only: myid
-        use actuatorRegionData, only: actuatorRegions, nActuatorRegions
         !
         !      Subroutine arguments.
         !
@@ -1475,27 +1474,6 @@ subroutine setBCData(bcDataNamesIn, bcDataIn, famLists, sps, &
             end do varLoop
         end do domainsLoop
 
-        ! Loop over any actuator regions since they also could have to set BCData
-        regionLoop: do iRegion = 1, nActuatorRegions
-            varLoop2: do iVar = 1, nVar
-                nFam = famLists(iVar, 1)
-                famInclude2: if (famInList(actuatorRegions(iRegion)%famID, famLists(iVar, 2:2 + nFam - 1))) then
-
-                    ! Extract the name
-                    varName = char2str(bcDataNamesIn(iVar, :), maxCGNSNameLen)
-
-                    if (trim(varName) == "Thrust") then
-                        actuatorRegions(iRegion)%force = actuatorRegions(iRegion)%axisVec * &
-                                                         bcDataIn(iVar)
-                    else if (trim(varName) == "Torque") then
-                        actuatorRegions(iRegion)%torque = bcDataIn(iVar)
-                    else if (trim(varName) == "Heat") then
-                        actuatorRegions(iRegion)%heat = bcDataIn(iVar)
-                    end if
-                end if famInclude2
-            end do varLoop2
-        end do regionLoop
-
     end subroutine setBCData
 
     subroutine setBCData_d(bcDataNamesIn, bcDataIn, bcDataInd, famLists, sps, &
@@ -1573,28 +1551,16 @@ subroutine setBCData_d(bcDataNamesIn, bcDataIn, bcDataInd, famLists, sps, &
             end do varLoop
         end do domainsLoop
 
-        ! Loop over any actuator regions since they also could have to set BCData
+        ! Loop over any actuator regions and set the seeds to 0 (this is a temporary fix until the BC seeds can be set in python)
         regionLoop: do iRegion = 1, nActuatorRegions
             varLoop2: do iVar = 1, nVar
                 nFam = famLists(iVar, 1)
-                famInclude2: if (famInList(actuatorRegions(iRegion)%famID, famLists(iVar, 2:2 + nFam - 1))) then
-
-                    ! Extract the name
-                    varName = char2str(bcDataNamesIn(iVar, :), maxCGNSNameLen)
-
-                    if (trim(varName) == "Thrust") then
-                        actuatorRegions(iRegion)%force = actuatorRegions(iRegion)%axisVec * &
-                                                         bcDataIn(iVar)
-                        actuatorRegionsd(iRegion)%force = actuatorRegions(iRegion)%axisVec * &
-                                                          bcDataInd(iVar)
-                    else if (trim(varName) == "Torque") then
-                        actuatorRegions(iRegion)%torque = bcDataIn(iVar)
-                        actuatorRegionsd(iRegion)%torque = bcDataInd(iVar)
-                    else if (trim(varName) == "Heat") then
-                        actuatorRegions(iRegion)%heat = bcDataIn(iVar)
-                        actuatorRegionsd(iRegion)%heat = bcDataInd(iVar)
-                    end if
-                end if famInclude2
+                if (.not. famInList(actuatorRegions(iRegion)%famID, famLists(iVar, 2:2 + nFam - 1))) then
+                    cycle
+                end if
+
+                actuatorRegionsd(iRegion)%force =  zero
+                actuatorRegionsd(iRegion)%heat =  zero
             end do varLoop2
         end do regionLoop
 
@@ -1678,24 +1644,20 @@ subroutine setBCData_b(bcDataNamesIn, bcDataIn, bcDataInd, famLists, sps, &
             end do varLoop
         end do domainsLoop
 
-        ! Loop over any actuator regions since they also could have to set BCData
+
+        ! Loop over any actuator regions and set the seeds to 0 (this is a temporary fix until the BC seeds can be set in python)
         regionLoop: do iRegion = 1, nActuatorRegions
             varLoop2: do iVar = 1, nVar
                 nFam = famLists(iVar, 1)
-                famInclude2: if (famInList(actuatorRegions(iRegion)%famID, famLists(iVar, 2:2 + nFam - 1))) then
-
-                    ! Extract the name
-                    varName = char2str(bcDataNamesIn(iVar, :), maxCGNSNameLen)
-
-                    if (trim(varName) == "Thrust") then
-                        bcDataInd(ivar) = &
-                            sum(actuatorRegions(iRegion)%axisVec * actuatorRegionsd(iRegion)%force)
-                    else if (trim(varName) == "Torque") then
-                        bcDataInd(ivar) = actuatorRegionsd(iRegion)%torque
-                    else if (trim(varName) == "Heat") then
-                        bcDataInd(ivar) = actuatorRegionsd(iRegion)%heat
-                    end if
-                end if famInclude2
+                if (.not. famInList(actuatorRegions(iRegion)%famID, famLists(iVar, 2:2 + nFam - 1))) then
+                    cycle
+                end if
+
+                ! Extract the name
+                varName = char2str(bcDataNamesIn(iVar, :), maxCGNSNameLen)
+
+                bcDataInd(ivar) = zero
+
             end do varLoop2
         end do regionLoop
 
diff --git a/src/f2py/adflow.pyf b/src/f2py/adflow.pyf
index 4eec9145d..3e4d147a9 100644
--- a/src/f2py/adflow.pyf
+++ b/src/f2py/adflow.pyf
@@ -685,22 +685,44 @@ python module libadflow
        end module usersurfaceintegrations
 
        module actuatorregion ! in :test:actuatorDiskRegion.F90
-         subroutine addactuatorregion(pts,conn,axis1,axis2,famname,famid,thrust,torque,heat,relaxstart,relaxend,npts,nconn) ! in :test:actuatorDiskRegion.F90:actuatorregion
-           real(kind=realtype) dimension(3,npts),intent(in) :: pts
-           integer(kind=inttype) dimension(4,nconn),intent(in) :: conn
-           real(kind=realtype), dimension(3),intent(in) :: axis1
-           real(kind=realtype), dimension(3),intent(in) :: axis2
-           character*(*) :: famname
+         subroutine computeinitialspatialmetrics(centerpoint, thrustvector, n, distance2plane, distance2axis, tangent) ! in :test:actuatorDiskRegion.F90:actuatorregion
+           real(kind=realtype) dimension(3),intent(in) :: centerpoint
+           real(kind=realtype) dimension(3),intent(in) :: thrustvector
+           integer(kind=inttype) intent(in) :: n
+           real(kind=realtype) dimension(n),intent(out),depend(n) :: distance2plane
+           real(kind=realtype) dimension(n),intent(out),depend(n) :: distance2axis
+           real(kind=realtype) dimension(3, n),intent(out),depend(n) :: tangent
+         end subroutine computeinitialspatialmetrics
+
+         subroutine addactuatorregion(flag,n,famname,famid,relaxstart,relaxend,iregion,nlocalcells) ! in :test:actuatorDiskRegion.F90:actuatorregion
+           integer(kind=inttype) dimension(n) :: flag
+           integer(kind=inttype), optional,intent(in),check(len(flag)>=n),depend(flag) :: n=len(flag)
+           character*(*) intent(in):: famname
            integer(kind=inttype) intent(in) :: famid
-           real(kind=realtype) :: thrust
-           real(kind=realtype) :: torque
-           real(kind=realtype) :: heat
-           real(kind=realtype) :: relaxstart
-           real(kind=realtype) :: relaxend
-           integer(kind=inttype), optional,intent(in),check(shape(pts,1)==npts),depend(pts) :: npts=shape(pts,1)
-           integer(kind=inttype), optional,intent(in),check(shape(conn,1)==nconn),depend(conn) :: nconn=shape(conn,1)
+           real(kind=realtype) intent(in):: relaxstart
+           real(kind=realtype) intent(in):: relaxend
+           integer(kind=inttype), intent(out) :: iregion
+           integer(kind=inttype), intent(out) :: nlocalcells
          end subroutine addactuatorregion
 
+         subroutine computespatialmetrics(iregion, nlocalcells, centerpoint, thrustvector, distance2plane, distance2axis, tangent, volume) ! in :test:actuatorDiskRegion.F90:actuatorregion
+           real(kind=realtype) dimension(3),intent(in) :: centerpoint
+           real(kind=realtype) dimension(3),intent(in) :: thrustvector
+           integer(kind=inttype), intent(in) :: iregion
+           integer(kind=inttype) intent(in) :: nlocalcells
+           real(kind=realtype) dimension(nlocalcells),intent(out),depend(nlocalcells) :: distance2plane
+           real(kind=realtype) dimension(nlocalcells),intent(out),depend(nlocalcells) :: distance2axis
+           real(kind=realtype) dimension(3, nlocalcells),intent(out),depend(nlocalcells) :: tangent
+           real(kind=realtype) dimension(nlocalcells),intent(out),depend(nlocalcells) :: volume
+         end subroutine addactuatorregion
+
+         subroutine populatebcvalues(iregion, nlocalcells, force, heat) ! in :test:actuatorDiskRegion.F90:actuatorregion
+           integer(kind=inttype), intent(in) :: iregion
+           integer(kind=inttype) intent(in) :: nlocalcells
+           real(kind=realtype) dimension(3, nlocalcells),intent(in),depend(nlocalcells) :: force
+           real(kind=realtype) dimension(nlocalcells),intent(in),depend(nlocalcells) :: heat
+         end subroutine populatebcvalues
+
          subroutine writeactuatorregions(filename)
            character*(*) intent(in) :: filename
          end subroutine writeactuatorregions
diff --git a/src/modules/actuatorRegionData.F90 b/src/modules/actuatorRegionData.F90
index b49d82c89..ffb811e79 100644
--- a/src/modules/actuatorRegionData.F90
+++ b/src/modules/actuatorRegionData.F90
@@ -11,22 +11,8 @@ module actuatorRegionData
         ! The total number of cells included this proc has
         integer(kind=intType) :: nCellIDs
 
-        ! the force vector to be applied on this region
-        ! this is equal to torque * axisVec
-        real(kind=realType) :: force(3)
-
-        ! magnitude of the total torque to be applied on this region
-        real(kind=realType) :: torque
-        ! vector that determines the direction of the applied torque
-        real(kind=realType), dimension(3) :: axisVec
-
-        ! total heat flux to be added on this regoin
-        real(kind=realType) :: heat
-
-        ! Volume is the total integrated volume of all cells (on all
-        ! procs) included in this region
-        real(kind=realType) :: volume
-        real(kind=realType) :: volLocal
+        real(kind=realType), dimension(:, :), pointer :: force
+        real(kind=realType), dimension(:), pointer :: heat
 
         integer(kind=intType), dimension(:), allocatable :: blkPtr
 
diff --git a/src/modules/diffSizes.f90 b/src/modules/diffSizes.f90
index b9e12494e..4641ab4dd 100644
--- a/src/modules/diffSizes.f90
+++ b/src/modules/diffSizes.f90
@@ -3,7 +3,7 @@ module diffSizes
     implicit none
     save
 
-    ! These are the diff sizes reqruied for the forward mode AD
+    ! These are the diff sizes required for the forward mode AD
     integer(kind=intType), parameter :: ISIZE3ofviscsubface = 3
     integer(kind=intType) :: ISIZE1OFDrfbcdata
 
@@ -13,7 +13,7 @@ module diffSizes
     integer(kind=intType) :: ISIZE3OFDrfflowdoms
     integer(kind=intType) :: ISIZE1OFDrfflowdoms_bcdata
 
-    ! These are the diff sizes reqruied for the reverse mode AD
+    ! These are the diff sizes required for the reverse mode AD
     integer(kind=intType) :: ISIZE3OFDrfrlv, ISIZE2OFDrfrlv, ISIZE1OFDrfrlv
     integer(kind=intType) :: ISIZE4OFDrfw, ISIZE3OFDrfw, ISIZE2OFDrfw, ISIZE1OFDrfw
     integer(kind=intType) :: ISIZE4OFDrffw, ISIZE3OFDrffw, ISIZE2OFDrffw, ISIZE1OFDrffw
@@ -68,7 +68,7 @@ module diffSizes
     integer(kind=intType) :: ISIZE2OFDrfbmtk1, ISIZE2OFDrfbvtk1
     integer(kind=intType) :: ISIZE1OFDrfbmtk1, ISIZE1OFDrfbvtk1
 
-    ! These are the diff sizes reqruied for the forward mode debug
+    ! These are the diff sizes required for the forward mode debug
     integer(kind=intType) :: ISIZE1OFDrfDrfbcdata_m, ISIZE2OFDrfDrfbcdata_m, ISIZE3OFDrfDrfbcdata_m
     integer(kind=intType) :: ISIZE1OFDrfDrfbcdata_fp, ISIZE2OFDrfDrfbcdata_fp, ISIZE3OFDrfDrfbcdata_fp
     integer(kind=intType) :: ISIZE1OFDrfDrfbcdata_fv, ISIZE2OFDrfDrfbcdata_fv, ISIZE3OFDrfDrfbcdata_fv
@@ -107,7 +107,9 @@ module diffSizes
     integer(kind=intType) :: ISIZE1OFDu1, ISIZE1OFDu2, ISIZE1OFDu3
     integer(kind=intType) :: ISIZE1OFLeft, ISIZE1OFRight, ISIZE1OFFlux
 
-    ! These are the diff sizes reqruied for the reverse mode
+    integer(kind=intType) :: ISIZE1OFtemp
+
+    ! These are the diff sizes required for the reverse mode
     integer(kind=intType) :: ISIZE1OFDRFCOEFTIME
     integer(kind=intType) :: ISIZE1OFDRFDTL, ISIZE2OFDRFDTL, ISIZE3OFDRFDTL
     integer(kind=intType) :: ISIZE1OFDRFDRFVISCSUBFACE_UTAU, ISIZE2OFDRFDRFVISCSUBFACE_UTAU
diff --git a/src/output/outputMod.F90 b/src/output/outputMod.F90
index 6214fae12..40c752127 100644
--- a/src/output/outputMod.F90
+++ b/src/output/outputMod.F90
@@ -1498,7 +1498,7 @@ subroutine storeSurfsolInBuffer(sps, buffer, nn, blockID, &
 
             select case (solName)
 
-             case (cgnsSkinFmag, cgnsStanton, cgnsYplus, &
+            case (cgnsSkinFmag, cgnsStanton, cgnsYplus, &
                   cgnsSkinFx, cgnsSkinFy, cgnsSkinFz, cgnsForceInDragDir, cgnsForceInLiftDir, &
                   cgnsSepSensor, cgnsSepSensorKs, cgnsSepSensorKsArea)
 
diff --git a/src/preprocessing/preprocessingAPI.F90 b/src/preprocessing/preprocessingAPI.F90
index 6516357ec..9991634ef 100644
--- a/src/preprocessing/preprocessingAPI.F90
+++ b/src/preprocessing/preprocessingAPI.F90
@@ -2732,7 +2732,6 @@ subroutine metric(level)
         use communication
         use inputTimeSpectral
         use checkVolBlock
-        use actuatorRegion, only: computeActuatorRegionVolume
         use actuatorRegionData, only: nActuatorRegions, actuatorRegions
         use inputIteration, only: printWarnings, printNegativeVolumes, &
                                   useSkewnessCheck, meshMaxSkewness, printBadlySkewedCells
@@ -2787,11 +2786,6 @@ subroutine metric(level)
         nBlockBad = 0
         nBadSkew = 0
 
-        ! Zero out the local volume pointers for the actuator zone
-        do iRegion = 1, nActuatorRegions
-            actuatorRegions(iRegion)%volLocal = zero
-        end do
-
         ! Loop over the number of spectral solutions and local blocks.
 
         spectral: do sps = 1, nTimeIntervalsSpectral
@@ -3063,11 +3057,6 @@ subroutine metric(level)
                     end do
                 end do
 
-                ! Compute the volume of the actuator regions on each domain
-                do iRegion = 1, nActuatorRegions
-                    call computeActuatorRegionVolume(nn, iRegion)
-                end do
-
                 ! Determine the orientation of the block. For the fine level
                 ! this is based on the number of positive and negative
                 ! volumes; on the coarse levels the corresponding fine level
@@ -3347,13 +3336,6 @@ subroutine metric(level)
             end do domains
         end do spectral
 
-        ! Loop over the actuator regions again to compute the total volumes
-        do iRegion = 1, nActuatorRegions
-            call mpi_allreduce(actuatorRegions(iRegion)%volLocal, actuatorRegions(iRegion)%volume, 1, &
-                               adflow_real, mpi_sum, adflow_comm_world, ierr)
-            call ECHK(ierr, __FILE__, __LINE__)
-        end do
-
         ! Determine the global number of bad blocks. The result must be
         ! known on all processors and thus an allreduce is needed.
 
diff --git a/src/solver/actuatorRegion.F90 b/src/solver/actuatorRegion.F90
index e372cceae..b8d394332 100644
--- a/src/solver/actuatorRegion.F90
+++ b/src/solver/actuatorRegion.F90
@@ -6,29 +6,51 @@ module actuatorRegion
     implicit none
 
 contains
-
-    subroutine computeActuatorRegionVolume(nn, iRegion)
-        use blockPointers, only: nDom, vol
+    subroutine computeCellSpatialMetrics(i, j, k, centerPoint, thrustVector, distance2plane, distance2axis, tangent)
+        use constants
+        use blockPointers, only: x
         implicit none
 
         ! Inputs
-        integer(kind=intType), intent(in) :: nn, iRegion
+        real(kind=realType), dimension(3), intent(in) :: centerPoint, thrustVector
+        integer(kind=intType), intent(in) :: i, j, k
+
+        ! Outputs
+        real(kind=realType), intent(out) :: distance2axis, distance2plane
+        real(kind=realType), dimension(3), intent(out) :: tangent
 
         ! Working
-        integer(kind=intType) :: iii
-        integer(kind=intType) :: i, j, k
+        real(kind=realType)  :: thrustVectorNorm, dotProduct, tangentNorm
+        real(kind=realType), dimension(3) ::  xCen, distance2center, rawTangent
 
-        ! Loop over the region for this block
-        do iii = actuatorRegions(iRegion)%blkPtr(nn - 1) + 1, actuatorRegions(iRegion)%blkPtr(nn)
-            i = actuatorRegions(iRegion)%cellIDs(1, iii)
-            j = actuatorRegions(iRegion)%cellIDs(2, iii)
-            k = actuatorRegions(iRegion)%cellIDs(3, iii)
+        ! Compute the cell center
+        xCen = eighth * (x(i - 1, j - 1, k - 1, :) + x(i, j - 1, k - 1, :) &
+                         + x(i - 1, j, k - 1, :) + x(i, j, k - 1, :) &
+                         + x(i - 1, j - 1, k, :) + x(i, j - 1, k, :) &
+                         + x(i - 1, j, k, :) + x(i, j, k, :))
 
-            ! Sum the volume of each cell within the region on this proc
-            actuatorRegions(iRegion)%volLocal = actuatorRegions(iRegion)%volLocal + vol(i, j, k)
-        end do
+        thrustVectorNorm = sqrt(thrustVector(1)**2 + thrustVector(2)**2 + thrustVector(3)**2)
+        distance2center = xCen - centerPoint
+
+        ! compute distance to plane
+        dotProduct = thrustVector(1)*distance2center(1) + thrustVector(2)*distance2center(2) + thrustVector(3)*distance2center(3)
+        distance2plane = dotProduct / thrustVectorNorm
+
+        ! compute distance to axis
+        rawTangent(1) = distance2center(2)*thrustVector(3) - distance2center(3)*thrustVector(2)
+        rawTangent(2) = distance2center(3)*thrustVector(1) - distance2center(1)*thrustVector(3)
+        rawTangent(3) = distance2center(1)*thrustVector(2) - distance2center(2)*thrustVector(1)
+
+        tangentNorm = sqrt(rawTangent(1)**2 + rawTangent(2)**2 + rawTangent(3)**2)
+
+        distance2axis = tangentNorm / thrustVectorNorm
+
+        ! compute tangential vector
+        tangent = rawTangent / tangentNorm
+
+
+    end subroutine computeCellSpatialMetrics
 
-    end subroutine computeActuatorRegionVolume
 
     ! ----------------------------------------------------------------------
     !                                                                      |
@@ -37,46 +59,83 @@ end subroutine computeActuatorRegionVolume
     ! ----------------------------------------------------------------------
 
 #ifndef USE_TAPENADE
-    subroutine addActuatorRegion(pts, conn, axis1, axis2, famName, famID, &
-                                 thrust, torque, heat, relaxStart, relaxEnd, nPts, nConn)
+    subroutine computeInitialSpatialMetrics(centerPoint, thrustVector, n, distance2plane, distance2axis, tangent)
+
+        use constants
+        use blockPointers, only: x, il, jl, kl, nDom, iBlank
+        use utils, only: setPointers
+        implicit none
+
+        ! Input variables
+        real(kind=realType), dimension(3), intent(in) :: centerPoint, thrustVector
+        integer(kind=intType), intent(in) :: n
+        real(kind=realType), dimension(n), intent(out) :: distance2plane, distance2axis
+        real(kind=realType), dimension(3, n), intent(out) :: tangent
+
+
+        ! Working variables
+        integer(kind=intType) :: ii, i, j, k, nn, iDim, level, sps
+        real(kind=realType), dimension(3) ::  xCen
+
+        ! Loop through all the cells and compute the distance of each cell center from the thrust axis (centerPoint + thrustVector)
+        ! also compute the distance from the plane formed by the centerPoint and the normal trhustVector
+        distance2axis = large
+        distance2plane = large
+
+        sps = 1
+        level = 1
+
+        ii = 0
+        do nn = 1, nDom
+            call setPointers(nn, level, sps)
+            do k = 2, kl
+                do j = 2, jl
+                    do i = 2, il
+                        ii = ii + 1
+
+                        ! cycle if it is not a compute cell
+                        if (iblank(i, j, k) /= 1) then
+                            cycle
+                        end if
+
+                        call computeCellSpatialMetrics(i, j, k, &
+                            centerPoint, thrustVector, &
+                            distance2plane(ii), distance2axis(ii), tangent(:, ii))
+                    end do
+                end do
+            end do
+        end do
+
+    end subroutine computeInitialSpatialMetrics
+
+
+    subroutine addActuatorRegion(flag, n, famName, famID, relaxStart, relaxEnd, iRegion, nLocalCells)
         ! Add a user-supplied integration surface.
 
         use communication, only: myID, adflow_comm_world
         use constants
-        use adtBuild, only: buildSerialQuad, destroySerialQuad
-        use adtLocalSearch, only: minDistanceTreeSearchSinglePoint
-        use ADTUtils, only: stack
-        use ADTData
         use blockPointers, only: x, il, jl, kl, nDom, iBlank, vol
         use adjointVars, only: nCellsLocal
         use utils, only: setPointers, EChk
         implicit none
 
         ! Input variables
-        integer(kind=intType), intent(in) :: nPts, nConn, famID
-        real(kind=realType), dimension(3, nPts), intent(in), target :: pts
-        integer(kind=intType), dimension(4, nConn), intent(in), target :: conn
-        real(kind=realType), intent(in), dimension(3) :: axis1, axis2
-        character(len=*) :: famName
-        real(kind=realType) :: thrust, torque, heat, relaxStart, relaxEnd
+        integer(kind=intType), intent(in) :: famID
+        character(len=*), intent(in) :: famName
+        real(kind=realType), intent(in) :: relaxStart, relaxEnd
+
+        integer(kind=intType), intent(in) :: n
+        integer(kind=intType), dimension(n), intent(in) :: flag
+
+        ! Output Variables
+        integer(kind=intType), intent(out) :: iRegion
+        integer(kind=intType), intent(out) :: nLocalCells
 
         ! Working variables
-        integer(kind=intType) :: i, j, k, nn, iDim, cellID, intInfo(3), sps, level, iii, ierr
-        real(kind=realType) :: dStar, frac, volLocal
+        integer(kind=intType) :: i, j, k, nn, sps, level, ii, iii, ierr
         type(actuatorRegionType), pointer :: region
-        real(kind=realType), dimension(3) :: minX, maxX, v1, v2, v3, xCen, axisVec
-        type(adtType) :: ADT
-        real(kind=realType) :: axisVecNorm
-        real(kind=realType), dimension(:, :), allocatable :: norm
-        integer(kind=intType), dimension(:), allocatable :: normCount
         integer(kind=intType), dimension(:, :), pointer :: tmp
 
-        ! ADT Type required data
-        integer(kind=intType), dimension(:), pointer :: frontLeaves, frontLeavesNew
-        type(adtBBoxTargetType), dimension(:), pointer :: BB
-        real(kind=realType) :: coor(4), uvw(5)
-        real(kind=realType) :: dummy(3, 2)
-
         nActuatorRegions = nActuatorRegions + 1
         if (nActuatorRegions > nActuatorRegionsMax) then
             print *, "Error: Exceeded the maximum number of actuatorDiskRegions.  "&
@@ -86,90 +145,15 @@ subroutine addActuatorRegion(pts, conn, axis1, axis2, famName, famID, &
 
         ! Save the input information
         region => actuatorRegions(nActuatorRegions)
+        iRegion = nActuatorRegions
         region%famName = famName
         region%famID = famID
-        region%torque = torque
-        region%heat = heat
         region%relaxStart = relaxStart
         region%relaxEnd = relaxEnd
-        ! We use the axis to define the direction of F. Since we are
-        ! dealing with rotating machinary, it is pretty good approximation
-        ! to assume that the thrust is going to be in the direction of the
-        ! axis.
-        axisVec = axis2 - axis1
-        axisVecNorm = sqrt((axisVec(1)**2 + axisvec(2)**2 + axisVec(3)**2))
-        if (axisVecNorm < 1e-12) then
-            print *, "Error: Axis cannot be determined by the supplied points. They are too close"
-            stop
-        end if
 
-        axisVec = axisVec / axisVecNorm
-
-        region%force = axisVec * thrust
-        region%axisVec = axisVec
 
         allocate (region%blkPtr(0:nDom))
         region%blkPtr(0) = 0
-
-        ! Next thing we need to do is to figure out if any of our cells
-        ! are inside the actuator disk region. If so we will save them in
-        ! the actuatorRegionType data structure
-
-        ! Since this is effectively a wall-distance calc it gets super
-        ! costly for the points far away. Luckly, we can do a fairly
-        ! simple shortcut: Just compute the bounding box of the region and
-        ! use that as the "already found" distance in the cloest point
-        ! search. This will eliminate all the points further away
-        ! immediately and this should be sufficiently fast.
-
-        ! So...compute that bounding box:
-        do iDim = 1, 3
-            minX(iDim) = minval(pts(iDim, :))
-            maxX(iDim) = maxval(pts(iDim, :))
-        end do
-
-        ! Get the max distance. This should be quite conservative.
-        dStar = (maxX(1) - minx(1))**2 + (maxX(2) - minX(2))**2 + (maxX(3) - minX(3))**2
-
-        ! Now build the tree.
-        call buildSerialQuad(size(conn, 2), size(pts, 2), pts, conn, ADT)
-
-        ! Compute the (averaged) unique nodal vectors:
-        allocate (norm(3, size(pts, 2)), normCount(size(pts, 2)))
-
-        norm = zero
-        normCount = 0
-
-        do i = 1, size(conn, 2)
-
-            ! Compute cross product normal and normalize
-            v1 = pts(:, conn(3, i)) - pts(:, conn(1, i))
-            v2 = pts(:, conn(4, i)) - pts(:, conn(2, i))
-
-            v3(1) = (v1(2) * v2(3) - v1(3) * v2(2))
-            v3(2) = (v1(3) * v2(1) - v1(1) * v2(3))
-            v3(3) = (v1(1) * v2(2) - v1(2) * v2(1))
-            v3 = v3 / sqrt(v3(1)**2 + v3(2)**2 + v3(3)**2)
-
-            ! Add to each of the four pts and increment the number added
-            do j = 1, 4
-                norm(:, conn(j, i)) = norm(:, conn(j, i)) + v3
-                normCount(conn(j, i)) = normCount(conn(j, i)) + 1
-            end do
-        end do
-
-        ! Now just divide by the norm count
-        do i = 1, size(pts, 2)
-            norm(:, i) = norm(:, i) / normCount(i)
-        end do
-
-        ! Norm count is no longer needed
-        deallocate (normCount)
-
-        ! Allocate the extra data the tree search requires.
-        allocate (stack(100), BB(20), frontLeaves(25), frontLeavesNew(25))
-
-        ! Allocate sufficient space for the maximum possible number of cellIDs
         allocate (region%cellIDs(3, nCellsLocal(1)))
 
         ! Now search for all the coordinate. Note that We have explictly
@@ -177,37 +161,26 @@ subroutine addActuatorRegion(pts, conn, axis1, axis2, famName, famID, &
         sps = 1
         level = 1
 
+        ii = 0
         do nn = 1, nDom
             call setPointers(nn, level, sps)
             do k = 2, kl
                 do j = 2, jl
                     do i = 2, il
-                        ! Only check real cells
-                        if (iblank(i, j, k) == 1) then
-                            ! Compute the cell center
-                            xCen = eighth * (x(i - 1, j - 1, k - 1, :) + x(i, j - 1, k - 1, :) &
-                                             + x(i - 1, j, k - 1, :) + x(i, j, k - 1, :) &
-                                             + x(i - 1, j - 1, k, :) + x(i, j - 1, k, :) &
-                                             + x(i - 1, j, k, :) + x(i, j, k, :))
-
-                            ! The current point to search for and continually
-                            ! reset the "closest point already found" variable.
-                            coor(1:3) = xCen
-                            coor(4) = dStar
-                            intInfo(3) = 0
-                            call minDistancetreeSearchSinglePoint(ADT, coor, intInfo, &
-                                                                  uvw, dummy, 0, BB, frontLeaves, frontLeavesNew)
-                            cellID = intInfo(3)
-                            if (cellID > 0) then
-                                ! Now check if this was successful or now:
-                                if (checkInside()) then
-                                    ! Whoohoo! We are inside the region. Add this cell
-                                    ! to the list.
-                                    region%nCellIDs = region%nCellIDs + 1
-                                    region%cellIDs(:, region%nCellIDs) = (/i, j, k/)
-                                end if
-                            end if
+                        ii = ii + 1
+
+                        ! cycle if it is not a compute cell
+                        if (iblank(i, j, k) /= 1) then
+                            cycle
+                        end if
+
+                        ! cycle if cell was not tagged
+                        if (flag(ii) /= 1) then
+                            cycle
                         end if
+
+                        region%nCellIDs = region%nCellIDs + 1
+                        region%cellIDs(:, region%nCellIDs) = (/i, j, k/)
                     end do
                 end do
             end do
@@ -218,69 +191,103 @@ subroutine addActuatorRegion(pts, conn, axis1, axis2, famName, famID, &
 
         end do
 
+
         ! Resize the cellIDs to the correct size now that we know the
         ! correct exact number.
         tmp => region%cellIDs
         allocate (region%cellIDs(3, region%nCellIDs))
         region%cellIDs = tmp(:, 1:region%nCellIDs)
         deallocate (tmp)
+        nLocalCells = region%nCellIDs
+
+
+        ! Allocate variables for the source terms
+        allocate (region%force(3, region%nCellIDs))
+        allocate (region%heat(region%nCellIDs))
+
+
+        ! same for derivative seeds
+        allocate (actuatorRegionsd(nActuatorRegions)%force(3, region%nCellIDs))
+        allocate (actuatorRegionsd(nActuatorRegions)%heat(region%nCellIDs))
+
+    end subroutine addActuatorRegion
+
+    subroutine computeSpatialMetrics(iRegion, nLocalCells, centerPoint, thrustVector, &
+            distance2plane, distance2axis, tangent, volume)
+        use constants
+        use blockPointers, only: nDom, vol
+        use utils, only: setPointers
+        implicit none
+
+        ! Inputs
+        real(kind=realType), dimension(3), intent(in) :: centerPoint, thrustVector
+        integer(kind=intType), intent(in) :: nLocalCells
+        integer(kind=intType), intent(in) :: iRegion
+
+        ! Outputs
+        real(kind=realType), intent(out), dimension(nLocalCells) :: distance2axis, distance2plane, volume
+        real(kind=realType), dimension(3, nLocalCells), intent(out) :: tangent
+
+        ! Working
+        type(actuatorRegionType), pointer :: region
+        integer(kind=intType) :: sps, level, nn, i, j, k, iii
 
-        ! Now go back and generate the total volume of the the cells we've flagged
-        volLocal = zero
+        region => actuatorRegions(iRegion)
 
+
+        sps = 1
+        level = 1
         do nn = 1, nDom
             call setPointers(nn, level, sps)
+            do iii = actuatorRegions(iRegion)%blkPtr(nn - 1) + 1, actuatorRegions(iRegion)%blkPtr(nn)
+                i = actuatorRegions(iRegion)%cellIDs(1, iii)
+                j = actuatorRegions(iRegion)%cellIDs(2, iii)
+                k = actuatorRegions(iRegion)%cellIDs(3, iii)
+
+
+                call computeCellSpatialMetrics(i, j, k, &
+                    centerPoint, thrustVector, &
+                    distance2plane(iii), distance2axis(iii), tangent(:, iii))
 
-            ! Loop over the region for this block
-            do iii = region%blkPtr(nn - 1) + 1, region%blkPtr(nn)
-                i = region%cellIDs(1, iii)
-                j = region%cellIDs(2, iii)
-                k = region%cellIDs(3, iii)
-                volLocal = volLocal + vol(i, j, k)
+                volume(iii) = vol(i, j, k)
             end do
         end do
 
-        call mpi_allreduce(volLocal, region%volume, 1, adflow_real, &
-                           MPI_SUM, adflow_comm_world, ierr)
-        call ECHK(ierr, __FILE__, __LINE__)
-
-        ! Final memory cleanup
-        deallocate (stack, norm, frontLeaves, frontLeavesNew, BB)
-        call destroySerialQuad(ADT)
-    contains
-
-        function checkInside()
-
-            implicit none
-            logical :: checkInside
-            integer(kind=intType) :: jj
-            real(kind=realType) :: shp(4), xp(3), normal(3), v1(3), dp
-
-            ! bi-linear shape functions (CCW ordering)
-            shp(1) = (one - uvw(1)) * (one - uvw(2))
-            shp(2) = (uvw(1)) * (one - uvw(2))
-            shp(3) = (uvw(1)) * (uvw(2))
-            shp(4) = (one - uvw(1)) * (uvw(2))
-
-            xp = zero
-            normal = zero
-            do jj = 1, 4
-                xp = xp + shp(jj) * pts(:, conn(jj, cellID))
-                normal = normal + shp(jj) * norm(:, conn(jj, cellID))
+    end subroutine computeSpatialMetrics
+
+    subroutine populateBCValues(iRegion, nLocalCells, force, heat)
+        use constants
+        use blockPointers, only: nDom
+        use utils, only: setPointers
+        implicit none
+
+        ! Inputs
+        integer(kind=intType), intent(in) :: iRegion
+        integer(kind=intType), intent(in) :: nLocalCells
+        real(kind=realType), dimension(3, nLocalCells), intent(in) :: force
+        real(kind=realType), dimension(nLocalCells), intent(in) :: heat
+
+
+        ! Working
+        type(actuatorRegionType), pointer :: region
+        integer(kind=intType) :: sps, level, nn, i, j, k, iii
+
+
+        region => actuatorRegions(iRegion)
+
+
+        sps = 1
+        level = 1
+        do nn = 1, nDom
+            call setPointers(nn, level, sps)
+            do iii = actuatorRegions(iRegion)%blkPtr(nn - 1) + 1, actuatorRegions(iRegion)%blkPtr(nn)
+                actuatorRegions(iRegion)%force(:, iii) = force(:, iii)
+                actuatorRegions(iRegion)%heat(iii) = heat(iii)
             end do
+        end do
 
-            ! Compute the dot product of normal with cell center
-            ! (stored in coor) with the point on the surface.
-            v1 = coor(1:3) - xp
-            dp = normal(1) * v1(1) + normal(2) * v1(2) + normal(3) * v1(3)
+    end subroutine populateBCValues
 
-            if (dp < zero) then
-                checkInside = .True.
-            else
-                checkInside = .False.
-            end if
-        end function checkInside
-    end subroutine addActuatorRegion
 
     subroutine writeActuatorRegions(fileName)
 
diff --git a/src/solver/residuals.F90 b/src/solver/residuals.F90
index 44d8a39a6..9e7dcc7ed 100644
--- a/src/solver/residuals.F90
+++ b/src/solver/residuals.F90
@@ -364,7 +364,7 @@ subroutine sourceTerms_block(nn, res, iRegion, pLocal)
 
         ! Working
         integer(kind=intType) :: i, j, k, ii, iStart, iEnd
-        real(kind=realType) :: Ftmp(3), Vx, Vy, Vz, F_fact(3), Q_fact, Qtmp, reDim, factor, oStart, oEnd
+        real(kind=realType) :: Vx, Vy, Vz, Fx, Fy, Fz, Q, reDim, factor, oStart, oEnd
 
         reDim = pRef * uRef
 
@@ -381,12 +381,6 @@ subroutine sourceTerms_block(nn, res, iRegion, pLocal)
             factor = (ordersConverged - oStart) / (oEnd - oStart)
         end if
 
-        ! Compute the constant force factor
-        F_fact = factor * actuatorRegions(iRegion)%force / actuatorRegions(iRegion)%volume / pRef
-
-        ! Heat factor. This is heat added per unit volume per unit time
-        Q_fact = factor * actuatorRegions(iRegion)%heat / actuatorRegions(iRegion)%volume / (pRef * uRef * LRef * LRef)
-
         ! Loop over the ranges for this block
         iStart = actuatorRegions(iRegion)%blkPtr(nn - 1) + 1
         iEnd = actuatorRegions(iRegion)%blkPtr(nn)
@@ -400,25 +394,29 @@ subroutine sourceTerms_block(nn, res, iRegion, pLocal)
             k = actuatorRegions(iRegion)%cellIDs(3, ii)
 
             ! This actually gets the force
-            FTmp = vol(i, j, k) * F_fact
+            Fx = factor * actuatorRegions(iRegion)%force(1, ii) / pRef
+            Fy = factor * actuatorRegions(iRegion)%force(2, ii) / pRef
+            Fz = factor * actuatorRegions(iRegion)%force(3, ii) / pRef
 
             Vx = w(i, j, k, iVx)
             Vy = w(i, j, k, iVy)
             Vz = w(i, j, k, iVz)
 
             ! this gets the heat addition rate
-            QTmp = vol(i, j, k) * Q_fact
+            Q = factor * actuatorRegions(iRegion)%heat(ii)  / (pRef * uRef * LRef * LRef)
 
             if (res) then
                 ! Momentum residuals
-                dw(i, j, k, imx:imz) = dw(i, j, k, imx:imz) - Ftmp
+                dw(i, j, k, imx) = dw(i, j, k, imx) - Fx
+                dw(i, j, k, imy) = dw(i, j, k, imy) - Fy
+                dw(i, j, k, imz) = dw(i, j, k, imz) - Fz
 
                 ! energy residuals
                 dw(i, j, k, iRhoE) = dw(i, j, k, iRhoE) - &
-                                     Ftmp(1) * Vx - Ftmp(2) * Vy - Ftmp(3) * Vz - Qtmp
+                                     Fx * Vx - Fy * Vy - Fz * Vz - Q
             else
                 ! Add in the local power contribution:
-                pLocal = pLocal + (Vx * Ftmp(1) + Vy * FTmp(2) + Vz * Ftmp(3)) * reDim
+                pLocal = pLocal + (Vx * Fx + Vy * Fy + Vz * Fz) * reDim
             end if
         end do
 
diff --git a/src/solver/surfaceIntegrations.F90 b/src/solver/surfaceIntegrations.F90
index 357bc18d4..e236eca7d 100644
--- a/src/solver/surfaceIntegrations.F90
+++ b/src/solver/surfaceIntegrations.F90
@@ -159,10 +159,10 @@ subroutine getCostFunctions(globalVals, funcValues)
                 funcValues(costFuncSepSensorKs) = funcValues(costFuncSepSensorKs) + ks_comp
 
                 funcValues(costFuncsepSensorKsArea) = funcValues(costFuncsepSensorKsArea) + &
-                                                    ovrNTS * globalVals(iSepSensorKsArea, sps) * ks_comp * &
-                                                    one / (one + exp(2 * sepSensorKsSharpness &
-                                                                     * (ks_comp + sepSensorKsOffset))) + &
-                                                    ovrNTS * globalVals(iSepSensorArea, sps)
+                                                      ovrNTS * globalVals(iSepSensorKsArea, sps) * ks_comp * &
+                                                      one / (one + exp(2 * sepSensorKsSharpness &
+                                                                       * (ks_comp + sepSensorKsOffset))) + &
+                                                      ovrNTS * globalVals(iSepSensorArea, sps)
 
             end if
 
diff --git a/tests/reg_tests/test_actuator.py b/tests/reg_tests/test_actuator.py
index 25364fa90..922bb2e54 100644
--- a/tests/reg_tests/test_actuator.py
+++ b/tests/reg_tests/test_actuator.py
@@ -8,7 +8,7 @@
 from mpi4py import MPI
 
 # MACH classes
-from adflow import ADFLOW
+from adflow import ADFLOW, UniformActuatorRegion
 from idwarp import USMesh
 from pygeo import DVGeometry
 from adflow import ADFLOW_C
@@ -91,16 +91,19 @@ def setUp(self):
         # this is imported from reg_aeroproblems utility script
         self.ap = ap_actuator_pipe
 
-        actuatorFile = os.path.join(baseDir, "../../input_files/actuator_test_disk.xyz")
+        actuator_region = UniformActuatorRegion(
+            centerPoint=np.array([-0.38, 0, 0.3]),
+            thrustVector=np.array([-1, 0, 0]),
+            innerDiameter=0.0,
+            outerDiameter=0.354,
+            regionDepth=0.12,
+            thrust=0,
+            heat=0,
+        )
+
         self.CFDSolver.addActuatorRegion(
-            actuatorFile,
-            np.array([0, 0, 0]),
-            np.array([1, 0, 0]),
+            actuator_region,
             "actuator_region",
-            # we will set these individually in the tests below
-            thrust=0.0,
-            torque=0.0,
-            heat=0.0,
         )
 
         # add thrust and heat as AP DVs
@@ -681,6 +684,7 @@ def setUp(self):
         actuatorFile = os.path.join(baseDir, "../../input_files/actuator_test_disk.xyz")
         self.CFDSolver.addActuatorRegion(
             actuatorFile,
+            "uniform",
             np.array([0, 0, 0]),
             np.array([1, 0, 0]),
             "actuator_region",