Functional QSat module. Table calculations are still incorrect due to…

… differences between Fortran and C interpretation of constants (see GEOS-ESM/SMT-Nebulae#88 for more information). Table values are instead read in from netCDF file as a temporary solution.

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/constants.py

@@ -11,7 +11,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 / DEGSUBS)
 TABLESIZE = np.int32(TMAXTBL - TMINTBL) * DEGSUBS + 1
 TMIX = Float(-20.0)
 ESFAC = MAPL_H2OMW / MAPL_AIRMW

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/qsat.py

@@ -1,10 +1,12 @@
+import os
 from ndsl import StencilFactory, QuantityFactory, orchestrate
 from ndsl.dsl.typing import Float, FloatField, Int, IntField
 from ndsl.constants import X_DIM, Y_DIM, Z_DIM
 from gt4py.cartesian.gtscript import computation, interval, PARALLEL, floor
 import gt4py.cartesian.gtscript as gtscript
 import copy
 from typing import Optional
+import xarray as xr
 from pyMoist.saturation.formulation import SaturationFormulation
 from pyMoist.saturation.table import get_table
 from pyMoist.saturation.constants import (
@@ -271,20 +273,39 @@ def __init__(
         stencil_factory: StencilFactory,
         quantity_factory: QuantityFactory,
         formulation: SaturationFormulation = SaturationFormulation.Staars,
+        table_method: Float = 0,
     ) -> None:
-
-        self.table = get_table(formulation)
-
         self.extra_dim_quantity_factory = self.make_extra_dim_quantity_factory(
             quantity_factory
         )
 
         self.ese = self.extra_dim_quantity_factory.zeros([Z_DIM, "table_axis"], "n/a")
         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")
-        self.ese.view[:] = self.table.ese
-        self.esw.view[:] = self.table.esw
-        self.esx.view[:] = self.table.esx
+
+        if table_method == 0:
+            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
+
+        if table_method == 1:
+            raise NotImplementedError(
+                """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
 
         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")

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/qsat_ice.py

@@ -54,7 +54,7 @@ def _saturation_formulation(
 ):
     if formulation == SaturationFormulation.Staars:
         TT = t - MAPL_TICE
-        if TT <= TSTARR1:
+        if TT < TSTARR1:
             LOC = 1.1
             EX = (
                 TT

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/saturation/table.py

@@ -29,16 +29,20 @@ def __init__(
         self._estimated_ese = np.empty(TABLESIZE, dtype=Float)
         self._estimated_esx = np.empty(TABLESIZE, dtype=Float)
         self._TI = np.empty(TABLESIZE, dtype=Float)
+        self._LOC = np.empty(TABLESIZE, dtype=Float)
 
         for i in range(TABLESIZE):
             t = i * DELTA_T + TMINTBL
-
-            self._estimated_esw[i], self._TI[i], _ = qsat_liquid_scalar_exact(t, formulation)
+            self._estimated_esw[i], self._TI[i], _, self._LOC[i] = (
+                qsat_liquid_scalar_exact(t, formulation)
+            )
 
             if t > MAPL_TICE:
                 self._estimated_ese[i] = self._estimated_esw[i]
             else:
-                self._estimated_ese[i], self._TI[i], _ = qsat_ice_scalar_exact(t, formulation)
+                self._estimated_ese[i], self._TI[i], _, self._LOC[i] = (
+                    qsat_ice_scalar_exact(t, formulation)
+                )
 
             t = t - MAPL_TICE
 
@@ -49,8 +53,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):

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/translate_qsat.py

@@ -41,11 +41,11 @@ def __init__(self, grid, namelist: Namelist, stencil_factory: StencilFactory):
             "ESTBLE_TEST": self.grid.compute_dict(),
             "ESTBLW_TEST": self.grid.compute_dict(),
             "ESTBLX_TEST": self.grid.compute_dict(),
-            "IT": self.grid.compute_dict(),
-            "IFELSE": self.grid.compute_dict(),
-            "QSAT_HALFWAY": self.grid.compute_dict(),
-            "TI": self.grid.compute_dict(),
-            "LOC": self.grid.compute_dict(),
+            # "IT": self.grid.compute_dict(),
+            # "IFELSE": self.grid.compute_dict(),
+            # "QSAT_HALFWAY": self.grid.compute_dict(),
+            # "TI": self.grid.compute_dict(),
+            # "LOC": self.grid.compute_dict(),
         }
 
     def make_ij_field(self, data) -> Quantity:
@@ -95,19 +95,19 @@ def compute(self, inputs):
         code(
             T,
             PL,
-            ese=ese,
-            esw=esw,
-            esx=esx,
+            # ese=ese,
+            # esw=esw,
+            # esx=esx,
         )
 
         return {
             "QSAT": code.QSat.view[:],
             "ESTBLE_TEST": code.ese.view[0],
             "ESTBLW_TEST": code.esw.view[0],
             "ESTBLX_TEST": code.esx.view[0],
-            "IT": code._IT.view[:],
-            "IFELSE": code._IFELSE.view[:],
-            "QSAT_HALFWAY": code._QSAT_HALFWAY.view[:],
-            "TI": code._TI.view[:],
-            "LOC": code.table._LOC,
+            # "IT": code._IT.view[:],
+            # "IFELSE": code._IFELSE.view[:],
+            # "QSAT_HALFWAY": code._QSAT_HALFWAY.view[:],
+            # "TI": code._TI.view[:],
+            # "LOC": code.table._LOC,
         }

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/scripts/run_tests.sh

@@ -1,11 +1,10 @@
 #!/bin/bash
-rm -rf ./.gt_cache_*
+rm -rf ../.gt_cache_*
 export PACE_FLOAT_PRECISION=32
-export FV3_DACEMODE=Python
-python -m pytest -v -s -x\
-    --pdb \
-    --data_path=../../test_data/11.5.2/Moist/TBC_C24_L72_Debug \
+python -m pytest -s --disable-warnings --multimodal_metric \
+    --data_path=/home/charleskrop/netcdfs \
     --backend=dace:cpu \
+    --which_rank=0 \
+    --which_modules=QSat \
     --grid=default \
-    --multimodal_metric \
-    ..
+    ..