diff --git a/flow/src/wAdjoint.cpp b/flow/src/wAdjoint.cpp
index a5d9be24ef513795d82f9b6fce75f29a7433ccca..11eebf9865a1b6cad58ff8ec79553f6f884c2c13 100644
--- a/flow/src/wAdjoint.cpp
+++ b/flow/src/wAdjoint.cpp
@@ -267,8 +267,8 @@ void Adjoint::buildGradientLoadsFlow(Eigen::VectorXd &dL, Eigen::VectorXd &dD)
             tbb::spin_mutex::scoped_lock lock(mutex);
             for (size_t i = 0; i < e->nodes.size(); ++i)
             {
-                Eigen::RowVectorXd AeL = -dirL * Ae.block(i * 3, 0, 3, eV->nodes.size());
-                Eigen::RowVectorXd AeD = -dirD * Ae.block(i * 3, 0, 3, eV->nodes.size());
+                Eigen::RowVectorXd AeL = dirL * Ae.block(i * 3, 0, 3, eV->nodes.size());
+                Eigen::RowVectorXd AeD = dirD * Ae.block(i * 3, 0, 3, eV->nodes.size());
                 for (size_t j = 0; j < eV->nodes.size(); ++j)
                 {
                     Node *nodj = eV->nodes[j];
@@ -342,7 +342,7 @@ void Adjoint::buildGradientLoadsAoa(double &dL, double &dD)
     Eigen::Vector3d Cf(0, 0, 0);
     for (auto bnd : sol->pbl->bnds)
         for (size_t i = 0; i < bnd->nodes.size(); ++i)
-            Cf -= bnd->nLoads[i];
+            Cf += bnd->nLoads[i];
     // Rotate to gradient of flow direction
     dD = Cf.dot(dirD);
     dL = Cf.dot(dirL);
@@ -533,10 +533,10 @@ void Adjoint::buildGradientLoadsMesh(Eigen::RowVectorXd &dL, Eigen::RowVectorXd
             for (size_t i = 0; i < e->nodes.size(); ++i)
             {
                 // rotate to flow direction
-                Eigen::RowVectorXd AesL = -dirL * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size());
-                Eigen::RowVectorXd AesD = -dirD * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size());
-                Eigen::RowVectorXd AevL = -dirL * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size());
-                Eigen::RowVectorXd AevD = -dirD * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size());
+                Eigen::RowVectorXd AesL = dirL * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size());
+                Eigen::RowVectorXd AesD = dirD * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size());
+                Eigen::RowVectorXd AevL = dirL * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size());
+                Eigen::RowVectorXd AevD = dirD * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size());
                 for (size_t m = 0; m < sol->pbl->nDim; ++m)
                 {
                     // surface
diff --git a/flow/src/wLoadFunctional.cpp b/flow/src/wLoadFunctional.cpp
index 6c0bbe17b7325badaadfe0f851fc64b87873e393..1113e3d8ed0702e7d9bae259ed983873fefdc6b1 100644
--- a/flow/src/wLoadFunctional.cpp
+++ b/flow/src/wLoadFunctional.cpp
@@ -36,7 +36,7 @@ Eigen::VectorXd LoadFunctional::build(Element const &e, Element const &eV, std::
     Gauss &gauss = cache.getVGauss();
 
     // Pressure force coefficient in associated volume
-    Eigen::Vector3d fcp = cp.eval(eV, phi, 0) * e.getNormal();
+    Eigen::Vector3d fcp = -cp.eval(eV, phi, 0) * e.getNormal(); // "-" because force on body not on fluid
 
     // Integrate cp on the surface
     Eigen::VectorXd b = Eigen::VectorXd::Zero(3 * e.nodes.size());
@@ -62,7 +62,7 @@ Eigen::MatrixXd LoadFunctional::buildGradientFlow(Element const &e, Element cons
     Gauss &gauss = cache.getVGauss();
 
     // Gradient of pressure force coefficient in associated volume
-    Eigen::MatrixXd dfcp = e.getNormal() * cp.evalD(eV, phi, 0) * eV.computeGradient(phi, 0).transpose() * eV.getJinv(0) * eV.getVCache().getDsf(0);
+    Eigen::MatrixXd dfcp = -e.getNormal() * cp.evalD(eV, phi, 0) * eV.computeGradient(phi, 0).transpose() * eV.getJinv(0) * eV.getVCache().getDsf(0); // "-" because force on body not on fluid
 
     // Integrate cp on the surface
     Eigen::MatrixXd A = Eigen::MatrixXd::Zero(3 * e.nodes.size(), eV.nodes.size());
@@ -91,10 +91,10 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> LoadFunctional::buildGradientMesh(E
     Gauss &gauss = cache.getVGauss();
 
     // Variables
-    double cP = cp.eval(eV, phi, 0);     // pressure coefficient
+    double cP = -cp.eval(eV, phi, 0);    // pressure coefficient ("-" because force on body not on fluid)
     Eigen::Vector3d nrm = e.getNormal(); // unit normal
     // Gradients
-    double dCp = cp.evalD(eV, phi, 0);                       // pressure coefficient
+    double dCp = -cp.evalD(eV, phi, 0);                      // pressure coefficient ("-" because force on body not on fluid)
     Eigen::VectorXd dPhi = eV.computeGradient(phi, 0);       // potential
     Eigen::MatrixXd const &iJ = eV.getJinv(0);               // inverse Jacobian
     std::vector<Eigen::MatrixXd> const &dJ = eV.getGradJ(0); // Jacobian
diff --git a/flow/src/wSolver.cpp b/flow/src/wSolver.cpp
index 1b15b0495f31c1a4316d4127ee95800d050293d5..7fae0493b0a2c12bdee06baf62efc0fe6f609eb6 100644
--- a/flow/src/wSolver.cpp
+++ b/flow/src/wSolver.cpp
@@ -197,8 +197,8 @@ void Solver::computeLoad()
         {
             Eigen::Vector3d l = sur->nodes[i]->pos - pbl->x_ref; // lever arm
             Eigen::Vector3d cf = sur->nLoads[i];                 // normalized force
-            Cf -= cf;
-            sur->Cm -= (l(pbl->nDim - 1) * cf(0) - l(0) * cf(pbl->nDim - 1)) / (pbl->S_ref * pbl->c_ref); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
+            Cf += cf;
+            sur->Cm += (l(pbl->nDim - 1) * cf(0) - l(0) * cf(pbl->nDim - 1)) / (pbl->S_ref * pbl->c_ref); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
         }
         // rotate to flow direction
         sur->Cd = Cf.dot(dirD);