diff --git a/.github/workflows/CPU_inferencce_validation.yml b/.github/workflows/CPU_inferencce_validation.yml
index de3113b8c..22ed1873e 100644
--- a/.github/workflows/CPU_inferencce_validation.yml
+++ b/.github/workflows/CPU_inferencce_validation.yml
@@ -12,6 +12,7 @@ jobs:
     if: github.repository_owner == 'deepmodeling'     
     container: 
       image: ubuntu:20.04
+      options: --shm-size=5g  # shared memory size = 5GB
     steps:
      - name: Checkout
        uses: actions/checkout@v3
@@ -46,11 +47,10 @@ jobs:
         OMPI_MCA_btl_vader_single_copy_mechanism: none
        run: 
         /bin/bash -c " echo $PWD && ls $PWD
-        && wget https://github.com/JX278/DeepFlame-flareFGM-table/raw/main/flareFGM_Table_Download.zip
-        && unzip flareFGM_Table_Download.zip
+        && wget --content-disposition https://aisquare.oss-us-east-1.aliyuncs.com/data/datasets/14b50df5-dbe9-4f1c-bf58-032b8bc40a20
+        && unzip flare_CH4_SandiaD_4D.zip
         && ls $PWD
-        && ls flareFGM_Table_Download
-        && cp -r flareFGM_Table_Download/SandiaD/flare.tbl examples/dfLowMachFoam/2DSandiaD_flareFGM/
+        && cp -r flare_CH4_SandiaD_4D.tbl examples/dfLowMachFoam/2DSandiaD_flareFGM/
         && git clone https://github.com/intelligent-algorithm-team/intelligent-combustion.git 
         && cp -r intelligent-combustion/DeePCK/Model/HE04_Hydrogen_ESH2_GMS_sub_20221101/ mechanisms/ && source ~/miniconda3/etc/profile.d/conda.sh && conda activate libcantera && source /opt/openfoam7/etc/bashrc 
         && . configure.sh --use_pytorch&& source ./bashrc && . install.sh 
diff --git a/applications/solvers/dfHighSpeedFoam/createFields.H b/applications/solvers/dfHighSpeedFoam/createFields.H
index 852122f45..4d5e8b4c4 100644
--- a/applications/solvers/dfHighSpeedFoam/createFields.H
+++ b/applications/solvers/dfHighSpeedFoam/createFields.H
@@ -186,8 +186,15 @@ volScalarField rhoE
     rho*(ea + 0.5*magSqr(U))
 );
 
