diff --git a/src/adjoint/DAObjFunc/DAObjFunc.C b/src/adjoint/DAObjFunc/DAObjFunc.C
index 89a8c5bc..22c7c999 100755
--- a/src/adjoint/DAObjFunc/DAObjFunc.C
+++ b/src/adjoint/DAObjFunc/DAObjFunc.C
@@ -59,7 +59,7 @@ DAObjFunc::DAObjFunc(
     calcRefDiffSquare_ = objFuncDict_.lookupOrDefault<label>("calcRefDiffSquare", 0);
     if (calcRefDiffSquare_)
     {
-        objFuncDict_.readEntry<scalar>("ref", ref_);
+        objFuncDict_.readEntry<scalarList>("ref", ref_);
     }
 }
 
@@ -318,6 +318,29 @@ scalar DAObjFunc::getObjFuncValue()
     return objFuncValue_;
 }
 
+void DAObjFunc::calcRefDiffSquare(scalar& objFuncValue)
+{
+    /*
+    Description:
+        Call the variable difference with respect to a given reference and take a square of it.
+        This can be used in FIML. This function is for calcRefDiffSquare_ == 1
+    */
+
+    if (calcRefDiffSquare_)
+    {
+        if (ref_.size() == 1)
+        {
+            objFuncValue = (objFuncValue - ref_[0]) * (objFuncValue - ref_[0]);
+        }
+        else
+        {
+            label idxI = mesh_.time().timeIndex() - 1;
+            objFuncValue = (objFuncValue - ref_[idxI]) * (objFuncValue - ref_[idxI]);
+        }
+    }
+
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/adjoint/DAObjFunc/DAObjFunc.H b/src/adjoint/DAObjFunc/DAObjFunc.H
index eb62b229..f6bf097c 100755
--- a/src/adjoint/DAObjFunc/DAObjFunc.H
+++ b/src/adjoint/DAObjFunc/DAObjFunc.H
@@ -95,9 +95,8 @@ protected:
     /// whether to calculate (obj-ref)^2
     label calcRefDiffSquare_;
 
-    /// if calcRefDiffSquare_ is True, set the reference value
-    scalar ref_;
-
+    /// if calcRefDiffSquare_ is True, set the reference value list
+    scalarList ref_;
 
 public:
     /// Runtime type information
@@ -243,6 +242,9 @@ public:
 
     /// whether to compute reference coefficients for special objFunc treatment such as totalPressureRatio
     label calcRefCoeffs = 1;
+
+    /// calculate the varDiffSquare
+    void calcRefDiffSquare(scalar& objFuncValue);
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/adjoint/DAObjFunc/DAObjFuncCenterOfPressure.C b/src/adjoint/DAObjFunc/DAObjFuncCenterOfPressure.C
index 62f97a0d..d3df0556 100644
--- a/src/adjoint/DAObjFunc/DAObjFuncCenterOfPressure.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncCenterOfPressure.C
@@ -128,10 +128,8 @@ void DAObjFuncCenterOfPressure::calcObjFunc(
 
     objFuncValue = scale_ * (weightedPressure / totalPressure + (center_ & axis_));
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }
diff --git a/src/adjoint/DAObjFunc/DAObjFuncForce.C b/src/adjoint/DAObjFunc/DAObjFuncForce.C
index 7d5beaf8..06b3cd30 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncForce.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncForce.C
@@ -191,10 +191,8 @@ void DAObjFuncForce::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }
diff --git a/src/adjoint/DAObjFunc/DAObjFuncLocation.C b/src/adjoint/DAObjFunc/DAObjFuncLocation.C
index db848099..5587ed4f 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncLocation.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncLocation.C
@@ -219,10 +219,8 @@ void DAObjFuncLocation::calcObjFunc(
                                           << abort(FatalError);
     }
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }
diff --git a/src/adjoint/DAObjFunc/DAObjFuncMass.C b/src/adjoint/DAObjFunc/DAObjFuncMass.C
index a2d34599..1a54be48 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncMass.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncMass.C
@@ -96,10 +96,8 @@ void DAObjFuncMass::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }
diff --git a/src/adjoint/DAObjFunc/DAObjFuncMassFlowRate.C b/src/adjoint/DAObjFunc/DAObjFuncMassFlowRate.C
index 4bb8b9a8..db9dbb3d 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncMassFlowRate.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncMassFlowRate.C
@@ -119,10 +119,8 @@ void DAObjFuncMassFlowRate::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }
diff --git a/src/adjoint/DAObjFunc/DAObjFuncMeshQualityKS.C b/src/adjoint/DAObjFunc/DAObjFuncMeshQualityKS.C
index 7606313a..9f58f497 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncMeshQualityKS.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncMeshQualityKS.C
@@ -219,10 +219,8 @@ void DAObjFuncMeshQualityKS::calcObjFunc(
 
     objFuncValue = log(objFuncValue) / coeffKS_;
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }
diff --git a/src/adjoint/DAObjFunc/DAObjFuncMoment.C b/src/adjoint/DAObjFunc/DAObjFuncMoment.C
index 640b0c21..fbfd062c 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncMoment.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncMoment.C
@@ -165,10 +165,8 @@ void DAObjFuncMoment::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }
diff --git a/src/adjoint/DAObjFunc/DAObjFuncPatchMean.C b/src/adjoint/DAObjFunc/DAObjFuncPatchMean.C
index ad166c43..2476f2b2 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncPatchMean.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncPatchMean.C
@@ -132,10 +132,8 @@ void DAObjFuncPatchMean::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }
diff --git a/src/adjoint/DAObjFunc/DAObjFuncPower.C b/src/adjoint/DAObjFunc/DAObjFuncPower.C
index 513306f2..1fdaf56f 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncPower.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncPower.C
@@ -168,10 +168,8 @@ void DAObjFuncPower::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }
diff --git a/src/adjoint/DAObjFunc/DAObjFuncTotalPressure.C b/src/adjoint/DAObjFunc/DAObjFuncTotalPressure.C
index b1566ab8..5f28fa1e 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncTotalPressure.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncTotalPressure.C
@@ -132,10 +132,8 @@ void DAObjFuncTotalPressure::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }
diff --git a/src/adjoint/DAObjFunc/DAObjFuncTotalPressureRatio.C b/src/adjoint/DAObjFunc/DAObjFuncTotalPressureRatio.C
index c2a6f8fe..34fe7b4b 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncTotalPressureRatio.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncTotalPressureRatio.C
@@ -217,10 +217,8 @@ void DAObjFuncTotalPressureRatio::calcObjFunc(
         }
     }
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }
diff --git a/src/adjoint/DAObjFunc/DAObjFuncTotalTemperatureRatio.C b/src/adjoint/DAObjFunc/DAObjFuncTotalTemperatureRatio.C
index c8c32fdd..9e0c1f5e 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncTotalTemperatureRatio.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncTotalTemperatureRatio.C
@@ -207,10 +207,8 @@ void DAObjFuncTotalTemperatureRatio::calcObjFunc(
         }
     }
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }
diff --git a/src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C b/src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C
index 0cd576a0..a95d5f23 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C
@@ -167,10 +167,8 @@ void DAObjFuncVariableVolSum::calcObjFunc(
 
     objFuncValue /= totalVol;
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }
diff --git a/src/adjoint/DAObjFunc/DAObjFuncVonMisesStressKS.C b/src/adjoint/DAObjFunc/DAObjFuncVonMisesStressKS.C
index 20771c81..5e3a0969 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncVonMisesStressKS.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncVonMisesStressKS.C
@@ -121,10 +121,8 @@ void DAObjFuncVonMisesStressKS::calcObjFunc(
 
     objFuncValue = log(objValTmp) / coeffKS_;
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }
diff --git a/src/adjoint/DAObjFunc/DAObjFuncWallHeatFlux.C b/src/adjoint/DAObjFunc/DAObjFuncWallHeatFlux.C
index c9a633f9..73108496 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncWallHeatFlux.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncWallHeatFlux.C
@@ -246,10 +246,8 @@ void DAObjFuncWallHeatFlux::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }
diff --git a/tests/refs/DAFoam_Test_DAPimpleFoamRef.txt b/tests/refs/DAFoam_Test_DAPimpleFoamRef.txt
index 961f0683..4f26c994 100755
--- a/tests/refs/DAFoam_Test_DAPimpleFoamRef.txt
+++ b/tests/refs/DAFoam_Test_DAPimpleFoamRef.txt
@@ -1,5 +1,7 @@
 Dictionary Key: CD
 @value    7.747035248407224 1e-08 1e-10
