diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp
index 45319be..a5f709a 100644
--- a/src/pmpo_MPMesh.cpp
+++ b/src/pmpo_MPMesh.cpp
@@ -301,6 +301,7 @@ void MPMesh::reconstructSlices() {
     if (reconstructSlice.size() == 0) return;
     Kokkos::Timer timer;
     calcBasis();
+    resetPreComputeFlag();
     for (auto const& [index, reconstruct] : reconstructSlice) {
         if (reconstruct) reconstruct();
     }
diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp
index 5a63ef9..33f0599 100644
--- a/src/pmpo_MPMesh.hpp
+++ b/src/pmpo_MPMesh.hpp
@@ -16,7 +16,17 @@ template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_RotLatLonIn
 #define maxMPsPerElm 8
 
 class MPMesh{
+  private:
+   
+    bool isPreComputed;
+  
   public:
+    
+    MPMesh() : isPreComputed(false){};
+    void computeMatricesAndSolve(); 
+    void resetPreComputeFlag();
+    Kokkos::View<double*[vec4d_nEntries]> precomputedVtxCoeffs;
+
     Mesh* p_mesh;
     MaterialPoints* p_MPs;
 
@@ -49,7 +59,7 @@ class MPMesh{
     void assemblyElm0();
     template <MeshFieldIndex meshFieldIndex>
     void assemblyVtx1();
-
+    
     template<MeshFieldIndex meshFieldIndex>
     void setReconstructSlice(int order, MeshFieldType type);
     void reconstructSlices();
diff --git a/src/pmpo_MPMesh_assembly.hpp b/src/pmpo_MPMesh_assembly.hpp
index e5dedb5..25899b9 100644
--- a/src/pmpo_MPMesh_assembly.hpp
+++ b/src/pmpo_MPMesh_assembly.hpp
@@ -98,35 +98,36 @@ void MPMesh::assemblyElm0() {
   pumipic::RecordTime("PolyMPO_Reconstruct_Elm0", timer.seconds());
 }
 
-template <MeshFieldIndex meshFieldIndex>
-void MPMesh::assemblyVtx1() {
+
+void MPMesh::resetPreComputeFlag(){
+  isPreComputed = false;
+}
+
+void MPMesh::computeMatricesAndSolve(){
+
   //Mesh Information
   auto elm2VtxConn = p_mesh->getElm2VtxConn();  
   int numVtx = p_mesh->getNumVertices();
   auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
 
-  //Mesh Field
-  constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice<meshFieldIndex>;
-  const int numEntries = mpSliceToNumEntries<mpfIndex>();
-  p_mesh->fillMeshField<meshFieldIndex>(numVtx, numEntries, 0.0);
-  auto meshField = p_mesh->getMeshField<meshFieldIndex>();
+  //Dual Element Area for Regularization
+  auto dual_triangle_area=p_mesh->getMeshField<MeshF_DualTriangleArea>();
 
   //Material Points
-  auto mpData = p_MPs->getData<mpfIndex>();
   auto weight = p_MPs->getData<MPF_Basis_Vals>();
   auto mpPositions = p_MPs->getData<MPF_Cur_Pos_XYZ>();
 
   //Matrix for each vertex
   Kokkos::View<double*[vec4d_nEntries][vec4d_nEntries]> VtxMatrices("VtxMatrices", p_mesh->getNumVertices());
 
-  //Reconstructed values
-  Kokkos::View<double**> reconVals("meshField", p_mesh->getNumVertices(), numEntries);
-  
   //Earth Radius
   double radius = 1.0;
   if(p_mesh->getGeomType() == geom_spherical_surf)
     radius=p_mesh->getSphereRadius();
 
+  bool scaling=true;
+  int reg_method = 2;
+  
   //Assemble matrix for each vertex
   auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
     if(mask) { //if material point is 'active'/'enabled'
@@ -134,37 +135,112 @@ void MPMesh::assemblyVtx1() {
       for(int i=0; i<nVtxE; i++){
         int vID = elm2VtxConn(elm,i+1)-1; //vID = vertex id
         double w_vtx=weight(mp,i);
-        
+        double mScale=1;
+        if(scaling)
+          mScale=sqrt(dual_triangle_area(vID,0))/radius;
         double CoordDiffs[vec4d_nEntries] = {1, (vtxCoords(vID,0) - mpPositions(mp,0))/radius, 
                                                 (vtxCoords(vID,1) - mpPositions(mp,1))/radius, 
                                                 (vtxCoords(vID,2) - mpPositions(mp,2))/radius};		 
-        for (int k=0; k<vec4d_nEntries; k++)
-          for (int l=0; l<vec4d_nEntries; l++)
+        //All entries except first row and column
+        for (int k=1; k<vec4d_nEntries; k++)
+          for (int l=1; l<vec4d_nEntries; l++)
             Kokkos::atomic_add(&VtxMatrices(vID,k,l), CoordDiffs[k] * CoordDiffs[l] * w_vtx);
+	//First entry
+        Kokkos::atomic_add(&VtxMatrices(vID,0,0), CoordDiffs[0] * CoordDiffs[0] * w_vtx*mScale*mScale); 
+	//First row and column except the first entry
+	for (int k=1; k<vec4d_nEntries; k++){
+          Kokkos::atomic_add(&VtxMatrices(vID,0,k), CoordDiffs[0] * CoordDiffs[k] * w_vtx*mScale);
+	  Kokkos::atomic_add(&VtxMatrices(vID,k,0), CoordDiffs[k] * CoordDiffs[0] * w_vtx*mScale);
+	}
       }
     }
   };
   p_MPs->parallel_for(assemble, "assembly");
   
-  //Solve Ax=b for each vertex  
+  //Assemble matrix for each vertex and apply regularization
   Kokkos::View<double*[vec4d_nEntries]> VtxCoeffs("VtxCoeffs", p_mesh->getNumVertices());
+
   Kokkos::parallel_for("solving Ax=b", numVtx, KOKKOS_LAMBDA(const int vtx){
     Vec4d v0 = {VtxMatrices(vtx,0,0), VtxMatrices(vtx,0,1), VtxMatrices(vtx,0,2), VtxMatrices(vtx,0,3)};
     Vec4d v1 = {VtxMatrices(vtx,1,0), VtxMatrices(vtx,1,1), VtxMatrices(vtx,1,2), VtxMatrices(vtx,1,3)};
     Vec4d v2 = {VtxMatrices(vtx,2,0), VtxMatrices(vtx,2,1), VtxMatrices(vtx,2,2), VtxMatrices(vtx,2,3)};
     Vec4d v3 = {VtxMatrices(vtx,3,0), VtxMatrices(vtx,3,1), VtxMatrices(vtx,3,2), VtxMatrices(vtx,3,3)};
-     
-    Matrix4d A = {v0,v1,v2,v3}; 
-    double A_trace = A.trace();
+    //Define the matrices
+    Matrix4d A = {v0,v1,v2,v3};
     Matrix4d A_regularized = {v0, v1, v2, v3};
-    A_regularized.addToDiag(A_trace*1e-8);
-   
+    //Regularization
+    switch(reg_method){
+      case 0:{
+        break;
+      }   
+      case 1:{
+        double A_trace = A.trace();
+        A_regularized.addToDiag(A_trace*1e-8);
+        break;
+      }
+      case 2:{
+        double regParam=sqrt(EPSILON)*(VtxMatrices(vtx,0,0)+VtxMatrices(vtx,1,1)+
+			               VtxMatrices(vtx,2,2)+VtxMatrices(vtx,3,3));
+        A_regularized.addToDiag(regParam);
+        break;
+      }
+      default:{
+        printf("Invalid regularization method \n");
+        break;	
+      }
+    }
+    //Solve Ax=b 
     double coeff[vec4d_nEntries]={0.0, 0.0, 0.0, 0.0};
     CholeskySolve4d_UnitRHS(A_regularized, coeff);
+    // Undo scaling
+    double mScale=1;
+    if(scaling)
+      mScale=sqrt(dual_triangle_area(vtx,0))/radius;
+        
+    coeff[0]=coeff[0]*mScale*mScale;
+    coeff[1]=coeff[1]*mScale;
+    coeff[2]=coeff[2]*mScale;
+    coeff[3]=coeff[3]*mScale;
     for (int i=0; i<vec4d_nEntries; i++) 
       VtxCoeffs(vtx,i)=coeff[i];
   });
+  this->precomputedVtxCoeffs = VtxCoeffs;
+}
+
+template <MeshFieldIndex meshFieldIndex>
+void MPMesh::assemblyVtx1() {
+   
+  //If no reconstruction till now calculate the coeffs
+  if (!isPreComputed) {
+    computeMatricesAndSolve();
+    isPreComputed=true;
+  }
   
+  auto VtxCoeffs=this->precomputedVtxCoeffs;
+  //Mesh Information
+  auto elm2VtxConn = p_mesh->getElm2VtxConn();  
+  int numVtx = p_mesh->getNumVertices();
+  auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
+
+  //Mesh Field
+  constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice<meshFieldIndex>;
+  const int numEntries = mpSliceToNumEntries<mpfIndex>();
+  p_mesh->fillMeshField<meshFieldIndex>(numVtx, numEntries, 0.0);
+  auto meshField = p_mesh->getMeshField<meshFieldIndex>();
+
+  //Material Points
+  auto mpData = p_MPs->getData<mpfIndex>();
+  auto weight = p_MPs->getData<MPF_Basis_Vals>();
+  auto mpPositions = p_MPs->getData<MPF_Cur_Pos_XYZ>();
+
+  //Earth Radius
+  double radius = 1.0;
+  if(p_mesh->getGeomType() == geom_spherical_surf)
+    radius=p_mesh->getSphereRadius();
+
+  //Reconstructed values
+  Kokkos::View<double**> reconVals("meshField", p_mesh->getNumVertices(), numEntries);
+
   //Reconstruct
   auto reconstruct = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
     if(mask) { //if material point is 'active'/'enabled'
@@ -175,7 +251,7 @@ void MPMesh::assemblyVtx1() {
         double CoordDiffs[vec4d_nEntries] = {1, (vtxCoords(vID,0) - mpPositions(mp,0))/radius,
                                                 (vtxCoords(vID,1) - mpPositions(mp,1))/radius, 
                                                 (vtxCoords(vID,2) - mpPositions(mp,2))/radius};
-        
+
         auto factor = w_vtx*(VtxCoeffs(vID,0) + VtxCoeffs(vID,1)*CoordDiffs[1] + 
                                                 VtxCoeffs(vID,2)*CoordDiffs[2] + 
                                                 VtxCoeffs(vID,3)*CoordDiffs[3]);
@@ -188,7 +264,8 @@ void MPMesh::assemblyVtx1() {
     }
   };
   p_MPs->parallel_for(reconstruct, "reconstruct");
- 
+
+  //Assign as a field 
   Kokkos::parallel_for("assigning", numVtx, KOKKOS_LAMBDA(const int vtx){
     for(int k=0; k<numEntries; k++)
       meshField(vtx, k) = reconVals(vtx,k);
diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp
index d5d0650..67653b5 100644
--- a/src/pmpo_c.cpp
+++ b/src/pmpo_c.cpp
@@ -743,6 +743,42 @@ void polympo_getMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, double* x
   }
 }
 
+void polympo_setMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* areaTriangle){
+
+  //chech validity
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+
+  //check the size
+  PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); 
+
+  //copy the host array to the device
+  auto dualArea = p_mesh->getMeshField<polyMPO::MeshF_DualTriangleArea>();
+  auto h_dualArea = Kokkos::create_mirror_view(dualArea);
+  for(int i=0; i<nVertices; i++)
+    h_dualArea(i,0) = areaTriangle[i];
+  Kokkos::deep_copy(dualArea, h_dualArea);
+
+}
+
+void polympo_getMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, double* areaTriangle){
+  
+  //chech validity
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+
+  //check the size
+  PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); 
+
+  //copy the device to host 
+  auto dualArea = p_mesh->getMeshField<polyMPO::MeshF_DualTriangleArea>();
+  auto h_dualArea = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), dualArea);
+  for(int i=0; i<nVertices; i++)
+    areaTriangle[i] = h_dualArea(i,0);
+
+}
+
+
 void polympo_setMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* uVelIn, const double* vVelIn){
   //check mpMesh is valid
   checkMPMeshValid(p_mpmesh);
diff --git a/src/pmpo_c.h b/src/pmpo_c.h
index 6c38106..174f862 100644
--- a/src/pmpo_c.h
+++ b/src/pmpo_c.h
@@ -77,6 +77,11 @@ void polympo_setMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, c
 void polympo_getMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array);//vec2d
 void polympo_setMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, const double* xArray, const double* yArray, const double* zArray);
 void polympo_getMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, double* xArray, double* yArray, double* zArray);