-chemistry->correctThermo();
-Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
+const word combModelName(mesh.objectRegistry::lookupObject<IOdictionary>("combustionProperties").lookup("combustionModel"));
+Info << "Combustion Model Name is confirmed as "<< combModelName << endl;
+
+
+if (combModelName != "flareFGM")
+{
+    chemistry->correctThermo();
+    Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
+}
 
 forAll(Y, i)
 {
diff --git a/applications/solvers/dfLowMachFoam/createFields.H b/applications/solvers/dfLowMachFoam/createFields.H
index 9e750c334..15413d462 100644
--- a/applications/solvers/dfLowMachFoam/createFields.H
+++ b/applications/solvers/dfLowMachFoam/createFields.H
@@ -93,8 +93,11 @@ chemistry->setEnergyName("ha");
 chemistry->updateEnergy();
 
 
-chemistry->correctThermo();
-Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
+if (combModelName != "flareFGM")
+{
+    chemistry->correctThermo();
+    Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
+}
 
 //for dpdt
 
diff --git a/applications/solvers/dfSprayFoam/createFields.H b/applications/solvers/dfSprayFoam/createFields.H
index 3c33fce70..7506f74b8 100644
--- a/applications/solvers/dfSprayFoam/createFields.H
+++ b/applications/solvers/dfSprayFoam/createFields.H
@@ -78,11 +78,30 @@ autoPtr<compressible::turbulenceModel> turbulence
 
 const word turbName(mesh.objectRegistry::lookupObject<IOdictionary>("turbulenceProperties").lookup("simulationType"));
 
+Info<< "Creating field dpdt\n" << endl;
+volScalarField dpdt
+(
+    IOobject
+    (
+        "dpdt",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar(p.dimensions()/dimTime, 0)
+);
+
 Info<< "Creating reaction model\n" << endl;
 autoPtr<CombustionModel<basicThermo>> combustion
 (
     CombustionModel<basicThermo>::New(thermo, turbulence())
 );
+
+const word combModelName(mesh.objectRegistry::lookupObject<IOdictionary>
+("combustionProperties").lookup("combustionModel"));
+Info << "Combustion Model Name is confirmed as "<< combModelName << endl;
+
+
 Info<< "end Creating reaction model\n" << endl;
 dfChemistryModel<basicThermo>* chemistry = combustion->chemistry();
 PtrList<volScalarField>& Y = chemistry->Y();
@@ -91,23 +110,12 @@ const label inertIndex(chemistry->species()[inertSpecie]);
 chemistry->setEnergyName("ha");
 chemistry->updateEnergy();
 
+if (combModelName != "flareFGM")
+{
+    chemistry->correctThermo();
+    Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
+}
 
-chemistry->correctThermo();
-Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
-
-
-Info<< "Creating field dpdt\n" << endl;
-volScalarField dpdt
-(
-    IOobject
-    (
-        "dpdt",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar(p.dimensions()/dimTime, 0)
-);
 
 Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
diff --git a/applications/solvers/dfSprayFoam/dfSprayFoam.C b/applications/solvers/dfSprayFoam/dfSprayFoam.C
index ff044db6d..363508475 100644
--- a/applications/solvers/dfSprayFoam/dfSprayFoam.C
+++ b/applications/solvers/dfSprayFoam/dfSprayFoam.C
@@ -153,9 +153,17 @@ int main(int argc, char *argv[])
         while (pimple.loop())
         {
             #include "UEqn.H"
-            #include "YEqn.H"
-            #include "EEqn.H"
-            chemistry->correctThermo();
+
+            if(combModelName!="ESF" && combModelName!="flareFGM" )
+            {
+                #include "YEqn.H"
+                #include "EEqn.H"
+                chemistry->correctThermo();
+            }
+            else
+            {
+                combustion->correct();
+            }
             Info<< "T gas min/max   " << min(T).value() << ", "
                 << max(T).value() << endl;
 
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/G.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/G.gz
deleted file mode 100644
index a075abd5e..000000000
Binary files a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/G.gz and /dev/null differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/He b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/Ha
similarity index 96%
rename from examples/dfLowMachFoam/2DSandiaD_flareFGM/0/He
rename to examples/dfLowMachFoam/2DSandiaD_flareFGM/0/Ha
index 88f68ca3a..c2db898b2 100644
--- a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/He
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/Ha
@@ -34,7 +34,7 @@ boundaryField
     inletPilot
     {
         type            fixedValue;
-        value           uniform -170863.69815950916;
+        value           uniform -294024.3725;
     }
 
     inletAir
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/T b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/T
index f90678557..b529e59df 100644
--- a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/T
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/T
@@ -26,39 +26,36 @@ boundaryField
         type            fixedValue;
         value           uniform 294;
     }
-
     wallOutside
     {
         type            zeroGradient;
     }
-
     wallTube
     {
         type            zeroGradient;
     }
-
     inletPilot
     {
         type            fixedValue;
         value           uniform 1880;
     }
-
     inletAir
     {
         type            fixedValue;
         value           uniform 291;
     }
-
     outlet
     {
         type            zeroGradient;
     }
-
+    axis
+    {
+        type            empty;
+    }
     frontAndBack_pos
     {
         type            wedge;
     }
-
     frontAndBack_neg
     {
         type            wedge;
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/U b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/U
new file mode 100644
index 000000000..0c21ab649
--- /dev/null
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/U
@@ -0,0 +1,70 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  7
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0.9);
+
+boundaryField
+{
+    wallTube
+    {
+        type            noSlip;
+    }
+
+    outlet
+    {
+        type            pressureInletOutletVelocity;
+        value           $internalField;
+    }
+
+    inletPilot
+    {
+        type            fixedValue;
+        value           uniform (0 0 11.4);
+    }
+
+    inletAir
+    {
+        type            fixedValue;
+        value           uniform (0 0 0.9);
+    }
+
+    wallOutside
+    {
+        type            zeroGradient;
+    }
+
+    inletCH4
+    {
+        type            fixedValue;
+        value           uniform (0 0 49.6);
+    }
+
+    frontAndBack_pos
+    {
+        type            wedge;
+    }
+
+    frontAndBack_neg
+    {
+        type            wedge;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/U.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/U.gz
deleted file mode 100644
index eb990d607..000000000
Binary files a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/U.gz and /dev/null differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/Ydefault b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/Ydefault
new file mode 100644
index 000000000..6d85f7b3d
--- /dev/null
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/Ydefault
@@ -0,0 +1,69 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  7
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      Ydefault;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    wallTube
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    inletPilot
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    inletAir
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    wallOutside
+    {
+        type            zeroGradient;
+    }
+
+    inletCH4
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    frontAndBack_pos
+    {
+        type            wedge;
+    }
+
+    frontAndBack_neg
+    {
+        type            wedge;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/Ydefault.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/Ydefault.gz
deleted file mode 100644
index 1982e123a..000000000
Binary files a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/Ydefault.gz and /dev/null differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/alphat b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/alphat
new file mode 100644
index 000000000..d2bf2613e
--- /dev/null
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/alphat
@@ -0,0 +1,74 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  7
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    wallTube
+    {
+        type            compressible::alphatWallFunction;
+        Prt             0.85;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    inletPilot
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    inletAir
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    wallOutside
+    {
+        type            compressible::alphatWallFunction;
+        Prt             0.85;
+        value           $internalField;
+    }
+
+    inletCH4
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    frontAndBack_pos
+    {
+        type            wedge;
+    }
+
+    frontAndBack_neg
+    {
+        type            wedge;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/alphat.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/alphat.gz
deleted file mode 100644
index 363e369b6..000000000
Binary files a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/alphat.gz and /dev/null differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/epsilon b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/epsilon
new file mode 100644
index 000000000..426806689
--- /dev/null
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/epsilon
@@ -0,0 +1,84 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  7
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 30000;
+
+boundaryField
+{
+    wallTube
+    {
+        type            epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 30000;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    inletPilot
+    {
+        type            turbulentMixingLengthDissipationRateInlet;
+        mixingLength    0.000735;
+        phi             phi;
+        k               k;
+        value           uniform 1;
+    }
+
+    inletAir
+    {
+        type            turbulentMixingLengthDissipationRateInlet;
+        mixingLength    0.019677;
+        value           uniform 1;
+    }
+
+    wallOutside
+    {
+        type            epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 30000;
+    }
+
+    inletCH4
+    {
+        type            turbulentMixingLengthDissipationRateInlet;
+        mixingLength    0.000504;
+        phi             phi;
+        k               k;
+        value           uniform 1;
+    }
+
+    frontAndBack_pos
+    {
+        type            wedge;
+    }
+
+    frontAndBack_neg
+    {
+        type            wedge;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/epsilon.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/epsilon.gz
deleted file mode 100644
index 8f25949fd..000000000
Binary files a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/epsilon.gz and /dev/null differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/k b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/k
new file mode 100644
index 000000000..25602f936
--- /dev/null
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/k
@@ -0,0 +1,74 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  7
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 30;
+
+boundaryField
+{
+    wallTube
+    {
+        type            kqRWallFunction;
+        value           uniform 30;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    inletPilot
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.0628;
+        value           uniform 1;
+    }
+
+    inletAir
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.0471;
+        value           uniform 1;
+    }
+
+    wallOutside
+    {
+        type            kqRWallFunction;
+        value           uniform 30;
+    }
+
+    inletCH4
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.0458;
+        value           uniform 1;
+    }
+
+    frontAndBack_pos
+    {
+        type            wedge;
+    }
+
+    frontAndBack_neg
+    {
+        type            wedge;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/k.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/k.gz
deleted file mode 100644
index c07786222..000000000
Binary files a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/k.gz and /dev/null differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/nut b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/nut
new file mode 100644
index 000000000..ee38ce327
--- /dev/null
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/nut
@@ -0,0 +1,78 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  7
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    wallTube
+    {
+        type            nutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    inletPilot
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    inletAir
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    wallOutside
+    {
+        type            nutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           $internalField;
+    }
+
+    inletCH4
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    frontAndBack_pos
+    {
+        type            wedge;
+    }
+
+    frontAndBack_neg
+    {
+        type            wedge;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/nut.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/nut.gz
deleted file mode 100644
index efdec743e..000000000
Binary files a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/nut.gz and /dev/null differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/p b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/p
new file mode 100644
index 000000000..0a4e7409d
--- /dev/null
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/p
@@ -0,0 +1,67 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  7
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    wallTube
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            totalPressure;
+        p0              $internalField;
+    }
+
+    inletPilot
+    {
+        type            zeroGradient;
+    }
+
+    inletAir
+    {
+        type            zeroGradient;
+    }
+
+    wallOutside
+    {
+        type            zeroGradient;
+    }
+
+    inletCH4
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack_pos
+    {
+        type            wedge;
+    }
+
+    frontAndBack_neg
+    {
+        type            wedge;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/p.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/p.gz
deleted file mode 100644
index 75cef9487..000000000
Binary files a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/p.gz and /dev/null differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/rho b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/rho
new file mode 100644
index 000000000..0e729c8e1
--- /dev/null
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/0/rho
@@ -0,0 +1,72 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  4.1                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      rho;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -3 0 0 0 0 0];
+
+internalField   uniform 1.177;
+
+boundaryField
+{
+    wallTube
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    inletPilot
+    {
+        type            zeroGradient;
+    }
+
+    inletAir
+    {
+        type            zeroGradient;
+    }
+
+    wallOutside
+    {
+        type            zeroGradient;
+    }
+
+    inletCH4
+    {
+        type            zeroGradient;
+    }
+
+    axis
+    {
+        type            empty;
+    }
+    
+    frontAndBack_pos
+    {
+        type            wedge;
+    }
+    frontAndBack_neg
+    {
+        type            wedge;
+    }
+
+
+}
+
+
+// ************************************************************************* //
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/AR.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/AR.gz
new file mode 100644
index 000000000..532d6a70e
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/AR.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C.gz
new file mode 100644
index 000000000..bf9f6acbd
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H.gz
new file mode 100644
index 000000000..72cd9af70
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H2.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H2.gz
new file mode 100644
index 000000000..18b517ef7
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H2.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H3.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H3.gz
new file mode 100644
index 000000000..37d316347
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H3.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H4.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H4.gz
new file mode 100644
index 000000000..2ae1c1fba
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H4.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H5.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H5.gz
new file mode 100644
index 000000000..01b7fd997
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H5.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H6.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H6.gz
new file mode 100644
index 000000000..0e597e171
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C2H6.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C3H7.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C3H7.gz
new file mode 100644
index 000000000..b3ee39589
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C3H7.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C3H8.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C3H8.gz
new file mode 100644
index 000000000..ad3b29f63
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/C3H8.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH.gz
new file mode 100644
index 000000000..61e050d55
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2(S).gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2(S).gz
new file mode 100644
index 000000000..998b2845c
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2(S).gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2.gz
new file mode 100644
index 000000000..5e77499a3
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2CHO.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2CHO.gz
new file mode 100644
index 000000000..40dffa261
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2CHO.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2CO.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2CO.gz
new file mode 100644
index 000000000..16b67dac0
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2CO.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2O.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2O.gz
new file mode 100644
index 000000000..424279d2a
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2O.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2OH.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2OH.gz
new file mode 100644
index 000000000..7938aa781
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH2OH.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH3.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH3.gz
new file mode 100644
index 000000000..d5e4deba8
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH3.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH3CHO.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH3CHO.gz
new file mode 100644
index 000000000..df0001c3d
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH3CHO.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH3O.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH3O.gz
new file mode 100644
index 000000000..37609deb7
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH3O.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH3OH.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH3OH.gz
new file mode 100644
index 000000000..8d6853acb
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH3OH.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH4.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH4.gz
new file mode 100644
index 000000000..c801edf22
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CH4.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CO.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CO.gz
new file mode 100644
index 000000000..1aaea37e4
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CO.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CO2.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CO2.gz
new file mode 100644
index 000000000..e0defed33
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/CO2.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Cp.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Cp.gz
new file mode 100644
index 000000000..ca2836928
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Cp.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/H.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/H.gz
new file mode 100644
index 000000000..af0270781
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/H.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/H2.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/H2.gz
new file mode 100644
index 000000000..558ab2002
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/H2.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/H2O.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/H2O.gz
new file mode 100644
index 000000000..0629effca
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/H2O.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/H2O2.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/H2O2.gz
new file mode 100644
index 000000000..f8712791f
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/H2O2.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/HCCO.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/HCCO.gz
new file mode 100644
index 000000000..0e93c4fae
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/HCCO.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/HCCOH.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/HCCOH.gz
new file mode 100644
index 000000000..d012eacec
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/HCCOH.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/HCO.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/HCO.gz
new file mode 100644
index 000000000..8ef9cb997
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/HCO.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/HO2.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/HO2.gz
new file mode 100644
index 000000000..c1976fa91
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/HO2.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Ha.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Ha.gz
new file mode 100644
index 000000000..bf2af4f0d
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Ha.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Ha_s.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Ha_s.gz
new file mode 100644
index 000000000..2b3487729
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Ha_s.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Hf.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Hf.gz
new file mode 100644
index 000000000..6aa3343c8
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Hf.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/N2.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/N2.gz
new file mode 100644
index 000000000..e3f9ec44e
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/N2.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/O.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/O.gz
new file mode 100644
index 000000000..f9deec678
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/O.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/O2.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/O2.gz
new file mode 100644
index 000000000..779bf46e4
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/O2.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/OH.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/OH.gz
new file mode 100644
index 000000000..714e7c0fa
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/OH.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Qdot.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Qdot.gz
new file mode 100644
index 000000000..a72f3ffd9
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Qdot.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/T.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/T.gz
new file mode 100644
index 000000000..680eaea0c
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/T.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/U.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/U.gz
new file mode 100644
index 000000000..518ea1fcc
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/U.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Wt.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Wt.gz
new file mode 100644
index 000000000..c6305af24
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Wt.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Z.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Z.gz
new file mode 100644
index 000000000..342687ba4
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Z.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Zcvar.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Zcvar.gz
new file mode 100644
index 000000000..3dc021c5e
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Zcvar.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Zcvar_s.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Zcvar_s.gz
new file mode 100644
index 000000000..b7279fb39
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Zcvar_s.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Zvar.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Zvar.gz
new file mode 100644
index 000000000..088819196
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Zvar.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Zvar_s.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Zvar_s.gz
new file mode 100644
index 000000000..7cb7252c6
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/Zvar_s.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/alphat.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/alphat.gz
new file mode 100644
index 000000000..957f9a88a
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/alphat.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/betac.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/betac.gz
new file mode 100644
index 000000000..ff06e67a7
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/betac.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/c.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/c.gz
new file mode 100644
index 000000000..2384621f6
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/c.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/c_s.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/c_s.gz
new file mode 100644
index 000000000..d7131ac3b
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/c_s.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/cellCpuTimes.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/cellCpuTimes.gz
new file mode 100644
index 000000000..e600d0054
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/cellCpuTimes.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/chi_Z.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/chi_Z.gz
new file mode 100644
index 000000000..f5ea0403f
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/chi_Z.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/chi_Zc.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/chi_Zc.gz
new file mode 100644
index 000000000..5bf167ac9
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/chi_Zc.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/chi_c.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/chi_c.gz
new file mode 100644
index 000000000..15f63a7bf
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/chi_c.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/cvar.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/cvar.gz
new file mode 100644
index 000000000..3be08fba7
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/cvar.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/cvar_s.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/cvar_s.gz
new file mode 100644
index 000000000..fbdce24de
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/cvar_s.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/epsilon.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/epsilon.gz
new file mode 100644
index 000000000..34f4e1b9c
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/epsilon.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/k.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/k.gz
new file mode 100644
index 000000000..ee3392c37
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/k.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/nut.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/nut.gz
new file mode 100644
index 000000000..f072d3b63
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/nut.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/omega_c.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/omega_c.gz
new file mode 100644
index 000000000..b8451327a
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/omega_c.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/p.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/p.gz
new file mode 100644
index 000000000..ab23eff7d
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/p.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/phi.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/phi.gz
new file mode 100644
index 000000000..df142c4c3
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/phi.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/rho.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/rho.gz
new file mode 100644
index 000000000..4fae20606
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/rho.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/selectDNN.gz b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/selectDNN.gz
new file mode 100644
index 000000000..b05845914
Binary files /dev/null and b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/selectDNN.gz differ
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/combustionProperties.org b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/uniform/time
similarity index 78%
rename from examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/combustionProperties.org
rename to examples/dfLowMachFoam/2DSandiaD_flareFGM/1/uniform/time
index 13bf4db65..c34718405 100644
--- a/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/combustionProperties.org
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/1/uniform/time
@@ -10,21 +10,20 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
-    object      combustionProperties;
+    location    "1/uniform";
+    object      time;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-combustionModel  EDC;
+value           1;
 
-EDCCoeffs
-{
-    version v2005;
-}
+name            "1";
+
+index           17469;
+
+deltaT          5.7162129e-05;
+
+deltaT0         5.7162129e-05;
 
-PaSRCoeffs
-{
-   Cmix  0.03;
-}
 
 // ************************************************************************* //
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/Allrun b/examples/dfLowMachFoam/2DSandiaD_flareFGM/Allrun
index 39b0cddcc..fcb5cef98 100755
--- a/examples/dfLowMachFoam/2DSandiaD_flareFGM/Allrun
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/Allrun
@@ -12,7 +12,7 @@ runApplication blockMesh
 
 runApplication decomposePar -force 
 
-mpirun --oversubscribe -np 8 --allow-run-as-root dfLowMachFoam -parallel
+mpirun --oversubscribe -np 4 --allow-run-as-root dfLowMachFoam -parallel
 
 ./postProcess
  
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/CanteraTorchProperties b/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/CanteraTorchProperties
index 817b43040..d7ffb0ba7 100644
--- a/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/CanteraTorchProperties
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/CanteraTorchProperties
@@ -19,7 +19,7 @@ chemistry       off;
 
 CanteraMechanismFile "gri30.yaml";
 
-transportModel       "Mix";
+transportModel       "mixture-averaged";
 
 odeCoeffs
 {
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/combustionProperties b/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/combustionProperties
index d40666475..b166acd70 100644
--- a/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/combustionProperties
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/combustionProperties
@@ -66,7 +66,24 @@ PaSRCoeffs
 
 flareFGMCoeffs
 {
-  speciesName ("CO");
-  combustion on;
+  scaledPV         true;
+  combustion       true;
+  solveEnthalpy    true;
+  flameletT        false;
+  relaxation       false;
+  DpDt             false;
+  buffer           false;
+  bufferTime       0.0;
+  ignition         false;
+  /*
+  ignBeginTime     0.0;
+  ignDurationTime  0.0;
+  x0               0.0;
+  y0               0.0;
+  z0               0.0;
+  R0               0.0;
+  */
+  Sct              0.7;
+  tablePath    flare_CH4_SandiaD_4D.tbl;
 }
 // ************************************************************************* //
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/g b/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/g
deleted file mode 100644
index 113310120..000000000
--- a/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/g
+++ /dev/null
@@ -1,21 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Version:  7
-     \\/     M anipulation  |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       uniformDimensionedVectorField;
-    location    "constant";
-    object      g;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 1 -2 0 0 0 0];
-value           (0 0 -9.81);
-
-// ************************************************************************* //
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/thermophysicalProperties b/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/thermophysicalProperties
index 5b0c5908f..3940ea956 100644
--- a/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/thermophysicalProperties
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/thermophysicalProperties
@@ -15,5 +15,33 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+thermoType
+{
+    type            hePsiThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          eConst;
+    equationOfState perfectGas;
+    specie          specie;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    specie
+    {
+        molWeight   28.9;
+    }
+    thermodynamics
+    {
+        Cv          712;
+        Hf          0;
+    }
+    transport
+    {
+        mu          1.8e-05;
+        Pr          0.7;
+    }
+}
 
 // ************************************************************************* //
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/turbulenceProperties b/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/turbulenceProperties
index 6efd85e88..44a2ddc9b 100644
--- a/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/turbulenceProperties
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/constant/turbulenceProperties
@@ -19,10 +19,12 @@ simulationType   RAS;
 
 RAS
 {
-    RASModel        kEpsilon;
+    RASModel        RNGkEpsilon;
 
     turbulence      on;
     printCoeffs     on;
+
+    C1                1.52;
 }
 
 // ************************************************************************* //
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/postProcess b/examples/dfLowMachFoam/2DSandiaD_flareFGM/postProcess
index 33bdeb38b..68041df97 100755
--- a/examples/dfLowMachFoam/2DSandiaD_flareFGM/postProcess
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/postProcess
@@ -1,6 +1,6 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # Run from this directory
 
-reconstructPar -time 0.2
-postProcess -func sample -time 0.2
+reconstructPar -time 1.1
+postProcess -func sample -time 1.1
 
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/controlDict b/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/controlDict
index 555aaf4cb..a1583249b 100644
--- a/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/controlDict
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/controlDict
@@ -22,21 +22,21 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         0.2;
+endTime         1.1;
 
 deltaT          1e-05;
 
-writeControl    runTime;
+writeControl    adjustableRunTime;
 
-writeInterval   0.2;
+writeInterval   0.1;
 
 purgeWrite      0;
 
-writeFormat     binary;
+writeFormat     ascii;
 
-writePrecision  10;
+writePrecision  8;
 
-writeCompression off;
+writeCompression on;
 
 timeFormat      general;
 
@@ -46,7 +46,13 @@ runTimeModifiable true;
 
 adjustTimeStep  no;
 
-maxCo           0.1;
+adjustTimeStep  yes;
+
+maxCo           0.8;
+
+maxAlphaCo      1;
+
+maxDeltaT       1;
 
 
 // ************************************************************************* //
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/decomposeParDict b/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/decomposeParDict
index ad6a27c5f..634f839e2 100644
--- a/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/decomposeParDict
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/decomposeParDict
@@ -14,7 +14,7 @@ FoamFile
     object          decomposeParDict;
 }
 
-numberOfSubdomains 8;
+numberOfSubdomains 4;
 
 method          scotch;
 
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/fvSchemes b/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/fvSchemes
index 9fdf7a0e1..63840d6ce 100644
--- a/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/fvSchemes
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/fvSchemes
@@ -33,7 +33,7 @@ divSchemes
 //    div(phi,U)       Gauss linearUpwind grad(U);
     div(phi,Ulinear)    Gauss linearUpwind grad(U);
     div(phi,scalarUW)   Gauss limitedLinear 1;
-    div(phi,H)      Gauss limitedLinear 1;
+    div(phi,Ha)      Gauss limitedLinear 1;
     div(phi,Zvar)   Gauss limitedLinear 1;
     div(phi,Z)      Gauss limitedLinear01 1;
     div((muEff*dev2(T(grad(U))))) Gauss linear;
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/fvSolution b/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/fvSolution
index add1de749..20a540dc8 100644
--- a/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/fvSolution
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/fvSolution
@@ -68,7 +68,6 @@ solvers
         relTol          0.001;
     }
 
-/*
     "(k|epsilon)"
     {
         solver          PBiCG;
@@ -76,7 +75,6 @@ solvers
         tolerance       1e-09;
         relTol          0.001;
     }
-*/
 
     "(k|epsilon)Final"
     {
@@ -86,7 +84,7 @@ solvers
         // minIter         1;
     }
 
-    "(Z|H|Zvar|c|cvar|Zcvar).*"
+    "(Z|Ha|Zvar|c|cvar|Zcvar).*"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -103,15 +101,6 @@ PIMPLE
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
     momentumPredictor yes;
-
-    // outerCorrectorResidualControl
-    // {
-    //     "(U)"
-    //     {
-    //         tolerance       1e-04;
-    //         relTol          0;
-    //     }
-    // }
 }
 
 /*
@@ -124,7 +113,7 @@ relaxationFactors
     equations
     {
        U               0.8;
-//       phiFilt         0.5;
+       phiFilt         0.5;
     }
 }
 */
diff --git a/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/sample b/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/sample
index 3d5d75d9b..53e2752d6 100644
--- a/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/sample
+++ b/examples/dfLowMachFoam/2DSandiaD_flareFGM/system/sample
@@ -26,10 +26,10 @@ sets
     data
     {
         type    lineUniform;
-        axis    x;
-        start   (0    0 0.1);
-        end     (0.03 0 0.1);
-        nPoints 7;
+        axis    z;
+        start   (0   0   0);
+        end     (0   0   0.5);
+        nPoints 11;
     }
 );
 
diff --git a/src/dfChemistryModel/dfChemistryModel.H b/src/dfChemistryModel/dfChemistryModel.H
index 7705dc77f..c10035380 100644
--- a/src/dfChemistryModel/dfChemistryModel.H
+++ b/src/dfChemistryModel/dfChemistryModel.H
@@ -313,6 +313,9 @@ public:
         const volScalarField& rhoD(const label i) const {return rhoD_[i];}
 
         const volScalarField& hai(const label i) {return hai_[i];}
+        
+        const scalar & hci(const label i) {return hc_[i];}
+
         // update T, psi, mu, alpha, rhoD, hai (if needed)
         void correctThermo();
 
diff --git a/src/dfCombustionModels/FGM/FGMLiquidEvaporationBoil/FGMLiquidEvaporationBoil.C b/src/dfCombustionModels/FGM/FGMLiquidEvaporationBoil/FGMLiquidEvaporationBoil.C
new file mode 100644
index 000000000..bedb12220
--- /dev/null
+++ b/src/dfCombustionModels/FGM/FGMLiquidEvaporationBoil/FGMLiquidEvaporationBoil.C
@@ -0,0 +1,228 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "FGMLiquidEvaporationBoil.H"
+// #include "specie.H"
+#include "mathematicalConstants.H"
+
+
+using namespace Foam::constant::mathematical;
+
+// // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+// template<class CloudType>
+// Foam::tmp<Foam::scalarField> Foam::LiquidEvaporationBoil<CloudType>::calcXc
+// (
+//     const label celli
+// ) const
+// {
+//     scalarField Xc(this->owner().thermo().carrier().Y().size());
+
+//     forAll(Xc, i)
+//     {
+//         Xc[i] =
+//             this->owner().thermo().carrier().Y()[i][celli]
+//            /this->owner().thermo().carrier().Wi(i);
+//     }
+
+//     return Xc/sum(Xc);
+// }
+
+
+// template<class CloudType>
+// Foam::scalar Foam::LiquidEvaporationBoil<CloudType>::Sh
+// (
+//     const scalar Re,
+//     const scalar Sc
+// ) const
+// {
+//     return 2.0 + 0.6*Foam::sqrt(Re)*cbrt(Sc);
+// }
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::FGMLiquidEvaporationBoil<CloudType>::FGMLiquidEvaporationBoil
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    LiquidEvaporationBoil<CloudType>(dict, owner),
+    W_(this->owner().mesh().objectRegistry::lookupObject<const volScalarField>("Wt"))
+{
+
+}
+
+
+template<class CloudType>
+Foam::FGMLiquidEvaporationBoil<CloudType>::FGMLiquidEvaporationBoil
+(
+    const FGMLiquidEvaporationBoil<CloudType>& pcm
+)
+:
+    LiquidEvaporationBoil<CloudType>(pcm),
+
+    W_(this->owner().mesh().objectRegistry::lookupObject<const volScalarField>("Wt"))
+{
+
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::FGMLiquidEvaporationBoil<CloudType>::~FGMLiquidEvaporationBoil()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::FGMLiquidEvaporationBoil<CloudType>::calculate
+(
+    const scalar dt,
+    const label celli,
+    const scalar Re,
+    const scalar Pr,
+    const scalar d,
+    const scalar nu,
+    const scalar T,
+    const scalar Ts,
+    const scalar pc,
+    const scalar Tc,
+    const scalarField& X,
+    scalarField& dMassPC
+) const
+{
+    // immediately evaporate mass that has reached critical condition
+    if ((this->liquids_.Tc(X) - T) < small)
+    {
+        if (debug)
+        {
+            WarningInFunction
+                << "Parcel reached critical conditions: "
+                << "evaporating all available mass" << endl;
+        }
+
+        forAll(this->activeLiquids_, i)
+        {
+            const label lid = this->liqToLiqMap_[i];
+            dMassPC[lid] = great;
+        }
+
+        return;
+    }
+
+    // droplet surface pressure assumed to surface vapour pressure
+    scalar ps = this->liquids_.pv(pc, Ts, X);
+
+    // vapour density at droplet surface [kg/m^3]
+    const scalar RR = 1000.0*constant::physicoChemical::R.value(); // J/(kmol·k)
+    // vapour density at droplet surface [kg/m3]
+    scalar rhos = ps*this->liquids_.W(X)/(RR*Ts);
+
+    // // construct carrier phase species volume fractions for cell, celli
+    // const scalarField XcMix(calcXc(celli));
+
+    // // carrier thermo properties
+    // scalar Hsc = 0.0;
+    // scalar Hc = 0.0;
+    // scalar Cpc = 0.0;
+    // scalar kappac = 0.0;
+    // this->owner().thermo().carrier().calcH(Tc, pc);
+    // forAll(this->owner().thermo().carrier().Y(), i)
+    // {
+    //     const scalar Yc = this->owner().thermo().carrier().Y()[i][celli];
+    //     Hc += Yc*this->owner().thermo().carrier().Ha(i, pc, Tc);
+    // }
+
+    // this->owner().thermo().carrier().calcH(Ts, ps);
+    // this->owner().thermo().carrier().calcCp(Ts, ps);
+    // this->owner().thermo().carrier().calcMu(Ts, ps);
+    // forAll(this->owner().thermo().carrier().Y(), i)
+    // {
+    //     const scalar Yc = this->owner().thermo().carrier().Y()[i][celli];
+    //     Hsc += Yc*this->owner().thermo().carrier().Ha(i, ps, Ts);
+    //     Cpc += Yc*this->owner().thermo().carrier().Cp(i, ps, Ts);
+    //     kappac += Yc*this->owner().thermo().carrier().kappa(i, ps, Ts);
+    // }
+
+    // calculate mass transfer of each specie in liquid
+    forAll(this->activeLiquids_, i)
+    {
+        const label gid = this->liqToCarrierMap_[i];
+        const label lid = this->liqToLiqMap_[i];
+
+        // boiling temperature at cell pressure for liquid species lid [K]
+        const scalar TBoil = this->liquids_.properties()[lid].pvInvert(pc);
+        // limit droplet temperature to boiling/critical temperature
+        const scalar Td = min(T, 0.999*TBoil);
+        // saturation pressure for liquid species lid [Pa]
+        const scalar pSat = this->liquids_.properties()[lid].pv(pc, Td);
+        // surface molar fraction - Raoult's Law
+        const scalar Xs = X[lid]*pSat/pc;
+
+
+        // if (Xc*pc > pSat)
+        // {
+        //     // saturated vapour - no phase change
+        // }
+        // else
+        {
+            // vapour diffusivity [m2/s]
+            const scalar Dab = this->liquids_.properties()[lid].D(ps, Ts);
+
+            // Schmidt number
+            const scalar Sc = nu/(Dab + rootVSmall);
+
+            // Sherwood number
+            const scalar Sh = this->Sh(Re, Sc);
+
+            // mixture molar weight
+            scalar W_gas = W_[celli];//FGMmodel_.getValue(celli, "W");
+
+            // carrier phase concentration-fuel
+            const scalar Yc = this->owner().thermo().carrier().Y()[gid][celli];
+
+            // fuel molar weight [kg/kmol]
+            const scalar W_fuel = this->liquids_.properties()[lid].W();
+            const scalar Ys = Xs*W_fuel/W_gas/(1. + Xs*W_fuel/W_gas - Xs); // ref Hu.2017, Spalding.1953
+            const scalar Bm = (Ys - Yc)/max(small, 1. - Ys);
+
+            if (Bm > 0)
+            {
+                // mass transfer [kg]
+                dMassPC[lid] += pi*d*Sh*Dab*rhos*log(1.0 + Bm)*dt;
+            }
+        }
+    }
+}
+
+
+
+
+// ************************************************************************* //
diff --git a/src/dfCombustionModels/FGM/FGMLiquidEvaporationBoil/FGMLiquidEvaporationBoil.H b/src/dfCombustionModels/FGM/FGMLiquidEvaporationBoil/FGMLiquidEvaporationBoil.H
new file mode 100644
index 000000000..699177630
--- /dev/null
+++ b/src/dfCombustionModels/FGM/FGMLiquidEvaporationBoil/FGMLiquidEvaporationBoil.H
@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::FGMLiquidEvaporationBoil
+
+Description
+    Liquid evaporation model
+    - uses ideal gas assumption
+    - includes boiling model based on:
+
+    \verbatim
+        "Studies of Superheated Fuel Spray Structures and Vaporization in
+        GDI Engines"
+
+        Zuo, B., Gomes, A. M. and Rutland C. J.
+
+        International Journal of Engine Research, 2000, Vol. 1(4), pp. 321-336
+    \endverbatim
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef FGMLiquidEvaporationBoil_H
+#define FGMLiquidEvaporationBoil_H
+
+#include "LiquidEvaporationBoil.H"
+#include "liquidMixtureProperties.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                   Class FGMLiquidEvaporationBoil Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class FGMLiquidEvaporationBoil
+:
+    public LiquidEvaporationBoil<CloudType>
+{
+protected:
+
+    // Protected data
+
+        // //- Global liquid properties data
+        // const liquidMixtureProperties& liquids_;
+
+        // //- List of active liquid names
+        // List<word> activeLiquids_;
+
+        // //- Mapping between liquid and carrier species
+        // List<label> liqToCarrierMap_;
+
+        // //- Mapping between local and global liquid species
+        // List<label> liqToLiqMap_;
+
+        //- molar weight of gas mixture
+        const volScalarField& W_;
+
+
+    // // Protected Member Functions
+
+    //     //- Sherwood number as a function of Reynolds and Schmidt numbers
+    //     scalar Sh(const scalar Re, const scalar Sc) const;
+
+    //     //- Calculate the carrier phase component volume fractions at celli
+    //     tmp<scalarField> calcXc(const label celli) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("FGMliquidEvaporationBoil");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        FGMLiquidEvaporationBoil(const dictionary& dict, CloudType& cloud);
+
+        //- Construct copy
+        FGMLiquidEvaporationBoil(const FGMLiquidEvaporationBoil<CloudType>& pcm);
+
+        //- Construct and return a clone
+        virtual autoPtr<PhaseChangeModel<CloudType>> clone() const
+        {
+            return autoPtr<PhaseChangeModel<CloudType>>
+            (
+                new FGMLiquidEvaporationBoil<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~FGMLiquidEvaporationBoil();
+
+
+    // Member Functions
+
+        //- Update model
+        virtual void calculate
+        (
+            const scalar dt,
+            const label celli,
+            const scalar Re,
+            const scalar Pr,
+            const scalar d,
+            const scalar nu,
+            const scalar T,
+            const scalar Ts,
+            const scalar pc,
+            const scalar Tc,
+            const scalarField& X,
+            scalarField& dMassPC
+        ) const;
+
+        // //- Return the enthalpy per unit mass
+        // virtual scalar dh
+        // (
+        //     const label idc,
+        //     const label idl,
+        //     const scalar p,
+        //     const scalar T
+        // ) const;
+
+        // //- Return vapourisation temperature
+        // virtual scalar Tvap(const scalarField& X) const;
+
+        // //- Return maximum/limiting temperature
+        // virtual scalar TMax(const scalar p, const scalarField& X) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "FGMLiquidEvaporationBoil.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/dfCombustionModels/FGM/FGMLiquidEvaporationBoil/FGMLiquidEvaporationBoils.C b/src/dfCombustionModels/FGM/FGMLiquidEvaporationBoil/FGMLiquidEvaporationBoils.C
new file mode 100644
index 000000000..8a888c5aa
--- /dev/null
+++ b/src/dfCombustionModels/FGM/FGMLiquidEvaporationBoil/FGMLiquidEvaporationBoils.C
@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "basicSprayCloud.H"
+#include "basicReactingCloud.H"
+
+
+#include "FGMLiquidEvaporationBoil.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+makePhaseChangeModelType(FGMLiquidEvaporationBoil, basicSprayCloud );
+
+
+// ************************************************************************* //
diff --git a/src/dfCombustionModels/FGM/baseFGM/baseFGM.C b/src/dfCombustionModels/FGM/baseFGM/baseFGM.C
index 850656a23..d2654eb84 100644
--- a/src/dfCombustionModels/FGM/baseFGM/baseFGM.C
+++ b/src/dfCombustionModels/FGM/baseFGM/baseFGM.C
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "baseFGM.H"
+#include "basicSprayCloud.H" 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -37,6 +38,7 @@ Foam::combustionModels::baseFGM<ReactionThermo>::baseFGM
 )
 :
     laminar<ReactionThermo>(modelType, thermo, turb, combustionProperties),
+    fvOptions(fv::options::New(this->mesh())),
     buffer_(this->coeffs().lookupOrDefault("buffer", false)),
     scaledPV_(this->coeffs().lookupOrDefault("scaledPV", false)),
     incompPref_(this->coeffs().lookupOrDefault("incompPref", -10.0)),
@@ -44,7 +46,7 @@ Foam::combustionModels::baseFGM<ReactionThermo>::baseFGM
     combustion_(this->coeffs().lookupOrDefault("combustion", false)),
     solveEnthalpy_(this->coeffs().lookupOrDefault("solveEnthalpy", false)),
     flameletT_(this->coeffs().lookupOrDefault("flameletT", false)),
-    // psi_(const_cast<volScalarField&>(dynamic_cast<psiThermo&>(thermo).psi())),
+    tablePath_(this->coeffs().lookup("tablePath")),
     psi_(const_cast<volScalarField&>(dynamic_cast<rhoThermo&>(thermo).psi())),
     Wt_ 
     (
@@ -100,14 +102,14 @@ Foam::combustionModels::baseFGM<ReactionThermo>::baseFGM
     (
         IOobject
         (
-            "He",
+            "Ha",
             this->mesh().time().timeName(),
             this->mesh(),
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         this->mesh(),
-        dimensionedScalar("He",dimensionSet(0,2,-2,0,0,0,0),0.0)
+        dimensionedScalar("Ha",dimensionSet(0,2,-2,0,0,0,0),0.0)
     ),
     Hf_
     (
@@ -210,47 +212,97 @@ Foam::combustionModels::baseFGM<ReactionThermo>::baseFGM
         this->mesh(),  
          dimensionedScalar("chi_Zc",dimensionSet(0,0,-1,0,0,0,0),0.0)
     ),
-    YCO2_  
+    chi_Zfltd_
     (
         IOobject
         (
-            "YCO2",
+            "chi_Zfltd",
+            this->mesh().time().timeName(),
+            this->mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        this->mesh(),  
+        dimensionedScalar("chi_Zfltd",dimensionSet(0,0,-1,0,0,0,0),0)
+    ),
+    betac_  
+    (
+        IOobject
+        (
+            "betac",
             this->mesh().time().timeName(),
             this->mesh(),
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
         this->mesh(),  
-         dimensionedScalar("YCO2",dimless,0.0)  
+        dimensionedScalar("betac",dimensionSet(0,0,0,0,0,0,0),4.0) 
+    ),  
+    He_s_
+    (
+        IOobject
+        (
+            "Ha_s",
+            this->mesh().time().timeName(),
+            this->mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh(),
+        dimensionedScalar("Ha_s",dimensionSet(0,2,-2,0,0,0,0),0.0)
     ),
-    YCO_  
+    Zvar_s_
     (
         IOobject
         (
-            "YCO",
+            "Zvar_s",
             this->mesh().time().timeName(),
             this->mesh(),
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
         this->mesh(),  
-        dimensionedScalar("YCO",dimless,0.0)  
+        dimensionedScalar("Zvar_s",dimensionSet(0,0,0,0,0,0,0),0.0) 
     ),
-    YH2O_  
+    c_s_
+    (
+        IOobject
+        (
+            "c_s",
+            this->mesh().time().timeName(),
+            this->mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh(),  
+        dimensionedScalar("c_s",dimensionSet(0,0,0,0,0,0,0),0.0) 
+    ), 
+    cvar_s_
+    (
+        IOobject
+        (
+            "cvar_s",
+            this->mesh().time().timeName(),
+            this->mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh(),  
+        dimensionedScalar("cvar_s",dimensionSet(0,0,0,0,0,0,0),0.0) 
+    ),    
+    Zcvar_s_
     (
         IOobject
         (
-            "YH2O",
+            "Zcvar_s",
             this->mesh().time().timeName(),
             this->mesh(),
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
         this->mesh(),  
-        dimensionedScalar("YH2O",dimless,0.0) 
+        dimensionedScalar("Zcvar_s",dimensionSet(0,0,0,0,0,0,0),0.0) 
     ),
-    speciesNames_(this->coeffs().lookup("speciesName")),
-    Y_(speciesNames_.size()),
     cOmega_c_(omega_c_),
     ZOmega_c_(omega_c_), 
     WtCells_ (Wt_.primitiveFieldRef()),
@@ -267,10 +319,9 @@ Foam::combustionModels::baseFGM<ReactionThermo>::baseFGM
     ZOmega_cCells_(ZOmega_c_.primitiveFieldRef()),
     chi_cCells_(chi_c_.primitiveFieldRef()),
     chi_ZCells_(chi_Z_.primitiveFieldRef()),
-    chi_ZcCells_(chi_Zc_.primitiveFieldRef()), 
-    YCO2Cells_(YCO2_.primitiveFieldRef()),     
-    YCOCells_(YCO_.primitiveFieldRef()),   
-    YH2OCells_(YH2O_.primitiveFieldRef()),        
+    chi_ZcCells_(chi_Zc_.primitiveFieldRef()),     
+    chi_ZfltdCells_(chi_Zfltd_.primitiveFieldRef()),  
+    betacCells_(betac_.primitiveFieldRef()),     
     ZMax_(1.0),
     ZMin_(0.0),
     ZvarMax_(0.25),
@@ -286,6 +337,7 @@ Foam::combustionModels::baseFGM<ReactionThermo>::baseFGM
     rho_inRhoThermo_(dynamic_cast<rhoThermo&>(thermo).rho()),
     p_(this->thermo().p()),
     T_(this->thermo().T()),
+    Y_(this->chemistryPtr_->Y()),
     U_(this->mesh().objectRegistry::lookupObject<volVectorField>("U")),
     dpdt_(this->mesh().objectRegistry::lookupObject<volScalarField>("dpdt")),        
     phi_(this->mesh().objectRegistry::lookupObject<surfaceScalarField>("phi")),
@@ -301,7 +353,10 @@ Foam::combustionModels::baseFGM<ReactionThermo>::baseFGM
     Sc_(this->coeffs().lookupOrDefault("Sc", 1.0)),   
     bufferTime_(this->coeffs().lookupOrDefault("bufferTime", 0.0)),
     relaxation_(this->coeffs().lookupOrDefault("relaxation", false)),
-    DpDt_(this->coeffs().lookupOrDefault("DpDt", false))
+    DpDt_(this->coeffs().lookupOrDefault("DpDt", false)),
+    magUPrime_(U_.component(vector::X)),
+    magUPrimeCells_(magUPrime_.primitiveFieldRef())
+    // isLES_(this->coeffs().lookupOrDefault("isLES", false)),  
 {
     if(incompPref_ < 0.0)  //local pressure used to calculate EOS
     {
@@ -313,29 +368,43 @@ Foam::combustionModels::baseFGM<ReactionThermo>::baseFGM
         << incompPref_ << " Pa" << nl << endl;
     }
 
-    //initialize species fields
-    forAll(Y_, speciesI)
+    //- LES
+    this->isLES_ = this->mesh().objectRegistry::foundObject<compressible::LESModel>(turbulenceModel::propertiesName);
+
+    if (this->isLES_)
     {
+        const Foam::compressible::LESModel& lesModel = this->mesh().objectRegistry::lookupObject<compressible::LESModel>(turbulenceModel::propertiesName);
+            
+        const scalarField& LESdeltaCells=lesModel.delta().internalField();
 
-     Y_.set
-     (
-      speciesI,
-        new volScalarField
-        (
-            IOobject
-            (
-                "Y_"  + speciesNames_[speciesI],
-                this->mesh().time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            this->mesh(),
-            dimensionedScalar(dimless, Zero)
-        )
-     );
+        this->deltaCells_ = LESdeltaCells;
+
+    }
 
-   }
+
+//     //initialize species fields
+//     forAll(Y_, speciesI)
+//     {
+
+//      Y_.set
+//      (
+//       speciesI,
+//         new volScalarField
+//         (
+//             IOobject
+//             (
+//                 "Y_"  + speciesNames_[speciesI],
+//                 this->mesh().time().timeName(),
+//                 this->mesh(),
+//                 IOobject::NO_READ,
+//                 IOobject::AUTO_WRITE
+//             ),
+//             this->mesh(),
+//             dimensionedScalar(dimless, Zero)
+//         )
+//      );
+
+//    }
 
    // add fields
      fields_.add(Z_);
@@ -381,68 +450,250 @@ void Foam::combustionModels::baseFGM<ReactionThermo>::transport()
     );
 
 
-    // Solve the mixture fraction transport equation
-    if(buffer_) 
+    if (this->mesh().objectRegistry::foundObject<sprayCloud>(Foam::sprayCloud::typeName))
+
     {
-        Info<< "UW schemes used for scalars" << endl;
+        Info << "Solve Spray Z" << endl;
+ 
+        //- create spray object 
+        basicSprayCloud& spray = dynamic_cast<basicSprayCloud&>
+        (this->mesh().objectRegistry::lookupObjectRef<sprayCloud>(Foam::sprayCloud::typeName));        
+
+        tmp<fvScalarMatrix> tZSource(new fvScalarMatrix(this->Z_, dimMass/dimTime));
+        fvScalarMatrix& SZ = tZSource.ref();
+        forAll(spray.rhoTrans(), i)
+        {
+            SZ.source() += spray.rhoTrans()[i]/this->mesh().time().deltaT();
+            spray.rhoTrans()[i].writeOpt() = IOobject::NO_WRITE;
+        }
 
-        fvScalarMatrix ZEqn
-        (
-            fvm::ddt(rho_,Z_)
-            +scalarUWConvection->fvmDiv(phi_, Z_)  
-            -fvm::laplacian( mut/Sct_ + mu/Sc_ , Z_)     
-        );
-        if(relaxation_)
+        // Solve the mixture fraction transport equation
+        if(buffer_) 
         {
-            ZEqn.relax();
+            Info<< "UW schemes used for scalars" << endl;
+
+            fvScalarMatrix ZEqn
+            (
+                fvm::ddt(rho_,Z_)
+                +scalarUWConvection->fvmDiv(phi_, Z_)  
+                -fvm::laplacian( mut/Sct_ + mu/Sc_ , Z_)   
+                ==
+                - SZ  
+            );
+            if(relaxation_)
+            {
+                ZEqn.relax();
+            }
+            ZEqn.solve();
         }
-        ZEqn.solve();
-    }
 
+        else
+        {
+            Info<< "TVD schemes used for scalars" << endl;
+
+            fvScalarMatrix ZEqn
+            (
+                fvm::ddt(rho_,Z_)
+                +fvm::div(phi_, Z_)
+                -fvm::laplacian(mut/Sct_ + mu/Sc_ , Z_) 
+                ==
+                - SZ                  
+            );
+
+            if(relaxation_)
+            {
+                ZEqn.relax();
+            }
+            ZEqn.solve();
+        }
+        Z_.min(ZMax_);  
+        Z_.max(ZMin_);
+
+
+        if (combustion_)
+        {
+
+            tmp<fvScalarMatrix> tcSource(new fvScalarMatrix(this->c_, dimMass/dimTime));
+            fvScalarMatrix& S_c = tcSource.ref();
+            forAll(spray.rhoTrans(), i)
+            {
+                S_c.source() += spray.rhoTrans()[i]/this->mesh().time().deltaT();
+                spray.rhoTrans()[i].writeOpt() = IOobject::NO_WRITE;
+            }
+
+            // At initial time, cMax is set as Ycmaxall when unscaled PV employed
+            if(this->mesh().time().timeIndex() == 1 && !scaledPV_) cMax_ = Ycmaxall_;  
+
+            // Solve the progress variable transport equation
+            fvScalarMatrix cEqn
+            (
+                fvm::ddt(rho_, c_)
+                +(
+                    buffer_
+                    ? scalarUWConvection->fvmDiv(phi_, c_)
+                    : fvm::div(phi_, c_)
+                )
+                -fvm::laplacian( mut/Sct_ + mu/Sc_, c_)
+                -omega_c_
+                + fvOptions(rho_, c_)
+                ==
+                c_ * S_c
+            );  
+            fvOptions.correct(c_);
+            
+            if(relaxation_)
+            {
+                cEqn.relax();
+            }
+            cEqn.solve();
+            c_.min(cMax_);
+            c_.max(cMin_); 
+        } // end if (combustion_)
+
+    }
     else
     {
-        Info<< "TVD schemes used for scalars" << endl;
 
-        fvScalarMatrix ZEqn
-        (
-            fvm::ddt(rho_,Z_)
-            +fvm::div(phi_, Z_)
-            -fvm::laplacian(mut/Sct_ + mu/Sc_ , Z_) 
-        );
+        // Solve the mixture fraction transport equation
+        if(buffer_) 
+        {
+            Info<< "UW schemes used for scalars" << endl;
 
-        if(relaxation_)
+            fvScalarMatrix ZEqn
+            (
+                fvm::ddt(rho_,Z_)
+                +scalarUWConvection->fvmDiv(phi_, Z_)  
+                -fvm::laplacian( mut/Sct_ + mu/Sc_ , Z_)     
+            );
+            if(relaxation_)
+            {
+                ZEqn.relax();
+            }
+            ZEqn.solve();
+        }
+
+        else
         {
-            ZEqn.relax();
+            Info<< "TVD schemes used for scalars" << endl;
+
+            fvScalarMatrix ZEqn
+            (
+                fvm::ddt(rho_,Z_)
+                +fvm::div(phi_, Z_)
+                -fvm::laplacian(mut/Sct_ + mu/Sc_ , Z_) 
+            );
+
+            if(relaxation_)
+            {
+                ZEqn.relax();
+            }
+            ZEqn.solve();
         }
-        ZEqn.solve();
+        Z_.min(ZMax_);  
+        Z_.max(ZMin_);
+
+        if (combustion_)
+        {
+
+            // At initial time, cMax is set as Ycmaxall when unscaled PV employed
+            if(this->mesh().time().timeIndex() == 1 && !scaledPV_) cMax_ = Ycmaxall_;  
+
+            // Solve the progress variable transport equation
+            fvScalarMatrix cEqn
+            (
+                fvm::ddt(rho_, c_)
+                +(
+                    buffer_
+                    ? scalarUWConvection->fvmDiv(phi_, c_)
+                    : fvm::div(phi_, c_)
+                )
+                -fvm::laplacian( mut/Sct_ + mu/Sc_, c_)
+                -omega_c_
+                + fvOptions(rho_, c_)
+            );  
+            fvOptions.correct(c_);
+            
+            if(relaxation_)
+            {
+                cEqn.relax();
+            }
+            cEqn.solve();
+            c_.min(cMax_);
+            c_.max(cMin_); 
+        } // end if (combustion_)
     }
-    Z_.min(ZMax_);  
-    Z_.max(ZMin_);
 
    // solve the total enthalpy transport equation
     if(solveEnthalpy_)
     {
-        fvScalarMatrix HEqn
-        (
-            fvm::ddt(rho_,He_)
-            +(
-                buffer_
-                ? scalarUWConvection->fvmDiv(phi_, He_)
-                : fvm::div(phi_, He_)
-            )
-            +(
-                DpDt_
-                ? - dpdt_ - ( U_ & fvc::grad(p_) ) - fvm::laplacian( mut/Sct_ + mu/Sc_, He_)
-                : - fvm::laplacian( mut/Sct_ + mu/Sc_, He_) 
-            )
-        ); 
+        if (this->mesh().objectRegistry::foundObject<sprayCloud>(Foam::sprayCloud::typeName))
+        {
 
-        if(relaxation_)
+            Info << "Solve Spray He" << endl;
+            //- create spray object 
+            basicSprayCloud& spray = dynamic_cast<basicSprayCloud&>
+            (this->mesh().objectRegistry::lookupObjectRef<sprayCloud>(Foam::sprayCloud::typeName));     
+                      
+            tmp<fvScalarMatrix> thSource(new fvScalarMatrix(this->He_, dimEnergy/dimTime));
+            fvScalarMatrix& hSource = thSource.ref();
+
+            forAll(this->chemistryPtr_->species(), i)
+            {
+                hSource.source() += spray.rhoTrans(i)*this->chemistryPtr_->hci(i)/this->mesh().time().deltaT();
+            }
+
+            hSource.source() += spray.hsTrans()/this->mesh().time().deltaT();
+
+
+            fvScalarMatrix HEqn
+            (
+                fvm::ddt(rho_,He_)
+                +(
+                    buffer_
+                    ? scalarUWConvection->fvmDiv(phi_, He_)
+                    : fvm::div(phi_, He_)
+                )
+                +(
+                    DpDt_
+                    ? - dpdt_ - ( U_ & fvc::grad(p_) ) - fvm::laplacian( mut/Sct_ + mu/Sc_, He_)
+                    : - fvm::laplacian( mut/Sct_ + mu/Sc_, He_) 
+                )
+                ==
+                - hSource
+            ); 
+
+            if(relaxation_)
+            {
+                HEqn.relax();
+            }
+            HEqn.solve();
+        }
+        else
         {
-            HEqn.relax();
+            fvScalarMatrix HEqn
+            (
+                fvm::ddt(rho_,He_)
+                +(
+                    buffer_
+                    ? scalarUWConvection->fvmDiv(phi_, He_)
+                    : fvm::div(phi_, He_)
+                )
+                +(
+                    DpDt_
+                    ? - dpdt_ - ( U_ & fvc::grad(p_) ) - fvm::laplacian( mut/Sct_ + mu/Sc_, He_)
+                    : - fvm::laplacian( mut/Sct_ + mu/Sc_, He_) 
+                )
+            ); 
+
+            if(relaxation_)
+            {
+                HEqn.relax();
+            }
+            HEqn.solve();
         }
-        HEqn.solve();
-    }
+ 
+    }    
+
 
 
     // Solve the mixture fraction variance transport equation
@@ -470,31 +721,15 @@ void Foam::combustionModels::baseFGM<ReactionThermo>::transport()
     Zvar_.min(ZvarMax_);
     Zvar_.max(ZvarMin_);    
 
+    if (this->isLES_)
+    {
+        chi_Zfltd_ = this->turbulence().mu()/rho_*(fvc::grad(Z_) & fvc::grad(Z_)); 
+    }
+
+
 
     if(combustion_ && reactFlowTime_ > 0.0)
     {
-        // At initial time, cMax is set as Ycmaxall when unscaled PV employed
-        if(this->mesh().time().timeIndex() == 1 && !scaledPV_) cMax_ = Ycmaxall_;  
-
-        // Solve the progress variable transport equation
-        fvScalarMatrix cEqn
-        (
-            fvm::ddt(rho_, c_)
-            +(
-                buffer_
-                ? scalarUWConvection->fvmDiv(phi_, c_)
-                : fvm::div(phi_, c_)
-            )
-            -fvm::laplacian( mut/Sct_ + mu/Sc_, c_)
-            -omega_c_
-        );  
-        if(relaxation_)
-        {
-            cEqn.relax();
-        }
-        cEqn.solve();
-        c_.min(cMax_);
-        c_.max(cMin_); 
 
         // Solve the progress variable variance transport equation
         fvScalarMatrix cvarEqn
@@ -530,8 +765,7 @@ void Foam::combustionModels::baseFGM<ReactionThermo>::transport()
                 :  fvm::div(phi_, Zcvar_)
             )
             -fvm::laplacian( mut/Sct_+mu/Sc_, Zcvar_)
-            -(1.0*mut/Sct_*(fvc::grad(c_) & fvc::grad(Z_)))
-            -(1.0*mut/Sct_*(fvc::grad(Z_) & fvc::grad(c_)))
+            -(2.0*mut/Sct_*(fvc::grad(c_) & fvc::grad(Z_)))
             +(2.0*rho_*chi_Zc_)  
             -1.0*(ZOmega_c_-omega_c_*Z_)  
         );
@@ -581,21 +815,14 @@ void Foam::combustionModels::baseFGM<ReactionThermo>::initialiseFalmeKernel()
 }
 
 template<class ReactionThermo>
-Foam::tmp<Foam::fvScalarMatrix>
-Foam::combustionModels::baseFGM<ReactionThermo>::R(volScalarField& Y) const
+void Foam::combustionModels::baseFGM<ReactionThermo>::magUPrime()
 {
-    return laminar<ReactionThermo>::R(Y);
-}
+    autoPtr<LESfilter> filterPtr_(LESfilter::New(this->mesh(), this->turbulence().coeffDict()));
+    LESfilter& filter_(filterPtr_());
 
-template<class ReactionThermo>
-Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::baseFGM<ReactionThermo>::Qdot() const
-{
-    return volScalarField::New
-    (
-        this->thermo().phasePropertyName("Qdot"),
-        laminar<ReactionThermo>::Qdot()
-    );
+    magUPrime_=mag(U_-filter_(U_));
+
+   // return magUPrime_;
 }
 
 
diff --git a/src/dfCombustionModels/FGM/baseFGM/baseFGM.H b/src/dfCombustionModels/FGM/baseFGM/baseFGM.H
index 59863fbd9..f263d0a31 100644
--- a/src/dfCombustionModels/FGM/baseFGM/baseFGM.H
+++ b/src/dfCombustionModels/FGM/baseFGM/baseFGM.H
@@ -42,6 +42,8 @@ SourceFiles
 #include "../laminar/laminar.H"
 #include "uniformDimensionedFields.H"
 #include "CanteraMixture.H"
+#include "fvOptions.H"
+#include "LESfilter.H"
 //#include "basicThermo.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -61,7 +63,7 @@ class baseFGM
 {
     // Protected Data
  protected:
-
+    fv::options& fvOptions; 
     //- buffer switch
     Switch buffer_;
 
@@ -82,7 +84,9 @@ class baseFGM
     
     //- T-lookup switch
     Switch flameletT_;
- 
+
+    word tablePath_;
+
     //- psi=\rho/p
     volScalarField& psi_;
 
@@ -125,16 +129,17 @@ class baseFGM
     //-covariance Z-c dissipation rate (s^{-1})
     volScalarField chi_Zc_;
 
-    //- default species needed to lookup
-    volScalarField YCO2_;
-    volScalarField YCO_;    
-    volScalarField YH2O_;
+    //- filtered SDR of mixture fraction
+    volScalarField chi_Zfltd_; 
 
-    //- self-defined species name
-    hashedWordList speciesNames_;
-    
-    //- self-defined speices mass fraction 
-    PtrList<volScalarField> Y_;
+    //- betac for LES
+    volScalarField betac_;
+
+    volScalarField He_s_;
+    volScalarField Zvar_s_; 
+    volScalarField c_s_;
+    volScalarField cvar_s_;
+    volScalarField Zcvar_s_;
 
     //- copy of omega_c_
     volScalarField cOmega_c_;
@@ -158,9 +163,12 @@ class baseFGM
     scalarField& chi_ZCells_;
     scalarField& chi_ZcCells_;
 
-    scalarField& YCO2Cells_; 
-    scalarField& YCOCells_; 
-    scalarField& YH2OCells_;    
+    //- filter size
+    scalarField deltaCells_; 
+
+    scalarField& chi_ZfltdCells_;
+
+    scalarField& betacCells_; 
 
     //- max and min Z
     scalar ZMax_;
@@ -183,6 +191,7 @@ class baseFGM
     scalar ZcvarMax_;
     scalar ZcvarMin_;   
 
+
     //- fluid density
      volScalarField& rho_;
     //- fluid density
@@ -194,6 +203,8 @@ class baseFGM
     //- fluid temperature
     volScalarField& T_;
 
+    PtrList<volScalarField>& Y_;
+
     //- fluid velocity
     const volVectorField& U_;     
 
@@ -238,6 +249,14 @@ class baseFGM
 
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields_;
 
+    //- LES
+    //- if in LES mode
+    Switch isLES_;
+
+    //SGS
+    volScalarField magUPrime_;
+    scalarField& magUPrimeCells_; 
+
  public:
     //- Runtime type information
     TypeName("baseFGM");
@@ -266,11 +285,8 @@ class baseFGM
         //- initialize flame kernel 
         virtual void initialiseFalmeKernel();        
 
-        //- Fuel consumption rate matrix.
-        virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
-
-        //- Heat release rate [kg/m/s^3]
-        virtual tmp<volScalarField> Qdot() const;
+        //- return mag SGS
+        virtual void magUPrime();
 
         //- reaturn buffer_
         virtual bool buffer();
diff --git a/src/dfCombustionModels/FGM/flameletTableSolver/readThermChemTables.H b/src/dfCombustionModels/FGM/flameletTableSolver/readThermChemTables.H
index dbea81fde..b625b6a8a 100644
--- a/src/dfCombustionModels/FGM/flameletTableSolver/readThermChemTables.H
+++ b/src/dfCombustionModels/FGM/flameletTableSolver/readThermChemTables.H
@@ -1,199 +1,47 @@
-    Info<< "Reading " << H_fuel << "\n" << endl;
-    Info<< "Reading " << H_ox << "\n" << endl;
-
-    z_Tb5={new double[NZL]{}};
-    sl_Tb5={new double[NZL]{}};
-    th_Tb5={new double[NZL]{}};
-    tau_Tb5={new double[NZL]{}};
-    kctau_Tb5={new double[NZL]{}};
-
-    Info<< "Reading laminar flame properties\n" << endl;
-
-    for(int jj = 0; jj < NZL; jj++)  
-    {
-         f1=fscanf( table, "%lf %lf %lf %lf %lf\n",
-                &z_Tb5[jj],&sl_Tb5[jj],&th_Tb5[jj],
-                &tau_Tb5[jj],&kctau_Tb5[jj] );
-    }
-
-    fNZ=fscanf(table, "%d",&NZ);
-    z_Tb3={ new double[NZ]{} };
-    for(int ii = 0; ii < NZ; ii++)
-    {
-         f2=fscanf( table, "%lf\n", &z_Tb3[ii] );   //load z to z_Tb3
-    }
-
-    fNC=fscanf(table, "%d",&NC);
-    c_Tb3={ new double[NC]{} };
-    for(int ii = 0; ii < NC; ii++)  
-    {
-         f3=fscanf( table, "%lf\n", &c_Tb3[ii] );  ///load c to c_Tb3
-    }
-
-    fNGZ=fscanf(table, "%d",&NGZ);
-    gz_Tb3={ new double[NGZ]{} };
-    for(int ii = 0; ii < NGZ; ii++)   
-    {
-         f4=fscanf( table, "%lf\n", &gz_Tb3[ii] );  //load z-varance to gz_Tb3
-    }
-
-    fNGC=fscanf(table, "%d",&NGC);
-    gc_Tb3={ new double[NGC]{} };
-    for(int ii = 0; ii < NGC; ii++)
-    {
-         f5=fscanf( table, "%lf\n", &gc_Tb3[ii] );   //load c-variance to gc_Tb3
-    }   
-
-    fNZC=fscanf(table, "%d",&NZC);
-    gzc_Tb3={ new double[NZC]{} };
-    for(int ii = 0; ii < NZC; ii++)
-    {
-         f6=fscanf( table, "%lf\n", &gzc_Tb3[ii] );   //load zc-vovariance to gzc_Tb3
-    }    
-
-    fNS_NY=fscanf(table, "%d %d",&NS,&NY);
-
-    omgc_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} };     
-    cOc_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} };      
-    ZOc_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} };     
-    cp_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} };       
-    hiyi_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} };     
-    mwt_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} };      
-    Tf_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} };       
-    nu_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} };       
-    Ycmax_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} };   
-
-    Yi01_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} },    
-    Yi02_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} },    
-    Yi03_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} };    
-
-    Info << "NS is read as: " << NS << endl;
-
-    if(NS == 8)
-    {
-        scaledPV_ = true;
-        Info<< "=============== Using scaled PV ==============="
-            << "\n" << endl;
-    }
-    else if(NS == 9)
-    {
-        scaledPV_ = false;
-        Info<< "=============== Using unscaled PV ==============="
-            << "\n" << endl;
-    }
-    else
-    {
-        WarningInFunction << "Number of columns wrong in flare.tbl !!!"
-                            << "\n" << endl;
-    }
-
-    const char* fmt8{"%lf %lf %lf %lf %lf %lf %lf %lf"};
-    const char* fmt9{"%lf %lf %lf %lf %lf %lf %lf %lf %lf"};
-
-    //Note：if nYis is changed，fmt_nYis is also needed to change！
-    const char* fmt_nYis{"%lf %lf %lf"};
-
-
-    Info<< "Reading turbulent flame properties\n" << endl;   
-
-    int count = 0;    
-    for(int ii = 0; ii < NZ; ii++)
-    {
-        for(int jj = 0; jj < NC; jj++)
-        {
-            for(int kk = 0; kk < NGZ; kk++)
+            value_temp = std::vector<std::vector<double>> ((NY+NS),std::vector<double>(singleTableSize_));
+            count = 0;    //counter
+            for (int hh=0; hh<NH; ++hh)
             {
-                for(int ll = 0; ll < NGC; ll++)
+                for(int ii = 0; ii < NZ; ++ii)
                 {
-                    for(int mm = 0; mm < NZC; mm++)
+                    for(int jj = 0; jj < NC; ++jj)
                     {
-                        if(scaledPV_)  
+                        for(int kk = 0; kk < NGZ; ++kk)
                         {
-                             f7=fscanf
-                                     (
-                                        table,fmt8,&omgc_Tb3[count],&cOc_Tb3[count],&ZOc_Tb3[count],&cp_Tb3[count],
-                                        &mwt_Tb3[count],&hiyi_Tb3[count],&Tf_Tb3[count],&nu_Tb3[count]
-                                     );   
-                        }
-                        else
-                        {
-                             f8=fscanf
-                                     (
-                                        table,fmt9,&omgc_Tb3[count],&cOc_Tb3[count],&ZOc_Tb3[count],&cp_Tb3[count],
-                                        &mwt_Tb3[count],&hiyi_Tb3[count],&Tf_Tb3[count],&nu_Tb3[count],&Ycmax_Tb3[count]
-                                     );
-                        }
-
-                        if(NY > 0)   
-                        {
-                        
-                         f9=fscanf(
-                                     table,fmt_nYis,&Yi01_Tb3[count],&Yi02_Tb3[count],&Yi03_Tb3[count]
-                                  );   
-                                  
-                         //- loop speciesNames
-                    /*   forAll(Yi_Tb3,speciesI)  
-                         {
-                            f9=fscanf(
-                                        table,fmt_nYis,&Yi_Tb3[speciesI][count]
-                                    );   //再读入'H2O','CO','CO2'的数据
-                         }  
-                    */
+                            for(int ll = 0; ll < NGC; ++ll)
+                            {
+                                for(int mm = 0; mm < NZC; ++mm)
+                                {
+                                    if (std::getline(table, line))
+                                    {
+                                        std::istringstream iss(line);
+                                        for (int ivar=0; ivar<NY+NS; ++ivar)
+                                        {
+                                            iss >> value_temp[ivar][count];
+                                        } 
+                                    }
+                                    else
+                                    {
+                                        Info << "Error: no more line in: h=" << hh << ", Z=" << ii << ", c=" << jj << ", gz=" << kk << ", gc=" << ll <<", zc=" << mm <<endl;
+                                    }
+                                    count++;
+                                }
+                            }
                         }
+                    }
 
-                        count++;
+                    if(ii%10 == 0 && NH == 1) 
+                    {
+                        prog = 100.0*ii / ((NZ-1));
+                        Info << " Progress -- " << prog << "%" << endl;
                     }
                 }
-            }
-        }
-        if(ii%10 == 0)  
-        {
-            prog = 100.0*ii/(NZ-1);
-            Info<< " Progress -- " << prog << "%" << endl;
-        }
-    }
-
-    //- findthe maxmum PV value from Ycmax_Tb3
-    const scalar Ycmaxall{*std::max_element(Ycmax_Tb3,Ycmax_Tb3+count+1)};   
-    cMaxAll_ = Ycmaxall;
-
-    if(!scaledPV_)
-    {
-        Info<< "\nunscaled PV -- Ycmaxall = "<<Ycmaxall<< endl;   
-    }
-
-    if(flameletT_)  
-    {
-        delete [] cp_Tb3;
-        delete [] hiyi_Tb3;
-    }
-
-    Info<< "\nReading non-premixed properties\n" << endl;   
-    d2Yeq_Tb2={ new double[NZ*NGZ]{} };
-    int countN = 0;
-    for(int ii = 0; ii < NZ; ii++)
-    {
-        for(int jj = 0; jj < NGZ; jj++)
-        {
-             f10=fscanf( table, "%lf\n", &d2Yeq_Tb2[countN] );   
-
-            countN++;
-        }
-    }
-
-    fclose(table);
-
-    Info<< "* * * * * * * * * * Complete FlaRe table * * * * * * * * * \n"
-        << endl;
 