+Dictionary Key: CMZVAR
+@value    41.05825112111668 1e-08 1e-10
 Dictionary Key: fail
 @value                    0 1e-08 1e-10
 Dictionary Key: CD
@@ -25,5 +27,29 @@ Dictionary Key: twist
 @value   0.5092926662504045 0.0001 1e-06
 Dictionary Key: uin
 @value   0.1887121251618032 0.0001 1e-06
+Dictionary Key: CMZVAR
+Dictionary Key: actuator
+@value   0.0129164206206973 0.0001 1e-06
+@value -0.06721779867591618 0.0001 1e-06
+@value 1.01218774787428e-06 0.0001 1e-06
+@value                    0 0.0001 1e-06
+@value -0.03989829709014971 0.0001 1e-06
+@value -1.174777939902339e-05 0.0001 1e-06
+@value -0.01318488544862298 0.0001 1e-06
+@value  0.02098447283371081 0.0001 1e-06
+@value 7.847578532670927e-05 0.0001 1e-06
+@value 1.135345283703712e-06 0.0001 1e-06
+@value -0.00490063870041851 0.0001 1e-06
+@value -0.007270668470157968 0.0001 1e-06
+@value                    0 0.0001 1e-06
+Dictionary Key: alpha
+@value   0.1694703610086252 0.0001 1e-06
+Dictionary Key: parameter
+@value  0.01322002733662093 0.0001 1e-06
+Dictionary Key: twist
+@value    2.800461713324844 0.0001 1e-06
+Dictionary Key: uin
+@value   0.5352121060453168 0.0001 1e-06
+
 Dictionary Key: fail
 @value                    0 0.0001 1e-06