+//Area Triangle
+void polympo_setMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* areaTriangle);
+void polympo_getMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, double* areaTriangle);
+
+
 
 // calculations
 void polympo_push_f(MPMesh_ptr p_mpmesh);
diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90
index f16f57a..cf013ca 100644
--- a/src/pmpo_fortran.f90
+++ b/src/pmpo_fortran.f90
@@ -513,6 +513,34 @@ subroutine polympo_getMeshElmCenter(mpMesh, nCells, xArray, yArray, zArray) &
     integer(c_int), value :: nCells
     type(c_ptr), value :: xArray, yArray, zArray
   end subroutine
+
+  !---------------------------------------------------------------------------
+  !> @brief set the polympo dual mesh triangle area
+  !> @param mpmesh(in/out) MPMesh object
+  !> @param nVertices(in) length of array in, use for assertion 
+  !> @param Array(in) 1D array of area of dual triangle elements
+  !---------------------------------------------------------------------------
+  subroutine polympo_setMeshDualTriangleArea(mpMesh, nVertices, areaTriangle) &
+             bind(C, NAME='polympo_setMeshDualTriangleArea_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nVertices
+    type(c_ptr), intent(in), value :: areaTriangle
+  end subroutine
+  !---------------------------------------------------------------------------
+  !> @brief get the polympo mesh dual mesh triangle area
+  !> @param mpmesh(in/out) MPMesh object
+  !> @param nVetices(in) length of array in, use for assertion
+  !> @param Array(in/out) 1D array of area of dual triangle elements
+  !---------------------------------------------------------------------------
+  subroutine polympo_getMeshDualTriangleArea(mpMesh, nVertices, areaTriangle) &
+             bind(C, NAME='polympo_getMeshDualTriangleArea_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nVertices
+    type(c_ptr), value :: areaTriangle
+  end subroutine
+
   !---------------------------------------------------------------------------
   !> @brief set the vertices velocity from a host array
   !> @param mpmesh(in/out) MPMesh object
diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp
index cdf9b28..9232ad1 100644
--- a/src/pmpo_mesh.cpp
+++ b/src/pmpo_mesh.cpp
@@ -41,6 +41,10 @@ namespace polyMPO{
         auto vtxRotLatLonIncrMapEntry = meshFields2TypeAndString.at(MeshF_RotLatLonIncr);
         PMT_ALWAYS_ASSERT(vtxRotLatLonIncrMapEntry.first == MeshFType_VtxBased);
         vtxRotLatLonIncr_ = MeshFView<MeshF_RotLatLonIncr>(vtxRotLatLonIncrMapEntry.second,numVtxs_);
+
+        auto dualTriangleAreaEntry = meshFields2TypeAndString.at(MeshF_DualTriangleArea);
+        PMT_ALWAYS_ASSERT(dualTriangleAreaEntry.first == MeshFType_VtxBased);
+        dualTriangleArea_ = MeshFView<MeshF_DualTriangleArea>(dualTriangleAreaEntry.second,numVtxs_);
     }
     
     void Mesh::setMeshElmBasedFieldSize(){
diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp
index 3b91dc3..2058e13 100644
--- a/src/pmpo_mesh.hpp
+++ b/src/pmpo_mesh.hpp
@@ -20,6 +20,7 @@ enum MeshFieldIndex{
     MeshF_VtxCoords,
     MeshF_VtxRotLat,
     MeshF_ElmCenterXYZ,
+    MeshF_DualTriangleArea,
     MeshF_Vel,
     MeshF_VtxMass,
     MeshF_ElmMass,
@@ -38,6 +39,7 @@ template <MeshFieldIndex> struct meshFieldToType;
 template <> struct meshFieldToType < MeshF_VtxCoords         > { using type = Kokkos::View<vec3d_t*>; };
 template <> struct meshFieldToType < MeshF_VtxRotLat         > { using type = DoubleView; };
 template <> struct meshFieldToType < MeshF_ElmCenterXYZ      > { using type = Kokkos::View<vec3d_t*>; };
+template <> struct meshFieldToType < MeshF_DualTriangleArea  > { using type = Kokkos::View<doubleSclr_t*>; };
 template <> struct meshFieldToType < MeshF_Vel               > { using type = Kokkos::View<vec2d_t*>; };
 template <> struct meshFieldToType < MeshF_VtxMass           > { using type = Kokkos::View<doubleSclr_t*>; };
 template <> struct meshFieldToType < MeshF_ElmMass           > { using type = Kokkos::View<doubleSclr_t*>; };
@@ -55,6 +57,7 @@ const std::map<MeshFieldIndex, std::pair<MeshFieldType,
                {MeshF_VtxCoords,        {MeshFType_VtxBased,"MeshField_VerticesCoords"}},
                {MeshF_VtxRotLat,        {MeshFType_VtxBased,"MeshField_VerticesLatitude"}},
                {MeshF_ElmCenterXYZ,     {MeshFType_ElmBased,"MeshField_ElementCenterXYZ"}},
+               {MeshF_DualTriangleArea, {MeshFType_VtxBased,"MeshField_DualTriangleArea"}},
                {MeshF_Vel,              {MeshFType_VtxBased,"MeshField_Velocity"}},
                {MeshF_VtxMass,          {MeshFType_VtxBased,"MeshField_VerticesMass"}},
                {MeshF_ElmMass,          {MeshFType_ElmBased,"MeshField_ElementsMass"}},
@@ -90,6 +93,7 @@ class Mesh {
     MeshFView<MeshF_VtxCoords> vtxCoords_;
     MeshFView<MeshF_VtxRotLat> vtxRotLat_;
     MeshFView<MeshF_ElmCenterXYZ> elmCenterXYZ_;
+    MeshFView<MeshF_DualTriangleArea> dualTriangleArea_;
     MeshFView<MeshF_Vel> vtxVel_;
     MeshFView<MeshF_VtxMass> vtxMass_;
     MeshFView<MeshF_ElmMass> elmMass_;
@@ -181,6 +185,9 @@ auto Mesh::getMeshField(){
     else if constexpr (index==MeshF_ElmCenterXYZ){
         return elmCenterXYZ_;
     }
+    else if constexpr (index==MeshF_DualTriangleArea){
+        return dualTriangleArea_;
+    }
     else if constexpr (index==MeshF_Vel){
         return vtxVel_;
     }
diff --git a/src/pmpo_utils.hpp b/src/pmpo_utils.hpp
index b9d5f7d..7a2fdcd 100644
--- a/src/pmpo_utils.hpp
+++ b/src/pmpo_utils.hpp
@@ -274,6 +274,18 @@ class Matrix4d {
         data_[i][i]=data_[i][i]+eps;
       }
     }
+
+    KOKKOS_INLINE_FUNCTION
+    void scaleFirstRowAndColumn(double factor) {
+      // Scale the first row
+      for (int j = 0; j < 4; j++) {
+        data_[0][j] *= factor;
+      }
+      // Scale the first column
+      for (int i = 0; i < 4; i++) {
+        data_[i][0] *= factor;
+      }
+    }
 				        
 };
 
@@ -358,27 +370,24 @@ void QRDecomp4d(Matrix4d& A, Matrix4d& Q, Matrix4d& R){
       Vec4d u_column;
       for (int l=0; l<m_size; l++){
          A_column[l] = A(l,i);
-	 u_column[l] = A(l,i);
+         u_column[l] = A(l,i);
       }
 
       for (int j=0; j<i; j++){
          Vec4d q_column;
-	 for (int l=0; l<m_size; l++)
-	   q_column[l] = Q(l,j);
-	 double r=q_column.dot(A_column);
-	 R(j,i) = r;
+         for (int l=0; l<m_size; l++)
+           q_column[l] = Q(l,j);
+         double r=q_column.dot(A_column);
+         R(j,i) = r;
          for (int k=0; k<m_size; k++) u_column[k] -= r*Q(k,j);
       }
       double uNorm = std::sqrt(u_column.dot(u_column));
       R(i,i)=uNorm;
       for (int k=0; k<m_size; k++)
-          Q(k,i)= u_column[k]/uNorm; 
-
+          Q(k,i)= u_column[k]/uNorm;
   }
-
 }
 
-
 KOKKOS_INLINE_FUNCTION
 void initArray(Vec3d* arr, int n, Vec3d fill){
     for(int i=0; i<n; i++){
diff --git a/test/readMPAS.f90 b/test/readMPAS.f90
index ac7549d..c5a4116 100644
--- a/test/readMPAS.f90
+++ b/test/readMPAS.f90
@@ -27,7 +27,7 @@ subroutine loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, &
                         xVertex, yVertex, zVertex, &
                         latVertex, &
                         xCell, yCell, zCell, &
-                        verticesOnCell, cellsOnCell)
+                        verticesOnCell, cellsOnCell, areaTriangle)
     use :: netcdf
     use :: iso_c_binding
     implicit none
@@ -42,6 +42,8 @@ subroutine loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, &
     real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex
     real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell
     integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell
+    real(kind=MPAS_RKIND), dimension(:), pointer, optional :: areaTriangle
+
 
     call polympo_checkPrecisionForRealKind(MPAS_RKIND)
     !check on maxEdges and vertexDegree
@@ -75,6 +77,11 @@ subroutine loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, &
     
     !set mesh element center
     call polympo_setMeshElmCenter(mpMesh,nCells,c_loc(xCell),c_loc(yCell),c_loc(zCell))
+
+    if (present(areaTriangle)) then
+      call polympo_setMeshDualTriangleArea(mpMesh, nVertices, c_loc(areaTriangle)) 
+    end if
+
 end subroutine
 
 subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
@@ -83,7 +90,8 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
                         xVertex, yVertex, zVertex, &
                         latVertex, lonVertex, &
                         xCell, yCell, zCell, &
-                        verticesOnCell, cellsOnCell)
+                        verticesOnCell, cellsOnCell, &
+                        areaTriangle)
     use :: netcdf
     use :: iso_c_binding
     implicit none
@@ -99,12 +107,14 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
     real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex
     real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell
     integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell
+    
+    real(kind=MPAS_RKIND), dimension(:), pointer, optional :: areaTriangle(:)
 
     integer :: ncid, status, nCellsID, nVerticesID, maxEdgesID, vertexDegreeID, &
                nEdgesOnCellID, xVertexID, yVertexID, zVertexID, &
                latVertexID, lonVertexID, &
                xCellID, yCellID, zCellID, &
-               verticesOnCellID, cellsOnCellID
+               verticesOnCellID, cellsOnCellID, areaTriangleID
     
     status = nf90_open(path=trim(filename), mode=nf90_nowrite, ncid=ncid)
     if (status /= nf90_noerr) then
@@ -181,6 +191,10 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
     allocate(verticesOnCell(maxEdges, nCells))
     allocate(cellsOnCell(maxEdges, nCells))
    
+    if (present(areaTriangle)) then
+      allocate(areaTriangle(nVertices))
+    end if
+   
     status = nf90_get_att(ncid, nf90_global, "on_a_sphere", onSphere)
     if (status /= nf90_noerr) then
         write(0, *) "readMPASMeshFromNCFile: Error on get attribute 'on_a_sphere'"
@@ -275,6 +289,16 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
         stop
     end if
 
+    if (present(areaTriangle)) then
+      status = nf90_inq_varid(ncid, 'areaTriangle', areaTriangleID)
+      if (status /= nf90_noerr) then
+        write(0, *) "readMPASMeshFromNCFile: Error on inquire varid of 'areaTriangleID'"
+        write(0, *) trim(nf90_strerror(status))
+        stop
+      end if
+    end if
+
+    !Reading the variables
     status = nf90_get_var(ncid, xVertexID, xVertex)
     if (status /= nf90_noerr) then
         write(0, *) "readMPASMeshFromNCFile: Error on get var of 'xVertex'"
@@ -352,6 +376,18 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
         stop
     end if
 
+    if (present(areaTriangle)) then
+      print *, "Reading area of dual triangles"
+      status = nf90_get_var(ncid, areaTriangleID, areaTriangle)
+      if (status /= nf90_noerr) then
+        write(0, *) "readMPASMeshFromNCFile: Error on get var of 'areaTriangle'"
+        write(0, *) trim(nf90_strerror(status))
+        stop
+      end if
+    else
+      print *, "Not reading area of dual cells"
+    end if
+
     status = nf90_close(ncid)   
 end subroutine readMPASMeshFromNCFile
 subroutine setWithMPASMeshByFortran(mpMesh, fileName, n) bind(C, name="setWithMPASMeshByFortran")
diff --git a/test/testFortranMPReconstruction.f90 b/test/testFortranMPReconstruction.f90
index e88720e..03a4165 100644
--- a/test/testFortranMPReconstruction.f90
+++ b/test/testFortranMPReconstruction.f90
@@ -19,7 +19,7 @@ program main
   real(kind=MPAS_RKIND) :: pi = 4.0_MPAS_RKIND*atan(1.0_MPAS_RKIND)
   real(kind=MPAS_RKIND) :: TEST_VAL = 1.1_MPAS_RKIND
   real(kind=MPAS_RKIND) :: TOLERANCE = 0.0001_MPAS_RKIND
-  real(kind=MPAS_RKIND) :: TOLERANCE1 = 0.0000000001_MPAS_RKIND
+  real(kind=MPAS_RKIND) :: TOLERANCE1 = 0.000000001_MPAS_RKIND
   character (len=2048) :: filename
   real(kind=MPAS_RKIND), dimension(:,:), pointer :: dispIncr
   character (len=64) :: onSphere
@@ -28,12 +28,15 @@ program main
   real(kind=MPAS_RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex
   real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex
   real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell
+  real(kind=MPAS_RKIND), dimension(:), pointer :: areaTriangle
+  
   integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell
   integer :: numMPs, vID
   integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive
   real(kind=MPAS_RKIND), dimension(:,:), pointer :: mpPosition, mpLatLon
   real(kind=MPAS_RKIND), dimension(:,:), pointer :: mpMass, mpVel
   real(kind=MPAS_RKIND), dimension(:), pointer :: meshVtxMass, meshElmMass, meshVtxMass1
+  real(kind=MPAS_RKIND), dimension(:), pointer :: meshVtxVelu, meshVtxVelv
   logical :: inBound
   integer, parameter :: MP_ACTIVE = 1
   integer, parameter :: MP_INACTIVE = 0
@@ -58,7 +61,7 @@ program main
                         xVertex, yVertex, zVertex, &
                         latVertex, lonVertex, &
                         xCell, yCell, zCell, &
-                        verticesOnCell, cellsOnCell)
+                        verticesOnCell, cellsOnCell, areaTriangle)
   if (onSphere .ne. 'YES') then
     write (*,*) "The mesh is not spherical!"
     call exit(1)