+                if(NH != 1) 
+                {
+                    prog = 100.0*hh / ((NH-1));
+                    Info << " Progress -- " << prog << "%" << endl;
+                }
+            }
 
-/*
-   forAll(Yi_Tb3, speciesI)
-   {
-        Yi_Tb3.set
-        (
-            speciesI,
-            new scalarField(NZ*NC*NGZ*NGC*NZC)
-        );
-   }
-*/ 
\ No newline at end of file
+    
diff --git a/src/dfCombustionModels/FGM/flameletTableSolver/tableSolver.C b/src/dfCombustionModels/FGM/flameletTableSolver/tableSolver.C
index 33f7292c5..600ff33b6 100644
--- a/src/dfCombustionModels/FGM/flameletTableSolver/tableSolver.C
+++ b/src/dfCombustionModels/FGM/flameletTableSolver/tableSolver.C
@@ -30,7 +30,7 @@ License
 #include "IFstream.H"
 #include "Pstream.H"
 #include "clockTime.H"
-
+#include "scalar.H"
 
 namespace Foam
 {
@@ -39,25 +39,307 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-//tableSolver::tableSolver(const wordList& tableNames, string suffix_)
-tableSolver::tableSolver(wordList speciesNames,Switch& scaledPV, Switch flameletT, scalar& cMaxAll)
+tableSolver::tableSolver(word tablePath)
 :
 small(1.0e-4),
 smaller(1.0e-6),
 smallest(1.0e-12),
 T0(298.15),
-table(fopen("./flare.tbl", "r")),
-speciesNames_(speciesNames),
-scaledPV_(scaledPV),
-flameletT_(flameletT),
-cMaxAll_(cMaxAll),
-fHox_fu(fscanf(table, "%lf %lf",&Hfu,&Hox)),
-fNZL(fscanf(table, "%d",&NZL)),
 H_fuel("H_fuel",dimensionSet(0,2,-2,0,0,0,0),Hfu),
 H_ox("H_ox",dimensionSet(0,2,-2,0,0,0,0),Hox)
 {
+    Info<< "* * * * * * * * * * Opening FlaRe table * * * * * * * * * * \n" << endl;
+    std::ifstream table(tablePath.c_str());
+    if (table.is_open())
+    {
+        std::string line;
+
+        if (std::getline(table, line))
+        {
+            std::istringstream iss(line);
+            iss >> NH >> NZ >> NC >> NGZ >> NGC >> NZC >> NS >> NYomega >> NY >> NZL;
+        }
+        Info<< "Reading: NH=" << NH << ", NZ=" << NZ << ", NC="<< NC << ", NGZ=" << NGZ << ", NGC=" 
+            << NGC << ", NZC="<< NZC << ", NS=" << NS << ", NYomega=" << NYomega << ", NY=" << NY << ", NZL=" << NZL << "\n" << endl;
+
+        tableNames_ = wordList({"omgc_Tb3", "cOc_Tb3", "ZOc_Tb3", "cp_Tb3", "mwt_Tb3", "hiyi_Tb3", "Tf_Tb3", "nu_Tb3"});
+
+        if(NS == 8+NYomega)
+        {
+            scaledPV_ = true;
+            Info<< "=============== Using scaled PV ==============="
+                << "\n" << endl;
+        }
+        else if(NS == 9+NYomega)
+        {
+            scaledPV_ = false;
+            tableNames_.append("Ycmax_Tb3");
+            Info<< "=============== Using unscaled PV ==============="
+                << "\n" << endl;
+        }
+        else
+        {
+            WarningInFunction << "Number of columns wrong in flare.tbl !!!"
+                                << "\n" << endl;
+        }
+
+
+        if (std::getline(table, line))
+        {
+            std::istringstream iss(line);
+            std::string spc_name;
+            for (int ii=0; ii<NYomega; ++ii)
+            {
+                iss >> spc_name;
+                spc_omegaNames_table_.append(spc_name);
+            }
+        }
+
+        if (std::getline(table, line))
+        {
+            std::istringstream iss(line);
+            std::string spc_name;
+            for (int ii=0; ii<NY; ++ii)
+            {
+                iss >> spc_name;
+                speciesNames_table_.append(spc_name);
+                tableNames_.append(spc_name);
+            }
+        }
+
+
+        Info << "Load omega of species: " << spc_omegaNames_table_ << endl;
+        Info << "Load species: " << speciesNames_table_ << endl;
+
+        h_Tb3 = { new double[NH]{} };  z_Tb3 = { new double[NZ]{} }; 
+        c_Tb3 = { new double[NC]{} };  gz_Tb3 = { new double[NGZ]{} };
+        gc_Tb3 = { new double[NGC]{} };  gzc_Tb3 = { new double[NZC]{} };  
+
+        for(int ii = 0; ii < NH; ++ii)
+        {
+            if (std::getline(table, line))
+            {
+                std::istringstream iss(line);
+                iss >> h_Tb3[ii];
+            }
+        }
+        for(int ii = 0; ii < NZ; ++ii)
+        {
+            if (std::getline(table, line))
+            {
+                std::istringstream iss(line);
+                iss >> z_Tb3[ii];
+            }
+        }
+        for(int ii = 0; ii < NC; ++ii)
+        {
+            if (std::getline(table, line))
+            {
+                std::istringstream iss(line);
+                iss >> c_Tb3[ii];
+            }
+        }
+        for(int ii = 0; ii < NGZ; ++ii)
+        {
+            if (std::getline(table, line))
+            {
+                std::istringstream iss(line);
+                iss >> gz_Tb3[ii];
+            }
+        }
+        for(int ii = 0; ii < NGC; ++ii)
+        {
+            if (std::getline(table, line))
+            {
+                std::istringstream iss(line);
+                iss >> gc_Tb3[ii];
+            }
+        }
+        for(int ii = 0; ii < NZC; ++ii)
+        {
+            if (std::getline(table, line))
+            {
+                std::istringstream iss(line);
+                iss >> gzc_Tb3[ii];
+            }
+        }
+        
+        if (std::getline(table, line))
+        {
+            std::istringstream iss(line);
+            iss >> Hfu >> Hox;
+        }
+        Info<< "Reading: H_fuel=" << Hfu << ", H_ox=" << Hox << "\n" << endl;
+
+        z_Tb5 = { new double[NH*NZL]{} }; sl_Tb5 = { new double[NH*NZL]{} };
+        th_Tb5 = { new double[NH*NZL]{} }; tau_Tb5 = { new double[NH*NZL]{} };
+        kctau_Tb5 = { new double[NH*NZL]{} };
+
+        Info<< "Reading laminar flame properties\n" << endl;
+        int count = 0;
+        for (int ii=0; ii<NH; ++ii)
+        {
+            for(int jj = 0; jj < NZL; ++jj)  
+            {
+                if (std::getline(table, line))
+                {
+                    std::istringstream iss(line);
+                    iss >> z_Tb5[count] >> sl_Tb5[count] >> th_Tb5[count]
+                        >> tau_Tb5[count] >> kctau_Tb5[count];
+                    count++;
+                }
+            }
+        }
+
+        Info<< "Reading turbulence flame properties\n" << endl;
+
+        singleTableSize_ = NH*NZ*NC*NGZ*NGC*NZC;
+
+        if (Pstream::parRun()) // parallel computing
+        {
+            // Create node-local communicator
+            MPI_Comm nodecomm;
+            MPI_Comm_split_type(MPI_COMM_WORLD, MPI_COMM_TYPE_SHARED, Pstream::myProcNo(), MPI_INFO_NULL, &nodecomm);
+
+            int noderank;
+            MPI_Comm_rank(nodecomm, &noderank);
+
+            // wins_.resize(wins_.size()+1);
+            // const int disp_unit = sizeof(double);
+
+            // scalar *table_ptr;
+
+            // Initialise table on node-master with appropriate synchronisation
+            // MPI_Win_fence(0, wins_.last());
+            // MPI_Barrier(MPI_COMM_WORLD);
+
+            // if (noderank == 0) // node-master
+            // {
+            //     #include "readThermChemTables.H"
+            // } // end if (noderank == 0)
+
+            // MPI_Barrier(MPI_COMM_WORLD);
+
+            std::vector<std::vector<double>> value_temp;
+
+            if (noderank == 0)
+            {
+                #include "readThermChemTables.H"
+            }
+            else
+            {
+                for (int ii=0;ii<singleTableSize_; ++ii) std::getline(table, line);
+            }
+            
+            MPI_Barrier(MPI_COMM_WORLD);
+            // MPI_Win_fence(0, wins_.last());
+
+            for (int ivar=0; ivar<NS+NY; ++ivar)
+            {
+                scalar *table_ptr;
+                wins_.resize(wins_.size()+1);
+                int disp_unit = sizeof(double);
+                
+                // Only rank 0 on a node actually allocates memory
+                MPI_Win_allocate_shared(noderank==0 ? singleTableSize_/1.0*disp_unit : 0, disp_unit, MPI_INFO_NULL, 
+                                        nodecomm, &table_ptr, &wins_.last());
+
+                if (noderank != 0) // node-slave
+                {
+                    MPI_Aint winsize;
+                    int windisp;
+                    // table on node-slave, let it points to the allocated address of table on node-master
+                    MPI_Win_shared_query(wins_.last(), 0, &winsize, &windisp, &table_ptr);
+                }
+
+                MPI_Win_fence(0, wins_.last());
+
+                if (noderank == 0)
+                {
+                    std::copy(value_temp[ivar].begin(), value_temp[ivar].end(), table_ptr);
+                    value_temp[ivar].clear();  //清空内部向量的元素
+                    std::vector<double>().swap(value_temp[ivar]);  //使用临时向量进行swap，收缩容量
+                }
+
+                Info<<"ivar = "<<ivar<<endl;
+                MPI_Win_fence(0, wins_.last());
+                tableValues_.append(table_ptr);
+                table_ptr = nullptr;
+            }
+                
+
+        }
+        else  // 1 core computing
+        {
+            std::vector<std::vector<double>> value_temp;
+
+            #include "readThermChemTables.H"
+
+            double *table_ptr;
+
+            // double *table_ptr = nullptr;
+            for (int ivar=0; ivar<NS+NY; ++ivar)
+            {
+                table_ptr = new double[singleTableSize_];
+
+                std::copy(value_temp[ivar].begin(), value_temp[ivar].end(), table_ptr);
+
+                tableValues_.append(table_ptr); // let tableValues_ take over that piece of memory
+
+                // table_ptr = nullptr; // de-pointer with the allocated memory
+            }
+        }
+
+        Info<<"end reading turbulence properties"<<endl;
+
+        //- findthe maxmum PV value from Ycmax_Tb3
+        if (this->scaledPV_)
+        {
+            cMaxAll_ = 1.0;
+        }
+        else
+        {
+            const scalar Ycmaxall{*std::max_element(tableValues_[NS-1],tableValues_[NS-1]+singleTableSize_)};   
+            cMaxAll_ = Ycmaxall;
+            Info<< "\nunscaled PV -- Ycmaxall = "<<Ycmaxall<< endl;
+        }
+        
+        if (this->scaledPV_)
+        {
+            // dYeq_Tb2 = std::vector<std::vector<double>>(2,std::vector<double> (NH*NZ*NGZ));
+            d2Yeq_Tb2 = { new double[NH*NZ*NGZ]{} }; d1Yeq_Tb2 = { new double[NH*NZ*NGZ]{} };
+            Info<< "\nReading non-premixed properties\n" << endl;   //屏幕提示：读取非预混的属性
+            count = 0;
+            for (int hh=0; hh<NH; ++hh)
+            {
+                for(int ii = 0; ii < NZ; ++ii)
+                {
+                    for(int jj = 0; jj < NGZ; ++jj)
+                    {
+                        if (std::getline(table, line))
+                        {
+                            std::istringstream iss(line);
+                            iss >> d2Yeq_Tb2[count] >> d1Yeq_Tb2[count];
+                            count++;
+                        }
+                    }
+                }
+            }
+        }
+        else
+        {
+            Info << "no need to Reading non-premixed properties"<<endl;
+        }
+
+
+        table.close();
+        Info<< "* * * * * * * * * * Complete FlaRe table * * * * * * * * * " << endl;
+    }// end if (table.is_open())
+    else
+    {
+        Info << "Error: Unable to open chemTab file!" << endl;
+    }
 
-  #include "readThermChemTables.H"
 
 }
 
@@ -65,7 +347,23 @@ H_ox("H_ox",dimensionSet(0,2,-2,0,0,0,0),Hox)
 
 tableSolver::~tableSolver()
 {
-
+    if (Pstream::parRun())
+    {
+        forAll(wins_, tableI)
+        {
+            MPI_Win_free(&wins_[tableI]);
+            // after free window object, we cannot delete variables located in this windows
+            // because the memory has already been released
+        }
+    }
+    else
+    {
+        forAll(tableValues_, tableI)
+        {
+            delete[] tableValues_[tableI];
+            tableValues_[tableI] = nullptr;
+        }
+    }
 }
 
 // * * * * * * * * * * * * * *  Member Functions * * * * * * * * * * * * * * //
@@ -109,17 +407,13 @@ double Foam::tableSolver::cal_gcor
 {
     double gcor;
 
-    if(abs(Zcvar) < 1e-12)
-    {
-        gcor=0.0;
-    }
-    else if(abs(cvar) < 1e-12 || abs(Zvar) < 1e-12)
+    if ( (cvar < 1.0e-4) or (Zvar < 1.0e-6) )
     {
         gcor=0.0;
     }
     else
     {
-        gcor = Zcvar/(Foam::sqrt(Zvar * cvar));
+        gcor = (Zcvar) / (Foam::sqrt(Zvar) * Foam::sqrt(cvar));
     }
 
     gcor = fmin(1.0,gcor);
@@ -137,6 +431,8 @@ int Foam::tableSolver::locate_lower
 {
         int i,loc;
 
+        if (n == 1) return 0;
+
         loc = 0;
         if(x <= array[loc]) return loc;
 
@@ -155,6 +451,35 @@ int Foam::tableSolver::locate_lower
         return loc;    
 }
 
+int Foam::tableSolver::locate_higher
+(
+    int n, 
+    double array[], 
+    double x
+)
+{
+        int i,loc;
+
+        if (n == 1) return 0;
+
+        loc = 0;
+        if(x >= array[loc]) return loc;
+
+        loc = n-2;
+        if(x <= array[loc+1]) return loc;
+
+        for(i=0; i<n-1; i++)
+        {
+            if(x <= array[i] && x > array[i+1])
+            {
+                loc = i;
+                return loc;
+            }
+        }
+
+        return loc;    
+}
+
 double Foam::tableSolver::intfac
 (
      double xx, 
@@ -166,11 +491,28 @@ double Foam::tableSolver::intfac
 
         if(xx <= low) fac = 0.0;
         else if(xx >= high) fac = 1.0;
+        else if((high-low) < 1e-30) fac = 0.0;
         else fac = (xx-low)/(high-low);
 
         return fac;    
 }
 
+double Foam::tableSolver::intfac_high
+(
+     double xx,
+     double high, 
+     double low
+)
+{
+        double fac;
+
+        if(xx >= high) fac = 0.0;
+        else if(xx <= low) fac = 1.0;
+        else fac = (high-xx)/(high-low);
+
+        return fac;    
+}
+
 double Foam::tableSolver::interp1d
 (
   int n1, 
@@ -254,19 +596,93 @@ double Foam::tableSolver::interp5d
         return result;    
 }
 
-double Foam::tableSolver::interp2d
+/*------------------------------------------------------------------------*\
+                        6D linear interpolation
+\*------------------------------------------------------------------------*/
+double Foam::tableSolver::interp6d
 (
-    int n1, int n2, int loc_z, int loc_gz,
-    double zfac,double gzfac, double table2d[]
+    int nh, int nz, int nc, int ngz, int ngc, int ngcz,
+    int loc_h, int loc_z, int loc_c, int loc_gz, int loc_gc, int loc_gcz,
+    double hfac, double zfac, double cfac, double gzfac, double gcfac,
+    double gczfac, double table_6d[]
 )
 {
-        int i1,i2,j1,j2,loc;
-        double factor,result;
+    int i0,i1,i2,i3,i4,i5,j0,j1,j2,j3,j4,j5,loc;
+    double factor, result;
+
+    result =0.0;
+    for (i0=0; i0<2; i0++)
+    {
+        if(i0 == 1 and nh>1) j0 = loc_h+1;
+        else j0 = loc_h;
 
-        result = 0.0;
         for(i1=0; i1<2; i1++)
         {
-            if(i1 == 1) j1 = loc_z+1;
+            if(i1 == 1 and nz>1) j1 = loc_z+1;
+            else j1 = loc_z;
+
+            for(i2=0; i2<2; i2++)
+            {
+                if(i2 == 1 and nc>1) j2 = loc_c+1;
+                else j2 = loc_c;
+
+                for(i3=0; i3<2; i3++)
+                {
+                    if(i3 == 1 and ngz>1) j3 = loc_gz+1;
+                    else j3 = loc_gz;
+
+                    for(i4=0; i4<2; i4++)
+                    {
+                        if(i4 == 1 and ngc>1) j4 = loc_gc+1;
+                        else j4 = loc_gc;
+
+                        for(i5=0; i5<2; i5++)
+                        {
+                            factor = (1.0-hfac+i0*(2.0*hfac-1.0))
+                                    *(1.0-zfac+i1*(2.0*zfac-1.0))
+                                    *(1.0-cfac+i2*(2.0*cfac-1.0))
+                                    *(1.0-gzfac+i3*(2.0*gzfac-1.0))
+                                    *(1.0-gcfac+i4*(2.0*gcfac-1.0))
+                                    *(1.0-gczfac+i5*(2.0*gczfac-1.0));
+
+                            if(i5 == 1 and ngcz>1) j5 = loc_gcz+1;
+                            else j5 = loc_gcz;
+
+                            loc = j0*nz*nc*ngz*ngc*ngcz
+                                +j1*nc*ngz*ngc*ngcz
+                                +j2*ngz*ngc*ngcz
+                                +j3*ngc*ngcz
+                                +j4*ngcz
+                                +j5;
+                            result = result + factor*table_6d[loc] ;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    return result;
+}
+
+double Foam::tableSolver::interp3d
+(
+    int n0, int n1, int n2, int loc_h, int loc_z, int loc_gz,
+    double hfac, double zfac,double gzfac, double table3d[]
+)
+{
+    int i0,i1,i2,j0,j1,j2,loc;
+    double factor,result;
+
+    result = 0.0;
+    for (i0=0; i0<2; i0++)
+    {
+        if(i0 == 1 and n0>1) j0 = loc_h+1;
+        else j0 = loc_h;
+
+        for(i1=0; i1<2; i1++)
+        {
+            if(i1 == 1 and n1>1) j1 = loc_z+1;
             else j1 = loc_z;
 
             for(i2=0; i2<2; i2++)
@@ -274,9 +690,41 @@ double Foam::tableSolver::interp2d
                 factor = (1.0-zfac+i1*(2.0*zfac-1.0))
                         *(1.0-gzfac+i2*(2.0*gzfac-1.0));
 
-                if(i2 == 1) j2 = loc_gz+1;
+                if(i2 == 1 and n2>1) j2 = loc_gz+1;
                 else j2 = loc_gz;
 
+                loc = j0*n1*n2+j1*n2+j2;
+                result = result + factor*table3d[loc];
+            }
+        }
+    }
+
+    return result;
+}
+
+double Foam::tableSolver::interp2d
+(
+    int n1, int n2, int loc_h, int loc_z,
+    double hfac,double zfac, double table2d[]
+)
+{
+        int i1,i2,j1,j2,loc;
+        double factor,result;
+
+        result = 0.0;
+        for(i1=0; i1<2; i1++)
+        {
+            if(i1 == 1 and n1>1) j1 = loc_h+1;
+            else j1 = loc_h;
+
+            for(i2=0; i2<2; i2++)
+            {
+                factor = (1.0-hfac+i1*(2.0*hfac-1.0))
+                        *(1.0-zfac+i2*(2.0*zfac-1.0));
+
+                if(i2 == 1 and n2>1) j2 = loc_z+1;
+                else j2 = loc_z;
+
                 loc = j1*n2+j2;
                 result = result + factor*table2d[loc];
             }
@@ -345,6 +793,25 @@ double Foam::tableSolver::lookup2d
         return interp2d(n1,n2,loc_1,loc_2,fac_1,fac_2,table_2d);    
 }
 
+
+double Foam::tableSolver::lookup3d
+(
+    int n1, double list_1[], double x1,
+    int n2, double list_2[], double x2,
+    int n3, double list_3[], double x3,
+    double table_3d[]    
+)
+{
+    int loc_1 = locate_lower(n1,list_1,x1);
+    double fac_1 = intfac(x1,list_1[loc_1],list_1[loc_1+1]);
+    int loc_2 = locate_lower(n2,list_2,x2);
+    double fac_2 = intfac(x2,list_2[loc_2],list_2[loc_2+1]);
+    int loc_3 = locate_lower(n3,list_3,x3);
+    double fac_3 = intfac(x3,list_3[loc_3],list_3[loc_3+1]);
+
+    return interp3d(n1,n2,n3,loc_1,loc_2,loc_3,fac_1,fac_2,fac_3,table_3d);
+}
+
 double Foam::tableSolver::lookup5d
 (
     int n1, double list_1[], double x1,
@@ -371,6 +838,35 @@ double Foam::tableSolver::lookup5d
 }
 
 
+double Foam::tableSolver::lookup6d
+(
+    int n1, double list_1[], double x1,
+    int n2, double list_2[], double x2,
+    int n3, double list_3[], double x3,
+    int n4, double list_4[], double x4,
+    int n5, double list_5[], double x5,
+    int n6, double list_6[], double x6,
+    double table_6d[]
+)
+{
+
+    int loc_1 = locate_lower(n1,list_1,x1);
+    double fac_1 = intfac(x1,list_1[loc_1],list_1[loc_1+1]);
+    int loc_2 = locate_lower(n2,list_2,x2);
+    double fac_2 = intfac(x2,list_2[loc_2],list_2[loc_2+1]);
+    int loc_3 = locate_lower(n3,list_3,x3);
+    double fac_3 = intfac(x3,list_3[loc_3],list_3[loc_3+1]);
+    int loc_4 = locate_lower(n4,list_4,x4);
+    double fac_4 = intfac(x4,list_4[loc_4],list_4[loc_4+1]);
+    int loc_5 = locate_lower(n5,list_5,x5);
+    double fac_5 = intfac(x5,list_5[loc_5],list_5[loc_5+1]);
+    int loc_6 = locate_lower(n6,list_6,x6);
+    double fac_6 = intfac(x6,list_6[loc_6],list_6[loc_6+1]);
+
+    return interp6d(n1,n2,n3,n4,n5,n6,loc_1,loc_2,loc_3,loc_4,loc_5,loc_6,
+        fac_1,fac_2,fac_3,fac_4,fac_5,fac_6,table_6d);
+}
+
 double Foam::tableSolver::RANSsdrFLRmodel
 (
     double cvar, double epsilon, double k, double nu,
@@ -379,12 +875,12 @@ double Foam::tableSolver::RANSsdrFLRmodel
 {
         double beta=6.7, C3, C4, Kc, Ka;
 
-        Ka=(dl/sl)/Foam::sqrt(nu/epsilon);
+        Ka=(dl/(sl+SMALL))/(Foam::sqrt(nu/epsilon)+SMALL);
         Kc=kc_s*tau;
         C3=1.5*Foam::sqrt(Ka)/(1+Foam::sqrt(Ka));
         C4=1.1/Foam::pow((1+Ka),0.4);
         
-        return rho/beta*cvar*( (2.0*Kc-tau*C4)*sl/dl + C3*epsilon/k );    
+        return rho/beta*cvar*( (2.0*Kc-tau*C4)*sl/(dl+SMALL) + C3*epsilon/(k+SMALL) );    
 }
 
 
diff --git a/src/dfCombustionModels/FGM/flameletTableSolver/tableSolver.H b/src/dfCombustionModels/FGM/flameletTableSolver/tableSolver.H
index 4133f9f8e..caf15a35e 100644
--- a/src/dfCombustionModels/FGM/flameletTableSolver/tableSolver.H
+++ b/src/dfCombustionModels/FGM/flameletTableSolver/tableSolver.H
@@ -39,9 +39,12 @@ SourceFiles
 #include "scalarList.H"
 #include "scalarField.H"
 #include "scalar.H"
+#include "volFields.H"
 #include "dimensionedScalar.H"
 #include "Switch.H"
 #include "PtrList.H"
+#include <mpi.h>
+#include <vector>
 
 namespace Foam
 {
@@ -49,40 +52,56 @@ namespace Foam
 class tableSolver
 {
  //- protected data
- public:
-   double small;
-   double smaller;
-   double smallest;
-   double T0;
+ protected:
+    double small;
+    double smaller;
+    double smallest;
+    double T0;
 
 
-  //- create table
-   FILE *table;
+    //- create table
+    FILE *table;
+    word tablePath_;
 
-  //- self-defined species name
-   hashedWordList speciesNames_;
+    //- list of MPI windows, used for shared memory
+    List<MPI_Win> wins_;
+
+    //- self-defined species name
+    hashedWordList speciesNames_table_;
+    hashedWordList spc_omegaNames_table_;
+
+    //- List of the used tables
+    hashedWordList tableNames_;
+
+    //- tables of all tabulated properties, by tableNames_'s order
+    List<scalar*> tableValues_;
+
+    double *d2Yeq_Tb2, *d1Yeq_Tb2;
+
+    //- size of one signle table
+    unsigned long long int singleTableSize_;
 
    //- enthalpy solve witch
-   Switch& scaledPV_;
+    Switch scaledPV_;
 
-   //- T-lookup switch
-   Switch flameletT_;
+    //- T-lookup switch
+    // Switch flameletT_;
 
-   //- unscaled max C
-   scalar& cMaxAll_;
-   //PtrList<scalarField> Yi_Tb3;
+    //- unscaled max C
+    double cMaxAll_;
+    //PtrList<scalarField> Yi_Tb3;
 
-   int NZL,NZ,NC,NGZ,NGC,NZC,NY,NS;
+    int NZL,NH,NZ,NC,NGZ,NGC,NZC,NY,NS,NYomega;
 
-   double Hox,Hfu,prog;
+    double Hox,Hfu,prog;
 
-   int fHox_fu, fNZL, fNZ, fNC, fNGZ, fNGC, fNZC, fNS_NY;
+    int fHox_fu, fNZL, fNH, fNZ, fNC, fNGZ, fNGC, fNZC, fNS_NY;
 
-   int f1, f2, f3, f4, f5, f6, f7, f8, f9, f10;
+    int f0, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10;
 
-   //- fuel and oxidizer enthalpy
-   const dimensionedScalar H_fuel;
-   const dimensionedScalar H_ox;
+    //- fuel and oxidizer enthalpy
+    const dimensionedScalar H_fuel;
+    const dimensionedScalar H_ox;
 
 
     double  *z_Tb5,     //- interpolation z_space
@@ -91,7 +110,8 @@ class tableSolver
             *tau_Tb5,   //- heat release parameter
             *kctau_Tb5; //- for calculating chi_c
 
-    double  *z_Tb3,     //- z space
+    double  *h_Tb3,     //- h loss space
+            *z_Tb3,     //- z space
             *c_Tb3,     //- c space
             *gz_Tb3,    //- z-variance space
             *gc_Tb3,    //- c-variance space
@@ -113,13 +133,13 @@ class tableSolver
            *Yi03_Tb3;    //- mass fraction of CO2
 
 
-    double *d2Yeq_Tb2;
+    // double *d2Yeq_Tb2;
 
 
   public:
 
     //- Constructor
-     tableSolver(wordList speciesNames, Switch& scaledPV_, Switch flameletT_, scalar& cMaxAll_);
+     tableSolver(word tablePath_);
 
     //- Member function 
 
@@ -132,6 +152,12 @@ class tableSolver
     //- Locate lower index for 1D linear interpolation
     int locate_lower(int n, double array[], double x); 
 
+    //- Locate higher index for 1D linear interpolation
+    int locate_higher(int n, double array[], double x); 
+
+    //- Compute 1D linear interpolation factor [for Hloss]
+    double intfac_high(double xx, double high, double low);
+
     //- Compute 1D linear interpolation factor
     double intfac(double xx, double low, double high);
 
@@ -144,6 +170,16 @@ class tableSolver
                     ,double zfac, double cfac, double gzfac, double gcfac
                     ,double gczfac, double table_5d[]);
 
+    //- 6D linear interpolation
+    double interp6d(int nh,int nz, int nc, int ngz, int ngc, int ngcz
+                    , int loc_h, int loc_z, int loc_c, int loc_gz, int loc_gc, int loc_gcz
+                    , double hfac, double zfac, double cfac, double gzfac, double gcfac
+                    ,double gczfac, double table_6d[]);
+
+    //- 3D linear interpolation
+    double interp3d(int n0, int n1, int n2, int loc_h, int loc_z, int loc_gz, double hfac,double zfac
+                    ,double gzfac, double table3d[]); 
+
    //- 2D linear interpolation
     double interp2d(int n1, int n2, int loc_z, int loc_gz, double zfac
                     ,double gzfac, double table2d[]); 
@@ -166,6 +202,9 @@ class tableSolver
     double lookup2d(int n1, double list_1[], double x1,int n2,
                     double list_2[], double x2, double table_2d[]);
 
+   //- 3D table look-up
+    double lookup3d(int n1, double list_1[], double x1,int n2,
+                    double list_2[], double x2, int n3, double list_3[], double x3, double table_3d[]);
 
    //- 5D table look-up
     double lookup5d(int n1, double list_1[], double x1,
@@ -175,6 +214,15 @@ class tableSolver
                     int n5, double list_5[], double x5,
                     double table_5d[]); 
 
+    //- 6D table look-up
+    double lookup6d(int n1, double list_1[], double x1,
+                    int n2, double list_2[], double x2,
+                    int n3, double list_3[], double x3,
+                    int n4, double list_4[], double x4,
+                    int n5, double list_5[], double x5,
+                    int n6, double list_6[], double x6,
+                    double table_5d[]); 
+
 
    //- RANS SDR model for progress variable
     double RANSsdrFLRmodel(double cvar, double epsilon, double k, double nu,
diff --git a/src/dfCombustionModels/FGM/flareFGM/flareFGM.C b/src/dfCombustionModels/FGM/flareFGM/flareFGM.C
index cafdf3357..95e136081 100644
--- a/src/dfCombustionModels/FGM/flareFGM/flareFGM.C
+++ b/src/dfCombustionModels/FGM/flareFGM/flareFGM.C
@@ -37,15 +37,18 @@ Foam::combustionModels::flareFGM<ReactionThermo>::flareFGM
 )
 :
     baseFGM<ReactionThermo>(modelType, thermo, turb, combustionProperties),
-    tableSolver(
-                 baseFGM<ReactionThermo>::speciesNames_,
-                 baseFGM<ReactionThermo>::scaledPV_,
-                 baseFGM<ReactionThermo>::flameletT_,
-                 baseFGM<ReactionThermo>::Ycmaxall_
-               )
+    tableSolver(baseFGM<ReactionThermo>::tablePath_)
 {
     //- retrieval data from table
     retrieval();
+
+    Info<< "At flareFGM, min/max(T) = " << min(this->T_).value() << ", " << max(this->T_).value() << endl;    
+
+    if (tableSolver::scaledPV_ != baseFGM<ReactionThermo>::scaledPV_)
+    {
+        Info << "Warning! -- scaledPV in FGM table: " << tableSolver::scaledPV_
+            << "are not equal to that in combustionProperties: " << baseFGM<ReactionThermo>::scaledPV_ << endl;
+    }
 }
 
 
@@ -64,10 +67,15 @@ void Foam::combustionModels::flareFGM<ReactionThermo>::correct()
     //- solve transport equation
     baseFGM<ReactionThermo>::transport();
 
+    if (this->isLES_)
+    {
+        baseFGM<ReactionThermo>::magUPrime();
+    }
+
     //update enthalpy using lookup data
     if(!(this->solveEnthalpy_))
     {
-        this->He_ = this->Z_*(H_fuel-H_ox) + H_ox;
+        this->He_ = this->Z_*(this->H_fuel-this->H_ox) + this->H_ox;
     }
   
     //- retrieval data from table
@@ -91,177 +99,252 @@ void Foam::combustionModels::flareFGM<ReactionThermo>::retrieval()
     volScalarField& mu = const_cast<volScalarField&>(tmu());
     scalarField& muCells = mu.primitiveFieldRef();     
 
+    tmp<volScalarField> tmut = this->turbulence().mut();  
+    volScalarField& mut = const_cast<volScalarField&>(tmut());
+    scalarField& mutCells = mut.primitiveFieldRef();   
 
-    //- calculate reacting flow solution
-    const scalar Zl{z_Tb5[0]};  
-    const scalar Zr{z_Tb5[NZL-1]};  
+    int ih = 0;
+    scalar Zl, Zr;
 
-    dimensionedScalar TMin("TMin",dimensionSet(0,0,0,1,0,0,0),200.0);    
-    dimensionedScalar TMax("TMax",dimensionSet(0,0,0,1,0,0,0),3000.0);  
+    // dimensionedScalar TMin("TMin",dimensionSet(0,0,0,1,0,0,0),270.0);    
+    // dimensionedScalar TMax("TMax",dimensionSet(0,0,0,1,0,0,0),5000.0);  
+    const double TMin = 270.0;
+    const double TMax = 5000.0;
 
     forAll(this->rho_, celli)  
     {
+        if(this->isLES_)
+        {
+             this->chi_ZCells_[celli] = this->sdrLRXmodel(2.0,mutCells[celli]
+                                   /this->rho_[celli],this->deltaCells_[celli],this->ZvarCells_[celli]); 
+             this->chi_ZcCells_[celli] = this->sdrLRXmodel(2.0,mutCells[celli]
+                                  /this->rho_[celli],this->deltaCells_[celli],this->ZcvarCells_[celli]); 
+        }     
+        else
+        {
+            this->chi_ZCells_[celli] = 1.0*epsilonCells[celli]/kCells[celli] *this->ZvarCells_[celli]; 
+
+            this->chi_ZcCells_[celli] = 1.0*epsilonCells[celli]/kCells[celli] *this->ZcvarCells_[celli];             
+        }
+
 
-        this->chi_ZCells_[celli] = 1.0*epsilonCells[celli]/kCells[celli] *this->ZvarCells_[celli]; 
+        double hLoss = ( this->ZCells_[celli]*(this->Hfu-this->Hox) + this->Hox ) - this->HCells_[celli];
+        hLoss = max(hLoss, this->h_Tb3[0]);
+        hLoss = min(hLoss, this->h_Tb3[this->NH - 1]);
 
-        this->chi_ZcCells_[celli] = 1.0*epsilonCells[celli]/kCells[celli] *this->ZcvarCells_[celli]; 
+        ih = this->locate_lower(this->NH,this->h_Tb3,hLoss);
+
+        Zl = this->z_Tb5[ih*this->NZL];  
+        Zr = this->z_Tb5[(ih+1)*this->NZL-1];
 
         if(this->ZCells_[celli] >= Zl && this->ZCells_[celli] <= Zr     
         && this->combustion_ && this->cCells_[celli] > this->small) 
         {
 
-            double kc_s = this->lookup1d(NZL,z_Tb5,this->ZCells_[celli],kctau_Tb5);      
-            double tau = this->lookup1d(NZL,z_Tb5,this->ZCells_[celli],tau_Tb5);    
-            double sl = this->lookup1d(NZL,z_Tb5,this->ZCells_[celli],sl_Tb5);     
-            double dl = this->lookup1d(NZL,z_Tb5,this->ZCells_[celli],th_Tb5);   
+            double kc_s = this->lookup2d(this->NH,this->h_Tb3,hLoss,this->NZL,this->z_Tb5,this->ZCells_[celli],this->kctau_Tb5);      
+            double tau = this->lookup2d(this->NH,this->h_Tb3,hLoss,this->NZL,this->z_Tb5,this->ZCells_[celli],this->tau_Tb5);    
+            double sl = this->lookup2d(this->NH,this->h_Tb3,hLoss,this->NZL,this->z_Tb5,this->ZCells_[celli],this->sl_Tb5);     
+            double dl = this->lookup2d(this->NH,this->h_Tb3,hLoss,this->NZL,this->z_Tb5,this->ZCells_[celli],this->th_Tb5);   
 
-            this->chi_cCells_[celli] =
-                this->RANSsdrFLRmodel(this->cvarCells_[celli],epsilonCells[celli],
-                    kCells[celli],muCells[celli]/this->rho_[celli],
-                    sl,dl,tau,kc_s,this->rho_[celli]);  
+            if(this->isLES_)
+            {
+                this->chi_cCells_[celli] = this->sdrFLRmodel(this->cvarCells_[celli],this->magUPrimeCells_[celli],
+                                            this->deltaCells_[celli],sl,dl,tau,kc_s,this->betacCells_[celli]);
+            } 
+            else
+            {
+                this->chi_cCells_[celli] =
+                    this->RANSsdrFLRmodel(this->cvarCells_[celli],epsilonCells[celli],
+                        kCells[celli],muCells[celli]/this->rho_[celli],
+                        sl,dl,tau,kc_s,this->rho_[celli]);  
+            }
 
         }
         else
         {
-            this->chi_cCells_[celli] = 1.0*epsilonCells[celli]/kCells[celli]*this->cvarCells_[celli]; 
+            if(this->isLES_)
+            {
+                this->chi_cCells_[celli] = this->sdrLRXmodel(2.0,mutCells[celli]
+                           /this->rho_[celli],this->deltaCells_[celli],this->cvarCells_[celli]);     
+            }
+            else
+            {
+                this->chi_cCells_[celli] = 1.0*epsilonCells[celli]/kCells[celli]*this->cvarCells_[celli]; 
+            }
         }
+  
 
-    // = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
-    
-    double gz{cal_gvar(this->ZCells_[celli],this->ZvarCells_[celli])};   
-    double gcz{cal_gcor(this->ZCells_[celli],this->cCells_[celli],this->ZvarCells_[celli],this->cvarCells_[celli],this->ZcvarCells_[celli])},
-            Ycmax{-1.0},cNorm{},gc{};    
-
-    if(scaledPV_)    
-    {
-        cNorm = this->cCells_[celli];  
-    }
-    else
-    {
-        Ycmax = this->lookup5d(NZ,z_Tb3,this->ZCells_[celli],
-                              NC,c_Tb3,0.0,
-                              NGZ,gz_Tb3,gz,
-                              NGC,gc_Tb3,0.0,
-                              NZC,gzc_Tb3,0.0,
-                              Ycmax_Tb3);    
-        Ycmax = max(this->smaller,Ycmax);    
-        cNorm = this->cCells_[celli]/Ycmax; 
-    }
+        // = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+        
+        double gz{this->cal_gvar(this->ZCells_[celli],this->ZvarCells_[celli])};   
+        double gcz{this->cal_gcor(this->ZCells_[celli],this->cCells_[celli],this->ZvarCells_[celli],this->cvarCells_[celli],this->ZcvarCells_[celli])},
+                Ycmax{-1.0},cNorm{},gc{};    
 
-    gc = cal_gvar(this->cCells_[celli],this->cvarCells_[celli],Ycmax);  
+        if(tableSolver::scaledPV_)    
+        {
+            cNorm = this->cCells_[celli];  
+        }
+        else
+        {
+            Ycmax = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                this->NZ,this->z_Tb3,this->ZCells_[celli],
+                                this->NC,this->c_Tb3,0.0,
+                                this->NGZ,this->gz_Tb3,gz,
+                                this->NGC,this->gc_Tb3,0.0,
+                                this->NZC,this->gzc_Tb3,0.0,
+                                this->tableValues_[this->NS-1]);    
+            Ycmax = max(this->smaller,Ycmax);    
+            cNorm = this->cCells_[celli]/Ycmax; 
+        }
 
-    this->WtCells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli],    
-                                   NC,c_Tb3,cNorm,
-                                   NGZ,gz_Tb3,gz,
-                                   NGC,gc_Tb3,gc,
-                                   NZC,gzc_Tb3,gcz,
-                                   mwt_Tb3);      
+        gc = this->cal_gvar(this->cCells_[celli],this->cvarCells_[celli],Ycmax);  
 
-    muCells[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli],
-                                   NC,c_Tb3,cNorm,
-                                   NGZ,gz_Tb3,gz,
-                                   NGC,gc_Tb3,gc,
-                                   NZC,gzc_Tb3,gcz,
-                                   nu_Tb3)*this->rho_[celli];   
+        if(this->ZCells_[celli] >= Zl && this->ZCells_[celli] <= Zr
+        && this->combustion_ && this->cCells_[celli] > this->small)  
+        {
+            if (this->isLES_)
+            {
+                this->omega_cCells_[celli] =
+                    this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                this->NZ,this->z_Tb3,this->ZCells_[celli],
+                                this->NC,this->c_Tb3,cNorm,
+                                this->NGZ,this->gz_Tb3,gz,
+                                this->NGC,this->gc_Tb3,gc,
+                                this->NZC,this->gzc_Tb3,gcz,
+                                this->tableValues_[0])
+                    + (
+                        tableSolver::scaledPV_
+                        ? (this->chi_ZCells_[celli] + this->chi_ZfltdCells_[celli])*this->cCells_[celli]
+                            *this->lookup3d(this->NH,this->h_Tb3,hLoss,this->NZ,this->z_Tb3,this->ZCells_[celli],
+                                        this->NGZ,this->gz_Tb3,gz,this->d2Yeq_Tb2)
+                            + 2*this->chi_ZcCells_[celli]*this->lookup3d(this->NH,this->h_Tb3,hLoss,
+                                                    this->NZ,this->z_Tb3,this->ZCells_[celli],
+                                                    this->NGZ,this->gz_Tb3,gz,this->d1Yeq_Tb2)
+                        : 0.0
+                    );   
+            }
+            else
+            {
+                this->omega_cCells_[celli] =
+                    this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                this->NZ,this->z_Tb3,this->ZCells_[celli],
+                                this->NC,this->c_Tb3,cNorm,
+                                this->NGZ,this->gz_Tb3,gz,
+                                this->NGC,this->gc_Tb3,gc,
+                                this->NZC,this->gzc_Tb3,gcz,
+                                this->tableValues_[0])
+                    + (
+                        tableSolver::scaledPV_
+                        ? (  this->chi_ZCells_[celli]*this->cCells_[celli]
+                            *this->lookup3d(this->NH,this->h_Tb3,hLoss,this->NZ,this->z_Tb3,this->ZCells_[celli],
+                                        this->NGZ,this->gz_Tb3,gz,this->d2Yeq_Tb2)
+                            + 2*this->chi_ZcCells_[celli]*this->lookup3d(this->NH,this->h_Tb3,hLoss,
+                                                        this->NZ,this->z_Tb3,this->ZCells_[celli],
+                                                        this->NGZ,this->gz_Tb3,gz,this->d1Yeq_Tb2)  )
+                        : 0.0
+                    );   
+            }
 
-   // -------------------- Yis begin ------------------------------
-    if(NY > 0)
-    {
-        this->YH2OCells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli],
-                                        NC,c_Tb3,cNorm,
-                                        NGZ,gz_Tb3,gz,
-                                        NGC,gc_Tb3,gc,
-                                        NZC,gzc_Tb3,gcz,
-                                        Yi01_Tb3);   
-        this->YCOCells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli],
-                                       NC,c_Tb3,cNorm,
-                                       NGZ,gz_Tb3,gz,
-                                       NGC,gc_Tb3,gc,
-                                       NZC,gzc_Tb3,gcz,
-                                       Yi02_Tb3);   
-        this->YCO2Cells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli],
-                                       NC,c_Tb3,cNorm,
-                                       NGZ,gz_Tb3,gz,
-                                       NGC,gc_Tb3,gc,
-                                       NZC,gzc_Tb3,gcz,
-                                       Yi03_Tb3);   
-    }
 