\ No newline at end of file
diff --git a/tests/runTests_DAPimpleFoam.py b/tests/runTests_DAPimpleFoam.py
index d2945075..23f8a2f6 100755
--- a/tests/runTests_DAPimpleFoam.py
+++ b/tests/runTests_DAPimpleFoam.py
@@ -95,6 +95,19 @@
                 "addToAdjoint": False,
             }
         },
+        "CMZVAR": {
+            "part1": {
+                "type": "moment",
+                "source": "patchToFace",
+                "patches": ["wing"],
+                "axis": [0.0, 0.0, 1.0],
+                "center": [0.25, 0.0, 0.05],
+                "scale": 1.0,
+                "addToAdjoint": True,
+                "calcRefDiffSquare": True,
+                "ref": [0.1, 0.05, 0.04, 0.02, 0.02, 0.01, 0.0, -0.01, -0.01, -0.02]
+            }
+        }
     },
     "adjStateOrdering": "cell",
     "adjEqnOption": {
@@ -267,6 +280,9 @@ def regModel(val, geo):
     parameterNormU = np.linalg.norm(funcsSens["CD"]["parameter"])
     funcsSens["CD"]["parameter"] = parameterNormU
 
+    parameterNormM = np.linalg.norm(funcsSens["CMZVAR"]["parameter"])
+    funcsSens["CMZVAR"]["parameter"] = parameterNormM
+
     if gcomm.rank == 0:
         reg_write_dict(funcs, 1e-8, 1e-10)
         reg_write_dict(funcsSens, 1e-4, 1e-6)
diff --git a/tests/runTests_DASimpleFoamAD.py b/tests/runTests_DASimpleFoamAD.py
index 760c9696..8dcd88bb 100755
--- a/tests/runTests_DASimpleFoamAD.py
+++ b/tests/runTests_DASimpleFoamAD.py
@@ -131,7 +131,7 @@
                 "scale": 1.0 / (0.5 * U0 * U0 * A0 * 1.0),
                 "addToAdjoint": True,
                 "calcRefDiffSquare": True,
-                "ref": 0.1
+                "ref": [0.1]
             }
         },
     },