diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/__init__.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/__init__.py
index 29a42185e..7fdd22809 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/__init__.py
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/__init__.py
@@ -1,2 +1,2 @@
 from .qsat import QSat, QSat_Float, QSat_Float_Ice, QSat_Float_Liquid, QSat_FloatField
-from .types import SaturationFormulation, TableMethod
+from .types import SaturationFormulation
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/constants.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/constants.py
index fb072a3d4..ce7c083ca 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/constants.py
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/constants.py
@@ -13,7 +13,7 @@
 TMAXTBL = Float(333.0)
 TMINLQU = MAPL_TICE - Float(40.0)
 DEGSUBS = np.int32(100)
-DELTA_T = Float(1.0 / DEGSUBS)
+DELTA_T = Float(1.0) / Float(DEGSUBS)
 TABLESIZE = np.int32(TMAXTBL - TMINTBL) * DEGSUBS + 1
 TMIX = Float(-20.0)
 ESFAC = MAPL_H2OMW / MAPL_AIRMW
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/netCDFs/QSat_Tables.nc b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/netCDFs/QSat_Tables.nc
deleted file mode 100644
index 9a98982fe..000000000
Binary files a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/netCDFs/QSat_Tables.nc and /dev/null differ
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/qsat.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/qsat.py
index 3c6375e6f..a646e6ae8 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/qsat.py
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/qsat.py
@@ -1,9 +1,7 @@
 import copy
-import os
 from typing import Optional
 
 import gt4py.cartesian.gtscript as gtscript
-import xarray as xr
 from gt4py.cartesian.gtscript import PARALLEL, computation, floor, interval
 
 from ndsl import QuantityFactory, StencilFactory, orchestrate
@@ -21,7 +19,8 @@
     TMINTBL,
     TMIX,
 )
-from pyMoist.saturation.types import SaturationFormulation, TableMethod
+from pyMoist.saturation.table import get_table
+from pyMoist.saturation.types import SaturationFormulation
 
 
 # FloatField with extra dimension initialized to handle table data
@@ -194,17 +193,16 @@ def QSat_FloatField(
     QSAT: FloatField,
     DQSAT: FloatField,
     RAMP: FloatField,
-    PASCALS_trigger: bool,
-    RAMP_trigger: bool,
-    DQSAT_trigger: bool,
 ):
+    from __externals__ import FILL_DQSAT, USE_PASCALS, USE_RAMP
+
     with computation(PARALLEL), interval(...):
-        if RAMP_trigger:
+        if USE_RAMP:
             URAMP = -abs(RAMP)
         else:
             URAMP = TMIX
 
-        if PASCALS_trigger:
+        if USE_PASCALS:
             PP = PL
         else:
             PP = PL * 100.0