-    // -------------------- Yis end ------------------------------
+            this->cOmega_cCells_[celli] = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                                this->NZ,this->z_Tb3,this->ZCells_[celli],
+                                                this->NC,this->c_Tb3,cNorm,
+                                                this->NGZ,this->gz_Tb3,gz,
+                                                this->NGC,this->gc_Tb3,gc,
+                                                this->NZC,this->gzc_Tb3,gcz,
+                                                this->tableValues_[1]);    
 
-    if(this->ZCells_[celli] >= Zl && this->ZCells_[celli] <= Zr
-       && this->combustion_ && this->cCells_[celli] > this->small)  
-    {
-        this->omega_cCells_[celli] =
-            lookup5d(NZ,z_Tb3,this->ZCells_[celli],
-                          NC,c_Tb3,cNorm,
-                          NGZ,gz_Tb3,gz,
-                          NGC,gc_Tb3,gc,
-                          NZC,gzc_Tb3,gcz,
-                          omgc_Tb3)
-            + (
-                  scaledPV_
-                  ? this->chi_ZCells_[celli]*this->cCells_[celli]
-                    *lookup2d(NZ,z_Tb3,this->ZCells_[celli],
-                                   NGZ,gz_Tb3,gz,d2Yeq_Tb2)
-                  : 0.0
-              );   
-
-         this->cOmega_cCells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli],
-                                             NC,c_Tb3,cNorm,
-                                             NGZ,gz_Tb3,gz,
-                                             NGC,gc_Tb3,gc,
-                                             NZC,gzc_Tb3,gcz,
-                                             cOc_Tb3);    
-
-         this->ZOmega_cCells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli],
-                                             NC,c_Tb3,cNorm,
-                                             NGZ,gz_Tb3,gz,
-                                             NGC,gc_Tb3,gc,
-                                             NZC,gzc_Tb3,gcz,
-                                             ZOc_Tb3);    
+            this->ZOmega_cCells_[celli] = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                                this->NZ,this->z_Tb3,this->ZCells_[celli],
+                                                this->NC,this->c_Tb3,cNorm,
+                                                this->NGZ,this->gz_Tb3,gz,
+                                                this->NGC,this->gc_Tb3,gc,
+                                                this->NZC,this->gzc_Tb3,gcz,
+                                                this->tableValues_[2]);    
 
