[NDSL] Feature/evap subl pdf (#1008)

* First bit of the evap_subl_pdf port. Initial calculation verifies, meltfrz does not.

* Updates to the evap_subl_pdf loop. meltfrz now verifies (excluding minor random errors on the order of 10^-8 which have yet to be investigated). evap and subl are believed to verify pending a fix to the cloud_effective_radius functions. Error in cloud_effective_radius functions appears to be a reappearance of the 32-bit precision errors previous noted when using these functions (errors are on the order of 10^-12). First attempt at hystpdf port, but still missing many pieces and does not execute.

* Additional minor changes to hystpdf, various other minor QOL changes

* SaturationVaporPressureTable code first pass
QSatIce missing before test can be done

* QSat_Ice_scalar_exact

* Fix calculation. Table generates OK.

* Fix Murphy and CAM

* Rename estimate_table -> table
Verbose the README

* Creation of QSat stencil, which still needs a lot of work. Need to push this code back to dsl/qsat branch as well at some point to have a working qsat function to pull on for all needs, but for now it is in evap_subl_pdf for development so that I can test incorporating it into the evap loop

* Functional QSat module, fails at majority of points but error pattern suggests numerics are as close as I can get them while remaining true to Fortran

* Updates to QSat and a few other minor QOL changes

* One more change that wasn't staged for some reason

* Functional QSat module. Can handle optional RAMP/PASCALS/DQSAT inputs. Fails to verify at majority of points using original metric, has not been tested with new metrics. Error patterns of the ES tables suggest arise during math within qsat_liquid and qsat_ice. Errors are small (10^-4). Further testing using find_klcl is in process to see if these errors propogate further or are locked within QSat.

* More clean up, first attempt at QSat function

* Working QSat function. QSat class can still be called outside of stencil to calculate QSat for entire field, or function can be called from within the stencil to calculate QSat at individual points. To use function, the QSat class must be initialized before stencil call and the tables must be passed down to the function as an input.

* Functional find_klcl port. Verifies at all but 2 points. Two failed cells are believed to be due to errors in QSat.

* Minor change for improved DQSAT functionality

* Minor change for improved DQSAT functionality

* Reorganization of functions/stencils from evap_subl_pdf loop to new design that (should) work for the full GFDL_1M module.

New stencil: hystpdf (in process)
New functions: pdffrac, pdfcondensate, bergeron_partition (in process)

Bergeron_partition calls GEOS_QsatLQU/ICE, which reaches into the non-table side of the QSAT world (not currently implemented, going to tackle that next).

All new additions will be tested as one once the entire hystpdf is ported.

* Two new functions: QSat_Float_Liquid and QSat_Float_Ice. These two functions are a port of QSATLQU0/QSATICE0 when UTBL is false.

Notable difference from QSat_Float: QSat_Float features a smooth transition from ice to liquid water (based on the TMIX parameter) while these two functions have a hard transition at 0C (should transition from calling one to the other).

UNTESTED FOR NOW, NO EASILY ACCESSABLE TEST CASE

* Completed (but untested) hystpdf. Lots of problems compiling the stencil.

* Tiny little bugfix fro QSat_Float_Liquid and QSat_Float_Ice

* hystpdf compiles, few fixes to deal with various gt4py/ndsl related errors/workarounds

* Fixed indexing error in QSat_Float_Liquid/Ice, and inserted a workaround for bad interaction between dace backend and cast to int

* hystpdf nearly verifies. lots of small errors, potentially rounding/32bit problems

* Code cleanup, ready to be pulled back into main branch

* Few more cleanups, modification to translate test to test sat tables

* QSat verifies if tables are read in from netcdf (this is the current state of the code). Generated tables do not yet verify.

* Fix run tests

* 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.

* Merged new QSat, manually fixed the qsat bugs that didn't make it pre-merge

* Debugging hystpdf, commiting now to pull new dsl/develop pr with QSat changes

* Functional hystpdf, verifies... when constants behave

* Functional, verified evap_subl_pdf loop.

Moved find_klcl into new misc_unused directory

Pre-commit run, fails in one spot b/c one comment is "too long"

* Update GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/misc_unused/find_klcl.py

Co-authored-by: Florian Deconinck <deconinck.florian@gmail.com>

* Bunch of changes per Florian's comments on the PR

* Linting (and one file rename)

* Update GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M.py

Co-authored-by: Florian Deconinck <deconinck.florian@gmail.com>

* More changes for the PR

* One last change

* Lint

* Changes per Tobias' comments less a proper variable naming scheme and implementation of new k-level feature. Will rename after k-level feature is implemented, want to remain true to fortran until development is finished for easier tracking.

* Couple more changes based on Florian & Tobias comments. find_klcl removed b/c it is unfinished and unused. will be recreated later on when it is needed

* Delete GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/scripts/run.out

* Lint & revert run_tests.sh

* ToDo for renaming hystpdf

---------

Co-authored-by: Florian Deconinck <deconinck.florian@gmail.com>

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M.py

@@ -0,0 +1,238 @@
+"""This module is the wrapper for the GFDL_1M microphysics scheme (in progress).
+I/O and errorhandling is performed here.
+Calculations can be found in deeper functions."""
+
+from ndsl import QuantityFactory, StencilFactory, orchestrate
+from ndsl.constants import X_DIM, Y_DIM, Z_DIM
+from ndsl.dsl.typing import Float, FloatField, FloatFieldIJ
+from pyMoist.GFDL_1M.evap_subl_pdf_core import (
+    evaporate,
+    hystpdf,
+    initial_calc,
+    melt_freeze,
+    sublimate,
+)
+from pyMoist.saturation.formulation import SaturationFormulation
+from pyMoist.saturation.qsat import QSat
+from pyMoist.shared_gt4py_workarounds import get_last, hybrid_index_2dout
+from pyMoist.shared_incloud_processes import fix_up_clouds
+
+
+class GFDL_1M:
+    """This class is the wrapper for the GFDL_1M microphysics scheme. I/O and error handling
+    are perfromed at this level, all calculations are performed within deeper functions.
+    """
+
+    def __init__(
+        self,
+        stencil_factory: StencilFactory,
+        quantity_factory: QuantityFactory,
+        formulation: SaturationFormulation = SaturationFormulation.Staars,
+        use_bergeron: bool = True,
+    ):
+        if use_bergeron is not True:
+            raise NotImplementedError(
+                "Untested option for use_bergeron. Code may be missing or incomplete. \
+                    Disable this error manually to continue."
+            )
+
+        self.stencil_factory = stencil_factory
+        self.quantity_factory = quantity_factory
+
+        # Initalize QSat tables for later calculations
+        self.qsat = QSat(
+            self.stencil_factory,
+            self.quantity_factory,
+            formulation=formulation,
+        )
+
+        orchestrate(obj=self, config=stencil_factory.config.dace_config)
+        self._get_last = self.stencil_factory.from_dims_halo(
+            func=get_last,
+            compute_dims=[X_DIM, Y_DIM, Z_DIM],
+        )
+        self._hybrid_index_2dout = self.stencil_factory.from_dims_halo(
+            func=hybrid_index_2dout,
+            compute_dims=[X_DIM, Y_DIM, Z_DIM],
+        )
+        self._initial_calc = self.stencil_factory.from_dims_halo(
+            func=initial_calc,
+            compute_dims=[X_DIM, Y_DIM, Z_DIM],
+        )
+        # TODO: rename: is this "Hydrostatic PDF"?
+        self._hystpdf = self.stencil_factory.from_dims_halo(
+            func=hystpdf,
+            compute_dims=[X_DIM, Y_DIM, Z_DIM],
+            externals={
+                "use_bergeron": use_bergeron,
+            },
+        )
+        self._meltfrz = self.stencil_factory.from_dims_halo(
+            func=melt_freeze,
+            compute_dims=[X_DIM, Y_DIM, Z_DIM],
+        )
+        self._evap = self.stencil_factory.from_dims_halo(
+            func=evaporate,
+            compute_dims=[X_DIM, Y_DIM, Z_DIM],
+        )
+        self._subl = self.stencil_factory.from_dims_halo(
+            func=sublimate,
+            compute_dims=[X_DIM, Y_DIM, Z_DIM],
+        )
+        self._fix_up_clouds = self.stencil_factory.from_dims_halo(
+            func=fix_up_clouds,
+            compute_dims=[X_DIM, Y_DIM, Z_DIM],
+        )
+        self._tmp = self.quantity_factory.zeros([X_DIM, Y_DIM], "n/a")
+        self._minrhcrit = self.quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a")
+        self._PLEmb_top = self.quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a")
+        self._PLmb_at_klcl = self.quantity_factory.zeros([X_DIM, Y_DIM], "n/a")
+        self._halo = self.stencil_factory.grid_indexing.n_halo
+        self._k_mask = self.quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a")
+        self._alpha = self.quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a")
+        self.rhx = self.quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a")
+        self.evapc = self.quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a")
+        self.sublc = self.quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a")
+        for i in range(0, self._k_mask.view[:].shape[0]):
+            for j in range(0, self._k_mask.view[:].shape[1]):
+                for k in range(0, self._k_mask.view[:].shape[2]):
+                    self._k_mask.view[i, j, k] = k + 1
+
+    def __call__(
+        self,
+        EIS: FloatFieldIJ,
+        dw_land: Float,
+        dw_ocean: Float,
+        PDFSHAPE: Float,
+        TURNRHCRIT_PARAM: Float,
+        PLmb: FloatField,
+        KLCL: FloatFieldIJ,
+        PLEmb: FloatField,
+        AREA: FloatFieldIJ,
+        DT_MOIST: Float,
+        CNV_FRC: FloatFieldIJ,
+        SRF_TYPE: FloatFieldIJ,
+        T: FloatField,
+        QLCN: FloatField,
+        QICN: FloatField,
+        QLLS: FloatField,
+        QILS: FloatField,
+        CCW_EVAP_EFF: Float,
+        CCI_EVAP_EFF: Float,
+        Q: FloatField,
+        CLLS: FloatField,
+        CLCN: FloatField,
+        NACTL: FloatField,
+        NACTI: FloatField,
+        QST: FloatField,
+        LMELTFRZ: bool = True,
+    ):
+        self.check_flags(LMELTFRZ, PDFSHAPE, CCW_EVAP_EFF, CCI_EVAP_EFF)
+
+        self._get_last(PLEmb, self._tmp, self._PLEmb_top)
+
+        self._hybrid_index_2dout(PLmb, self._k_mask, KLCL, self._PLmb_at_klcl)
+
+        self._initial_calc(
+            EIS,
+            dw_land,
+            dw_ocean,
+            TURNRHCRIT_PARAM,
+            self._minrhcrit,
+            self._PLmb_at_klcl,
+            PLmb,
+            self._PLEmb_top,
+            AREA,
+            self._alpha,
+        )
+
+        self._hystpdf(
+            DT_MOIST,
+            self._alpha,
+            PDFSHAPE,
+            CNV_FRC,
+            SRF_TYPE,
+            PLmb,
+            Q,
+            QLLS,
+            QLCN,
+            QILS,
+            QICN,
+            T,
+            CLLS,
+            CLCN,
+            NACTL,
+            NACTI,
+            self.rhx,
+            self.qsat.ese,
+            self.qsat.esw,
+            self.qsat.esx,
+            self.qsat.esw.view[0][12316],
+            self.qsat.esw.view[0][8316],
+        )
+
+        if LMELTFRZ:
+            self._meltfrz(DT_MOIST, CNV_FRC, SRF_TYPE, T, QLCN, QICN)
+            self._meltfrz(DT_MOIST, CNV_FRC, SRF_TYPE, T, QLLS, QILS)
+
+        if CCW_EVAP_EFF > 0.0:
+            self._evap(
+                DT_MOIST,
+                CCW_EVAP_EFF,
+                PLmb,
+                T,
+                Q,
+                QLCN,
+                QICN,
+                CLCN,
+                NACTL,
+                NACTI,
+                QST,
+                self.evapc,
+            )
+
+        if CCI_EVAP_EFF > 0.0:
+            self._subl(
+                DT_MOIST,
+                CCW_EVAP_EFF,
+                PLmb,
+                T,
+                Q,
+                QLCN,
+                QICN,
+                CLCN,
+                NACTL,
+                NACTI,
+                QST,
+                self.sublc,
+            )
+
+        self._fix_up_clouds(Q, T, QLLS, QILS, CLLS, QLCN, QICN, CLCN)
+
+    def check_flags(
+        self,
+        LMELTFRZ: bool,
+        PDFSHAPE: Float,
+        CCW_EVAP_EFF: Float,
+        CCI_EVAP_EFF: Float,
+    ):
+        if LMELTFRZ is not True:
+            raise NotImplementedError(
+                "Untested option for LMELTFRZ. Code may be missing or incomplete. \
+                    Disable this error manually to continue."
+            )
+        if PDFSHAPE != 1:
+            raise NotImplementedError(
+                "Untested option for PDFSHAPE. Code may be missing or incomplete. \
+                    Disable this error manually to continue."
+            )
+        if CCW_EVAP_EFF <= 0:
+            raise NotImplementedError(
+                "Untested option for CCW_EVAP_EFF. Code may be missing or incomplete. \
+                    Disable this error manually to continue."
+            )
+        if CCI_EVAP_EFF <= 0:
+            raise NotImplementedError(
+                "Untested option for CCI_EVAP_EFF. Code may be missing or incomplete. \
+                    Disable this error manually to continue."
+            )