@@ -74,7 +77,7 @@ program main
                         xVertex, yVertex, zVertex, &
                         latVertex, &
                         xCell, yCell, zCell, &
-                        verticesOnCell, cellsOnCell)
+                        verticesOnCell, cellsOnCell, areaTriangle)
  
   nCompsDisp = 2
   allocate(dispIncr(nCompsDisp,nVertices))
@@ -90,6 +93,8 @@ program main
   allocate(meshVtxMass(nVertices))
   allocate(meshVtxMass1(nVertices))
   allocate(meshElmMass(nCells))
+  allocate(meshVtxVelu(nVertices))
+  allocate(meshVtxVelv(nVertices))
 
   isMPActive = MP_ACTIVE !all active MPs and some changed below
   mpsPerElm = 1 !all elements have 1 MP and some changed below
@@ -120,14 +125,20 @@ program main
 
   do i = 1, nVertices
     call assert(meshVtxMass(i) < TEST_VAL+TOLERANCE .and. meshVtxMass(i) > TEST_VAL-TOLERANCE, "Error: wrong vtx mass")
+    !write(*, *) 'The value of LRV is:', meshVtxMass(i)
   end do
   
   !Test vtx order 1 reconstruction
   call polympo_setReconstructionOfMass(mpMesh,1,polympo_getMeshFVtxType())