-    }
-    else   
-    {
-        this->omega_cCells_[celli] = 0.0;
-        this->cOmega_cCells_[celli] = 0.0;  
-        this->ZOmega_cCells_[celli] = 0.0; 
-    }
+        }
+        else   
+        {
+            this->omega_cCells_[celli] = 0.0;
+            this->cOmega_cCells_[celli] = 0.0;  
+            this->ZOmega_cCells_[celli] = 0.0; 
+        }
+        this->omega_cCells_[celli] = this->omega_cCells_[celli]*this->rho_[celli];   
+        this->cOmega_cCells_[celli] = this->cOmega_cCells_[celli]*this->rho_[celli];
+        this->ZOmega_cCells_[celli] = this->ZOmega_cCells_[celli]*this->rho_[celli];
+
+        this->He_s_[celli] = hLoss;
+        this->c_s_[celli] = cNorm;
+        this->Zvar_s_[celli] = gz;
+        this->cvar_s_[celli] = gc;
+        this->Zcvar_s_[celli] = gcz;
+
+        this->WtCells_[celli] = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                    this->NZ,this->z_Tb3,this->ZCells_[celli],    
+                                    this->NC,this->c_Tb3,cNorm,
+                                    this->NGZ,this->gz_Tb3,gz,
+                                    this->NGC,this->gc_Tb3,gc,
+                                    this->NZC,this->gzc_Tb3,gcz,
+                                    this->tableValues_[4]);      
+
+        muCells[celli] = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                    this->NZ,this->z_Tb3,this->ZCells_[celli],
+                                    this->NC,this->c_Tb3,cNorm,
+                                    this->NGZ,this->gz_Tb3,gz,
+                                    this->NGC,this->gc_Tb3,gc,
+                                    this->NZC,this->gzc_Tb3,gcz,
+                                    this->tableValues_[7])*this->rho_[celli];   
+
+        // // -------------------- Yis begin ------------------------------
+        for (int yi=0; yi<this->NY; yi++)
+        {
+            word specieName2Update = this->speciesNames_table_[yi];
+            const label& specieLabel2Update = this->chemistryPtr_->species()[specieName2Update];
+            this->Y_[specieLabel2Update].primitiveFieldRef()[celli] = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                        this->NZ,this->z_Tb3,this->ZCells_[celli],
+                                        this->NC,this->c_Tb3,cNorm,
+                                        this->NGZ,this->gz_Tb3,gz,
+                                        this->NGC,this->gc_Tb3,gc,
+                                        this->NZC,this->gzc_Tb3,gcz,
+                                        this->tableValues_[this->NS+yi]);  
+        }
 