@@ -230,12 +228,12 @@ def QSat_FloatField(
             QSAT = (TI - IT) * DQ + ese[0][IT_MINUS_1]  # type: ignore
         if PP <= QSAT:
             QSAT = MAX_MIXING_RATIO
-            if DQSAT_trigger:
+            if FILL_DQSAT:
                 DQSAT = 0.0
         else:
             DD = 1.0 / (PP - (1.0 - ESFAC) * QSAT)
             QSAT = ESFAC * QSAT * DD
-            if DQSAT_trigger:
+            if FILL_DQSAT:
                 DQSAT = ESFAC * DQ * DEGSUBS * PP * (DD * DD)
 
 
@@ -269,7 +267,9 @@ def __init__(
         stencil_factory: StencilFactory,
         quantity_factory: QuantityFactory,
         formulation: SaturationFormulation = SaturationFormulation.Staars,
-        table_method: TableMethod = TableMethod.NetCDF,
+        use_ramp: bool = False,
+        use_pascals: bool = False,
+        fill_dqsat: bool = False,
     ) -> None:
         self.extra_dim_quantity_factory = self.make_extra_dim_quantity_factory(
             quantity_factory
@@ -279,34 +279,10 @@ def __init__(
         self.esw = self.extra_dim_quantity_factory.zeros([Z_DIM, "table_axis"], "n/a")
         self.esx = self.extra_dim_quantity_factory.zeros([Z_DIM, "table_axis"], "n/a")
 
-        if table_method == TableMethod.NetCDF:
-            with xr.open_dataset(
-                os.path.join(os.path.dirname(__file__), "netCDFs", "QSat_Tables.nc")
-            ) as ds:
-                ese_array = ds.data_vars["ese"].values[0, 0, :]
-                esw_array = ds.data_vars["esw"].values[0, 0, :]
-                esx_array = ds.data_vars["esx"].values[0, 0, :]
-
-                self.ese.view[:] = ese_array
-                self.esw.view[:] = esw_array
-                self.esx.view[:] = esx_array
-        elif table_method == TableMethod.Computation:
-            raise NotImplementedError(
-                """[QSat]  Calculated tables are incorrect due to differences between
-                C and Fortran calculations.\n
-                For now the tables are being read in from a netCDF file.\n
-                See https://github.com/GEOS-ESM/SMT-Nebulae/issues/88
-                for more information."""
-            )
-            # self.table = get_table(formulation)
-
-            # self.ese.view[:] = self.table.ese
-            # self.esw.view[:] = self.table.esw
-            # self.esx.view[:] = self.table.esx
-        else:
-            raise NotImplementedError(
-                f"[QSat] Unknown {table_method} to create estimation table"
-            )
+        self.table = get_table(formulation)
+        self.ese.view[:] = self.table.ese
+        self.esw.view[:] = self.table.esw
+        self.esx.view[:] = self.table.esx
 
         self._RAMP = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a")
         self._DQSAT = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a")
@@ -317,6 +293,11 @@ def __init__(
         self._QSat_FloatField = stencil_factory.from_dims_halo(
             func=QSat_FloatField,
             compute_dims=[X_DIM, Y_DIM, Z_DIM],
+            externals={
+                "USE_RAMP": use_ramp,
+                "USE_PASCALS": use_pascals,
+                "FILL_DQSAT": fill_dqsat,
+            },
         )
 
     @staticmethod
@@ -333,16 +314,13 @@ def __call__(
         self,
         T: FloatField,
         PL: FloatField,
-        RAMP: Optional[FloatField] = None,
-        PASCALS: bool = False,
-        DQSAT: bool = False,
+        DQSat: Optional[FloatField] = None,
         use_table_lookup: bool = True,
     ):
-        if RAMP is None:
-            RAMP = self._RAMP
-            RAMP_trigger = False
+        if DQSat:
+            dqsat = DQSat
         else:
-            RAMP_trigger = True
+            dqsat = self._DQSAT
 
         if use_table_lookup:
             self._QSat_FloatField(
@@ -351,11 +329,8 @@ def __call__(
                 T,
                 PL,
                 self.QSat,
-                self._DQSAT,
-                RAMP,
-                PASCALS,
-                RAMP_trigger,
-                DQSAT,
+                dqsat,
+                self._RAMP,
             )
 
         if not use_table_lookup:
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/qsat_ice.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/qsat_ice.py
index 4aa51e008..0d4b81dcd 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/qsat_ice.py
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/qsat_ice.py
@@ -13,10 +13,11 @@
 from pyMoist.saturation.types import SaturationFormulation
 
 
+f64 = np.float64
+
 TMINSTR = Float(-95.0)
 TMINICE = MAPL_TICE + TMINSTR
 
-
 TMINSTR = Float(-95.0)
 TSTARR1 = Float(-75.0)
 TSTARR2 = Float(-65.0)
@@ -24,31 +25,41 @@
 TSTARR4 = Float(-40.0)
 TMAXSTR = Float(+60.0)
 
-DI = [57518.5606e08, 2.01889049, 3.56654, 20.947031]
-CI = [9.550426, -5723.265, 3.53068, -0.00728332]
+DI = [
+    f64(Float(57518.5606e08)),
+    f64(Float(2.01889049)),
+    f64(Float(3.56654)),
+    f64(Float(20.947031)),
+]
+CI = [
+    f64(Float(9.550426)),
+    f64(Float(-5723.265)),
+    f64(Float(3.53068)),
+    f64(Float(-0.00728332)),
+]
 
 # 64-bit float in Fortran
-S16 = 0.516000335e-11 * 100.0
-S15 = 0.276961083e-8 * 100.0
-S14 = 0.623439266e-6 * 100.0
-S13 = 0.754129933e-4 * 100.0
-S12 = 0.517609116e-2 * 100.0
-S11 = 0.191372282e0 * 100.0
-S10 = 0.298152339e1 * 100.0
-S26 = 0.314296723e-10 * 100.0
-S25 = 0.132243858e-7 * 100.0
-S24 = 0.236279781e-5 * 100.0
-S23 = 0.230325039e-3 * 100.0
-S22 = 0.129690326e-1 * 100.0
-S21 = 0.401390832e0 * 100.0
-S20 = 0.535098336e1 * 100.0
-BI6 = 1.838826904e-10 * 100.0
-BI5 = 4.838803174e-8 * 100.0
-BI4 = 5.824720280e-6 * 100.0
-BI3 = 4.176223716e-4 * 100.0
-BI2 = 1.886013408e-2 * 100.0
-BI1 = 5.034698970e-1 * 100.0
-BI0 = 6.109177956e0 * 100.0
+S16 = f64(Float(0.516000335e-11) * Float(100.0))
+S15 = f64(Float(0.276961083e-8) * Float(100.0))
+S14 = f64(Float(0.623439266e-6) * Float(100.0))
+S13 = f64(Float(0.754129933e-4) * Float(100.0))
+S12 = f64(Float(0.517609116e-2) * Float(100.0))
+S11 = f64(Float(0.191372282e0) * Float(100.0))
+S10 = f64(Float(0.298152339e1) * Float(100.0))
+S26 = f64(Float(0.314296723e-10) * Float(100.0))
+S25 = f64(Float(0.132243858e-7) * Float(100.0))
+S24 = f64(Float(0.236279781e-5) * Float(100.0))
+S23 = f64(Float(0.230325039e-3) * Float(100.0))
+S22 = f64(Float(0.129690326e-1) * Float(100.0))
+S21 = f64(Float(0.401390832e0) * Float(100.0))
+S20 = f64(Float(0.535098336e1) * Float(100.0))
+BI6 = f64(Float(1.838826904e-10) * Float(100.0))
+BI5 = f64(Float(4.838803174e-8) * Float(100.0))
+BI4 = f64(Float(5.824720280e-6) * Float(100.0))
+BI3 = f64(Float(4.176223716e-4) * Float(100.0))
+BI2 = f64(Float(1.886013408e-2) * Float(100.0))
+BI1 = f64(Float(5.034698970e-1) * Float(100.0))
+BI0 = f64(Float(6.109177956e0) * Float(100.0))
 
 
 def _saturation_formulation(
@@ -69,7 +80,7 @@ def _saturation_formulation(
                 TT
                 * (TT * (TT * (TT * (TT * (TT * S16 + S15) + S14) + S13) + S12) + S11)
                 + S10
-            ) + (1.0 - W) * (
+            ) + (Float(1.0) - W) * (
                 TT
                 * (TT * (TT * (TT * (TT * (TT * S26 + S25) + S24) + S23) + S22) + S21)
                 + S20
@@ -86,7 +97,7 @@ def _saturation_formulation(
                 TT
                 * (TT * (TT * (TT * (TT * (TT * S26 + S25) + S24) + S23) + S22) + S21)
                 + S20
-            ) + (1.0 - W) * (
+            ) + (Float(1.0) - W) * (
                 TT
                 * (TT * (TT * (TT * (TT * (TT * BI6 + BI5) + BI4) + BI3) + BI2) + BI1)
                 + BI0
@@ -119,14 +130,14 @@ def qsat_ice_scalar_exact(
     else:
         TI = temperature
 
-    DX = 0.0  # only calulcated when DQ is not none
+    DX = Float(0.0)  # only calulcated when DQ is not none
     EX = _saturation_formulation(formulation, TI)
 
     if DQ is not None:
         if temperature < TMINICE:
-            DDQ = 0.0
+            DDQ = Float(0.0)
         elif temperature > MAPL_TICE:
-            DDQ = 0.0
+            DDQ = Float(0.0)
         else:
             if PL > EX:
                 DD = EX
@@ -136,16 +147,16 @@ def qsat_ice_scalar_exact(
                 EX = DD
     if PL is not None:
         if PL > EX:
-            DD = ESFAC / (PL - (1.0 - ESFAC) * EX)
+            DD = ESFAC / (PL - (Float(1.0) - ESFAC) * EX)
             EX = EX * DD
             if DQ is not None:
                 DX = DDQ * ERFAC * PL * DD * DD
         else:
             EX = MAX_MIXING_RATIO
             if DQ is not None:
-                DX = 0.0
+                DX = Float(0.0)
     else:
         if DQ is not None:
-            DX = DDQ * (1.0 / DELTA_T)
+            DX = DDQ * (Float(1.0) / DELTA_T)
 
     return EX, TI, DX
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/qsat_liquid.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/qsat_liquid.py
index 64da52861..b83bdecd8 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/qsat_liquid.py
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/qsat_liquid.py
@@ -15,50 +15,59 @@
 from pyMoist.saturation.types import SaturationFormulation
 
 
+f64 = np.float64
+
 # Below are actual 64-bit float in Fortran
-B6 = 6.136820929e-11 * 100.0
-B5 = 2.034080948e-8 * 100.0
-B4 = 3.031240396e-6 * 100.0
-B3 = 2.650648471e-4 * 100.0
-B2 = 1.428945805e-2 * 100.0
-B1 = 4.436518521e-1 * 100.0
-B0 = 6.107799961e0 * 100.0
-BI6 = 1.838826904e-10 * 100.0
-BI5 = 4.838803174e-8 * 100.0
-BI4 = 5.824720280e-6 * 100.0
-BI3 = 4.176223716e-4 * 100.0
-BI2 = 1.886013408e-2 * 100.0
-BI1 = 5.034698970e-1 * 100.0
-BI0 = 6.109177956e0 * 100.0
-S16 = 0.516000335e-11 * 100.0
-S15 = 0.276961083e-8 * 100.0
-S14 = 0.623439266e-6 * 100.0
-S13 = 0.754129933e-4 * 100.0
-S12 = 0.517609116e-2 * 100.0
-S11 = 0.191372282e0 * 100.0
-S10 = 0.298152339e1 * 100.0
-S26 = 0.314296723e-10 * 100.0
-S25 = 0.132243858e-7 * 100.0
-S24 = 0.236279781e-5 * 100.0
-S23 = 0.230325039e-3 * 100.0
-S22 = 0.129690326e-1 * 100.0
-S21 = 0.401390832e0 * 100.0
-S20 = 0.535098336e1 * 100.0
+B6 = f64(Float(6.136820929e-11) * Float(100.0))
+B5 = f64(Float(2.034080948e-8) * Float(100.0))
+B4 = f64(Float(3.031240396e-6) * Float(100.0))
+B3 = f64(Float(2.650648471e-4) * Float(100.0))
+B2 = f64(Float(1.428945805e-2) * Float(100.0))
+B1 = f64(Float(4.436518521e-1) * Float(100.0))
+B0 = f64(Float(6.107799961e0) * Float(100.0))
+BI6 = f64(Float(1.838826904e-10) * Float(100.0))
+BI5 = f64(Float(4.838803174e-8) * Float(100.0))
+BI4 = f64(Float(5.824720280e-6) * Float(100.0))
+BI3 = f64(Float(4.176223716e-4) * Float(100.0))
+BI2 = f64(Float(1.886013408e-2) * Float(100.0))
+BI1 = f64(Float(5.034698970e-1) * Float(100.0))
+BI0 = f64(Float(6.109177956e0) * Float(100.0))
+S16 = f64(Float(0.516000335e-11) * Float(100.0))
+S15 = f64(Float(0.276961083e-8) * Float(100.0))
+S14 = f64(Float(0.623439266e-6) * Float(100.0))
+S13 = f64(Float(0.754129933e-4) * Float(100.0))
+S12 = f64(Float(0.517609116e-2) * Float(100.0))
+S11 = f64(Float(0.191372282e0) * Float(100.0))
+S10 = f64(Float(0.298152339e1) * Float(100.0))
+S26 = f64(Float(0.314296723e-10) * Float(100.0))
+S25 = f64(Float(0.132243858e-7) * Float(100.0))
+S24 = f64(Float(0.236279781e-5) * Float(100.0))
+S23 = f64(Float(0.230325039e-3) * Float(100.0))
+S22 = f64(Float(0.129690326e-1) * Float(100.0))
+S21 = f64(Float(0.401390832e0) * Float(100.0))
+S20 = f64(Float(0.535098336e1) * Float(100.0))
 
-DL = [-7.902980, 5.02808, -1.3816, 11.344, 8.1328, -3.49149]
-TS = 373.16
-LOGPS = 3.005714898  # log10(1013.246)
+DL = [
+    f64(Float(-7.902980)),
+    f64(Float(5.02808)),
+    f64(Float(-1.3816)),
+    f64(Float(11.344)),
+    f64(Float(8.1328)),
+    f64(Float(-3.49149)),
+]
+TS = f64(Float(373.16))
+LOGPS = f64(Float(3.005714898))  # log10(1013.246)
 CL = [
-    54.842763,
-    -6763.22,
-    -4.21000,
-    0.000367,
-    0.0415,
-    218.8,
-    53.878000,
-    -1331.22,
-    -9.44523,
-    0.014025,
+    f64(Float(54.842763)),
+    f64(Float(-6763.22)),
+    f64(Float(-4.21000)),
+    f64(Float(0.000367)),
+    f64(Float(0.0415)),
+    f64(Float(218.8)),
+    f64(Float(53.878000)),
+    f64(Float(-1331.22)),
+    f64(Float(-9.44523)),
+    f64(Float(0.014025)),
 ]
 
 
@@ -68,13 +77,17 @@ def _saturation_formulation(formulation: SaturationFormulation, t: Float):
         EX = TT * (TT * (TT * (TT * (TT * (TT * B6 + B5) + B4) + B3) + B2) + B1) + B0
     elif formulation == SaturationFormulation.CAM:
         TT = TS / t
-        EX = 10.0 ** (
-            DL[0] * (TT - 1.0)
+        EX = Float(10.0) ** (
+            DL[0] * (TT - Float(1.0))
             + DL[1] * np.log10(TT)
-            + DL[2] * (10.0 ** (DL[3] * (1.0 - (1.0 / TT))) - 1.0) / 10000000.0
-            + DL[4] * (10.0 ** (DL[5] * (TT - 1.0)) - 1.0) / 1000.0
+            + DL[2]
+            * (Float(10.0) ** (DL[3] * (Float(1.0) - (Float(1.0) / TT))) - Float(1.0))
+            / Float(10000000.0)
+            + DL[4]
+            * (Float(10.0) ** (DL[5] * (TT - Float(1.0))) - Float(1.0))
+            / Float(1000.0)
             + LOGPS
-            + 2.0
+            + Float(2.0)
         )
     elif formulation == SaturationFormulation.MurphyAndKoop:
         EX = np.exp(
@@ -100,14 +113,14 @@ def qsat_liquid_scalar_exact(
     else:
         TI = temperature
 
-    DX = 0.0  # DX only calculated when DQ is present
+    DX = Float(0.0)  # DX only calculated when DQ is present
     EX = _saturation_formulation(formulation, TI)
 
     if DQ is not None:
         if temperature < TMINLQU:
-            DDQ = 0.0
+            DDQ = Float(0.0)
         elif temperature > TMAXTBL:
-            DDQ = 0.0
+            DDQ = Float(0.0)
         else:
             if PL > EX:
                 DD = EX
@@ -118,15 +131,15 @@ def qsat_liquid_scalar_exact(
 
     if PL is not None:
         if PL > EX:
-            DD = ESFAC / (PL - (1.0 - ESFAC) * EX)
+            DD = ESFAC / (PL - (Float(1.0) - ESFAC) * EX)
             EX = EX * DD
             if DQ is not None:
                 DX = DDQ * ERFAC * PL * DD * DD
         else:
             EX = MAX_MIXING_RATIO
             if DQ is not None:
-                DX = 0.0
+                DX = Float(0.0)
     elif DQ is not None:
-        DX = DDQ * (1.0 / DELTA_T)
+        DX = DDQ * (Float(1.0) / DELTA_T)
 
     return EX, TI, DX
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/table.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/table.py
index a146524cf..71fb33f27 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/table.py
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/table.py
@@ -35,12 +35,11 @@ def __init__(
         self._LOC = np.empty(TABLESIZE, dtype=Float)
 
         for i in range(TABLESIZE):
-            t = i * DELTA_T + TMINTBL
+            t = Float(i * DELTA_T) + TMINTBL
             (
                 self._estimated_esw[i],
                 self._TI[i],
                 _,
-                self._LOC[i],
             ) = qsat_liquid_scalar_exact(t, formulation)
 
             if t > MAPL_TICE:
@@ -50,7 +49,6 @@ def __init__(
                     self._estimated_ese[i],
                     self._TI[i],
                     _,
-                    self._LOC[i],
                 ) = qsat_ice_scalar_exact(t, formulation)
 
             t = t - MAPL_TICE
@@ -62,8 +60,8 @@ def __init__(
             else:
                 self._estimated_esx[i] = self._estimated_ese[i]
 
-        self._estimated_frz, _, _, _ = qsat_liquid_scalar_exact(MAPL_TICE, formulation)
-        self._estimated_lqu, _, _, _ = qsat_liquid_scalar_exact(TMINLQU, formulation)
+        self._estimated_frz, _, _ = qsat_liquid_scalar_exact(MAPL_TICE, formulation)
+        self._estimated_lqu, _, _ = qsat_liquid_scalar_exact(TMINLQU, formulation)
 
     @property
     def ese(self):
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/types.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/types.py
index 2d2044427..c59413cf7 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/types.py
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/types.py
@@ -11,14 +11,3 @@ class SaturationFormulation(enum.Enum):
     Staars = 1
     CAM = 2
     MurphyAndKoop = 3
-
-
-class TableMethod(enum.Enum):
-    """
-    The table creation methods:
-        - load from netCDF
-        - compute form exact formulation
-    """
-
-    NetCDF = 1
-    Computation = 2
diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/translate_qsat.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/translate_qsat.py
index 4e10cb258..ff38999d7 100644
--- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/translate_qsat.py
+++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/translate_qsat.py
@@ -1,7 +1,7 @@
 from ndsl import Namelist, Quantity, StencilFactory
 from ndsl.constants import X_DIM, Y_DIM, Z_DIM
 from ndsl.stencils.testing.translate import TranslateFortranData2Py
-from pyMoist.saturation import QSat, TableMethod
+from pyMoist.saturation import QSat
 
 
 class TranslateQSat(TranslateFortranData2Py):
@@ -58,7 +58,6 @@ def compute(self, inputs):
         qsat = QSat(
             self.stencil_factory,
             self.quantity_factory,
-            table_method=TableMethod.NetCDF,
         )
 
         # FloatField Variables