+  call polympo_setReconstructionOfVel(mpMesh, 1, polympo_getMeshFVtxType())
   call polympo_applyReconstruction(mpMesh)
   call polympo_getMeshVtxMass(mpMesh,nVertices,c_loc(meshVtxMass1))
+  call polympo_getMeshVtxVel(mpMesh, nVertices, c_loc(meshVtxVelu), c_loc(meshVtxVelv))
   do i = 1, nVertices
-    call assert(meshVtxMass1(i) < TEST_VAL+TOLERANCE1 .and. meshVtxMass1(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx mass LINEAR")
+    call assert(meshVtxMass1(i) < TEST_VAL+TOLERANCE1 .and. meshVtxMass1(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx mass order 1")
+    call assert(meshVtxVelu(i) < TEST_VAL+TOLERANCE1 .and. meshVtxVelu(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx velU order 1")
+    call assert(meshVtxVelv(i) < TEST_VAL+TOLERANCE1 .and. meshVtxVelv(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx velV order 1")
+    write(*, *) 'The value of LRV is:', meshVtxVelu(i)-TEST_VAL, meshVtxVelv(i)-TEST_VAL
   end do
 
 
@@ -143,7 +154,7 @@ program main
 
     do i = 1, numMPs
       vID = verticesOnCell(1,mp2Elm(i))
-      call assert(meshVtxMass(vID) < TEST_VAL+TOLERANCE .and. meshVtxMass(vID) > TEST_VAL-TOLERANCE, "Error: wrong vtx mass")
+      call assert(meshVtxMass(vID) < TEST_VAL+TOLERANCE1 .and. meshVtxMass(vID) > TEST_VAL-TOLERANCE1, "Error: wrong vtx mass push")
 
       call assert(meshElmMass(mp2Elm(i)) < TEST_VAL+TOLERANCE &
             .and. meshElmMass(mp2Elm(i)) > TEST_VAL-TOLERANCE, "Error: wrong elm mass")
@@ -175,6 +186,8 @@ program main
   deallocate(meshVtxMass)
   deallocate(meshVtxMass1)
   deallocate(meshElmMass)
+  deallocate(meshVtxVelu)
+  deallocate(meshVtxVelv)
 
   stop