-    if(flameletT_)   
-    {
-        this->TCells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli],
-                                      NC,c_Tb3,cNorm,
-                                      NGZ,gz_Tb3,gz,
-                                      NGC,gc_Tb3,gc,
-                                      NZC,gzc_Tb3,gcz,
-                                      Tf_Tb3);   
-    }
-    else
-    {
-        this->CpCells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli],
-                                       NC,c_Tb3,cNorm,
-                                       NGZ,gz_Tb3,gz,
-                                       NGC,gc_Tb3,gc,
-                                       NZC,gzc_Tb3,gcz,
-                                       cp_Tb3);   
-
-         this->HfCells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli],
-                                       NC,c_Tb3,cNorm,
-                                       NGZ,gz_Tb3,gz,
-                                       NGC,gc_Tb3,gc,
-                                       NZC,gzc_Tb3,gcz,
-                                       hiyi_Tb3);   
-
-        this->TCells_[celli] = (this->HCells_[celli]-this->HfCells_[celli])/this->CpCells_[celli]
-                        + this->T0;   
-    }
+        // -------------------- Yis end ------------------------------
 
-    this->omega_cCells_[celli] = this->omega_cCells_[celli]*this->rho_[celli];   
-    this->cOmega_cCells_[celli] = this->cOmega_cCells_[celli]*this->rho_[celli];
-    this->ZOmega_cCells_[celli] = this->ZOmega_cCells_[celli]*this->rho_[celli];
+        if(baseFGM<ReactionThermo>::flameletT_)   
+        {
+            this->TCells_[celli] = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                        this->NZ,this->z_Tb3,this->ZCells_[celli],
+                                        this->NC,this->c_Tb3,cNorm,
+                                        this->NGZ,this->gz_Tb3,gz,
+                                        this->NGC,this->gc_Tb3,gc,
+                                        this->NZC,this->gzc_Tb3,gcz,
+                                        this->tableValues_[6]);   
+        }
+        else
+        {
+            this->CpCells_[celli] = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                        this->NZ,this->z_Tb3,this->ZCells_[celli],
+                                        this->NC,this->c_Tb3,cNorm,
+                                        this->NGZ,this->gz_Tb3,gz,
+                                        this->NGC,this->gc_Tb3,gc,
+                                        this->NZC,this->gzc_Tb3,gcz,
+                                        this->tableValues_[3]);   
+
+            this->HfCells_[celli] = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                        this->NZ,this->z_Tb3,this->ZCells_[celli],
+                                        this->NC,this->c_Tb3,cNorm,
+                                        this->NGZ,this->gz_Tb3,gz,
+                                        this->NGC,this->gc_Tb3,gc,
+                                        this->NZC,this->gzc_Tb3,gcz,
+                                        this->tableValues_[5]);   
+
+            this->TCells_[celli] = (this->HCells_[celli]-this->HfCells_[celli])/this->CpCells_[celli]
+                            + this->T0;   
+        }
 
+        this->TCells_[celli] = max(this->TCells_[celli], TMin);
+        this->TCells_[celli] = min(this->TCells_[celli], TMax);
     }
 
 
@@ -290,104 +373,170 @@ void Foam::combustionModels::flareFGM<ReactionThermo>::retrieval()
         tmp<scalarField> tmuw = this->turbulence().mu(patchi);
         scalarField& pmu = const_cast<scalarField&>(tmuw());    
 
+        const tmp<scalarField> tmutw = this->turbulence().mut(patchi);      
+        const scalarField& pmut = tmutw();  
+
+        fvPatchScalarField& pmagUPrime = this->magUPrime_.boundaryFieldRef()[patchi];    
+        fvPatchScalarField& pchi_Zfltd = this->chi_Zfltd_.boundaryFieldRef()[patchi];    
+
         fvPatchScalarField& pk = k.boundaryFieldRef()[patchi];
         fvPatchScalarField& pepsilon = epsilon.boundaryFieldRef()[patchi];
 
 
         forAll(prho_, facei)   
         {
+            label cellID = this->mesh().boundary()[patchi].faceCells()[facei]; 
 
-            pchi_Z[facei] = 1.0*pepsilon[facei]/pk[facei] *pZvar[facei]; 
+            double hLoss = ( pZ[facei]*(this->Hfu-this->Hox) + this->Hox ) - pH[facei];
+            hLoss = max(hLoss, this->h_Tb3[0]);
+            hLoss = min(hLoss, this->h_Tb3[this->NH - 1]);
+
+            ih = this->locate_lower(this->NH,this->h_Tb3,hLoss);
+
+            Zl = this->z_Tb5[ih*this->NZL];  
+            Zr = this->z_Tb5[(ih+1)*this->NZL-1];
+
+            if(this->isLES_)
+            {
+
+                pchi_Z[facei] = sdrLRXmodel(2.0,pmut[facei]
+                    /this->rho_[facei],this->deltaCells_[cellID],pZvar[facei]);
+                pchi_Zc[facei]= sdrLRXmodel(2.0,pmut[facei]
+                    /this->rho_[facei],this->deltaCells_[cellID],pZcvar[facei]); 
+            }
+            else
+            {
+                pchi_Z[facei] = 1.0*pepsilon[facei]/pk[facei] *pZvar[facei]; 
+
+                pchi_Zc[facei]= 1.0*pepsilon[facei]/pk[facei] *pZcvar[facei]; 
+            }
 
-            pchi_Zc[facei]= 1.0*pepsilon[facei]/pk[facei] *pZcvar[facei]; 
 
             if(pZ[facei] >= Zl && pZ[facei] <= Zr
             && this->combustion_ && pc[facei] > this->small) 
             {
 
-                double kc_s = lookup1d(NZL,z_Tb5,pZ[facei],kctau_Tb5);     
-                double tau = lookup1d(NZL,z_Tb5,pZ[facei],tau_Tb5);    
-                double sl = lookup1d(NZL,z_Tb5,pZ[facei],sl_Tb5);      
-                double dl = lookup1d(NZL,z_Tb5,pZ[facei],th_Tb5);      
-
-                pchi_c[facei] =
-                    RANSsdrFLRmodel(pcvar[facei],pepsilon[facei],
-                        pk[facei],pmu[facei]/prho_[facei],
-                        sl,dl,tau,kc_s,prho_[facei]);
+                double kc_s = this->lookup2d(this->NH,this->h_Tb3,hLoss,this->NZL,this->z_Tb5,pZ[facei],this->kctau_Tb5);     
+                double tau = this->lookup2d(this->NH,this->h_Tb3,hLoss,this->NZL,this->z_Tb5,pZ[facei],this->tau_Tb5);    
+                double sl = this->lookup2d(this->NH,this->h_Tb3,hLoss,this->NZL,this->z_Tb5,pZ[facei],this->sl_Tb5);      
+                double dl = this->lookup2d(this->NH,this->h_Tb3,hLoss,this->NZL,this->z_Tb5,pZ[facei],this->th_Tb5);      
+
+                if(this->isLES_)
+                {
+                    pchi_c[facei] =  sdrFLRmodel(pcvar[facei],pmagUPrime[facei],
+                                      this->deltaCells_[cellID],sl,dl,tau,kc_s,this->betacCells_[cellID]);   
+                }
+                else
+                {
+                    pchi_c[facei] =
+                        this->RANSsdrFLRmodel(pcvar[facei],pepsilon[facei],
+                            pk[facei],pmu[facei]/prho_[facei],
+                            sl,dl,tau,kc_s,prho_[facei]);
+                }
             }
             else
             {      
-                pchi_c[facei] = 1.0*pepsilon[facei]/pk[facei] *pcvar[facei]; 
+                if(this->isLES_)
+                {
+                    pchi_c[facei] = sdrLRXmodel(2.0,pmut[facei]
+                        /this->rho_[facei],this->deltaCells_[cellID],pcvar[facei]);                    
+                }
+                else
+                {
+                    pchi_c[facei] = 1.0*pepsilon[facei]/(pk[facei]+SMALL)*pcvar[facei];    
+                }
 
             }
 
+
          // = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
 
-            double gz{cal_gvar(pZ[facei],pZvar[facei])};  
-            double gcz{cal_gcor(pZ[facei],pc[facei],pZvar[facei],pcvar[facei],pZcvar[facei])},
+            double gz{this->cal_gvar(pZ[facei],pZvar[facei])};  
+            double gcz{this->cal_gcor(pZ[facei],pc[facei],pZvar[facei],pcvar[facei],pZcvar[facei])},
                     Ycmax{-1.0},cNorm{},gc{};     
 
-            if(scaledPV_)
+            if(tableSolver::scaledPV_)
             {
                 cNorm = pc[facei];   
             }
             else
             {
-                Ycmax = lookup5d(NZ,z_Tb3,pZ[facei],
-                                        NC,c_Tb3,0.0,
-                                        NGZ,gz_Tb3,gz,
-                                        NGC,gc_Tb3,0.0,
-                                        NZC,gzc_Tb3,0.0,
-                                        Ycmax_Tb3);    
+                Ycmax = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                        this->NZ,this->z_Tb3,pZ[facei],
+                                        this->NC,this->c_Tb3,0.0,
+                                        this->NGZ,this->gz_Tb3,gz,
+                                        this->NGC,this->gc_Tb3,0.0,
+                                        this->NZC,this->gzc_Tb3,0.0,
+                                        this->tableValues_[this->NS-1]);    
                 Ycmax = max(this->smaller,Ycmax);  
                 cNorm = pc[facei]/Ycmax;   
             }
 
-            gc = cal_gvar(pc[facei],pcvar[facei],Ycmax);   
-
-            pWt[facei] = lookup5d(NZ,z_Tb3,pZ[facei],
-                                        NC,c_Tb3,cNorm,
-                                        NGZ,gz_Tb3,gz,
-                                        NGC,gc_Tb3,gc,
-                                        NZC,gzc_Tb3,gcz,
-                                        mwt_Tb3);    
+            gc = this->cal_gvar(pc[facei],pcvar[facei],Ycmax);   
 
-            pmu[facei] = lookup5d(NZ,z_Tb3,pZ[facei],
-                                        NC,c_Tb3,cNorm,
-                                        NGZ,gz_Tb3,gz,
-                                        NGC,gc_Tb3,gc,
-                                        NZC,gzc_Tb3,gcz,
-                                        nu_Tb3)*prho_[facei];  
+            this->He_s_.boundaryFieldRef()[patchi][facei] = hLoss;
+            this->c_s_.boundaryFieldRef()[patchi][facei] = cNorm;
+            this->Zvar_s_.boundaryFieldRef()[patchi][facei] = gz;
+            this->cvar_s_.boundaryFieldRef()[patchi][facei] = gc;
+            this->Zcvar_s_.boundaryFieldRef()[patchi][facei] = gcz;
 
             if(pZ[facei] >= Zl && pZ[facei] <= Zr
                 && this->combustion_ && pc[facei] > this->small) 
             {
-                pomega_c[facei] =
-                        lookup5d(NZ,z_Tb3,pZ[facei],
-                                    NC,c_Tb3,cNorm,
-                                    NGZ,gz_Tb3,gz,
-                                    NGC,gc_Tb3,gc,
-                                    NZC,gzc_Tb3,gcz,
-                                    omgc_Tb3)   
-                    + (
-                            scaledPV_
-                            ? pchi_Z[facei]*pc[facei]
-                            *lookup2d(NZ,z_Tb3,pZ[facei],
-                                            NGZ,gz_Tb3,gz,d2Yeq_Tb2)
-                            : 0.0
-                        );  
-
-                pcOmega_c[facei] = lookup5d(NZ,z_Tb3,pZ[facei],
-                                                    NC,c_Tb3,cNorm,
-                                                    NGZ,gz_Tb3,gz,
-                                                    NGC,gc_Tb3,gc,
-                                                    NZC,gzc_Tb3,gcz,cOc_Tb3);
-
-                pZOmega_c[facei] = lookup5d(NZ,z_Tb3,pZ[facei],
-                                            NC,c_Tb3,cNorm,
-                                            NGZ,gz_Tb3,gz,
-                                            NGC,gc_Tb3,gc,
-                                            NZC,gzc_Tb3,gcz,ZOc_Tb3);
+                if (this->isLES_)
+                {
+                    pomega_c[facei] =
+                            this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                        this->NZ,this->z_Tb3,pZ[facei],
+                                        this->NC,this->c_Tb3,cNorm,
+                                        this->NGZ,this->gz_Tb3,gz,
+                                        this->NGC,this->gc_Tb3,gc,
+                                        this->NZC,this->gzc_Tb3,gcz,
+                                        this->tableValues_[0])    
+                        + (
+                                tableSolver::scaledPV_
+                                ? (pchi_Z[facei]+ pchi_Zfltd[facei])*pc[facei]
+                                    *this->lookup3d(this->NH,this->h_Tb3,hLoss,this->NZ,this->z_Tb3,pZ[facei],
+                                                    this->NGZ,this->gz_Tb3,gz,this->d2Yeq_Tb2)
+                                    + 2.0*pchi_Zc[facei]*this->lookup3d(this->NH,this->h_Tb3,hLoss,this->NZ,this->z_Tb3,pZ[facei],
+                                    this->NGZ,this->gz_Tb3,gz,this->d1Yeq_Tb2)
+                                : 0.0
+                            );  
+                }
+                else
+                {
+                    pomega_c[facei] =
+                            this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                        this->NZ,this->z_Tb3,pZ[facei],
+                                        this->NC,this->c_Tb3,cNorm,
+                                        this->NGZ,this->gz_Tb3,gz,
+                                        this->NGC,this->gc_Tb3,gc,
+                                        this->NZC,this->gzc_Tb3,gcz,
+                                        this->tableValues_[0])   
+                        + (
+                                tableSolver::scaledPV_
+                                ? ( pchi_Z[facei]*pc[facei]
+                                *this->lookup3d(this->NH,this->h_Tb3,hLoss,this->NZ,this->z_Tb3,pZ[facei],
+                                                this->NGZ,this->gz_Tb3,gz,this->d2Yeq_Tb2)
+                                + 2.0*pchi_Zc[facei]*this->lookup3d(this->NH,this->h_Tb3,hLoss,this->NZ,this->z_Tb3,pZ[facei],
+                                    this->NGZ,this->gz_Tb3,gz,this->d1Yeq_Tb2)  )
+                                : 0.0
+                            );  
+                }
+
+                pcOmega_c[facei] = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                                    this->NZ,this->z_Tb3,pZ[facei],
+                                                    this->NC,this->c_Tb3,cNorm,
+                                                    this->NGZ,this->gz_Tb3,gz,
+                                                    this->NGC,this->gc_Tb3,gc,
+                                                    this->NZC,this->gzc_Tb3,gcz,this->tableValues_[1]);
+
+                pZOmega_c[facei] = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                            this->NZ,this->z_Tb3,pZ[facei],
+                                            this->NC,this->c_Tb3,cNorm,
+                                            this->NGZ,this->gz_Tb3,gz,
+                                            this->NGC,this->gc_Tb3,gc,
+                                            this->NZC,this->gzc_Tb3,gcz,this->tableValues_[2]);
             }    
             else
             {
@@ -395,46 +544,109 @@ void Foam::combustionModels::flareFGM<ReactionThermo>::retrieval()
                 pcOmega_c[facei] = 0.0;
                 pZOmega_c[facei] = 0.0;
             }
+            pomega_c[facei] = pomega_c[facei]*prho_[facei];    
+            pcOmega_c[facei] = pcOmega_c[facei]*prho_[facei];
+            pZOmega_c[facei] = pZOmega_c[facei]*prho_[facei];
 
-            if(flameletT_)  
+            pWt[facei] = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                        this->NZ,this->z_Tb3,pZ[facei],
+                                        this->NC,this->c_Tb3,cNorm,
+                                        this->NGZ,this->gz_Tb3,gz,
+                                        this->NGC,this->gc_Tb3,gc,
+                                        this->NZC,this->gzc_Tb3,gcz,
+                                        this->tableValues_[4]);    
+
+            pmu[facei] = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                        this->NZ,this->z_Tb3,pZ[facei],
+                                        this->NC,this->c_Tb3,cNorm,
+                                        this->NGZ,this->gz_Tb3,gz,
+                                        this->NGC,this->gc_Tb3,gc,
+                                        this->NZC,this->gzc_Tb3,gcz,
+                                        this->tableValues_[7])*prho_[facei];  
+
+            if(baseFGM<ReactionThermo>::flameletT_)  
             {
-                pT[facei] = lookup5d(NZ,z_Tb3,pZ[facei],
-                                        NC,c_Tb3,cNorm,
-                                        NGZ,gz_Tb3,gz,
-                                        NGC,gc_Tb3,gc,
-                                        NZC,gzc_Tb3,gcz,
-                                        Tf_Tb3);   
+                pT[facei] = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                        this->NZ,this->z_Tb3,pZ[facei],
+                                        this->NC,this->c_Tb3,cNorm,
+                                        this->NGZ,this->gz_Tb3,gz,
+                                        this->NGC,this->gc_Tb3,gc,
+                                        this->NZC,this->gzc_Tb3,gcz,
+                                        this->tableValues_[6]);   
             }
             else
             {
-                pCp[facei] = lookup5d(NZ,z_Tb3,pZ[facei],
-                                        NC,c_Tb3,cNorm,
-                                        NGZ,gz_Tb3,gz,
-                                        NGC,gc_Tb3,gc,
-                                        NZC,gzc_Tb3,gcz,
-                                        cp_Tb3);   
-
-                pHf[facei] = lookup5d(NZ,z_Tb3,pZ[facei],
-                                        NC,c_Tb3,cNorm,
-                                        NGZ,gz_Tb3,gz,
-                                        NGC,gc_Tb3,gc,
-                                        NZC,gzc_Tb3,gcz,
-                                        hiyi_Tb3);  
+                pCp[facei] = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                        this->NZ,this->z_Tb3,pZ[facei],
+                                        this->NC,this->c_Tb3,cNorm,
+                                        this->NGZ,this->gz_Tb3,gz,
+                                        this->NGC,this->gc_Tb3,gc,
+                                        this->NZC,this->gzc_Tb3,gcz,
+                                        this->tableValues_[3]);   
+
+                pHf[facei] = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                        this->NZ,this->z_Tb3,pZ[facei],
+                                        this->NC,this->c_Tb3,cNorm,
+                                        this->NGZ,this->gz_Tb3,gz,
+                                        this->NGC,this->gc_Tb3,gc,
+                                        this->NZC,this->gzc_Tb3,gcz,
+                                        this->tableValues_[5]);  
 
                 pT[facei] = (pH[facei]-pHf[facei])/pCp[facei]
                             + this->T0;  
             }
 
-            pomega_c[facei] = pomega_c[facei]*prho_[facei];    
-            pcOmega_c[facei] = pcOmega_c[facei]*prho_[facei];
-            pZOmega_c[facei] = pZOmega_c[facei]*prho_[facei];
+            pT[facei] = max(pT[facei], TMin);
+            pT[facei] = min(pT[facei], TMax);
 
-        } 
 
+            // // -------------------- Yis begin ------------------------------
+
+            for (int yi=0; yi<this->NY; yi++)
+            {
+                word specieName2Update = this->speciesNames_table_[yi];
+                const label& specieLabel2Update = this->chemistryPtr_->species()[specieName2Update];
+                this->Y_[specieLabel2Update].boundaryFieldRef()[patchi][facei]  = this->lookup6d(this->NH,this->h_Tb3,hLoss,
+                                            this->NZ,this->z_Tb3,pZ[facei],
+                                            this->NC,this->c_Tb3,cNorm,
+                                            this->NGZ,this->gz_Tb3,gz,
+                                            this->NGC,this->gc_Tb3,gc,
+                                            this->NZC,this->gzc_Tb3,gcz,
+                                            this->tableValues_[this->NS+yi]);  
+            }
+
+            // -------------------- Yis end ------------------------------
+        } 
     }
 
-    this->T_.max(TMin);  
-    this->T_.min(TMax);  
+    // this->T_.max(TMin);  
+    // this->T_.min(TMax);  
+
+    double R_uniGas_ = 8.314e3;
+    double p_operateDim_ = this->coeffs().lookupOrDefault("p_operateDim", this->incompPref_);
+
+    forAll(this->rho_.boundaryFieldRef(), patchi)   
+    {
+        fvPatchScalarField& ppsi_ = this->psi_.boundaryFieldRef()[patchi];
+        fvPatchScalarField& prho_ = this->rho_.boundaryFieldRef()[patchi];
+        fvPatchScalarField& pWt = this->Wt_.boundaryFieldRef()[patchi];
+        fvPatchScalarField& pT = this->T_.boundaryFieldRef()[patchi];
+        fvPatchScalarField pp_ = this->p_.boundaryField()[patchi];
+
+        forAll(prho_, facei)   
+        {
+            ppsi_[facei] = pWt[facei] / (R_uniGas_*pT[facei]);
+
+            if(this->incompPref_ > 0.0) 
+            {  
+                prho_[facei] = p_operateDim_*ppsi_[facei]; 
+            }
+            else 
+            {
+                prho_[facei] = pp_[facei]*ppsi_[facei];
+            }
+        }
+    }
 
     dimensionedScalar R_uniGas("R_uniGas",dimensionSet(1,2,-2,-1,-1,0,0),8.314e3);
     this->psi_ = this->Wt_/(R_uniGas*this->T_);
diff --git a/src/dfCombustionModels/Make/files b/src/dfCombustionModels/Make/files
index 50951f153..fb57b6ac0 100644
--- a/src/dfCombustionModels/Make/files
+++ b/src/dfCombustionModels/Make/files
@@ -1,6 +1,8 @@
 combustionModel/combustionModel.C
 CombustionModel/CombustionModel/CombustionModels.C
 
+FGM/FGMLiquidEvaporationBoil/FGMLiquidEvaporationBoils.C
+
 /*diffusion/diffusions.C*/
 /*infinitelyFastChemistry/infinitelyFastChemistrys.C*/
 
@@ -10,9 +12,11 @@ EDC/EDCs.C
 FGM/flameletTableSolver/tableSolver.C
 FGM/baseFGM/baseFGMs.C
 FGM/flareFGM/flareFGMs.C
+
 #if defined USE_PYTORCH
 FGM/DeePFGM/DeePFGMs.C
 #endif
+
 laminar/laminars.C
 
 noCombustion/noCombustions.C
diff --git a/src/dfCombustionModels/Make/options b/src/dfCombustionModels/Make/options
index d909e6048..bd7e2982e 100644
--- a/src/dfCombustionModels/Make/options
+++ b/src/dfCombustionModels/Make/options
@@ -11,8 +11,20 @@ EXE_INC = -std=c++14 \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(DF_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(DF_SRC)/lagrangian/spray/lnInclude \
+    -I$(LIB_SRC)/lagrangian/spray/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(DF_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(DF_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(DF_SRC)/dfCanteraMixture/lnInclude \
     -I$(DF_SRC)/dfChemistryModel/lnInclude \
+    -I$(LIB_SRC)/functionObjects/field/lnInclude \
     -I$(LIB_SRC)/Pstream/mpi \
     -I$(CANTERA_ROOT)/include \
     $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include,) \
@@ -22,11 +34,19 @@ EXE_INC = -std=c++14 \
 LIB_LIBS = \
     -lcompressibleTransportModels \
     -lturbulenceModels \
-    -llagrangian \    
+    -llagrangian \
     -lfiniteVolume \
+    -lfvOptions \
     -lmeshTools \
     -L$(DF_LIBBIN) \
     -ldfCompressibleTurbulenceModels \
+    -ldfLagrangianIntermediate \
+    -ldfLagrangianTurbulence \
+    -ldfLagrangianSpray \
+    -ldfFluidThermophysicalModels \
+    -ldfThermophysicalProperties \
+    -ldfSLGThermo \
+    -lregionModels \
     -ldfCanteraMixture \
     -ldfChemistryModel \
     $(CANTERA_ROOT)/lib/libcantera.so \
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 67ea36ed3..48e00e9f3 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -21,7 +21,7 @@ file(COPY ./dfLowMachFoam/twoD_reactingTGV/H2/cvodeSolver/postProcessing/sample/
 file(COPY ./dfLowMachFoam/twoD_reactingTGV/H2/cvodeSolver/postProcessing/sample/0.0002/data_T.xy DESTINATION 2DTGV/2)
 file(COPY ./dfLowMachFoam/twoD_reactingTGV/H2/cvodeSolver/postProcessing/sample/0.0001/data_T.xy DESTINATION 2DTGV/1)
 
-file(COPY ./dfLowMachFoam/2DSandiaD_flareFGM/postProcessing/sample/0.2/data_T.xy DESTINATION 2DSandia)
+file(COPY ./dfLowMachFoam/2DSandiaD_flareFGM/postProcessing/sample/1.1/data_T.xy DESTINATION 2DSandia)
 
 
 enable_testing()
diff --git a/test/corrtest.cpp b/test/corrtest.cpp
index a6472d520..5ce84f6cc 100644
--- a/test/corrtest.cpp
+++ b/test/corrtest.cpp
@@ -31,6 +31,11 @@ float T3 = readSandia(3,"2DSandia/data_T.xy");
 float T4 = readSandia(4,"2DSandia/data_T.xy");
 float T5 = readSandia(5,"2DSandia/data_T.xy");
 float T6 = readSandia(6,"2DSandia/data_T.xy");
+float T7 = readSandia(7,"2DSandia/data_T.xy");
+float T8 = readSandia(8,"2DSandia/data_T.xy");
+float T9 = readSandia(9,"2DSandia/data_T.xy");
+float T10 = readSandia(10,"2DSandia/data_T.xy");
+float T11 = readSandia(11,"2DSandia/data_T.xy");
 
 // TEST(corrtest,df0DFoam_H2){
 //     EXPECT_FLOAT_EQ(H2maxT,2588.48);   // compare the maximum temperature of H2 case 
@@ -64,13 +69,17 @@ TEST(corrtest,dfLowMachFoam_TGV){
 }
 
 TEST(corrtest,2DSandia){
-    EXPECT_FLOAT_EQ(T1,802.1133601);   
-    EXPECT_FLOAT_EQ(T2,1146.975511);  
-    EXPECT_FLOAT_EQ(T3,1569.239903);
-    EXPECT_FLOAT_EQ(T4,1376.352674);
-    EXPECT_FLOAT_EQ(T5,720.0857963);
-    EXPECT_FLOAT_EQ(T6,474.650212);
-    
+    EXPECT_FLOAT_EQ(T1,307.93593);   
+    EXPECT_FLOAT_EQ(T2,311.36824);  
+    EXPECT_FLOAT_EQ(T3,378.9158);
+    EXPECT_FLOAT_EQ(T4,661.57222);
+    EXPECT_FLOAT_EQ(T5,1109.6588);
+    EXPECT_FLOAT_EQ(T6,1521.5545);
+    EXPECT_FLOAT_EQ(T7,1865.9998);
+    EXPECT_FLOAT_EQ(T8,1963.6765);
+    EXPECT_FLOAT_EQ(T9,1783.7844);
+    EXPECT_FLOAT_EQ(T10,1474.4433);
+    EXPECT_FLOAT_EQ(T11,1069.4404);
 }
 
 float readmaxTH2(){