diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/GFDL_1M_driver.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/GFDL_1M_driver.py deleted file mode 100644 index 452f7524b..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/GFDL_1M_driver.py +++ /dev/null @@ -1,1257 +0,0 @@ -"""GFDL_1M driver""" - -import numpy as np -from gt4py.cartesian.gtscript import i32 - -import pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver_constants as driver_constants -from ndsl import QuantityFactory, StencilFactory, orchestrate -from ndsl.constants import X_DIM, Y_DIM, Z_DIM, Z_INTERFACE_DIM -from ndsl.dsl.typing import Float, FloatField, FloatFieldIJ, Int -from pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver_core import ( - fall_speed, - fix_negative_values, - icloud_update, - init_temporaries, - terminal_fall_update, - update_tendencies, - warm_rain_update, -) -from pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver_tables import get_tables -from pyMoist.GFDL_1M.GFDL_1M_driver.icloud import icloud -from pyMoist.GFDL_1M.GFDL_1M_driver.terminal_fall import terminal_fall -from pyMoist.GFDL_1M.GFDL_1M_driver.warm_rain import warm_rain - - -class GFDL_1M_driver: - def __init__( - self, - stencil_factory: StencilFactory, - quantity_factory: QuantityFactory, - phys_hydrostatic: bool, - hydrostatic: bool, - dt_moist: Float, - mp_time: Float, - t_min: Float, - t_sub: Float, - tau_r2g: Float, - tau_smlt: Float, - tau_g2r: Float, - dw_land: Float, - dw_ocean: Float, - vi_fac: Float, - vr_fac: Float, - vs_fac: Float, - vg_fac: Float, - ql_mlt: Float, - do_qa: bool, - fix_negative: bool, - vi_max: Float, - vs_max: Float, - vg_max: Float, - vr_max: Float, - qs_mlt: Float, - qs0_crt: Float, - qi_gen: Float, - ql0_max: Float, - qi0_max: Float, - qi0_crt: Float, - qr0_crt: Float, - fast_sat_adj: bool, - rh_inc: Float, - rh_ins: Float, - rh_inr: Float, - const_vi: bool, - const_vs: bool, - const_vg: bool, - const_vr: bool, - use_ccn: bool, - rthreshu: Float, - rthreshs: Float, - ccn_l: Float, - ccn_o: Float, - qc_crt: Float, - tau_g2v: Float, - tau_v2g: Float, - tau_s2v: Float, - tau_v2s: Float, - tau_revp: Float, - tau_frz: Float, - do_bigg: bool, - do_evap: bool, - do_subl: bool, - sat_adj0: Float, - c_piacr: Float, - tau_imlt: Float, - tau_v2l: Float, - tau_l2v: Float, - tau_i2v: Float, - tau_i2s: Float, - tau_l2r: Float, - qi_lim: Float, - ql_gen: Float, - c_paut: Float, - c_psaci: Float, - c_pgacs: Float, - c_pgaci: Float, - z_slope_liq: bool, - z_slope_ice: bool, - prog_ccn: bool, - c_cracw: Float, - alin: Float, - clin: Float, - preciprad: bool, - cld_min: Float, - use_ppm: bool, - mono_prof: bool, - do_sedi_heat: bool, - sedi_transport: bool, - do_sedi_w: bool, - de_ice: bool, - icloud_f: Float, - irain_f: Float, - mp_print: bool, - ): - self.dt_moist = dt_moist - self.fix_negative = fix_negative - self.sedi_transport = sedi_transport - - self.namelist_constants = namelist_setup( - phys_hydrostatic, - hydrostatic, - dt_moist, - mp_time, - t_min, - t_sub, - tau_r2g, - tau_smlt, - tau_g2r, - dw_land, - dw_ocean, - vi_fac, - vr_fac, - vs_fac, - vg_fac, - ql_mlt, - do_qa, - fix_negative, - vi_max, - vs_max, - vg_max, - vr_max, - qs_mlt, - qs0_crt, - qi_gen, - ql0_max, - qi0_max, - qi0_crt, - qr0_crt, - fast_sat_adj, - rh_inc, - rh_ins, - rh_inr, - const_vi, - const_vs, - const_vg, - const_vr, - use_ccn, - rthreshu, - rthreshs, - ccn_l, - ccn_o, - qc_crt, - tau_g2v, - tau_v2g, - tau_s2v, - tau_v2s, - tau_revp, - tau_frz, - do_bigg, - do_evap, - do_subl, - sat_adj0, - c_piacr, - tau_imlt, - tau_v2l, - tau_l2v, - tau_i2v, - tau_i2s, - tau_l2r, - qi_lim, - ql_gen, - c_paut, - c_psaci, - c_pgacs, - c_pgaci, - z_slope_liq, - z_slope_ice, - prog_ccn, - c_cracw, - alin, - clin, - preciprad, - cld_min, - use_ppm, - mono_prof, - do_sedi_heat, - sedi_transport, - do_sedi_w, - de_ice, - icloud_f, - irain_f, - mp_print, - ) - - # Check values for untested code paths - self.check_flags( - phys_hydrostatic, - hydrostatic, - const_vi, - const_vs, - const_vg, - const_vr, - use_ppm, - use_ccn, - do_qa, - fix_negative, - fast_sat_adj, - do_bigg, - do_evap, - do_subl, - z_slope_liq, - z_slope_ice, - prog_ccn, - preciprad, - mono_prof, - do_sedi_heat, - sedi_transport, - do_sedi_w, - de_ice, - mp_print, - self.namelist_constants.DTS, - ) - - # ----------------------------------------------------------------------- - # initialize precipitation outputs - # ----------------------------------------------------------------------- - - self.rain = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") - self.snow = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") - self.ice = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") - self.graupel = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") - self.m2_rain = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.m2_sol = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.revap = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.isubl = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - - # ----------------------------------------------------------------------- - # initialize temporaries - # ----------------------------------------------------------------------- - - self.t1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.dp1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.omq = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.qv0 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.ql0 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.qr0 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.qi0 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.qs0 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.qg0 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.qa0 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.qv1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.ql1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.qr1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.qi1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.qs1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.qg1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.qa1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.dz1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.den = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.den1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.denfac = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.p_dry = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.m1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.u1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.v1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.w1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.onemsig = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") - self.ccn = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.c_praut = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.rh_limited = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.ze = quantity_factory.zeros([X_DIM, Y_DIM, Z_INTERFACE_DIM], "n/a") - self.zt = quantity_factory.zeros([X_DIM, Y_DIM, Z_INTERFACE_DIM], "n/a") - self.lhi = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.icpk = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.hold_data = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.vti = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.vts = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.vtg = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.vtr = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.m1_sol = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.m1_rain = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.rain1 = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") - self.graupel1 = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") - self.snow1 = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") - self.ice1 = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") - self.evap1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.subl1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - - # ----------------------------------------------------------------------- - # initialize masks - # ----------------------------------------------------------------------- - self.is_frozen = quantity_factory.zeros( - [X_DIM, Y_DIM, Z_DIM], "n/a", dtype=bool - ) - self.precip_fall = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") - self.melting_mask_1 = quantity_factory.zeros( - [X_DIM, Y_DIM, Z_DIM], "n/a", dtype=bool - ) - self.melting_mask_2 = quantity_factory.zeros( - [X_DIM, Y_DIM, Z_DIM], "n/a", dtype=bool - ) - self.current_k_level = quantity_factory.zeros( - [X_DIM, Y_DIM, Z_DIM], "n/a", dtype=Int - ) - for k in range(self.current_k_level.view[:].shape[2]): - self.current_k_level.view[:, :, k] = k - - # ----------------------------------------------------------------------- - # generate saturation specific humidity tables - # ----------------------------------------------------------------------- - - self.sat_tables = get_tables(stencil_factory.backend) - - # ----------------------------------------------------------------------- - # initialize stencils - # ----------------------------------------------------------------------- - - orchestrate(obj=self, config=stencil_factory.config.dace_config) - self._init_temporaries = stencil_factory.from_dims_halo( - func=init_temporaries, - compute_dims=[X_DIM, Y_DIM, Z_DIM], - externals={ - "cpaut": self.namelist_constants.CPAUT, - }, - ) - - self._gfdl_1m_driver_preloop = stencil_factory.from_dims_halo( - func=fix_negative_values, - compute_dims=[X_DIM, Y_DIM, Z_DIM], - externals={ - "c_air": self.namelist_constants.C_AIR, - "c_vap": self.namelist_constants.C_VAP, - "p_nonhydro": self.namelist_constants.P_NONHYDRO, - "d0_vap": self.namelist_constants.D0_VAP, - "lv00": self.namelist_constants.LV00, - "latv": self.namelist_constants.LATV, - "lati": self.namelist_constants.LATI, - "lats": self.namelist_constants.LATS, - "lat2": self.namelist_constants.LAT2, - "lcp": self.namelist_constants.LCP, - "icp": self.namelist_constants.ICP, - "tcp": self.namelist_constants.TCP, - "mpdt": self.namelist_constants.MPDT, - "rdt": self.namelist_constants.RDT, - "ntimes": self.namelist_constants.NTIMES, - "dts": self.namelist_constants.DTS, - "do_sedi_w": do_sedi_w, - "cpaut": self.namelist_constants.CPAUT, - "hydrostatic": hydrostatic, - "phys_hydrostatic": phys_hydrostatic, - "fix_negative": fix_negative, - "sedi_transport": sedi_transport, - "const_vi": const_vi, - "const_vs": const_vs, - "const_vg": const_vg, - }, - ) - - self._fall_speed = stencil_factory.from_dims_halo( - func=fall_speed, - compute_dims=[X_DIM, Y_DIM, Z_DIM], - externals={ - "p_nonhydro": self.namelist_constants.P_NONHYDRO, - "const_vi": const_vi, - "const_vs": const_vs, - "const_vg": const_vg, - "vi_fac": vi_fac, - "vi_max": vi_max, - "vs_fac": vs_fac, - "vs_max": vs_max, - "vg_fac": vg_fac, - "vg_max": vg_max, - }, - ) - - self._terminal_fall = stencil_factory.from_dims_halo( - func=terminal_fall, - compute_dims=[X_DIM, Y_DIM, Z_DIM], - externals={ - "dts": self.namelist_constants.DTS, - "tau_imlt": tau_imlt, - "ql_mlt": ql_mlt, - "vi_fac": vi_fac, - "do_sedi_w": do_sedi_w, - "use_ppm": use_ppm, - "tau_smlt": tau_smlt, - "tau_g2r": tau_g2r, - "c_air": self.namelist_constants.C_AIR, - "c_vap": self.namelist_constants.C_VAP, - "d0_vap": self.namelist_constants.D0_VAP, - "lv00": self.namelist_constants.LV00, - }, - ) - - self._terminal_fall_update = stencil_factory.from_dims_halo( - func=terminal_fall_update, - compute_dims=[X_DIM, Y_DIM, Z_DIM], - ) - - self._warm_rain = stencil_factory.from_dims_halo( - func=warm_rain, - compute_dims=[X_DIM, Y_DIM, Z_DIM], - externals={ - "dts": self.namelist_constants.DTS, - "do_qa": do_qa, - "rthreshs": rthreshs, - "rthreshu": rthreshu, - "irain_f": irain_f, - "ql0_max": ql0_max, - "z_slope_liq": z_slope_liq, - "vr_fac": vr_fac, - "const_vr": const_vr, - "vr_max": vr_max, - "tau_revp": tau_revp, - "lv00": self.namelist_constants.LV00, - "d0_vap": self.namelist_constants.D0_VAP, - "c_air": self.namelist_constants.C_AIR, - "c_vap": self.namelist_constants.C_VAP, - "crevp_0": self.namelist_constants.CREVP_0, - "crevp_1": self.namelist_constants.CREVP_1, - "crevp_2": self.namelist_constants.CREVP_2, - "crevp_3": self.namelist_constants.CREVP_3, - "crevp_4": self.namelist_constants.CREVP_4, - "cracw": self.namelist_constants.CRACW, - "do_sedi_w": do_sedi_w, - "use_ppm": use_ppm, - }, - ) - - self._warm_rain_update = stencil_factory.from_dims_halo( - func=warm_rain_update, - compute_dims=[X_DIM, Y_DIM, Z_DIM], - ) - - self._icloud = stencil_factory.from_dims_halo( - func=icloud, - compute_dims=[X_DIM, Y_DIM, Z_DIM], - externals={ - "c_air": self.namelist_constants.C_AIR, - "c_vap": self.namelist_constants.C_VAP, - "dts": self.namelist_constants.DTS, - "rdts": self.namelist_constants.RDTS, - "const_vi": const_vi, - "fac_g2v": self.namelist_constants.FAC_G2V, - "fac_i2s": self.namelist_constants.FAC_I2S, - "fac_imlt": self.namelist_constants.FAC_IMLT, - "fac_frz": self.namelist_constants.FAC_FRZ, - "fac_l2v": self.namelist_constants.FAC_L2V, - "fac_s2v": self.namelist_constants.FAC_S2V, - "fac_v2s": self.namelist_constants.FAC_V2S, - "fac_v2g": self.namelist_constants.FAC_V2G, - "cgacs": self.namelist_constants.CGACS, - "csacw": self.namelist_constants.CSACW, - "csaci": self.namelist_constants.CSACI, - "cgacw": self.namelist_constants.CGACW, - "cgaci": self.namelist_constants.CGACI, - "cgfr_0": self.namelist_constants.CGFR_0, - "cgfr_1": self.namelist_constants.CGFR_1, - "csmlt_0": self.namelist_constants.CSMLT_0, - "csmlt_1": self.namelist_constants.CSMLT_1, - "csmlt_2": self.namelist_constants.CSMLT_2, - "csmlt_3": self.namelist_constants.CSMLT_3, - "csmlt_4": self.namelist_constants.CSMLT_4, - "cgmlt_0": self.namelist_constants.CGMLT_0, - "cgmlt_1": self.namelist_constants.CGMLT_1, - "cgmlt_2": self.namelist_constants.CGMLT_2, - "cgmlt_3": self.namelist_constants.CGMLT_3, - "cgmlt_4": self.namelist_constants.CGMLT_4, - "cssub_0": self.namelist_constants.CSSUB_0, - "cssub_1": self.namelist_constants.CSSUB_1, - "cssub_2": self.namelist_constants.CSSUB_2, - "cssub_3": self.namelist_constants.CSSUB_3, - "cssub_4": self.namelist_constants.CSSUB_4, - "qi0_crt": qi0_crt, - "qs0_crt": qs0_crt, - "qs_mlt": qs_mlt, - "ql_mlt": ql_mlt, - "z_slope_ice": z_slope_ice, - "lv00": self.namelist_constants.LV00, - "d0_vap": self.namelist_constants.D0_VAP, - "lat2": self.namelist_constants.LAT2, - "do_qa": do_qa, - "do_evap": do_evap, - "do_bigg": do_bigg, - "qc_crt": qc_crt, - "qi_lim": qi_lim, - "rh_inc": rh_inc, - "rh_inr": rh_inr, - "t_min": t_min, - "t_sub": t_sub, - "preciprad": preciprad, - "icloud_f": icloud_f, - }, - ) - - self._icloud_update = stencil_factory.from_dims_halo( - func=icloud_update, - compute_dims=[X_DIM, Y_DIM, Z_DIM], - ) - - self._update_tendencies = stencil_factory.from_dims_halo( - func=update_tendencies, - compute_dims=[X_DIM, Y_DIM, Z_DIM], - externals={ - "c_air": self.namelist_constants.C_AIR, - "c_vap": self.namelist_constants.C_VAP, - "rdt": self.namelist_constants.RDT, - "do_sedi_w": do_sedi_w, - "sedi_transport": sedi_transport, - "do_qa": do_qa, - }, - ) - - def __call__( - self, - qv: FloatField, - ql: FloatField, - qr: FloatField, - qi: FloatField, - qs: FloatField, - qg: FloatField, - qa: FloatField, - qn: FloatField, # NACTL + NACTI - qv_dt: FloatField, - ql_dt: FloatField, - qr_dt: FloatField, - qi_dt: FloatField, - qs_dt: FloatField, - qg_dt: FloatField, - qa_dt: FloatField, - t_dt: FloatField, - t: FloatField, - w: FloatField, - uin: FloatField, - vin: FloatField, - udt: FloatField, - vdt: FloatField, - dz: FloatField, - dp: FloatField, - area: FloatFieldIJ, - fr_land: FloatFieldIJ, - cnv_frc: FloatFieldIJ, - srf_type: FloatFieldIJ, - eis: FloatFieldIJ, - rhcrit3d: FloatField, - anv_icefall: Float, - ls_icefall: Float, - ): - # The driver modifies a number of variables (t, p, qX) but does not pass - # the changes back to the rest of the model. To replicate this behavior, - # temporary copies of these variables are used throughout the driver. - self._init_temporaries( - t, - dp, - rhcrit3d, - qv, - ql, - qi, - qr, - qs, - qg, - qa, - qn, - self.qv0, - self.ql0, - self.qr0, - self.qi0, - self.qs0, - self.qg0, - self.qa0, - self.qv1, - self.ql1, - self.qr1, - self.qi1, - self.qs1, - self.qg1, - self.qa1, - dz, - uin, - vin, - w, - area, - self.t1, - self.dp1, - self.omq, - self.den, - self.p_dry, - self.m1, - self.u1, - self.v1, - self.w1, - self.onemsig, - self.ccn, - self.c_praut, - self.rh_limited, - self.rain, - self.snow, - self.graupel, - self.ice, - self.m2_rain, - self.m2_sol, - self.revap, - self.isubl, - ) - - self._gfdl_1m_driver_preloop( - self.t1, - self.qv1, - self.ql1, - self.qr1, - self.qi1, - self.qs1, - self.qg1, - self.dp1, - ) - - for n in range(self.namelist_constants.NTIMES): - self._fall_speed( - self.ql1, - self.qi1, - self.qs1, - self.qg1, - t, - self.t1, - dz, - self.dz1, - self.den, - self.den1, - self.denfac, - self.p_dry, - self.vti, - self.vts, - self.vtg, - cnv_frc, - anv_icefall, - ls_icefall, - ) - - self._terminal_fall( - self.t1, - self.qv1, - self.ql1, - self.qr1, - self.qg1, - self.qs1, - self.qi1, - self.dz1, - self.dp1, - self.den1, - self.vtg, - self.vts, - self.vti, - self.rain1, - self.graupel1, - self.snow1, - self.ice1, - self.m1_sol, - self.w1, - self.ze, - self.zt, - self.is_frozen, - self.precip_fall, - self.melting_mask_1, - self.melting_mask_2, - self.current_k_level, - ) - - self._terminal_fall_update( - self.rain, - self.graupel, - self.snow, - self.ice, - self.rain1, - self.graupel1, - self.snow1, - self.ice1, - ) - - self._warm_rain( - self.dp1, - self.dz1, - self.t1, - self.qv1, - self.ql1, - self.qr1, - self.qi1, - self.qs1, - self.qg1, - self.qa1, - self.ccn, - self.den, - self.denfac, - self.c_praut, - self.vtr, - self.evap1, - self.m1_rain, - self.w1, - self.rh_limited, - eis, - self.onemsig, - self.rain1, - self.ze, - self.zt, - self.precip_fall, - self.sat_tables.table1, - self.sat_tables.table2, - self.sat_tables.table3, - self.sat_tables.table4, - self.sat_tables.des1, - self.sat_tables.des2, - self.sat_tables.des3, - self.sat_tables.des4, - ) - - self._warm_rain_update( - self.rain, - self.rain1, - self.evap1, - self.revap, - self.m1_rain, - self.m2_rain, - self.m1_sol, - self.m2_sol, - self.m1, - ) - - self._icloud( - self.t1, - self.p_dry, - self.dp1, - self.qv1, - self.ql1, - self.qr1, - self.qi1, - self.qs1, - self.qg1, - self.qa1, - self.den1, - self.denfac, - self.vts, - self.vtg, - self.vtr, - self.subl1, - self.rh_limited, - self.ccn, - cnv_frc, - srf_type, - self.sat_tables.table1, - self.sat_tables.table2, - self.sat_tables.table3, - self.sat_tables.table4, - self.sat_tables.des1, - self.sat_tables.des2, - self.sat_tables.des3, - self.sat_tables.des4, - ) - - self._icloud_update( - self.isubl, - self.subl1, - ) - - self._update_tendencies( - self.qv0, - self.ql0, - self.qr0, - self.qi0, - self.qs0, - self.qg0, - self.qa0, - self.qv1, - self.ql1, - self.qr1, - self.qi1, - self.qs1, - self.qg1, - self.qa1, - self.ccn, - qv_dt, - ql_dt, - qr_dt, - qi_dt, - qs_dt, - qg_dt, - qa_dt, - t, - self.t1, - t_dt, - w, - self.w1, - uin, - self.u1, - udt, - vin, - self.v1, - vdt, - dz, - dp, - self.dp1, - self.den, - self.p_dry, - area, - self.dt_moist, - fr_land, - cnv_frc, - srf_type, - eis, - self.rh_limited, - self.m1, - anv_icefall, - ls_icefall, - self.revap, - self.isubl, - self.rain, - self.snow, - self.ice, - self.graupel, - self.m2_rain, - self.m2_sol, - ) - - def check_flags( - self, - phys_hydrostatic: bool, - hydrostatic: bool, - const_vi: bool, - const_vs: bool, - const_vg: bool, - const_vr: bool, - use_ppm: bool, - use_ccn: bool, - do_qa: bool, - fix_negative: bool, - fast_set_adj: bool, - do_bigg: bool, - do_evap: bool, - do_subl: bool, - z_slope_liq: bool, - z_slope_ice: bool, - prog_ccn: bool, - preciprad: bool, - mono_prof: bool, - do_sedi_heat: bool, - sedi_transport: bool, - do_sedi_w: bool, - de_ice: bool, - mp_print: bool, - dts: Float, - ): - failed_keywords = [] - if not phys_hydrostatic: - failed_keywords.append("phys_hydrostatic") - if hydrostatic: - failed_keywords.append("hydrostatic") - if const_vi: - failed_keywords.append("const_vi") - if const_vs: - failed_keywords.append("const_vs") - if const_vg: - failed_keywords.append("const_vg") - if const_vr: - failed_keywords.append("const_vr") - if use_ppm: - failed_keywords.append("use_ppm") - if not use_ccn: - failed_keywords.append("use_ccn") - if do_qa: - failed_keywords.append("do_qa") - if not fix_negative: - failed_keywords.append("fix_negative") - if fast_set_adj: - failed_keywords.append("fast_set_adj") - if do_bigg: - failed_keywords.append("do_bigg") - if do_evap: - failed_keywords.append("do_evap") - if do_subl: - failed_keywords.append("do_subl") - if not z_slope_liq: - failed_keywords.append("z_slope_liq") - if not z_slope_ice: - failed_keywords.append("z_slope_ice") - if not prog_ccn: - failed_keywords.append("prog_ccn") - if not preciprad: - failed_keywords.append("preciprad") - if not mono_prof: - failed_keywords.append("mono_prof") - if do_sedi_heat: - failed_keywords.append("do_sedi_heat") - if not sedi_transport: - failed_keywords.append("sedi_transport") - if do_sedi_w: - failed_keywords.append("do_sedi_w") - if de_ice: - failed_keywords.append("de_ice") - if mp_print: - failed_keywords.append("mp_print") - if dts >= 300: - failed_keywords.append("dts") - - if len(failed_keywords) > 0: - raise ValueError( - "One or more namelist parameters do not meet \ - expected values. Failing parameters: ", - failed_keywords, - ) - - -class namelist_setup: - def __init__( - self, - phys_hydrostatic: bool, - hydrostatic: bool, - dt_moist: Float, - mp_time: Float, - t_min: Float, - t_sub: Float, - tau_r2g: Float, - tau_smlt: Float, - tau_g2r: Float, - dw_land: Float, - dw_ocean: Float, - vi_fac: Float, - vr_fac: Float, - vs_fac: Float, - vg_fac: Float, - ql_mlt: Float, - do_qa: bool, - fix_negative: bool, - vi_max: Float, - vs_max: Float, - vg_max: Float, - vr_max: Float, - qs_mlt: Float, - qs0_crt: Float, - qi_gen: Float, - ql0_max: Float, - qi0_max: Float, - qi0_crt: Float, - qr0_crt: Float, - fast_sat_adj: bool, - rh_inc: Float, - rh_ins: Float, - rh_inr: Float, - const_vi: bool, - const_vs: bool, - const_vg: bool, - const_vr: bool, - use_ccn: bool, - rthreshu: Float, - rthreshs: Float, - ccn_l: Float, - ccn_o: Float, - qc_crt: Float, - tau_g2v: Float, - tau_v2g: Float, - tau_s2v: Float, - tau_v2s: Float, - tau_revp: Float, - tau_frz: Float, - do_bigg: bool, - do_evap: bool, - do_subl: bool, - sat_adj0: Float, - c_piacr: Float, - tau_imlt: Float, - tau_v2l: Float, - tau_l2v: Float, - tau_i2v: Float, - tau_i2s: Float, - tau_l2r: Float, - qi_lim: Float, - ql_gen: Float, - c_paut: Float, - c_psaci: Float, - c_pgacs: Float, - c_pgaci: Float, - z_slope_liq: bool, - z_slope_ice: bool, - prog_ccn: bool, - c_cracw: Float, - alin: Float, - clin: Float, - preciprad: bool, - cld_min: Float, - use_ppm: bool, - mono_prof: bool, - do_sedi_heat: bool, - sedi_transport: bool, - do_sedi_w: bool, - de_ice: bool, - icloud_f: Float, - irain_f: Float, - mp_print: bool, - ): - # ----------------------------------------------------------------------- - # define heat capacity of dry air and water vap based on hydrostatical property - # ----------------------------------------------------------------------- - - if phys_hydrostatic or hydrostatic: - self.C_AIR = driver_constants.CP_AIR - self.C_VAP = driver_constants.CP_VAP - self.P_NONHYDRO = False - else: - self.C_AIR = driver_constants.CV_AIR - self.C_VAP = driver_constants.CV_VAP - self.P_NONHYDRO = True - self.D0_VAP = self.C_VAP - driver_constants.C_LIQ - self.LV00 = driver_constants.HLV0 - self.D0_VAP * driver_constants.T_ICE - - if hydrostatic: - self.DO_SEDI_W = False - - # ----------------------------------------------------------------------- - # define latent heat coefficient used in wet bulb and bigg mechanism - # ----------------------------------------------------------------------- - - self.LATV = driver_constants.HLV - self.LATI = driver_constants.HLF - self.LATS = self.LATV + self.LATI - self.LAT2 = self.LATS * self.LATS - - self.LCP = self.LATV / driver_constants.CP_AIR - self.ICP = self.LATI / driver_constants.CP_AIR - self.TCP = (self.LATV + self.LATI) / driver_constants.CP_AIR - - # ----------------------------------------------------------------------- - # define cloud microphysics sub time step - # ----------------------------------------------------------------------- - - self.MPDT = min(dt_moist, mp_time) - self.RDT = Float(1.0) / dt_moist - self.NTIMES = i32(np.round(dt_moist / self.MPDT)) - - # small time step: - self.DTS = dt_moist / Float(self.NTIMES) - - # ----------------------------------------------------------------------- - # calculate cloud condensation nuclei (ccn) - # the following is based on klein eq. 15 - # ----------------------------------------------------------------------- - - self.CPAUT = c_paut * Float(0.104) * driver_constants.GRAV / Float(1.717e-5) - - # ----------------------------------------------------------------------- - # define conversion scalar / factor for icloud - # ----------------------------------------------------------------------- - self.RDTS = Float(1.0) / self.DTS - self.FAC_IMLT = Float(1.0) - np.exp(-self.DTS / tau_imlt, dtype=Float) - self.FAC_I2S = Float(1.0) - np.exp(-self.DTS / tau_i2s, dtype=Float) - self.FAC_V2L = Float(1.0) - np.exp(-self.DTS / tau_v2l, dtype=Float) - self.FAC_L2V = Float(1.0) - np.exp(-self.DTS / tau_l2v, dtype=Float) - self.FAC_I2V = Float(1.0) - np.exp(-self.DTS / tau_i2v, dtype=Float) - self.FAC_S2V = Float(1.0) - np.exp(-self.DTS / tau_s2v, dtype=Float) - self.FAC_V2S = Float(1.0) - np.exp(-self.DTS / tau_v2s, dtype=Float) - self.FAC_G2V = Float(1.0) - np.exp(-self.DTS / tau_g2v, dtype=Float) - self.FAC_V2G = Float(1.0) - np.exp(-self.DTS / tau_v2g, dtype=Float) - self.FAC_FRZ = Float(1.0) - np.exp(-self.DTS / tau_frz, dtype=Float) - - # ----------------------------------------------------------------------- - # constatns from setupm - # ----------------------------------------------------------------------- - - self.CGACS = ( - driver_constants.PISQ - * driver_constants.RNZG - * driver_constants.RNZS - * driver_constants.RHOS - ) - self.CGACS = self.CGACS * c_pgacs - - self.CSACW = ( - driver_constants.PIE - * driver_constants.RNZS - * clin - * driver_constants.GAM325 - / (Float(4.0) * np.power(driver_constants.ACT[0], 0.8125, dtype=Float)) - ) - # decreasing csacw to reduce cloud water --- > snow - - self.CRACI = ( - driver_constants.PIE - * driver_constants.RNZR - * alin - * driver_constants.GAM380 - / (Float(4.0) * np.power(driver_constants.ACT[1], 0.95, dtype=Float)) - ) - self.CSACI = self.CSACW * c_psaci - - self.CGACW = ( - driver_constants.PIE - * driver_constants.RNZG - * driver_constants.GAM350 - * driver_constants.GCON - / (Float(4.0) * np.power(driver_constants.ACT[5], 0.875, dtype=Float)) - ) - - self.CGACI = self.CGACW * c_pgaci - - self.CRACW = self.CRACI # cracw = 3.27206196043822 - self.CRACW = c_cracw * self.CRACW - - self.CSSUB = np.zeros(5) - self.CGSUB = np.zeros(5) - self.CREVP = np.zeros(5) - - self.CSSUB[0] = ( - Float(2.0) - * driver_constants.PIE - * driver_constants.VDIFU - * driver_constants.TCOND - * driver_constants.RVGAS - * driver_constants.RNZS - ) - self.CGSUB[0] = ( - Float(2.0) - * driver_constants.PIE - * driver_constants.VDIFU - * driver_constants.TCOND - * driver_constants.RVGAS - * driver_constants.RNZG - ) - self.CREVP[0] = ( - Float(2.0) - * driver_constants.PIE - * driver_constants.VDIFU - * driver_constants.TCOND - * driver_constants.RVGAS - * driver_constants.RNZR - ) - self.CSSUB[1] = Float(0.78) / np.sqrt(driver_constants.ACT[0], dtype=Float) - self.CGSUB[1] = Float(0.78) / np.sqrt(driver_constants.ACT[5], dtype=Float) - self.CREVP[1] = Float(0.78) / np.sqrt(driver_constants.ACT[1], dtype=Float) - self.CSSUB[2] = ( - Float(0.31) - * driver_constants.SCM3 - * driver_constants.GAM263 - * np.sqrt(clin / driver_constants.VISK, dtype=Float) - / np.power(driver_constants.ACT[0], Float(0.65625), dtype=Float) - ) - self.CGSUB[2] = ( - Float(0.31) - * driver_constants.SCM3 - * driver_constants.GAM275 - * np.sqrt(driver_constants.GCON / driver_constants.VISK, dtype=Float) - / np.power(driver_constants.ACT[5], Float(0.6875), dtype=Float) - ) - self.CREVP[2] = ( - Float(0.31) - * driver_constants.SCM3 - * driver_constants.GAM209 - * np.sqrt(alin / driver_constants.VISK, dtype=Float) - / np.power(driver_constants.ACT[1], Float(0.725), dtype=Float) - ) - self.CSSUB[3] = driver_constants.TCOND * driver_constants.RVGAS - self.CSSUB[4] = ( - np.power(driver_constants.HLTS, 2, dtype=Float) * driver_constants.VDIFU - ) - self.CGSUB[3] = self.CSSUB[3] - self.CREVP[3] = self.CSSUB[3] - self.CGSUB[4] = self.CSSUB[4] - self.CREVP[4] = ( - np.power(driver_constants.HLTC, 2, dtype=Float) * driver_constants.VDIFU - ) - - self.CGFR_0 = ( - Float(20.0e2) - * driver_constants.PISQ - * driver_constants.RNZR - * driver_constants.RHOR - / np.power(driver_constants.ACT[1], Float(1.75), dtype=Float) - ) - self.CGFR_1 = Float(0.66) - - self.CSSUB_0 = self.CSSUB[0] - self.CSSUB_1 = self.CSSUB[1] - self.CSSUB_2 = self.CSSUB[2] - self.CSSUB_3 = self.CSSUB[3] - self.CSSUB_4 = self.CSSUB[4] - - self.CGSUB_0 = self.CGSUB[0] - self.CGSUB_1 = self.CGSUB[1] - self.CGSUB_2 = self.CGSUB[2] - self.CGSUB_3 = self.CGSUB[3] - self.CGSUB_4 = self.CGSUB[4] - - self.CREVP_0 = self.CREVP[0] - self.CREVP_1 = self.CREVP[1] - self.CREVP_2 = self.CREVP[2] - self.CREVP_3 = self.CREVP[3] - self.CREVP_4 = self.CREVP[4] - - # smlt: five constants (lin et al. 1983) - - self.CSMLT = np.zeros(5) - self.CSMLT[0] = ( - Float(2.0) - * driver_constants.PIE - * driver_constants.TCOND - * driver_constants.RNZS - / driver_constants.HLTF - ) - self.CSMLT[1] = ( - Float(2.0) - * driver_constants.PIE - * driver_constants.VDIFU - * driver_constants.RNZS - * driver_constants.HLTC - / driver_constants.HLTF - ) - self.CSMLT[2] = self.CSSUB[1] - self.CSMLT[3] = self.CSSUB[2] - self.CSMLT[4] = driver_constants.CH2O / driver_constants.HLTF - - self.CSMLT_0 = self.CSMLT[0] - self.CSMLT_1 = self.CSMLT[1] - self.CSMLT_2 = self.CSMLT[2] - self.CSMLT_3 = self.CSMLT[3] - self.CSMLT_4 = self.CSMLT[4] - - # gmlt: five constants - - self.CGMLT = np.zeros(5) - self.CGMLT[0] = ( - Float(2.0) - * driver_constants.PIE - * driver_constants.TCOND - * driver_constants.RNZG - / driver_constants.HLTF - ) - self.CGMLT[1] = ( - Float(2.0) - * driver_constants.PIE - * driver_constants.VDIFU - * driver_constants.RNZG - * driver_constants.HLTC - / driver_constants.HLTF - ) - self.CGMLT[2] = self.CGSUB[1] - self.CGMLT[3] = self.CGSUB[2] - self.CGMLT[4] = driver_constants.CH2O / driver_constants.HLTF - - self.CGMLT_0 = self.CGMLT[0] - self.CGMLT_1 = self.CGMLT[1] - self.CGMLT_2 = self.CGMLT[2] - self.CGMLT_3 = self.CGMLT[3] - self.CGMLT_4 = self.CGMLT[4] diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/config.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/config.py new file mode 100644 index 000000000..6ec3157eb --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/config.py @@ -0,0 +1,175 @@ +import numpy as np +from gt4py.cartesian.gtscript import i32 +from ndsl.dsl.typing import Float + + +class config: + def __init__( + self, + phys_hydrostatic: bool, + hydrostatic: bool, + dt_moist: Float, + mp_time: Float, + t_min: Float, + t_sub: Float, + tau_r2g: Float, + tau_smlt: Float, + tau_g2r: Float, + dw_land: Float, + dw_ocean: Float, + vi_fac: Float, + vr_fac: Float, + vs_fac: Float, + vg_fac: Float, + ql_mlt: Float, + do_qa: bool, + fix_negative: bool, + vi_max: Float, + vs_max: Float, + vg_max: Float, + vr_max: Float, + qs_mlt: Float, + qs0_crt: Float, + qi_gen: Float, + ql0_max: Float, + qi0_max: Float, + qi0_crt: Float, + qr0_crt: Float, + fast_sat_adj: bool, + rh_inc: Float, + rh_ins: Float, + rh_inr: Float, + const_vi: bool, + const_vs: bool, + const_vg: bool, + const_vr: bool, + use_ccn: bool, + rthreshu: Float, + rthreshs: Float, + ccn_l: Float, + ccn_o: Float, + qc_crt: Float, + tau_g2v: Float, + tau_v2g: Float, + tau_s2v: Float, + tau_v2s: Float, + tau_revp: Float, + tau_frz: Float, + do_bigg: bool, + do_evap: bool, + do_subl: bool, + sat_adj0: Float, + c_piacr: Float, + tau_imlt: Float, + tau_v2l: Float, + tau_l2v: Float, + tau_i2v: Float, + tau_i2s: Float, + tau_l2r: Float, + qi_lim: Float, + ql_gen: Float, + c_paut: Float, + c_psaci: Float, + c_pgacs: Float, + c_pgaci: Float, + z_slope_liq: bool, + z_slope_ice: bool, + prog_ccn: bool, + c_cracw: Float, + alin: Float, + clin: Float, + preciprad: bool, + cld_min: Float, + use_ppm: bool, + mono_prof: bool, + do_sedi_heat: bool, + sedi_transport: bool, + do_sedi_w: bool, + de_ice: bool, + icloud_f: Float, + irain_f: Float, + mp_print: bool, + ): + self.phys_hydrostatic = phys_hydrostatic + self.hydrostatic = hydrostatic + self.dt_moist = dt_moist + self.mp_time = mp_time + self.t_min = t_min + self.t_sub = t_sub + self.tau_r2g = tau_r2g + self.tau_smlt = tau_smlt + self.tau_g2r = tau_g2r + self.dw_land = dw_land + self.dw_ocean = dw_ocean + self.vi_fac = vi_fac + self.vr_fac = vr_fac + self.vs_fac = vs_fac + self.vg_fac = vg_fac + self.ql_mlt = ql_mlt + self.do_qa = do_qa + self.fix_negative = fix_negative + self.vi_max = vi_max + self.vs_max = vs_max + self.vg_max = vg_max + self.vr_max = vr_max + self.qs_mlt = qs_mlt + self.qs0_crt = qs0_crt + self.qi_gen = qi_gen + self.ql0_max = ql0_max + self.qi0_max = qi0_max + self.qi0_crt = qi0_crt + self.qr0_crt = qr0_crt + self.fast_sat_adj = fast_sat_adj + self.rh_inc = rh_inc + self.rh_ins = rh_ins + self.rh_inr = rh_inr + self.const_vi = const_vi + self.const_vs = const_vs + self.const_vg = const_vg + self.const_vr = const_vr + self.use_ccn = use_ccn + self.rthreshu = rthreshu + self.rthreshs = rthreshs + self.ccn_l = ccn_l + self.ccn_o = ccn_o + self.qc_crt = qc_crt + self.tau_g2v = tau_g2v + self.tau_v2g = tau_v2g + self.tau_s2v = tau_s2v + self.tau_v2s = tau_v2s + self.tau_revp = tau_revp + self.tau_frz = tau_frz + self.do_bigg = do_bigg + self.do_evap = do_evap + self.do_subl = do_subl + self.sat_adj0 = sat_adj0 + self.c_piacr = c_piacr + self.tau_imlt = tau_imlt + self.tau_v2l = tau_v2l + self.tau_l2v = tau_l2v + self.tau_i2v = tau_i2v + self.tau_i2s = tau_i2s + self.tau_l2r = tau_l2r + self.qi_lim = qi_lim + self.ql_gen = ql_gen + self.c_paut = c_paut + self.c_psaci = c_psaci + self.c_pgacs = c_pgacs + self.c_pgaci = c_pgaci + self.z_slope_liq = z_slope_liq + self.z_slope_ice = z_slope_ice + self.prog_ccn = prog_ccn + self.c_cracw = c_cracw + self.alin = alin + self.clin = clin + self.preciprad = preciprad + self.cld_min = cld_min + self.use_ppm = use_ppm + self.mono_prof = mono_prof + self.do_sedi_heat = do_sedi_heat + self.sedi_transport = sedi_transport + self.do_sedi_w = do_sedi_w + self.de_ice = de_ice + self.icloud_f = icloud_f + self.irain_f = irain_f + self.mp_print = mp_print diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/GFDL_1M_driver_core.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/connections.py similarity index 87% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/GFDL_1M_driver_core.py rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/connections.py index 772fe6b2f..afd6bcc48 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/GFDL_1M_driver_core.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/connections.py @@ -12,12 +12,12 @@ sqrt, ) -import pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver_constants as driver_constants +import pyMoist.GFDL_1M.GFDL_1M_driver.constants as constants from ndsl.dsl.typing import Float, FloatField, FloatFieldIJ from pyMoist.shared_generic_math import sigma -GlobalTable_driver_qsat = gtscript.GlobalTable[(Float, (driver_constants.LENGTH))] +GlobalTable_driver_qsat = gtscript.GlobalTable[(Float, (constants.LENGTH))] def init_temporaries( @@ -119,8 +119,8 @@ def init_temporaries( qa0 = qa qa1 = qa - den = -dp1 / (driver_constants.GRAV * dz) # density of dry air - p_dry = den * driver_constants.RDGAS * t # dry air pressure + den = -dp1 / (constants.GRAV * dz) # density of dry air + p_dry = den * constants.RDGAS * t # dry air pressure # ----------------------------------------------------------------------- # for sedi_momentum @@ -133,7 +133,7 @@ def init_temporaries( # ccn needs units #/m^3 ccn = qn - c_praut = cpaut * (ccn * driver_constants.RHOR) ** (-1.0 / 3.0) + c_praut = cpaut * (ccn * constants.RHOR) ** (-1.0 / 3.0) # Reset precipitation aggregates to zero m2_rain = 0 @@ -179,11 +179,11 @@ def fix_negative_core( cvm = ( c_air + qv * c_vap - + (qr + ql) * driver_constants.C_LIQ - + (qi + qs + qg) * driver_constants.C_ICE + + (qr + ql) * constants.C_LIQ + + (qi + qs + qg) * constants.C_ICE ) lcpk = (lv00 + d0_vap * t) / cvm - icpk = (driver_constants.LI00 + driver_constants.DC_ICE * t) / cvm + icpk = (constants.LI00 + constants.DC_ICE * t) / cvm # ----------------------------------------------------------------------- # ice phase: @@ -293,19 +293,19 @@ def fall_speed_core( vs_max, ) - rhof = sqrt(min(10.0, driver_constants.SFCRHO / den)) + rhof = sqrt(min(10.0, constants.SFCRHO / den)) if const_vi == True: # noqa vti = vi_fac else: - if qi < driver_constants.THI: - vti = driver_constants.VF_MIN + if qi < constants.THI: + vti = constants.VF_MIN else: # ----------------------------------------------------------------------- # ice: # ----------------------------------------------------------------------- vi1 = 0.01 * vi_fac - tc = t - driver_constants.TICE # deg C + tc = t - constants.TICE # deg C IWC = qi * den * 1.0e3 # Units are g/m3 # ----------------------------------------------------------------------- # use deng and mace (2008, grl) @@ -314,27 +314,27 @@ def fall_speed_core( viLSC = lsc_icefall * 10.0 ** ( log10(IWC) * ( - tc * (driver_constants.AAL * tc + driver_constants.BBL) - + driver_constants.CCL + tc * (constants.AAL * tc + constants.BBL) + + constants.CCL ) - + driver_constants.DDL * tc - + driver_constants.EEL + + constants.DDL * tc + + constants.EEL ) viCNV = anv_icefall * 10.0 ** ( log10(IWC) * ( - tc * (driver_constants.AAC * tc + driver_constants.BBC) - + driver_constants.CCC + tc * (constants.AAC * tc + constants.BBC) + + constants.CCC ) - + driver_constants.DDC * tc - + driver_constants.EEC + + constants.DDC * tc + + constants.EEC ) # Combine vti = viLSC * (1.0 - cnv_frc) + viCNV * (cnv_frc) # Update units from cm/s to m/s vti = vi1 * vti # Limits - vti = min(vi_max, max(driver_constants.VF_MIN, vti)) + vti = min(vi_max, max(constants.VF_MIN, vti)) # ----------------------------------------------------------------------- # snow: @@ -343,16 +343,16 @@ def fall_speed_core( if const_vs == True: # noqa vts = vs_fac # 1. ifs_2016 else: - if qs < driver_constants.THS: - vts = driver_constants.VF_MIN + if qs < constants.THS: + vts = constants.VF_MIN else: vts = ( vs_fac - * driver_constants.VCONS + * constants.VCONS * rhof - * exp(0.0625 * log(qs * den / driver_constants.NORMS)) + * exp(0.0625 * log(qs * den / constants.NORMS)) ) - vts = min(vs_max, max(driver_constants.VF_MIN, vts)) + vts = min(vs_max, max(constants.VF_MIN, vts)) # ----------------------------------------------------------------------- # graupel: @@ -361,16 +361,16 @@ def fall_speed_core( if const_vg == True: # noqa vtg = vg_fac # 2. else: - if qg < driver_constants.THG: - vtg = driver_constants.VF_MIN + if qg < constants.THG: + vtg = constants.VF_MIN else: vtg = ( vg_fac - * driver_constants.VCONG + * constants.VCONG * rhof - * sqrt(sqrt(sqrt(qg * den / driver_constants.NORMG))) + * sqrt(sqrt(sqrt(qg * den / constants.NORMG))) ) - vtg = min(vg_max, max(driver_constants.VF_MIN, vtg)) + vtg = min(vg_max, max(constants.VF_MIN, vtg)) return vti, vts, vtg @@ -408,11 +408,11 @@ def fall_speed( if p_nonhydro: dz1 = dz den1 = den # dry air density remains the same - denfac = sqrt(driver_constants.SFCRHO / den1) + denfac = sqrt(constants.SFCRHO / den1) else: dz1 = dz * t1 / t # hydrostatic balance den1 = den * dz / dz1 - denfac = sqrt(driver_constants.SFCRHO / den1) + denfac = sqrt(constants.SFCRHO / den1) vti, vts, vtg = fall_speed_core( p_dry, cnv_frc, anv_icefall, lsc_icefall, den1, qs1, qi1, qg1, ql1, t1 @@ -591,22 +591,22 @@ def update_tendencies( cvm = ( c_air + qv1 * c_vap - + (qr1 + ql1) * driver_constants.C_LIQ - + (qi1 + qs1 + qg1) * driver_constants.C_ICE + + (qr1 + ql1) * constants.C_LIQ + + (qi1 + qs1 + qg1) * constants.C_ICE ) - t_dt = t_dt + rdt * (t1 - t) * cvm / driver_constants.CP_AIR + t_dt = t_dt + rdt * (t1 - t) * cvm / constants.CP_AIR # ----------------------------------------------------------------------- # update cloud fraction tendency # ----------------------------------------------------------------------- if do_qa == False: # noqa qa_dt = qa_dt + rdt * ( - qa0 * sqrt((qi1 + ql1) / max(qi0 + ql0, driver_constants.QCMIN)) - qa0 + qa0 * sqrt((qi1 + ql1) / max(qi0 + ql0, constants.QCMIN)) - qa0 ) # New Cloud - Old CloudCloud with computation(FORWARD), interval(0, 1): # convert to mm / day - conversion_factor = 86400.0 * rdt * driver_constants.RGRAV + conversion_factor = 86400.0 * rdt * constants.RGRAV rain = rain * conversion_factor snow = snow * conversion_factor ice = ice * conversion_factor diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/GFDL_1M_driver_constants.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/constants.py similarity index 98% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/GFDL_1M_driver_constants.py rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/constants.py index 5c0defb55..eca33ab22 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/GFDL_1M_driver_constants.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/constants.py @@ -6,9 +6,8 @@ """ import numpy as np -from gt4py.cartesian.gtscript import i32 -from ndsl.dsl.typing import Float +from ndsl.dsl.typing import Float, Int # constants from gfdl_1m_driver module @@ -230,4 +229,4 @@ # q table constants -LENGTH = i32(2621) +LENGTH = Int(2621) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/driver.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/driver.py new file mode 100644 index 000000000..71ebc0732 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/driver.py @@ -0,0 +1,565 @@ +"""GFDL_1M driver""" + +import numpy as np +from gt4py.cartesian.gtscript import i32 + +import pyMoist.GFDL_1M.GFDL_1M_driver.constants as constants +from ndsl import QuantityFactory, StencilFactory, orchestrate +from ndsl.constants import X_DIM, Y_DIM, Z_DIM, Z_INTERFACE_DIM +from ndsl.dsl.typing import Float, FloatField, FloatFieldIJ, Int +from pyMoist.GFDL_1M.GFDL_1M_driver.connections import ( + fall_speed, + fix_negative_values, + icloud_update, + init_temporaries, + terminal_fall_update, + update_tendencies, + warm_rain_update, +) +from pyMoist.GFDL_1M.GFDL_1M_driver.sat_tables import get_tables +from pyMoist.GFDL_1M.GFDL_1M_driver.icloud import icloud +from pyMoist.GFDL_1M.GFDL_1M_driver.terminal_fall.terminal_fall import terminal_fall +from pyMoist.GFDL_1M.GFDL_1M_driver.warm_rain import warm_rain +from pyMoist.GFDL_1M.GFDL_1M_driver.config import config +from pyMoist.GFDL_1M.GFDL_1M_driver.support import ( + check_flags, + config_constants, + outputs, + temporaries, + masks, +) + + +class driver: + def __init__( + self, + stencil_factory: StencilFactory, + quantity_factory: QuantityFactory, + GFDL_1M_config: config, + ): + + self.config_dependent_constants = config_constants(GFDL_1M_config) + + # Check values for untested code paths + check_flags( + GFDL_1M_config, + self.config_dependent_constants.DTS, + ) + + # ----------------------------------------------------------------------- + # initialize precipitation outputs + # ----------------------------------------------------------------------- + + self.outputs = outputs(quantity_factory) + + # ----------------------------------------------------------------------- + # initialize temporaries + # ----------------------------------------------------------------------- + + self.temporaries = temporaries(quantity_factory) + + # ----------------------------------------------------------------------- + # initialize masks + # ----------------------------------------------------------------------- + + self.masks = masks(quantity_factory) + + # ----------------------------------------------------------------------- + # generate saturation specific humidity tables + # ----------------------------------------------------------------------- + + self.sat_tables = get_tables(stencil_factory.backend) + + # ----------------------------------------------------------------------- + # initialize stencils + # ----------------------------------------------------------------------- + + orchestrate(obj=self, config=stencil_factory.config.dace_config) + self._init_temporaries = stencil_factory.from_dims_halo( + func=init_temporaries, + compute_dims=[X_DIM, Y_DIM, Z_DIM], + externals={ + "cpaut": self.config_dependent_constants.CPAUT, + }, + ) + + self._gfdl_1m_driver_preloop = stencil_factory.from_dims_halo( + func=fix_negative_values, + compute_dims=[X_DIM, Y_DIM, Z_DIM], + externals={ + "c_air": self.config_dependent_constants.C_AIR, + "c_vap": self.config_dependent_constants.C_VAP, + "p_nonhydro": self.config_dependent_constants.P_NONHYDRO, + "d0_vap": self.config_dependent_constants.D0_VAP, + "lv00": self.config_dependent_constants.LV00, + "latv": self.config_dependent_constants.LATV, + "lati": self.config_dependent_constants.LATI, + "lats": self.config_dependent_constants.LATS, + "lat2": self.config_dependent_constants.LAT2, + "lcp": self.config_dependent_constants.LCP, + "icp": self.config_dependent_constants.ICP, + "tcp": self.config_dependent_constants.TCP, + "mpdt": self.config_dependent_constants.MPDT, + "rdt": self.config_dependent_constants.RDT, + "ntimes": self.config_dependent_constants.NTIMES, + "dts": self.config_dependent_constants.DTS, + "do_sedi_w": GFDL_1M_config.do_sedi_w, + "cpaut": self.config_dependent_constants.CPAUT, + "hydrostatic": GFDL_1M_config.hydrostatic, + "phys_hydrostatic": GFDL_1M_config.phys_hydrostatic, + "fix_negative": GFDL_1M_config.fix_negative, + "sedi_transport": GFDL_1M_config.sedi_transport, + "const_vi": GFDL_1M_config.const_vi, + "const_vs": GFDL_1M_config.const_vs, + "const_vg": GFDL_1M_config.const_vg, + }, + ) + + self._fall_speed = stencil_factory.from_dims_halo( + func=fall_speed, + compute_dims=[X_DIM, Y_DIM, Z_DIM], + externals={ + "p_nonhydro": self.config_dependent_constants.P_NONHYDRO, + "const_vi": GFDL_1M_config.const_vi, + "const_vs": GFDL_1M_config.const_vs, + "const_vg": GFDL_1M_config.const_vg, + "vi_fac": GFDL_1M_config.vi_fac, + "vi_max": GFDL_1M_config.vi_max, + "vs_fac": GFDL_1M_config.vs_fac, + "vs_max": GFDL_1M_config.vs_max, + "vg_fac": GFDL_1M_config.vg_fac, + "vg_max": GFDL_1M_config.vg_max, + }, + ) + + self._terminal_fall = terminal_fall( + stencil_factory, + quantity_factory, + GFDL_1M_config, + self.config_dependent_constants, + ) + + self._terminal_fall_update = stencil_factory.from_dims_halo( + func=terminal_fall_update, + compute_dims=[X_DIM, Y_DIM, Z_DIM], + ) + + self._warm_rain = stencil_factory.from_dims_halo( + func=warm_rain, + compute_dims=[X_DIM, Y_DIM, Z_DIM], + externals={ + "dts": self.config_dependent_constants.DTS, + "do_qa": GFDL_1M_config.do_qa, + "rthreshs": GFDL_1M_config.rthreshs, + "rthreshu": GFDL_1M_config.rthreshu, + "irain_f": GFDL_1M_config.irain_f, + "ql0_max": GFDL_1M_config.ql0_max, + "z_slope_liq": GFDL_1M_config.z_slope_liq, + "vr_fac": GFDL_1M_config.vr_fac, + "const_vr": GFDL_1M_config.const_vr, + "vr_max": GFDL_1M_config.vr_max, + "tau_revp": GFDL_1M_config.tau_revp, + "lv00": self.config_dependent_constants.LV00, + "d0_vap": self.config_dependent_constants.D0_VAP, + "c_air": self.config_dependent_constants.C_AIR, + "c_vap": self.config_dependent_constants.C_VAP, + "crevp_0": self.config_dependent_constants.CREVP_0, + "crevp_1": self.config_dependent_constants.CREVP_1, + "crevp_2": self.config_dependent_constants.CREVP_2, + "crevp_3": self.config_dependent_constants.CREVP_3, + "crevp_4": self.config_dependent_constants.CREVP_4, + "cracw": self.config_dependent_constants.CRACW, + "do_sedi_w": GFDL_1M_config.do_sedi_w, + "use_ppm": GFDL_1M_config.use_ppm, + }, + ) + + self._warm_rain_update = stencil_factory.from_dims_halo( + func=warm_rain_update, + compute_dims=[X_DIM, Y_DIM, Z_DIM], + ) + + self._icloud = stencil_factory.from_dims_halo( + func=icloud, + compute_dims=[X_DIM, Y_DIM, Z_DIM], + externals={ + "c_air": self.config_dependent_constants.C_AIR, + "c_vap": self.config_dependent_constants.C_VAP, + "dts": self.config_dependent_constants.DTS, + "rdts": self.config_dependent_constants.RDTS, + "const_vi": GFDL_1M_config.const_vi, + "fac_g2v": self.config_dependent_constants.FAC_G2V, + "fac_i2s": self.config_dependent_constants.FAC_I2S, + "fac_imlt": self.config_dependent_constants.FAC_IMLT, + "fac_frz": self.config_dependent_constants.FAC_FRZ, + "fac_l2v": self.config_dependent_constants.FAC_L2V, + "fac_s2v": self.config_dependent_constants.FAC_S2V, + "fac_v2s": self.config_dependent_constants.FAC_V2S, + "fac_v2g": self.config_dependent_constants.FAC_V2G, + "cgacs": self.config_dependent_constants.CGACS, + "csacw": self.config_dependent_constants.CSACW, + "csaci": self.config_dependent_constants.CSACI, + "cgacw": self.config_dependent_constants.CGACW, + "cgaci": self.config_dependent_constants.CGACI, + "cgfr_0": self.config_dependent_constants.CGFR_0, + "cgfr_1": self.config_dependent_constants.CGFR_1, + "csmlt_0": self.config_dependent_constants.CSMLT_0, + "csmlt_1": self.config_dependent_constants.CSMLT_1, + "csmlt_2": self.config_dependent_constants.CSMLT_2, + "csmlt_3": self.config_dependent_constants.CSMLT_3, + "csmlt_4": self.config_dependent_constants.CSMLT_4, + "cgmlt_0": self.config_dependent_constants.CGMLT_0, + "cgmlt_1": self.config_dependent_constants.CGMLT_1, + "cgmlt_2": self.config_dependent_constants.CGMLT_2, + "cgmlt_3": self.config_dependent_constants.CGMLT_3, + "cgmlt_4": self.config_dependent_constants.CGMLT_4, + "cssub_0": self.config_dependent_constants.CSSUB_0, + "cssub_1": self.config_dependent_constants.CSSUB_1, + "cssub_2": self.config_dependent_constants.CSSUB_2, + "cssub_3": self.config_dependent_constants.CSSUB_3, + "cssub_4": self.config_dependent_constants.CSSUB_4, + "qi0_crt": GFDL_1M_config.qi0_crt, + "qs0_crt": GFDL_1M_config.qs0_crt, + "qs_mlt": GFDL_1M_config.qs_mlt, + "ql_mlt": GFDL_1M_config.ql_mlt, + "z_slope_ice": GFDL_1M_config.z_slope_ice, + "lv00": self.config_dependent_constants.LV00, + "d0_vap": self.config_dependent_constants.D0_VAP, + "lat2": self.config_dependent_constants.LAT2, + "do_qa": GFDL_1M_config.do_qa, + "do_evap": GFDL_1M_config.do_evap, + "do_bigg": GFDL_1M_config.do_bigg, + "qc_crt": GFDL_1M_config.qc_crt, + "qi_lim": GFDL_1M_config.qi_lim, + "rh_inc": GFDL_1M_config.rh_inc, + "rh_inr": GFDL_1M_config.rh_inr, + "t_min": GFDL_1M_config.t_min, + "t_sub": GFDL_1M_config.t_sub, + "preciprad": GFDL_1M_config.preciprad, + "icloud_f": GFDL_1M_config.icloud_f, + }, + ) + + self._icloud_update = stencil_factory.from_dims_halo( + func=icloud_update, + compute_dims=[X_DIM, Y_DIM, Z_DIM], + ) + + self._update_tendencies = stencil_factory.from_dims_halo( + func=update_tendencies, + compute_dims=[X_DIM, Y_DIM, Z_DIM], + externals={ + "c_air": self.config_dependent_constants.C_AIR, + "c_vap": self.config_dependent_constants.C_VAP, + "rdt": self.config_dependent_constants.RDT, + "do_sedi_w": GFDL_1M_config.do_sedi_w, + "sedi_transport": GFDL_1M_config.sedi_transport, + "do_qa": GFDL_1M_config.do_qa, + }, + ) + + def __call__( + self, + GFDL_1M_config: config, + qv: FloatField, + ql: FloatField, + qr: FloatField, + qi: FloatField, + qs: FloatField, + qg: FloatField, + qa: FloatField, + qn: FloatField, # NACTL + NACTI + qv_dt: FloatField, + ql_dt: FloatField, + qr_dt: FloatField, + qi_dt: FloatField, + qs_dt: FloatField, + qg_dt: FloatField, + qa_dt: FloatField, + t_dt: FloatField, + t: FloatField, + w: FloatField, + uin: FloatField, + vin: FloatField, + udt: FloatField, + vdt: FloatField, + dz: FloatField, + dp: FloatField, + area: FloatFieldIJ, + fr_land: FloatFieldIJ, + cnv_frc: FloatFieldIJ, + srf_type: FloatFieldIJ, + eis: FloatFieldIJ, + rhcrit3d: FloatField, + anv_icefall: Float, + ls_icefall: Float, + ): + # The driver modifies a number of variables (t, p, qX) but does not pass + # the changes back to the rest of the model. To replicate this behavior, + # temporary copies of these variables are used throughout the driver. + self._init_temporaries( + t, + dp, + rhcrit3d, + qv, + ql, + qi, + qr, + qs, + qg, + qa, + qn, + self.temporaries.qv0, + self.temporaries.ql0, + self.temporaries.qr0, + self.temporaries.qi0, + self.temporaries.qs0, + self.temporaries.qg0, + self.temporaries.qa0, + self.temporaries.qv1, + self.temporaries.ql1, + self.temporaries.qr1, + self.temporaries.qi1, + self.temporaries.qs1, + self.temporaries.qg1, + self.temporaries.qa1, + dz, + uin, + vin, + w, + area, + self.temporaries.t1, + self.temporaries.dp1, + self.temporaries.omq, + self.temporaries.den, + self.temporaries.p_dry, + self.temporaries.m1, + self.temporaries.u1, + self.temporaries.v1, + self.temporaries.w1, + self.temporaries.onemsig, + self.temporaries.ccn, + self.temporaries.c_praut, + self.temporaries.rh_limited, + self.outputs.rain, + self.outputs.snow, + self.outputs.graupel, + self.outputs.ice, + self.outputs.m2_rain, + self.outputs.m2_sol, + self.outputs.revap, + self.outputs.isubl, + ) + + self._gfdl_1m_driver_preloop( + self.temporaries.t1, + self.temporaries.qv1, + self.temporaries.ql1, + self.temporaries.qr1, + self.temporaries.qi1, + self.temporaries.qs1, + self.temporaries.qg1, + self.temporaries.dp1, + ) + + for n in range(self.config_dependent_constants.NTIMES): + self._fall_speed( + self.temporaries.ql1, + self.temporaries.qi1, + self.temporaries.qs1, + self.temporaries.qg1, + t, + self.temporaries.t1, + dz, + self.temporaries.dz1, + self.temporaries.den, + self.temporaries.den1, + self.temporaries.denfac, + self.temporaries.p_dry, + self.temporaries.vti, + self.temporaries.vts, + self.temporaries.vtg, + cnv_frc, + anv_icefall, + ls_icefall, + ) + + self._terminal_fall( + self.temporaries.t1, + self.temporaries.qv1, + self.temporaries.ql1, + self.temporaries.qr1, + self.temporaries.qg1, + self.temporaries.qs1, + self.temporaries.qi1, + self.temporaries.dz1, + self.temporaries.dp1, + self.temporaries.den1, + self.temporaries.vtg, + self.temporaries.vts, + self.temporaries.vti, + self.temporaries.rain1, + self.temporaries.graupel1, + self.temporaries.snow1, + self.temporaries.ice1, + self.temporaries.m1_sol, + self.temporaries.w1, + self.temporaries.ze, + self.temporaries.zt, + self.masks.is_frozen, + self.masks.precip_fall, + ) + + self._terminal_fall_update( + self.outputs.rain, + self.outputs.graupel, + self.outputs.snow, + self.outputs.ice, + self.temporaries.rain1, + self.temporaries.graupel1, + self.temporaries.snow1, + self.temporaries.ice1, + ) + + self._warm_rain( + self.temporaries.dp1, + self.temporaries.dz1, + self.temporaries.t1, + self.temporaries.qv1, + self.temporaries.ql1, + self.temporaries.qr1, + self.temporaries.qi1, + self.temporaries.qs1, + self.temporaries.qg1, + self.temporaries.qa1, + self.temporaries.ccn, + self.temporaries.den, + self.temporaries.denfac, + self.temporaries.c_praut, + self.temporaries.vtr, + self.temporaries.evap1, + self.temporaries.m1_rain, + self.temporaries.w1, + self.temporaries.rh_limited, + eis, + self.temporaries.onemsig, + self.temporaries.rain1, + self.temporaries.ze, + self.temporaries.zt, + self.masks.precip_fall, + self.sat_tables.table1, + self.sat_tables.table2, + self.sat_tables.table3, + self.sat_tables.table4, + self.sat_tables.des1, + self.sat_tables.des2, + self.sat_tables.des3, + self.sat_tables.des4, + ) + + self._warm_rain_update( + self.outputs.rain, + self.temporaries.rain1, + self.temporaries.evap1, + self.outputs.revap, + self.temporaries.m1_rain, + self.outputs.m2_rain, + self.temporaries.m1_sol, + self.outputs.m2_sol, + self.temporaries.m1, + ) + + self._icloud( + self.temporaries.t1, + self.temporaries.p_dry, + self.temporaries.dp1, + self.temporaries.qv1, + self.temporaries.ql1, + self.temporaries.qr1, + self.temporaries.qi1, + self.temporaries.qs1, + self.temporaries.qg1, + self.temporaries.qa1, + self.temporaries.den1, + self.temporaries.denfac, + self.temporaries.vts, + self.temporaries.vtg, + self.temporaries.vtr, + self.temporaries.subl1, + self.temporaries.rh_limited, + self.temporaries.ccn, + cnv_frc, + srf_type, + self.sat_tables.table1, + self.sat_tables.table2, + self.sat_tables.table3, + self.sat_tables.table4, + self.sat_tables.des1, + self.sat_tables.des2, + self.sat_tables.des3, + self.sat_tables.des4, + ) + + self._icloud_update( + self.outputs.isubl, + self.temporaries.subl1, + ) + + self._update_tendencies( + self.temporaries.qv0, + self.temporaries.ql0, + self.temporaries.qr0, + self.temporaries.qi0, + self.temporaries.qs0, + self.temporaries.qg0, + self.temporaries.qa0, + self.temporaries.qv1, + self.temporaries.ql1, + self.temporaries.qr1, + self.temporaries.qi1, + self.temporaries.qs1, + self.temporaries.qg1, + self.temporaries.qa1, + self.temporaries.ccn, + qv_dt, + ql_dt, + qr_dt, + qi_dt, + qs_dt, + qg_dt, + qa_dt, + t, + self.temporaries.t1, + t_dt, + w, + self.temporaries.w1, + uin, + self.temporaries.u1, + udt, + vin, + self.temporaries.v1, + vdt, + dz, + dp, + self.temporaries.dp1, + self.temporaries.den, + self.temporaries.p_dry, + area, + GFDL_1M_config.dt_moist, + fr_land, + cnv_frc, + srf_type, + eis, + self.temporaries.rh_limited, + self.temporaries.m1, + anv_icefall, + ls_icefall, + self.outputs.revap, + self.outputs.isubl, + self.outputs.rain, + self.outputs.snow, + self.outputs.ice, + self.outputs.graupel, + self.outputs.m2_rain, + self.outputs.m2_sol, + ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/icloud.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/icloud.py index 16a5ac1df..086dd06b0 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/icloud.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/icloud.py @@ -12,12 +12,12 @@ trunc, ) -import pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver_constants as driver_constants +import pyMoist.GFDL_1M.GFDL_1M_driver.constants as constants from ndsl.dsl.typing import Float, FloatField, FloatFieldIJ from pyMoist.shared_incloud_processes import ice_fraction -GlobalTable_driver_qsat = gtscript.GlobalTable[(Float, (int(driver_constants.LENGTH)))] +GlobalTable_driver_qsat = gtscript.GlobalTable[(Float, (int(constants.LENGTH)))] @gtscript.function @@ -58,14 +58,14 @@ def icloud_melt_freeze( # define heat capacity and latent heat coefficient # ----------------------------------------------------------------------- - lhi = driver_constants.LI00 + driver_constants.DC_ICE * t + lhi = constants.LI00 + constants.DC_ICE * t q_liq = ql + qr q_sol = qi + qs + qg cvm = ( c_air + qv * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) icpk = lhi / cvm @@ -75,7 +75,7 @@ def icloud_melt_freeze( melt = fac_imlt * max(0.0, newliq - ql) frez = fac_frz * max(0.0, newice - qi) - if melt > 0.0 and t > driver_constants.TICE and qi > driver_constants.QCMIN: + if melt > 0.0 and t > constants.TICE and qi > constants.QCMIN: # ----------------------------------------------------------------------- # pimlt: melting of cloud ice # ----------------------------------------------------------------------- @@ -92,7 +92,7 @@ def icloud_melt_freeze( 1.0, qa * max(qi + ql - melt + tmp, 0.0) - / max(qi + ql, driver_constants.QCMIN), + / max(qi + ql, constants.QCMIN), ), ) @@ -104,11 +104,11 @@ def icloud_melt_freeze( cvm = ( c_air + qv * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) t = t - melt * lhi / cvm - elif frez > 0.0 and t <= driver_constants.TICE and ql > driver_constants.QCMIN: + elif frez > 0.0 and t <= constants.TICE and ql > constants.QCMIN: # ----------------------------------------------------------------------- # pihom: homogeneous freezing of cloud water into cloud ice # this is the 1st occurance of liquid water freezing in the split mp process @@ -127,7 +127,7 @@ def icloud_melt_freeze( 1.0, qa * max(qi + ql - frez + tmp, 0.0) - / max(qi + ql, driver_constants.QCMIN), + / max(qi + ql, constants.QCMIN), ), ) @@ -139,8 +139,8 @@ def icloud_melt_freeze( cvm = ( c_air + qv * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) t = t + frez * lhi / cvm @@ -302,23 +302,23 @@ def snow_graupel_coldrain( pgacr = 0.0 pgacw = 0.0 - tc = t2 - driver_constants.TICE + tc = t2 - constants.TICE if tc >= 0.0: # ----------------------------------------------------------------------- # melting of snow # ----------------------------------------------------------------------- - dqs0 = driver_constants.CES0 / p_dry - qv2 + dqs0 = constants.CES0 / p_dry - qv2 - if qs2 > driver_constants.QPMIN: + if qs2 > constants.QPMIN: # ----------------------------------------------------------------------- # psacw: accretion of cloud water by snow # only rate is used (for snow melt) since tc > 0. # ----------------------------------------------------------------------- - if ql2 > driver_constants.QCMIN: + if ql2 > constants.QCMIN: factor = denfac * csacw * exp(0.8125 * log(qs2 * den1)) psacw = factor / (1.0 + dts * factor) * ql2 # rate else: @@ -329,17 +329,17 @@ def snow_graupel_coldrain( # pracs: accretion of snow by rain # ----------------------------------------------------------------------- - if qr2 > driver_constants.QPMIN: + if qr2 > constants.QPMIN: psacr = min( acr3d( vts, vtr, qr2, qs2, - driver_constants.CSARC, - driver_constants.ACCO_01, - driver_constants.ACCO_11, - driver_constants.ACCO_21, + constants.CSARC, + constants.ACCO_01, + constants.ACCO_11, + constants.ACCO_21, den1, ), qr2 * rdts, @@ -349,10 +349,10 @@ def snow_graupel_coldrain( vts, qs2, qr2, - driver_constants.CRACS, - driver_constants.ACCO_00, - driver_constants.ACCO_10, - driver_constants.ACCO_20, + constants.CRACS, + constants.ACCO_00, + constants.ACCO_10, + constants.ACCO_20, den1, ) else: @@ -397,7 +397,7 @@ def snow_graupel_coldrain( 1.0, qa1 * max(qi2 + ql2 + tmp, 0.0) - / max(qi2 + ql2, driver_constants.QCMIN), + / max(qi2 + ql2, constants.QCMIN), ), ) @@ -409,40 +409,40 @@ def snow_graupel_coldrain( cvm = ( c_air + qv2 * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) t2 = t2 - sink * lhi / cvm - tc = t2 - driver_constants.TICE + tc = t2 - constants.TICE # ----------------------------------------------------------------------- # update capacity heat and latent heat coefficient # ----------------------------------------------------------------------- - lhi = driver_constants.LI00 + driver_constants.DC_ICE * t2 + lhi = constants.LI00 + constants.DC_ICE * t2 icpk = lhi / cvm # ----------------------------------------------------------------------- # melting of graupel # ----------------------------------------------------------------------- - if qg2 > driver_constants.QPMIN and tc > 0.0: + if qg2 > constants.QPMIN and tc > 0.0: # ----------------------------------------------------------------------- # pgacr: accretion of rain by graupel # ----------------------------------------------------------------------- - if qr2 > driver_constants.QPMIN: + if qr2 > constants.QPMIN: pgacr = min( acr3d( vtg, vtr, qr2, qg2, - driver_constants.CGARC, - driver_constants.ACCO_02, - driver_constants.ACCO_12, - driver_constants.ACCO_22, + constants.CGARC, + constants.ACCO_02, + constants.ACCO_12, + constants.ACCO_22, den1, ), rdts * qr2, @@ -453,7 +453,7 @@ def snow_graupel_coldrain( # ----------------------------------------------------------------------- qden = qg2 * den1 - if ql2 > driver_constants.QCMIN: + if ql2 > constants.QCMIN: factor = cgacw * qden / sqrt(den1 * sqrt(sqrt(qden))) pgacw = factor / (1.0 + dts * factor) * ql2 # rate @@ -482,8 +482,8 @@ def snow_graupel_coldrain( cvm = ( c_air + qv2 * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) t2 = t2 - pgmlt * lhi / cvm @@ -498,7 +498,7 @@ def snow_graupel_coldrain( # ----------------------------------------------------------------------- if qi2 > 3.0e-7: # cloud ice sink terms - if qs2 > driver_constants.QPMIN: + if qs2 > constants.QPMIN: # ----------------------------------------------------------------------- # sjl added (following lin eq. 23) the temperature dependency # to reduce accretion, use esi = exp (0.05 * tc) as in hong et al 2004 @@ -531,9 +531,9 @@ def snow_graupel_coldrain( else: tmp = fac_i2s * exp(0.025 * tc) - di = max(di, driver_constants.QCMIN) + di = max(di, constants.QCMIN) q_plus = qi2 + di - if q_plus > (qim + driver_constants.QCMIN): + if q_plus > (qim + constants.QCMIN): if qim > (qi2 - di): dq = (0.25 * (q_plus - qim) ** 2) / di else: @@ -550,7 +550,7 @@ def snow_graupel_coldrain( 1.0, qa1 * max(qi2 + ql2 - sink + tmp, 0.0) - / max(qi2 + ql2, driver_constants.QCMIN), + / max(qi2 + ql2, constants.QCMIN), ), ) @@ -561,7 +561,7 @@ def snow_graupel_coldrain( # pgaci: accretion of cloud ice by graupel # ----------------------------------------------------------------------- - if qg2 > driver_constants.QPMIN: + if qg2 > constants.QPMIN: # ----------------------------------------------------------------------- # factor = dts * cgaci / sqrt (den (k)) * # exp (0.05 * tc + 0.875 * log (qg * den (k))) @@ -580,9 +580,9 @@ def snow_graupel_coldrain( # rain to ice, snow, graupel processes: # ----------------------------------------------------------------------- - tc = t2 - driver_constants.TICE + tc = t2 - constants.TICE - if qr2 > driver_constants.QPMIN and tc < 0.0: + if qr2 > constants.QPMIN and tc < 0.0: # ----------------------------------------------------------------------- # * sink * terms to qr: psacr + pgfr @@ -594,16 +594,16 @@ def snow_graupel_coldrain( # psacr accretion of rain by snow # ----------------------------------------------------------------------- - if qs2 > driver_constants.QPMIN: # if snow exists + if qs2 > constants.QPMIN: # if snow exists psacr = dts * acr3d( vts, vtr, qr2, qs2, - driver_constants.CSARC, - driver_constants.ACCO_01, - driver_constants.ACCO_11, - driver_constants.ACCO_21, + constants.CSARC, + constants.ACCO_01, + constants.ACCO_11, + constants.ACCO_21, den1, ) else: @@ -626,7 +626,7 @@ def snow_graupel_coldrain( # ----------------------------------------------------------------------- sink = psacr + pgfr - factor = min(sink, min(qr2, -tc / icpk)) / max(sink, driver_constants.QPMIN) + factor = min(sink, min(qr2, -tc / icpk)) / max(sink, constants.QPMIN) psacr = factor * psacr pgfr = factor * pgfr @@ -640,8 +640,8 @@ def snow_graupel_coldrain( cvm = ( c_air + qv2 * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) t2 = t2 + sink * lhi / cvm @@ -649,29 +649,29 @@ def snow_graupel_coldrain( # # update capacity heat and latent heat coefficient # # ----------------------------------------------------------------------- - lhi = driver_constants.LI00 + driver_constants.DC_ICE * t2 + lhi = constants.LI00 + constants.DC_ICE * t2 icpk = lhi / cvm # # ----------------------------------------------------------------------- # # graupel production terms: # # ----------------------------------------------------------------------- - if qs2 > driver_constants.QPMIN: + if qs2 > constants.QPMIN: # ----------------------------------------------------------------------- # accretion: snow -- > graupel # ----------------------------------------------------------------------- - if qg2 > driver_constants.QPMIN: + if qg2 > constants.QPMIN: sink = dts * acr3d( vtg, vts, qs2, qs2, cgacs, - driver_constants.ACCO_03, - driver_constants.ACCO_13, - driver_constants.ACCO_23, + constants.ACCO_03, + constants.ACCO_13, + constants.ACCO_23, den1, ) else: @@ -683,7 +683,7 @@ def snow_graupel_coldrain( qsm = qs0_crt / den1 if qs2 > qsm: - factor = dts * 1.0e-3 * exp(0.09 * (t2 - driver_constants.TICE)) + factor = dts * 1.0e-3 * exp(0.09 * (t2 - constants.TICE)) sink = sink + factor / (1.0 + factor) * (qs2 - qsm) sink = min(qs2, sink) @@ -691,13 +691,13 @@ def snow_graupel_coldrain( qs2 = qs2 - sink qg2 = qg2 + sink - if qg2 > driver_constants.QPMIN and t2 < (driver_constants.TICE - 0.01): + if qg2 > constants.QPMIN and t2 < (constants.TICE - 0.01): # ----------------------------------------------------------------------- # pgacw: accretion of cloud water by graupel # ----------------------------------------------------------------------- - if ql2 > driver_constants.QCMIN: + if ql2 > constants.QCMIN: qden = qg2 * den1 factor = dts * cgacw * qden / sqrt(den1 * sqrt(sqrt(qden))) pgacw = factor / (1.0 + factor) * ql2 @@ -708,7 +708,7 @@ def snow_graupel_coldrain( # pgacr: accretion of rain by graupel # ----------------------------------------------------------------------- - if qr2 > driver_constants.QPMIN: + if qr2 > constants.QPMIN: pgacr = min( dts * acr3d( @@ -716,10 +716,10 @@ def snow_graupel_coldrain( vtr, qr2, qg2, - driver_constants.CGARC, - driver_constants.ACCO_02, - driver_constants.ACCO_12, - driver_constants.ACCO_22, + constants.CGARC, + constants.ACCO_02, + constants.ACCO_12, + constants.ACCO_22, den1, ), qr2, @@ -728,11 +728,11 @@ def snow_graupel_coldrain( pgacr = 0.0 sink = pgacr + pgacw - if driver_constants.TICE - t2 > 0: - ans = driver_constants.TICE - t2 + if constants.TICE - t2 > 0: + ans = constants.TICE - t2 else: ans = 0 - factor = min(sink, ans / icpk) / max(sink, driver_constants.QPMIN) + factor = min(sink, ans / icpk) / max(sink, constants.QPMIN) pgacr = factor * pgacr pgacw = factor * pgacw @@ -745,8 +745,8 @@ def snow_graupel_coldrain( cvm = ( c_air + qv2 * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) t2 = t2 + sink * lhi / cvm @@ -797,7 +797,7 @@ def iqs1( reference Fortran: gfdl_cloud_microphys.F90: function iqs1 """ - tmin = driver_constants.TABLE_ICE - 160.0 + tmin = constants.TABLE_ICE - 160.0 if ta - tmin > 0: ans = ta - tmin else: @@ -806,7 +806,7 @@ def iqs1( ap1 = min(2621.0, ap1) it = i32(trunc(ap1)) es = table3.A[it - 1] + (ap1 - it) * des3.A[it - 1] - iqs1 = es / (driver_constants.RVGAS * ta * den) + iqs1 = es / (constants.RVGAS * ta * den) return iqs1 @@ -827,7 +827,7 @@ def iqs2( reference Fortran: gfdl_cloud_microphys.F90: function iqs2 """ - tmin = driver_constants.TABLE_ICE - 160.0 + tmin = constants.TABLE_ICE - 160.0 if ta - tmin > 0: ans = ta - tmin else: @@ -836,12 +836,12 @@ def iqs2( ap1 = min(2621.0, ap1) it = i32(trunc(ap1)) es = table3.A[it - 1] + (ap1 - it) * des3.A[it - 1] - iqs2 = es / (driver_constants.RVGAS * ta * den) + iqs2 = es / (constants.RVGAS * ta * den) it = i32(trunc(ap1 - 0.5)) dqdt = ( 10.0 * (des3.A[it - 1] + (ap1 - it) * (des3.A[it] - des3.A[it - 1])) - / (driver_constants.RVGAS * ta * den) + / (constants.RVGAS * ta * den) ) return iqs2, dqdt @@ -862,7 +862,7 @@ def wqs1( reference Fortran: gfdl_cloud_microphys.F90: function wqs1 """ - tmin = driver_constants.TABLE_ICE - 160.0 + tmin = constants.TABLE_ICE - 160.0 if ta - tmin > 0: ans = ta - tmin else: @@ -871,7 +871,7 @@ def wqs1( ap1 = min(2621.0, ap1) it = i32(trunc(ap1)) es = table2.A[it - 1] + (ap1 - it) * des2.A[it - 1] - wqs1 = es / (driver_constants.RVGAS * ta * den) + wqs1 = es / (constants.RVGAS * ta * den) return wqs1 @@ -892,7 +892,7 @@ def wqs2( reference Fortran: gfdl_cloud_microphys.F90: function wqs2 """ - tmin = driver_constants.TABLE_ICE - 160.0 + tmin = constants.TABLE_ICE - 160.0 if ta - tmin > 0: ans = ta - tmin @@ -902,13 +902,13 @@ def wqs2( ap1 = min(2621.0, ap1) it = i32(trunc(ap1)) es = table2.A[it - 1] + (ap1 - it) * des2.A[it - 1] - qsat = es / (driver_constants.RVGAS * ta * den) + qsat = es / (constants.RVGAS * ta * den) it = i32(trunc(ap1 - 0.5)) # finite diff, del_t = 0.1: dqdt = ( 10.0 * (des2.A[it - 1] + (ap1 - it) * (des2.A[it] - des2.A[it - 1])) - / (driver_constants.RVGAS * ta * den) + / (constants.RVGAS * ta * den) ) return qsat, dqdt @@ -982,24 +982,24 @@ def subgrid_z_proc( # ----------------------------------------------------------------------- lhl = lv00 + d0_vap * t1 - lhi = driver_constants.LI00 + driver_constants.DC_ICE * t1 + lhi = constants.LI00 + constants.DC_ICE * t1 q_liq = ql1 + qr1 q_sol = qi1 + qs1 + qg1 cvm = ( c_air + qv1 * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) lcpk = lhl / cvm icpk = lhi / cvm tcpk = lcpk + icpk - if driver_constants.TICE - t1 > 0: - ans = driver_constants.TICE - t1 + if constants.TICE - t1 > 0: + ans = constants.TICE - t1 else: ans = 0 tcp3 = lcpk + icpk * min( - 1.0, ans / (driver_constants.TICE - driver_constants.T_WFR) + 1.0, ans / (constants.TICE - constants.T_WFR) ) rh_adj = 1.0 - rh_limited - rh_inc @@ -1008,7 +1008,7 @@ def subgrid_z_proc( subl1 = 0.0 cycle = False - if p_dry < driver_constants.P_MIN: + if p_dry < constants.P_MIN: cycle = True # ----------------------------------------------------------------------- @@ -1016,8 +1016,8 @@ def subgrid_z_proc( # ----------------------------------------------------------------------- if t1 < t_min and cycle == False: # noqa - if qv1 - driver_constants.QVMIN > 0: - sink = qv1 - driver_constants.QVMIN + if qv1 - constants.QVMIN > 0: + sink = qv1 - constants.QVMIN else: sink = 0 qv1 = qv1 - sink @@ -1026,8 +1026,8 @@ def subgrid_z_proc( cvm = ( c_air + qv1 * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) t1 = t1 + sink * (lhl + lhi) / cvm if do_qa == True: # noqa @@ -1039,16 +1039,16 @@ def subgrid_z_proc( # update heat capacity and latent heat coefficient # ----------------------------------------------------------------------- lhl = lv00 + d0_vap * t1 - lhi = driver_constants.LI00 + driver_constants.DC_ICE * t1 + lhi = constants.LI00 + constants.DC_ICE * t1 lcpk = lhl / cvm icpk = lhi / cvm tcpk = lcpk + icpk - if driver_constants.TICE - t1 > 0: - ans = driver_constants.TICE - t1 + if constants.TICE - t1 > 0: + ans = constants.TICE - t1 else: ans = 0 tcp3 = lcpk + icpk * min( - 1.0, ans / (driver_constants.TICE - driver_constants.T_WFR) + 1.0, ans / (constants.TICE - constants.T_WFR) ) # ----------------------------------------------------------------------- @@ -1058,8 +1058,8 @@ def subgrid_z_proc( tin = t1 - (lhl * (ql1 + qi1) + lhi * qi1) / ( c_air + qpz * c_vap - + qr1 * driver_constants.C_LIQ - + (qs1 + qg1) * driver_constants.C_ICE + + qr1 * constants.C_LIQ + + (qs1 + qg1) * constants.C_ICE ) if tin > t_sub + 6.0: rh = qpz / iqs1(tin, den1, table3, des3) @@ -1079,7 +1079,7 @@ def subgrid_z_proc( if do_evap == True: # noqa qsw, dwsdt = wqs2(t1, den1, table2, des2) dq0 = qsw - qv1 - if dq0 > driver_constants.QVMIN: + if dq0 > constants.QVMIN: factor = min(1.0, fac_l2v * (10.0 * dq0 / qsw)) evap = min(ql1, factor * ql1 / (1.0 + tcp3 * dwsdt)) else: @@ -1090,8 +1090,8 @@ def subgrid_z_proc( cvm = ( c_air + qv1 * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) t1 = t1 - evap * lhl / cvm @@ -1099,7 +1099,7 @@ def subgrid_z_proc( # update heat capacity and latent heat coefficient # ----------------------------------------------------------------------- - lhi = driver_constants.LI00 + driver_constants.DC_ICE * t1 + lhi = constants.LI00 + constants.DC_ICE * t1 icpk = lhi / cvm # ----------------------------------------------------------------------- @@ -1107,7 +1107,7 @@ def subgrid_z_proc( # ----------------------------------------------------------------------- ifrac = ice_fraction(t1, cnv_frc, srf_type) - if ifrac == 1.0 and ql1 > driver_constants.QCMIN: + if ifrac == 1.0 and ql1 > constants.QCMIN: sink = ql1 ql1 = ql1 - sink qi1 = qi1 + sink @@ -1116,8 +1116,8 @@ def subgrid_z_proc( cvm = ( c_air + qv1 * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) t1 = t1 + sink * lhi / cvm @@ -1125,18 +1125,18 @@ def subgrid_z_proc( # update heat capacity and latent heat coefficient # ----------------------------------------------------------------------- - lhi = driver_constants.LI00 + driver_constants.DC_ICE * t1 + lhi = constants.LI00 + constants.DC_ICE * t1 icpk = lhi / cvm # ----------------------------------------------------------------------- # bigg mechanism heterogeneous freezing on existing cloud nuclei # ----------------------------------------------------------------------- - tc = driver_constants.TICE - t1 - if do_bigg == True and ql1 > driver_constants.QCMIN and tc > 0.0: # noqa + tc = constants.TICE - t1 + if do_bigg == True and ql1 > constants.QCMIN and tc > 0.0: # noqa sink = ( fac_frz - * (100.0 / driver_constants.RHOR / ccn) + * (100.0 / constants.RHOR / ccn) * dts * (exp(0.66 * tc) - 1.0) * den1 @@ -1151,8 +1151,8 @@ def subgrid_z_proc( cvm = ( c_air + qv1 * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) t1 = t1 + sink * lhi / cvm # significant ql existed @@ -1162,7 +1162,7 @@ def subgrid_z_proc( # ----------------------------------------------------------------------- lhl = lv00 + d0_vap * t1 - lhi = driver_constants.LI00 + driver_constants.DC_ICE * t1 + lhi = constants.LI00 + constants.DC_ICE * t1 lcpk = lhl / cvm icpk = lhi / cvm tcpk = lcpk + icpk @@ -1171,11 +1171,11 @@ def subgrid_z_proc( # sublimation / deposition of LS ice # ----------------------------------------------------------------------- - if t1 < driver_constants.TICE: + if t1 < constants.TICE: qsi, dqsdt = iqs2(t1, den1, table3, des3) dq = qv1 - qsi sink = min(qi1, dq / (1.0 + tcpk * dqsdt)) - if qi1 > driver_constants.QCMIN: + if qi1 > constants.QCMIN: # eq 9, hong et al. 2004, mwr # for a and b, see dudhia 1989: page 3103 eq (b7) and (b8) pidep = ( @@ -1184,7 +1184,7 @@ def subgrid_z_proc( * 349138.78 * exp(0.875 * log(qi1 * den1)) / ( - qsi * den1 * lat2 / (0.0243 * driver_constants.RVGAS * t1**2) + qsi * den1 * lat2 / (0.0243 * constants.RVGAS * t1**2) + 4.42478e4 ) ) @@ -1193,7 +1193,7 @@ def subgrid_z_proc( if dq > 0.0: # vapor - > ice # deposition ifrac = ice_fraction(t1, cnv_frc, srf_type) - tmp = driver_constants.TICE - t1 + tmp = constants.TICE - t1 qi_crt = 4.92e-11 * exp(1.33 * log(1.0e3 * exp(0.1 * tmp))) qi_crt = max(qi_crt, 1.82e-6) * qi_lim * ifrac / den1 sink = min(sink, min(max(qi_crt - qi1, pidep), tmp / tcpk)) @@ -1220,8 +1220,8 @@ def subgrid_z_proc( cvm = ( c_air + qv1 * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) t1 = t1 + sink * (lhl + lhi) / cvm @@ -1230,7 +1230,7 @@ def subgrid_z_proc( # ----------------------------------------------------------------------- lhl = lv00 + d0_vap * t1 - lhi = driver_constants.LI00 + driver_constants.DC_ICE * t1 + lhi = constants.LI00 + constants.DC_ICE * t1 lcpk = lhl / cvm icpk = lhi / cvm tcpk = lcpk + icpk @@ -1240,7 +1240,7 @@ def subgrid_z_proc( # this process happens for all temp rage # ----------------------------------------------------------------------- - if qs1 > driver_constants.QPMIN: + if qs1 > constants.QPMIN: qsi, dqsdt = iqs2(t1, den1, table3, des3) qden = qs1 * den1 tmp = exp(0.65625 * log(qden)) @@ -1261,11 +1261,11 @@ def subgrid_z_proc( pssub = min(fac_s2v * pssub * min(1.0, ans * 0.2), qs1) subl1 = subl1 + pssub / dts else: - if t1 > driver_constants.TICE: + if t1 > constants.TICE: pssub = 0.0 # no deposition else: pssub = max( - fac_v2s * pssub, max(dq, (t1 - driver_constants.TICE) / tcpk) + fac_v2s * pssub, max(dq, (t1 - constants.TICE) / tcpk) ) qs1 = qs1 - pssub qv1 = qv1 + pssub @@ -1273,8 +1273,8 @@ def subgrid_z_proc( cvm = ( c_air + qv1 * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) t1 = t1 - pssub * (lhl + lhi) / cvm @@ -1283,7 +1283,7 @@ def subgrid_z_proc( # ----------------------------------------------------------------------- lhl = lv00 + d0_vap * t1 - lhi = driver_constants.LI00 + driver_constants.DC_ICE * t1 + lhi = constants.LI00 + constants.DC_ICE * t1 lcpk = lhl / cvm icpk = lhi / cvm tcpk = lcpk + icpk @@ -1292,19 +1292,19 @@ def subgrid_z_proc( # simplified 2 - way grapuel sublimation - deposition mechanism # ----------------------------------------------------------------------- - if qg1 > driver_constants.QPMIN: + if qg1 > constants.QPMIN: qsi, dqsdt = iqs2(t1, den1, table3, des3) dq = (qv1 - qsi) / (1.0 + tcpk * dqsdt) pgsub = (qv1 / qsi - 1.0) * qg1 if pgsub > 0.0: # deposition - if t1 > driver_constants.TICE: + if t1 > constants.TICE: pgsub = 0.0 # no deposition else: pgsub = min( fac_v2g * pgsub, min( 0.2 * dq, - min(ql1 + qr1, (driver_constants.TICE - t1) / tcpk), + min(ql1 + qr1, (constants.TICE - t1) / tcpk), ), ) else: # submilation @@ -1320,8 +1320,8 @@ def subgrid_z_proc( cvm = ( c_air + qv1 * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) t1 = t1 + pgsub * (lhl + lhi) / cvm @@ -1367,10 +1367,10 @@ def subgrid_z_proc( # determine saturated specific humidity # ----------------------------------------------------------------------- - if tin <= driver_constants.T_WFR: + if tin <= constants.T_WFR: # ice phase: qstar = iqs1(tin, den1, table3, des3) - elif tin >= driver_constants.TICE: + elif tin >= constants.TICE: # liquid phase: qstar = wqs1(tin, den1, table2, des2) else: @@ -1388,9 +1388,9 @@ def subgrid_z_proc( # assuming subgrid linear distribution in horizontal; # this is effectively a smoother for the binary cloud scheme # ----------------------------------------------------------------------- - if qpz > driver_constants.QCMIN: + if qpz > constants.QCMIN: # partial cloudiness by pdf: - dq = max(driver_constants.QCMIN, rh_limited * qpz) + dq = max(constants.QCMIN, rh_limited * qpz) q_plus = qpz + dq # cloud free if qstar > q_plus q_minus = qpz - dq if icloud_f == 3: @@ -1400,7 +1400,7 @@ def subgrid_z_proc( do_nothing = True elif qpz <= qstar and qstar < q_plus: # partial cloud cover qa1 = max( - driver_constants.QCMIN, + constants.QCMIN, min( 1.0, qa1 @@ -1411,7 +1411,7 @@ def subgrid_z_proc( ) elif q_minus <= qstar and qstar < qpz: # partial cloud cover qa1 = max( - driver_constants.QCMIN, + constants.QCMIN, min( 1.0, qa1 @@ -1432,7 +1432,7 @@ def subgrid_z_proc( do_nothing = True elif qstar < q_plus and q_cond > qc_crt: qa1 = max( - driver_constants.QCMIN, + constants.QCMIN, min(1.0, qa1 + (q_plus - qstar) / (dq + dq)), ) # partial cloud cover elif qstar <= q_minus: @@ -1558,11 +1558,11 @@ def icloud( if z_slope_ice == True: # noqa dm_linear_prof = max( dm_linear_prof, - max(driver_constants.QVMIN, h_var_linear_prof * q_linear_prof), + max(constants.QVMIN, h_var_linear_prof * q_linear_prof), ) if z_slope_ice == False: # noqa dm_linear_prof = max( - driver_constants.QVMIN, h_var_linear_prof * q_linear_prof + constants.QVMIN, h_var_linear_prof * q_linear_prof ) # handle outputs of "function" @@ -1576,7 +1576,7 @@ def icloud( # update capacity heat and latent heat coefficient # ----------------------------------------------------------------------- lhl = lv00 + d0_vap * t1 - lhi = driver_constants.LI00 + driver_constants.DC_ICE * t1 + lhi = constants.LI00 + constants.DC_ICE * t1 lcpk = lhl / cvm icpk = lhi / cvm tcpk = lcpk + icpk @@ -1584,7 +1584,7 @@ def icloud( # ----------------------------------------------------------------------- # do nothing above p_min # ----------------------------------------------------------------------- - if p_dry >= driver_constants.P_MIN: + if p_dry >= constants.P_MIN: ( t1, qv1, diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/GFDL_1M_driver_tables.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/sat_tables.py similarity index 81% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/GFDL_1M_driver_tables.py rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/sat_tables.py index 80ebd68b8..5104bf36a 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/GFDL_1M_driver_tables.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/sat_tables.py @@ -9,7 +9,7 @@ log10, ) -import pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver_constants as driver_constants +import pyMoist.GFDL_1M.GFDL_1M_driver.constants as constants from ndsl.boilerplate import get_factories_single_tile from ndsl.constants import X_DIM, Y_DIM, Z_DIM from ndsl.dsl.typing import Float @@ -17,9 +17,9 @@ _FloatField_data_dim = gtscript.Field[ - gtscript.IJK, (Float, (int(driver_constants.LENGTH))) + gtscript.IJK, (Float, (int(constants.LENGTH))) ] -GlobalTable_driver_qsat = gtscript.GlobalTable[(Float, (driver_constants.LENGTH))] +GlobalTable_driver_qsat = gtscript.GlobalTable[(Float, (constants.LENGTH))] def qs_table_1(length: i32, table1: _FloatField_data_dim, esupc: _FloatField_data_dim): @@ -31,7 +31,7 @@ def qs_table_1(length: i32, table1: _FloatField_data_dim, esupc: _FloatField_dat """ with computation(PARALLEL), interval(...): delt = 0.1 - tmin = driver_constants.TABLE_ICE - 160.0 + tmin = constants.TABLE_ICE - 160.0 # ----------------------------------------------------------------------- # compute es over ice between - 160 deg c and 0 deg c. @@ -40,12 +40,12 @@ def qs_table_1(length: i32, table1: _FloatField_data_dim, esupc: _FloatField_dat i = 0 while i < 1600: tem = tmin + delt * i - fac0 = (tem - driver_constants.T_ICE) / (tem * driver_constants.T_ICE) - fac1 = fac0 * driver_constants.LI2 + fac0 = (tem - constants.T_ICE) / (tem * constants.T_ICE) + fac1 = fac0 * constants.LI2 fac2 = ( - driver_constants.D2ICE * log(tem / driver_constants.T_ICE) + fac1 - ) / driver_constants.RVGAS - table1[0, 0, 0][i] = driver_constants.E_00 * exp(fac2) + constants.D2ICE * log(tem / constants.T_ICE) + fac1 + ) / constants.RVGAS + table1[0, 0, 0][i] = constants.E_00 * exp(fac2) i = i + 1 # ----------------------------------------------------------------------- @@ -55,12 +55,12 @@ def qs_table_1(length: i32, table1: _FloatField_data_dim, esupc: _FloatField_dat i = 0 while i < 1421: tem = 233.16 + delt * i - fac0 = (tem - driver_constants.T_ICE) / (tem * driver_constants.T_ICE) - fac1 = fac0 * driver_constants.LV0 + fac0 = (tem - constants.T_ICE) / (tem * constants.T_ICE) + fac1 = fac0 * constants.LV0 fac2 = ( - driver_constants.DC_VAP * log(tem / driver_constants.T_ICE) + fac1 - ) / driver_constants.RVGAS - esh40 = driver_constants.E_00 * exp(fac2) + constants.DC_VAP * log(tem / constants.T_ICE) + fac1 + ) / constants.RVGAS + esh40 = constants.E_00 * exp(fac2) if i < 400: esupc[0, 0, 0][i] = esh40 else: @@ -93,17 +93,17 @@ def qs_table_2(length: i32, table2: _FloatField_data_dim): """ with computation(PARALLEL), interval(...): delt = 0.1 - tmin = driver_constants.TABLE_ICE - 160.0 + tmin = constants.TABLE_ICE - 160.0 i = 0 while i < length: tem = tmin + delt * i - fac0 = (tem - driver_constants.T_ICE) / (tem * driver_constants.T_ICE) - fac1 = fac0 * driver_constants.LV0 + fac0 = (tem - constants.T_ICE) / (tem * constants.T_ICE) + fac1 = fac0 * constants.LV0 fac2 = ( - driver_constants.DC_VAP * log(tem / driver_constants.T_ICE) + fac1 - ) / driver_constants.RVGAS - table2[0, 0, 0][i] = driver_constants.E_00 * exp(fac2) + constants.DC_VAP * log(tem / constants.T_ICE) + fac1 + ) / constants.RVGAS + table2[0, 0, 0][i] = constants.E_00 * exp(fac2) i = i + 1 @@ -116,29 +116,29 @@ def qs_table_3(length: i32, table3: _FloatField_data_dim, table1: _FloatField_da """ with computation(PARALLEL), interval(...): delt = 0.1 - tmin = driver_constants.TABLE_ICE - 160.0 + tmin = constants.TABLE_ICE - 160.0 i = 0 while i < length: tem0 = tmin + delt * i - fac0 = (tem0 - driver_constants.T_ICE) / (tem0 * driver_constants.T_ICE) + fac0 = (tem0 - constants.T_ICE) / (tem0 * constants.T_ICE) if i < 1600: # ----------------------------------------------------------------------- # compute es over ice between - 160 deg c and 0 deg c. # ----------------------------------------------------------------------- - fac1 = fac0 * driver_constants.LI2 + fac1 = fac0 * constants.LI2 fac2 = ( - driver_constants.D2ICE * log(tem0 / driver_constants.T_ICE) + fac1 - ) / driver_constants.RVGAS + constants.D2ICE * log(tem0 / constants.T_ICE) + fac1 + ) / constants.RVGAS else: # ----------------------------------------------------------------------- # compute es over water between 0 deg c and 102 deg c. # ----------------------------------------------------------------------- - fac1 = fac0 * driver_constants.LV0 + fac1 = fac0 * constants.LV0 fac2 = ( - driver_constants.DC_VAP * log(tem0 / driver_constants.T_ICE) + fac1 - ) / driver_constants.RVGAS - table3[0, 0, 0][i] = driver_constants.E_00 * exp(fac2) + constants.DC_VAP * log(tem0 / constants.T_ICE) + fac1 + ) / constants.RVGAS + table3[0, 0, 0][i] = constants.E_00 * exp(fac2) i = i + 1 # ----------------------------------------------------------------------- @@ -171,9 +171,9 @@ def qs_table_4(length: i32, table4: _FloatField_data_dim, table1: _FloatField_da with computation(PARALLEL), interval(...): delt = 0.1 esbasw = 1013246.0 - tbasw = driver_constants.TABLE_ICE + 100.0 + tbasw = constants.TABLE_ICE + 100.0 esbasi = 6107.1 - tmin = driver_constants.TABLE_ICE - 160.0 + tmin = constants.TABLE_ICE - 160.0 i = 0 while i < length: @@ -183,9 +183,9 @@ def qs_table_4(length: i32, table4: _FloatField_data_dim, table1: _FloatField_da # compute es over ice between - 160 deg c and 0 deg c. # see smithsonian meteorological tables page 350. # ----------------------------------------------------------------------- - aa = -9.09718 * (driver_constants.TABLE_ICE / tem - 1.0) - b = -3.56654 * log10(driver_constants.TABLE_ICE / tem) - c = 0.876793 * (1.0 - tem / driver_constants.TABLE_ICE) + aa = -9.09718 * (constants.TABLE_ICE / tem - 1.0) + b = -3.56654 * log10(constants.TABLE_ICE / tem) + c = 0.876793 * (1.0 - tem / constants.TABLE_ICE) e = log10(esbasi) table4[0, 0, 0][i] = 0.1 * 10 ** (aa + b + c + e) else: @@ -264,7 +264,7 @@ def __init__(self, backend): ) quantity_factory_data_dim.set_extra_dim_lengths( **{ - "table_axis": driver_constants.LENGTH, + "table_axis": constants.LENGTH, } ) @@ -322,12 +322,12 @@ def __init__(self, backend): domain=qsat_domain, ) - compute_qs_table_1(driver_constants.LENGTH, self._table1, self._esupc) - compute_qs_table_2(driver_constants.LENGTH, self._table2) - compute_qs_table_3(driver_constants.LENGTH, self._table3, self._table1) - compute_qs_table_4(driver_constants.LENGTH, self._table4, self._table1) + compute_qs_table_1(constants.LENGTH, self._table1, self._esupc) + compute_qs_table_2(constants.LENGTH, self._table2) + compute_qs_table_3(constants.LENGTH, self._table3, self._table1) + compute_qs_table_4(constants.LENGTH, self._table4, self._table1) compute_des_tables( - driver_constants.LENGTH, + constants.LENGTH, self._des1, self._des2, self._des3, diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/support.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/support.py new file mode 100644 index 000000000..f7d7d7122 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/support.py @@ -0,0 +1,428 @@ +import numpy as np +from gt4py.cartesian.gtscript import i32 +from ndsl.dsl.typing import Float +from ndsl.constants import X_DIM, Y_DIM, Z_DIM, Z_INTERFACE_DIM +import pyMoist.GFDL_1M.GFDL_1M_driver.constants as constants +from pyMoist.GFDL_1M.GFDL_1M_driver.config import config +from ndsl.dsl.typing import Int + + +def check_flags( + GFDL_1M_config: config, + dts: Float, +): + failed_keywords = [] + if not GFDL_1M_config.phys_hydrostatic: + failed_keywords.append("phys_hydrostatic") + if GFDL_1M_config.hydrostatic: + failed_keywords.append("hydrostatic") + if GFDL_1M_config.const_vi: + failed_keywords.append("const_vi") + if GFDL_1M_config.const_vs: + failed_keywords.append("const_vs") + if GFDL_1M_config.const_vg: + failed_keywords.append("const_vg") + if GFDL_1M_config.const_vr: + failed_keywords.append("const_vr") + if GFDL_1M_config.use_ppm: + failed_keywords.append("use_ppm") + if not GFDL_1M_config.use_ccn: + failed_keywords.append("use_ccn") + if GFDL_1M_config.do_qa: + failed_keywords.append("do_qa") + if not GFDL_1M_config.fix_negative: + failed_keywords.append("fix_negative") + if GFDL_1M_config.fast_sat_adj: + failed_keywords.append("fast_sat_adj") + if GFDL_1M_config.do_bigg: + failed_keywords.append("do_bigg") + if GFDL_1M_config.do_evap: + failed_keywords.append("do_evap") + if GFDL_1M_config.do_subl: + failed_keywords.append("do_subl") + if not GFDL_1M_config.z_slope_liq: + failed_keywords.append("z_slope_liq") + if not GFDL_1M_config.z_slope_ice: + failed_keywords.append("z_slope_ice") + if not GFDL_1M_config.prog_ccn: + failed_keywords.append("prog_ccn") + if not GFDL_1M_config.preciprad: + failed_keywords.append("preciprad") + if not GFDL_1M_config.mono_prof: + failed_keywords.append("mono_prof") + if GFDL_1M_config.do_sedi_heat: + failed_keywords.append("do_sedi_heat") + if not GFDL_1M_config.sedi_transport: + failed_keywords.append("sedi_transport") + if GFDL_1M_config.do_sedi_w: + failed_keywords.append("do_sedi_w") + if GFDL_1M_config.de_ice: + failed_keywords.append("de_ice") + if GFDL_1M_config.mp_print: + failed_keywords.append("mp_print") + if dts >= 300: + failed_keywords.append("dts") + + if len(failed_keywords) > 0: + raise ValueError( + "One or more namelist parameters do not meet \ + expected values. Failing parameters: ", + failed_keywords, + ) + + +class config_constants: + def __init__(self, GFDL_1M_config: config): + # ----------------------------------------------------------------------- + # define heat capacity of dry air and water vap based on hydrostatical property + # ----------------------------------------------------------------------- + + if GFDL_1M_config.phys_hydrostatic or GFDL_1M_config.hydrostatic: + self.C_AIR = constants.CP_AIR + self.C_VAP = constants.CP_VAP + self.P_NONHYDRO = False + else: + self.C_AIR = constants.CV_AIR + self.C_VAP = constants.CV_VAP + self.P_NONHYDRO = True + self.D0_VAP = self.C_VAP - constants.C_LIQ + self.LV00 = constants.HLV0 - self.D0_VAP * constants.T_ICE + + if GFDL_1M_config.hydrostatic: + self.DO_SEDI_W = False + + # ----------------------------------------------------------------------- + # define latent heat coefficient used in wet bulb and bigg mechanism + # ----------------------------------------------------------------------- + + self.LATV = constants.HLV + self.LATI = constants.HLF + self.LATS = self.LATV + self.LATI + self.LAT2 = self.LATS * self.LATS + + self.LCP = self.LATV / constants.CP_AIR + self.ICP = self.LATI / constants.CP_AIR + self.TCP = (self.LATV + self.LATI) / constants.CP_AIR + + # ----------------------------------------------------------------------- + # define cloud microphysics sub time step + # ----------------------------------------------------------------------- + + self.MPDT = min(GFDL_1M_config.dt_moist, GFDL_1M_config.mp_time) + self.RDT = Float(1.0) / GFDL_1M_config.dt_moist + self.NTIMES = i32(np.round(GFDL_1M_config.dt_moist / self.MPDT)) + + # small time step: + self.DTS = GFDL_1M_config.dt_moist / Float(self.NTIMES) + + # ----------------------------------------------------------------------- + # calculate cloud condensation nuclei (ccn) + # the following is based on klein eq. 15 + # ----------------------------------------------------------------------- + + self.CPAUT = ( + GFDL_1M_config.c_paut * Float(0.104) * constants.GRAV / Float(1.717e-5) + ) + + # ----------------------------------------------------------------------- + # define conversion scalar / factor for icloud + # ----------------------------------------------------------------------- + self.RDTS = Float(1.0) / self.DTS + self.FAC_IMLT = Float(1.0) - np.exp( + -self.DTS / GFDL_1M_config.tau_imlt, dtype=Float + ) + self.FAC_I2S = Float(1.0) - np.exp( + -self.DTS / GFDL_1M_config.tau_i2s, dtype=Float + ) + self.FAC_V2L = Float(1.0) - np.exp( + -self.DTS / GFDL_1M_config.tau_v2l, dtype=Float + ) + self.FAC_L2V = Float(1.0) - np.exp( + -self.DTS / GFDL_1M_config.tau_l2v, dtype=Float + ) + self.FAC_I2V = Float(1.0) - np.exp( + -self.DTS / GFDL_1M_config.tau_i2v, dtype=Float + ) + self.FAC_S2V = Float(1.0) - np.exp( + -self.DTS / GFDL_1M_config.tau_s2v, dtype=Float + ) + self.FAC_V2S = Float(1.0) - np.exp( + -self.DTS / GFDL_1M_config.tau_v2s, dtype=Float + ) + self.FAC_G2V = Float(1.0) - np.exp( + -self.DTS / GFDL_1M_config.tau_g2v, dtype=Float + ) + self.FAC_V2G = Float(1.0) - np.exp( + -self.DTS / GFDL_1M_config.tau_v2g, dtype=Float + ) + self.FAC_FRZ = Float(1.0) - np.exp( + -self.DTS / GFDL_1M_config.tau_frz, dtype=Float + ) + + # ----------------------------------------------------------------------- + # constatns from setupm + # ----------------------------------------------------------------------- + + self.CGACS = constants.PISQ * constants.RNZG * constants.RNZS * constants.RHOS + self.CGACS = self.CGACS * GFDL_1M_config.c_pgacs + + self.CSACW = ( + constants.PIE + * constants.RNZS + * GFDL_1M_config.clin + * constants.GAM325 + / (Float(4.0) * np.power(constants.ACT[0], 0.8125, dtype=Float)) + ) + # decreasing csacw to reduce cloud water --- > snow + + self.CRACI = ( + constants.PIE + * constants.RNZR + * GFDL_1M_config.alin + * constants.GAM380 + / (Float(4.0) * np.power(constants.ACT[1], 0.95, dtype=Float)) + ) + self.CSACI = self.CSACW * GFDL_1M_config.c_psaci + + self.CGACW = ( + constants.PIE + * constants.RNZG + * constants.GAM350 + * constants.GCON + / (Float(4.0) * np.power(constants.ACT[5], 0.875, dtype=Float)) + ) + + self.CGACI = self.CGACW * GFDL_1M_config.c_pgaci + + self.CRACW = self.CRACI # cracw = 3.27206196043822 + self.CRACW = GFDL_1M_config.c_cracw * self.CRACW + + self.CSSUB = np.zeros(5) + self.CGSUB = np.zeros(5) + self.CREVP = np.zeros(5) + + self.CSSUB[0] = ( + Float(2.0) + * constants.PIE + * constants.VDIFU + * constants.TCOND + * constants.RVGAS + * constants.RNZS + ) + self.CGSUB[0] = ( + Float(2.0) + * constants.PIE + * constants.VDIFU + * constants.TCOND + * constants.RVGAS + * constants.RNZG + ) + self.CREVP[0] = ( + Float(2.0) + * constants.PIE + * constants.VDIFU + * constants.TCOND + * constants.RVGAS + * constants.RNZR + ) + self.CSSUB[1] = Float(0.78) / np.sqrt(constants.ACT[0], dtype=Float) + self.CGSUB[1] = Float(0.78) / np.sqrt(constants.ACT[5], dtype=Float) + self.CREVP[1] = Float(0.78) / np.sqrt(constants.ACT[1], dtype=Float) + self.CSSUB[2] = ( + Float(0.31) + * constants.SCM3 + * constants.GAM263 + * np.sqrt(GFDL_1M_config.clin / constants.VISK, dtype=Float) + / np.power(constants.ACT[0], Float(0.65625), dtype=Float) + ) + self.CGSUB[2] = ( + Float(0.31) + * constants.SCM3 + * constants.GAM275 + * np.sqrt(constants.GCON / constants.VISK, dtype=Float) + / np.power(constants.ACT[5], Float(0.6875), dtype=Float) + ) + self.CREVP[2] = ( + Float(0.31) + * constants.SCM3 + * constants.GAM209 + * np.sqrt(GFDL_1M_config.alin / constants.VISK, dtype=Float) + / np.power(constants.ACT[1], Float(0.725), dtype=Float) + ) + self.CSSUB[3] = constants.TCOND * constants.RVGAS + self.CSSUB[4] = np.power(constants.HLTS, 2, dtype=Float) * constants.VDIFU + self.CGSUB[3] = self.CSSUB[3] + self.CREVP[3] = self.CSSUB[3] + self.CGSUB[4] = self.CSSUB[4] + self.CREVP[4] = np.power(constants.HLTC, 2, dtype=Float) * constants.VDIFU + + self.CGFR_0 = ( + Float(20.0e2) + * constants.PISQ + * constants.RNZR + * constants.RHOR + / np.power(constants.ACT[1], Float(1.75), dtype=Float) + ) + self.CGFR_1 = Float(0.66) + + self.CSSUB_0 = self.CSSUB[0] + self.CSSUB_1 = self.CSSUB[1] + self.CSSUB_2 = self.CSSUB[2] + self.CSSUB_3 = self.CSSUB[3] + self.CSSUB_4 = self.CSSUB[4] + + self.CGSUB_0 = self.CGSUB[0] + self.CGSUB_1 = self.CGSUB[1] + self.CGSUB_2 = self.CGSUB[2] + self.CGSUB_3 = self.CGSUB[3] + self.CGSUB_4 = self.CGSUB[4] + + self.CREVP_0 = self.CREVP[0] + self.CREVP_1 = self.CREVP[1] + self.CREVP_2 = self.CREVP[2] + self.CREVP_3 = self.CREVP[3] + self.CREVP_4 = self.CREVP[4] + + # smlt: five constants (lin et al. 1983) + + self.CSMLT = np.zeros(5) + self.CSMLT[0] = ( + Float(2.0) + * constants.PIE + * constants.TCOND + * constants.RNZS + / constants.HLTF + ) + self.CSMLT[1] = ( + Float(2.0) + * constants.PIE + * constants.VDIFU + * constants.RNZS + * constants.HLTC + / constants.HLTF + ) + self.CSMLT[2] = self.CSSUB[1] + self.CSMLT[3] = self.CSSUB[2] + self.CSMLT[4] = constants.CH2O / constants.HLTF + + self.CSMLT_0 = self.CSMLT[0] + self.CSMLT_1 = self.CSMLT[1] + self.CSMLT_2 = self.CSMLT[2] + self.CSMLT_3 = self.CSMLT[3] + self.CSMLT_4 = self.CSMLT[4] + + # gmlt: five constants + + self.CGMLT = np.zeros(5) + self.CGMLT[0] = ( + Float(2.0) + * constants.PIE + * constants.TCOND + * constants.RNZG + / constants.HLTF + ) + self.CGMLT[1] = ( + Float(2.0) + * constants.PIE + * constants.VDIFU + * constants.RNZG + * constants.HLTC + / constants.HLTF + ) + self.CGMLT[2] = self.CGSUB[1] + self.CGMLT[3] = self.CGSUB[2] + self.CGMLT[4] = constants.CH2O / constants.HLTF + + self.CGMLT_0 = self.CGMLT[0] + self.CGMLT_1 = self.CGMLT[1] + self.CGMLT_2 = self.CGMLT[2] + self.CGMLT_3 = self.CGMLT[3] + self.CGMLT_4 = self.CGMLT[4] + + +class outputs: + def __init__(self, quantity_factory): + # ----------------------------------------------------------------------- + # initialize precipitation outputs + # ----------------------------------------------------------------------- + + self.rain = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") + self.snow = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") + self.ice = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") + self.graupel = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") + self.m2_rain = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.m2_sol = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.revap = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.isubl = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + + +class temporaries: + def __init__(self, quantity_factory): + self.t1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.dp1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.omq = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.qv0 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.ql0 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.qr0 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.qi0 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.qs0 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.qg0 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.qa0 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.qv1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.ql1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.qr1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.qi1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.qs1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.qg1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.qa1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.dz1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.den = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.den1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.denfac = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.p_dry = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.m1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.u1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.v1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.w1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.onemsig = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") + self.ccn = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.c_praut = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.rh_limited = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.ze = quantity_factory.zeros([X_DIM, Y_DIM, Z_INTERFACE_DIM], "n/a") + self.zt = quantity_factory.zeros([X_DIM, Y_DIM, Z_INTERFACE_DIM], "n/a") + self.lhi = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.icpk = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.hold_data = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.vti = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.vts = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.vtg = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.vtr = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.m1_sol = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.m1_rain = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.rain1 = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") + self.graupel1 = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") + self.snow1 = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") + self.ice1 = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") + self.evap1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.subl1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + + +class masks: + def __init__(self, quantity_factory): + self.is_frozen = quantity_factory.zeros( + [X_DIM, Y_DIM, Z_DIM], "n/a", dtype=bool + ) + self.precip_fall = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") + # TODO: temporary code requires mask used within a double k loop + # will be removed once a proper feature for double k loop is introduces + self.melting_mask_1 = quantity_factory.zeros( + [X_DIM, Y_DIM, Z_DIM], "n/a", dtype=bool + ) + self.melting_mask_2 = quantity_factory.zeros( + [X_DIM, Y_DIM, Z_DIM], "n/a", dtype=bool + ) + self.current_k_level = quantity_factory.zeros( + [X_DIM, Y_DIM, Z_DIM], "n/a", dtype=Int + ) + for k in range(self.current_k_level.view[:].shape[2]): + self.current_k_level.view[:, :, k] = k diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/terminal_fall.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/terminal_fall.py deleted file mode 100644 index 58e04ebc9..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/terminal_fall.py +++ /dev/null @@ -1,784 +0,0 @@ -import gt4py.cartesian.gtscript as gtscript -from gt4py.cartesian.gtscript import ( - BACKWARD, - FORWARD, - PARALLEL, - computation, - exp, - i32, - interval, -) - -import pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver_constants as driver_constants -from ndsl.dsl.typing import BoolField, Float, FloatField, FloatFieldIJ, IntField - - -@gtscript.function -def prefall_melting( - t: Float, - qv: Float, - ql: Float, - qr: Float, - qg: Float, - qs: Float, - qi: Float, - m1_sol: Float, - is_frozen: bool, -): - from __externals__ import ql_mlt, tau_imlt - - """ - melt cloud ice before fall - - reference Fortran: gfdl_cloud_microphys.F90: subroutine terminal_fall - """ - from __externals__ import c_air, c_vap, d0_vap, dts, lv00 - - fac_imlt = 1.0 - exp(-dts / tau_imlt) - - # ----------------------------------------------------------------------- - # define heat capacity and latent heat coefficient - # ----------------------------------------------------------------------- - - m1_sol = 0.0 - lhl = lv00 + d0_vap * t - lhi = driver_constants.LI00 + driver_constants.DC_ICE * t - q_liq = ql + qr - q_sol = qi + qs + qg - cvm = ( - c_air - + qv * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE - ) - lcpk = lhl / cvm - icpk = lhi / cvm - - # ----------------------------------------------------------------------- - # melting of cloud_ice (before fall) : - # ----------------------------------------------------------------------- - - tc = t - driver_constants.TICE - if is_frozen == False and (qi > driver_constants.QCMIN and tc > 0.0): # noqa - sink = min(qi, fac_imlt * tc / icpk) - if ql_mlt - ql > 0: - ans = ql_mlt - ql - else: - ans = 0 - tmp = min(sink, ans) - ql = ql + tmp - qr = qr + sink - tmp - qi = qi - sink - q_liq = q_liq + sink - q_sol = q_sol - sink - cvm = ( - c_air - + qv * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE - ) - t = t - sink * lhi / cvm - - return t, ql, qr, qi, cvm, lhi, icpk, m1_sol - - -def terminal_fall( - t1: FloatField, - qv1: FloatField, - ql1: FloatField, - qr1: FloatField, - qg1: FloatField, - qs1: FloatField, - qi1: FloatField, - dz1: FloatField, - dp1: FloatField, - den1: FloatField, - vtg: FloatField, - vts: FloatField, - vti: FloatField, - precip_rain: FloatFieldIJ, - precip_graupel: FloatFieldIJ, - precip_snow: FloatFieldIJ, - precip_ice: FloatFieldIJ, - m1_sol: FloatField, - w1: FloatField, - ze: FloatField, - zt: FloatField, - is_frozen: BoolField, - precip_fall: FloatFieldIJ, - melting_mask_1: BoolField, - melting_mask_2: BoolField, - current_k_level: IntField, -): - """ - calculate terminal fall speed, accounting for - melting of ice, snow, and graupel during fall - - reference Fortran: gfdl_cloud_microphys.F90: subroutine terminal_fall - """ - from __externals__ import do_sedi_w, dts, k_end, use_ppm, vi_fac - - # determine frozen levels - # later operations will only be executed if frozen/melted - # initalized to is_frozen = False, True = frozen, False = melted - with computation(PARALLEL), interval(...): - if t1 <= driver_constants.TICE: - is_frozen = True - - # we only want the melting layer closest to the surface - with computation(BACKWARD), interval(0, -1): - if is_frozen[0, 0, 1] == True and is_frozen[0, 0, 0] == False: # noqa - is_frozen = True - - # force surface to "melt" for later calculations - with computation(PARALLEL), interval(-1, None): - is_frozen = False - - with computation(PARALLEL), interval(...): - t1, ql1, qr1, qi1, cvm, lhi, icpk, m1_sol = prefall_melting( - t1, - qv1, - ql1, - qr1, - qg1, - qs1, - qi1, - m1_sol, - is_frozen, - ) - - # if timestep is too small turn off melting - with computation(PARALLEL), interval(0, -1): - if dts < 300.0: - is_frozen = True - disable_melt = True - else: - disable_melt = False - - with computation(PARALLEL), interval(-1, None): - if dts < 300.0: - is_frozen = False - disable_melt = True - else: - disable_melt = False - - with computation(BACKWARD), interval(...): - ze = ze[0, 0, 1] - dz1 # dz < 0 - - with computation(FORWARD), interval(0, 1): - zt = ze - - # ----------------------------------------------------------------------- - # update capacity heat and latend heat coefficient - # ----------------------------------------------------------------------- - - with computation(PARALLEL), interval(...): - if is_frozen == False: # noqa - lhi = driver_constants.LI00 + driver_constants.DC_ICE * t1 - icpk = lhi / cvm - - with computation(FORWARD), interval(-1, None): - lhi = driver_constants.LI00 + driver_constants.DC_ICE * t1 - icpk = lhi / cvm - - # ----------------------------------------------------------------------- - # ----------------------------------------------------------------------- - # melting of falling cloud ice into rain - # ----------------------------------------------------------------------- - # ----------------------------------------------------------------------- - - # reference Fortran: gfdl_cloud_microphys.F90: subroutine check_column - # determine if any precip falls in the column - # if it falls anywhere in the column, the entire column becomes true - # initalized to 0 (false), potentially changed to 1 (true) - with computation(FORWARD), interval(...): - if qi1 > driver_constants.QPMIN: - precip_fall = 1 - # end reference Fortran: gfdl_cloud_microphys.F90: subroutine check_column - - with computation(FORWARD), interval(0, 1): - if vi_fac < 1.0e-5 or precip_fall == 0: - precip_ice = 0 - - with computation(FORWARD), interval(1, None): - if vi_fac >= 1.0e-5 and precip_fall == 1: - zt = ze - dts * (vti[0, 0, -1] + vti) / 2.0 - - with computation(FORWARD), interval(-1, None): - if vi_fac >= 1.0e-5 and precip_fall == 1: - zt[0, 0, 1] = driver_constants.ZS - dts * vti - - with computation(FORWARD), interval(...): - if vi_fac >= 1.0e-5 and precip_fall == 1: - if zt[0, 0, 1] >= zt: - zt[0, 0, 1] = zt - driver_constants.DZ_MIN - - with computation(PARALLEL), interval(...): - if vi_fac >= 1.0e-5 and precip_fall == 1 and disable_melt == False: # noqa - # copy frozen mask to make modifyable version - melting_mask_1 = is_frozen - # flip logic for clarity - if melting_mask_1 == True: # noqa - melting_mask_1 = False - else: - melting_mask_1 = True - - with computation(BACKWARD), interval(-1, None): - if vi_fac >= 1.0e-5 and precip_fall == 1 and disable_melt == False: # noqa - # ensure no operations are performed on the surface - melting_mask_1 = False - - with computation(BACKWARD), interval(...): - # THIS HAS NEVER BEEN TESTED B/C DTS WAS LESS THAN 300 IN THE TEST CASE - if vi_fac >= 1.0e-5 and precip_fall == 1 and disable_melt == False: # noqa - # only operate on melted layers - if melting_mask_1 == True and qi1 > driver_constants.QCMIN: # noqa - # initalize exit trigger - stop_melting = False - m: i32 = 0 - while m < k_end and stop_melting == False: # noqa - mplus1: i32 = ( - m + 1 - ) # TODO remove this line only, replace with better solution - # only opterate on previously iterated k-levels - # if melting_mask_2 == True: #noqa - # if zt[0, 0, 1] >= ze.at(K=m): - # stop_melting = ( - # True # if true exit early for ONLY this k level - # ) - # if stop_melting == False: #noqa - # dtime = min( - # 1.0, - # (ze.at(K=m) - ze.at(K=mplus1)) - # / (max(driver_constants.vr_min, vti) * tau_imlt), - # ) - # sink = min( - # qi1 * dp1 / dp1.at(K=m), - # dtime - # * (t1.at(K=m) - driver_constants.tice) - # / icpk.at(K=m), - # ) - # if ql_mlt - ql1.at(K=m) > 0: - # ans = ql_mlt - ql1.at(K=m) - # else: - # ans = 0 - # tmp = min(sink, ans) - # offset = i32(m - current_k_level) - # hold_ql1 = ql1.at(K=m) - # ql1[0, 0, offset] = hold_ql1 + tmp - # hold_qr1 = qr1.at(K=m) - # qr1[0, 0, offset] = hold_qr1 - tmp + sink - # hold_t1 = t1.at(K=m) - # t1[0, 0, offset] = hold_t1 - sink * icpk.at(K=m) - # qi1 = qi1 - sink * dp1.at(K=m) / dp1 - m = m + 1 - # set current layer as iterated layer - melting_mask_2 = True - - with computation(PARALLEL), interval(...): - if vi_fac >= 1.0e-5 and precip_fall == 1: - if do_sedi_w == True: # noqa - dm = dp1 * (1.0 + qv1 + ql1 + qr1 + qi1 + qs1 + qg1) - - # begin reference Fortran: gfdl_cloud_microphys.F90: subroutine implicit_fall - # this code computes the time-implicit monotonic scheme - # Fortran author: Shian-Jiann Lin, 2016 - # set up inputs to the "function" - with computation(PARALLEL), interval(...): - if vi_fac >= 1.0e-5 and precip_fall == 1: - if use_ppm == False: # noqa - q_implicit_fall = qi1 - vt_implicit_fall = vti - - with computation(PARALLEL), interval(...): - if vi_fac >= 1.0e-5 and precip_fall == 1: - if use_ppm == False: # noqa - hold_data = ze - ze[0, 0, 1] - dd = dts * vt_implicit_fall - q_implicit_fall = q_implicit_fall * dp1 - - # ----------------------------------------------------------------------- - # sedimentation: non - vectorizable loop - # ----------------------------------------------------------------------- - with computation(FORWARD), interval(0, 1): - if vi_fac >= 1.0e-5 and precip_fall == 1: - if use_ppm == False: # noqa - qm = q_implicit_fall / (hold_data + dd) - - with computation(FORWARD), interval(1, None): - if vi_fac >= 1.0e-5 and precip_fall == 1: - if use_ppm == False: # noqa - qm = (q_implicit_fall + dd[0, 0, -1] * qm[0, 0, -1]) / (hold_data + dd) - - # ----------------------------------------------------------------------- - # qm is density at this stage - # ----------------------------------------------------------------------- - with computation(PARALLEL), interval(...): - if vi_fac >= 1.0e-5 and precip_fall == 1: - if use_ppm == False: # noqa - qm = qm * hold_data - - # ----------------------------------------------------------------------- - # output mass fluxes: non - vectorizable loop - # ----------------------------------------------------------------------- - with computation(FORWARD), interval(0, 1): - if vi_fac >= 1.0e-5 and precip_fall == 1: - if use_ppm == False: # noqa - m1 = q_implicit_fall - qm - - with computation(FORWARD), interval(1, None): - if vi_fac >= 1.0e-5 and precip_fall == 1: - if use_ppm == False: # noqa - m1 = m1[0, 0, -1] + q_implicit_fall - qm - - with computation(FORWARD), interval(-1, None): - if vi_fac >= 1.0e-5 and precip_fall == 1: - if use_ppm == False: # noqa - precip = m1 - - # ----------------------------------------------------------------------- - # update: - # ----------------------------------------------------------------------- - with computation(PARALLEL), interval(...): - if vi_fac >= 1.0e-5 and precip_fall == 1: - if use_ppm == False: # noqa - q_implicit_fall = qm / dp1 - - # update "outputs" after "function" - with computation(PARALLEL), interval(...): - if vi_fac >= 1.0e-5 and precip_fall == 1: - if use_ppm == False: # noqa - qi1 = q_implicit_fall - m1_sol = ( - m1_sol + m1 - ) # NOTE: setting this to just m1_sol = m1 gives WILD values (1e31) - - with computation(FORWARD), interval(-1, None): - if vi_fac >= 1.0e-5 and precip_fall == 1: - if use_ppm == False: # noqa - precip_ice = precip - # end reference Fortran: gfdl_cloud_microphys.F90: subroutine implicit_fall - - with computation(FORWARD), interval(0, 1): - if vi_fac >= 1.0e-5 and precip_fall == 1: - if do_sedi_w: - w1 = (dm * w1 + m1_sol * vti) / (dm - m1_sol) - - with computation(FORWARD), interval(1, None): - if vi_fac >= 1.0e-5 and precip_fall == 1: - if do_sedi_w: - w1 = (dm * w1 - m1_sol[0, 0, -1] * vti[0, 0, -1] + m1_sol * vti) / ( - dm + m1_sol[0, 0, -1] - m1_sol - ) - - # reset masks - with computation(FORWARD), interval(0, 1): - precip_fall = 0 - - with computation(PARALLEL), interval(...): - melting_mask_1 = False - melting_mask_2 = False - m1 = 0.0 - - # ----------------------------------------------------------------------- - # ----------------------------------------------------------------------- - # melting of falling snow into rain - # ----------------------------------------------------------------------- - # ----------------------------------------------------------------------- - - # reference Fortran: gfdl_cloud_microphys.F90: subroutine check_column - # determine if any precip falls in the column - # if it falls anywhere in the column, the entire column becomes true - # initalized to 0 (false), potentially changed to 1 (true) - with computation(FORWARD), interval(...): - if qs1 > driver_constants.QPMIN: - precip_fall = 1 - # end reference Fortran: gfdl_cloud_microphys.F90: subroutine check_column - - with computation(FORWARD), interval(0, 1): - if precip_fall == 0: - precip_snow = 0 - - with computation(FORWARD), interval(1, None): - if precip_fall == 1: - zt = ze - dts * (vts[0, 0, -1] + vts) / 2.0 - - with computation(FORWARD), interval(-1, None): - if precip_fall == 1: - zt[0, 0, 1] = driver_constants.ZS - dts * vts - - with computation(FORWARD), interval(...): - if precip_fall == 1: - if zt[0, 0, 1] >= zt: - zt[0, 0, 1] = zt - driver_constants.DZ_MIN - - with computation(PARALLEL), interval(...): - if precip_fall == 1 and disable_melt == False: # noqa - # copy frozen mask to make modifyable version - melting_mask_1 = is_frozen - # flip logic for clarity - if melting_mask_1 == True: # noqa - melting_mask_1 = False - else: - melting_mask_1 = True - - with computation(BACKWARD), interval(-1, None): - if precip_fall == 1 and disable_melt == False: # noqa - # ensure no operations are performed on the surface - melting_mask_1 = False - - with computation(BACKWARD), interval(...): - # THIS HAS NEVER BEEN TESTED B/C DTS WAS LESS THAN 300 IN THE TEST CASE - if precip_fall == 1 and disable_melt == False: # noqa - # only operate on melted layers - if melting_mask_1 == True and qs1 > driver_constants.QPMIN: # noqa - # initalize exit trigger - stop_melting = False - m: i32 = 0 - while m < k_end and stop_melting == False: # noqa - mplus1: i32 = ( - m + 1 - ) # TODO remove this line only, replace with better solution - # only opterate on previously iterated k-levels - # if melting_mask_2 == True: #noqa - # if zt[0, 0, 1] >= ze.at(K=m): - # stop_melting = ( - # True # if true exit early for ONLY this k level - # ) - # if stop_melting == False: #noqa - # dtime = min( - # dts, - # (ze.at(K=m) - ze.at(K=mplus1)) - # / (driver_constants.vr_min + vts), - # ) - # if ( - # zt < ze.at(K=mplus1) - # and t1.at(K=m) > driver_constants.tice - # ): - # dtime = min(1.0, dtime / tau_smlt) - # sink = min( - # qs1 * dp1 / dp1.at(K=m), - # dtime - # * (t1.at(K=m) - driver_constants.tice) - # / icpk.at(K=m), - # ) - # offset = i32(m - current_k_level) - # hold_t1 = t1.at(K=m) - # t1[0, 0, offset] = hold_t1 - sink * icpk.at(K=m) - # qs1 = qs1 - sink * dp1.at(K=m) / dp1 - # if zt < 0: - # precip_rain = precip_rain + sink * dp1.at( - # K=m - # ) # precip as rain - # else: - # # qr source here will fall next time step - # # (therefore, can evap) - # hold_qr1 = qr1.at(K=m) - # qr1[0, 0, offset] = hold_qr1 + sink - # if stop_melting == False: #noqa - # if qs1 < driver_constants.qpmin: - # stop_melting = True - m = m + 1 - # set current layer as iterated layer - melting_mask_2 = True - - with computation(PARALLEL), interval(...): - if precip_fall == 1: - if do_sedi_w == True: # noqa - dm = dp1 * (1.0 + qv1 + ql1 + qr1 + qi1 + qs1 + qg1) - - # begin reference Fortran: gfdl_cloud_microphys.F90: subroutine implicit_fall - # this code computes the time-implicit monotonic scheme - # Fortran author: Shian-Jiann Lin, 2016 - # set up inputs to the "function" - with computation(PARALLEL), interval(...): - if precip_fall == 1: - if use_ppm == False: # noqa - q_implicit_fall = qs1 - vt_implicit_fall = vts - - with computation(PARALLEL), interval(...): - if precip_fall == 1: - if use_ppm == False: # noqa - hold_data = ze - ze[0, 0, 1] - dd = dts * vt_implicit_fall - q_implicit_fall = q_implicit_fall * dp1 - - # ----------------------------------------------------------------------- - # sedimentation: non - vectorizable loop - # ----------------------------------------------------------------------- - with computation(FORWARD), interval(0, 1): - if precip_fall == 1: - if use_ppm == False: # noqa - qm = q_implicit_fall / (hold_data + dd) - - with computation(FORWARD), interval(1, None): - if precip_fall == 1: - if use_ppm == False: # noqa - qm = (q_implicit_fall + dd[0, 0, -1] * qm[0, 0, -1]) / (hold_data + dd) - - # ----------------------------------------------------------------------- - # qm is density at this stage - # ----------------------------------------------------------------------- - with computation(PARALLEL), interval(...): - if precip_fall == 1: - if use_ppm == False: # noqa - qm = qm * hold_data - - # ----------------------------------------------------------------------- - # output mass fluxes: non - vectorizable loop - # ----------------------------------------------------------------------- - with computation(FORWARD), interval(0, 1): - if precip_fall == 1: - if use_ppm == False: # noqa - m1 = q_implicit_fall - qm - - with computation(FORWARD), interval(1, None): - if precip_fall == 1: - if use_ppm == False: # noqa - m1 = m1[0, 0, -1] + q_implicit_fall - qm - - with computation(FORWARD), interval(-1, None): - if precip_fall == 1: - if use_ppm == False: # noqa - precip = m1 - - # ----------------------------------------------------------------------- - # update: - # ----------------------------------------------------------------------- - with computation(PARALLEL), interval(...): - if precip_fall == 1: - if use_ppm == False: # noqa - q_implicit_fall = qm / dp1 - - # update "outputs" after "function" - with computation(PARALLEL), interval(...): - if precip_fall == 1: - if use_ppm == False: # noqa - qs1 = q_implicit_fall - m1_sol = ( - m1_sol + m1 - ) # NOTE: setting this to just m1_sol = m1 gives WILD values (1e31) - - with computation(FORWARD), interval(-1, None): - if precip_fall == 1: - if use_ppm == False: # noqa - precip_snow = precip - # end reference Fortran: gfdl_cloud_microphys.F90: subroutine implicit_fall - - with computation(FORWARD), interval(0, 1): - if precip_fall == 1: - if do_sedi_w: - w1 = (dm * w1 + m1 * vts) / (dm - m1) - - with computation(FORWARD), interval(1, None): - if precip_fall == 1: - if do_sedi_w: - w1 = (dm * w1 - m1[0, 0, -1] * vts[0, 0, -1] + m1 * vts) / ( - dm + m1[0, 0, -1] - m1 - ) - - # reset masks and temporaries - with computation(FORWARD), interval(0, 1): - precip_fall = 0 - - with computation(PARALLEL), interval(...): - melting_mask_1 = False - melting_mask_2 = False - # must be reset to zero everywhere in case there is no graupel - # in a column and no calculations are performed - m1 = 0.0 - - # ----------------------------------------------------------------------- - # ----------------------------------------------------------------------- - # melting of falling graupel into rain - # ----------------------------------------------------------------------- - # ----------------------------------------------------------------------- - - # reference Fortran: gfdl_cloud_microphys.F90: subroutine check_column - # determine if any precip falls in the column - # if it falls anywhere in the column, the entire column becomes true - # initalized to 0 (false), potentially changed to 1 (true) - with computation(FORWARD), interval(...): - if qg1 > driver_constants.QPMIN: - precip_fall = 1 - # end reference Fortran: gfdl_cloud_microphys.F90: subroutine check_column - - with computation(FORWARD), interval(0, 1): - if precip_fall == 0: - precip_graupel = 0 - - with computation(FORWARD), interval(1, None): - if precip_fall == 1: - zt = ze - dts * (vtg[0, 0, -1] + vtg) / 2.0 - - with computation(FORWARD), interval(-1, None): - if precip_fall == 1: - zt[0, 0, 1] = driver_constants.ZS - dts * vtg - - with computation(FORWARD), interval(...): - if precip_fall == 1: - if zt[0, 0, 1] >= zt: - zt[0, 0, 1] = zt - driver_constants.DZ_MIN - - with computation(PARALLEL), interval(...): - if precip_fall == 1 and disable_melt == False: # noqa - # copy frozen mask to make modifyable version - melting_mask_1 = is_frozen - # flip logic for clarity - if melting_mask_1 == True: # noqa - melting_mask_1 = False - else: - melting_mask_1 = True - - with computation(BACKWARD), interval(-1, None): - if precip_fall == 1 and disable_melt == False: # noqa - # ensure no operations are performed on the surface - melting_mask_1 = False - - with computation(BACKWARD), interval(...): - # THIS HAS NEVER BEEN TESTED B/C DTS WAS LESS THAN 300 IN THE TEST CASE - if precip_fall == 1 and disable_melt == False: # noqa - # only operate on melted layers - if melting_mask_1 == True and qg1 > driver_constants.QPMIN: # noqa - # initalize exit trigger - stop_melting = False - m: i32 = 0 - while m < k_end and stop_melting == False: # noqa - mplus1: i32 = ( - m + 1 - ) # TODO remove this line only, replace with better solution - # only opterate on previously iterated k-levels - # if melting_mask_2 == True: #noqa - # if zt[0, 0, 1] >= ze.at(K=m): - # stop_melting = ( - # True # if true exit early for ONLY this k level - # ) - # if stop_melting == False: #noqa - # dtime = min(dts, (ze.at(K=m) - ze.at(K=mplus1)) / vtg) - # if ( - # zt < ze.at(K=mplus1) - # and t1.at(K=m) > driver_constants.tice - # ): - # dtime = min(1.0, dtime / tau_g2r) - # sink = min( - # qg1 * dp1 / dp1.at(K=m), - # dtime - # * (t1.at(K=m) - driver_constants.tice) - # / icpk.at(K=m), - # ) - # offset = i32(m - current_k_level) - # hold_t1 = t1.at(K=m) - # t1[0, 0, offset] = hold_t1 - sink * icpk.at(K=m) - # qg1 = qg1 - sink * dp1.at(K=m) / dp1 - # if zt < 0: - # precip_rain = precip_rain + sink * dp1.at( - # K=m - # ) # precip as rain - # else: - # # qr source here will fall next time step - # # (therefore, can evap) - # hold_qr1 = qr1.at(K=m) - # qr1[0, 0, offset] = hold_qr1 + sink - # if stop_melting == False: #noqa - # if qg1 < driver_constants.qpmin: - # stop_melting = True - m = m + 1 - # set current layer as iterated layer - melting_mask_2 = True - - with computation(PARALLEL), interval(...): - if precip_fall == 1: - if do_sedi_w == True: # noqa - dm = dp1 * (1.0 + qv1 + ql1 + qr1 + qi1 + qs1 + qg1) - - # begin reference Fortran: gfdl_cloud_microphys.F90: subroutine implicit_fall - # this code computes the time-implicit monotonic scheme - # Fortran author: Shian-Jiann Lin, 2016 - # set up inputs to the "function" - with computation(PARALLEL), interval(...): - if precip_fall == 1: - if use_ppm == False: # noqa - q_implicit_fall = qg1 - vt_implicit_fall = vtg - - with computation(PARALLEL), interval(...): - if precip_fall == 1: - if use_ppm == False: # noqa - hold_data = ze - ze[0, 0, 1] - dd = dts * vt_implicit_fall - q_implicit_fall = q_implicit_fall * dp1 - - # ----------------------------------------------------------------------- - # sedimentation: non - vectorizable loop - # ----------------------------------------------------------------------- - with computation(FORWARD), interval(0, 1): - if precip_fall == 1: - if use_ppm == False: # noqa - qm = q_implicit_fall / (hold_data + dd) - - with computation(FORWARD), interval(1, None): - if precip_fall == 1: - if use_ppm == False: # noqa - qm = (q_implicit_fall + dd[0, 0, -1] * qm[0, 0, -1]) / (hold_data + dd) - - # ----------------------------------------------------------------------- - # qm is density at this stage - # ----------------------------------------------------------------------- - with computation(PARALLEL), interval(...): - if precip_fall == 1: - if use_ppm == False: # noqa - qm = qm * hold_data - - # ----------------------------------------------------------------------- - # output mass fluxes: non - vectorizable loop - # ----------------------------------------------------------------------- - with computation(FORWARD), interval(0, 1): - if precip_fall == 1: - if use_ppm == False: # noqa - m1 = q_implicit_fall - qm - - with computation(FORWARD), interval(1, None): - if precip_fall == 1: - if use_ppm == False: # noqa - m1 = m1[0, 0, -1] + q_implicit_fall - qm - - with computation(FORWARD), interval(-1, None): - if precip_fall == 1: - if use_ppm == False: # noqa - precip = m1 - - # ----------------------------------------------------------------------- - # update: - # ----------------------------------------------------------------------- - with computation(PARALLEL), interval(...): - if precip_fall == 1: - if use_ppm == False: # noqa - q_implicit_fall = qm / dp1 - - # update "outputs" after "function" - with computation(PARALLEL), interval(...): - if precip_fall == 1: - if use_ppm == False: # noqa - qg1 = q_implicit_fall - m1_sol = ( - m1_sol + m1 - ) # NOTE: setting this to just m1_sol = m1 gives WILD values (1e31) - - with computation(FORWARD), interval(-1, None): - if precip_fall == 1: - if use_ppm == False: # noqa - precip_graupel = precip - # end reference Fortran: gfdl_cloud_microphys.F90: subroutine implicit_fall - - with computation(FORWARD), interval(0, 1): - if precip_fall == 1: - if do_sedi_w: - w1 = (dm * w1 + m1 * vtg) / (dm - m1) - - with computation(FORWARD), interval(1, None): - if precip_fall == 1: - if do_sedi_w: - w1 = (dm * w1 - m1[0, 0, -1] * vtg[0, 0, -1] + m1 * vtg) / ( - dm + m1[0, 0, -1] - m1 - ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/terminal_fall/stencils.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/terminal_fall/stencils.py new file mode 100644 index 000000000..f639ac80b --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/terminal_fall/stencils.py @@ -0,0 +1,416 @@ +import gt4py.cartesian.gtscript as gtscript +from gt4py.cartesian.gtscript import ( + BACKWARD, + FORWARD, + PARALLEL, + computation, + exp, + i32, + interval, +) + +import pyMoist.GFDL_1M.GFDL_1M_driver.constants as constants +from ndsl.dsl.typing import BoolField, Float, FloatField, FloatFieldIJ, IntField + + +def check_precip_get_zt( + q: FloatField, + vt: FloatField, + ze: FloatField, + zt: FloatField, + precip: FloatFieldIJ, + precip_fall: FloatFieldIJ, +): + from __externals__ import dts + + # ----------------------------------------------------------------------- + # ----------------------------------------------------------------------- + # melting of falling snow into rain + # ----------------------------------------------------------------------- + # ----------------------------------------------------------------------- + + # reference Fortran: gfdl_cloud_microphys.F90: subroutine check_column + # determine if any precip falls in the column + # if it falls anywhere in the column, the entire column becomes true + # initalized to 0 (false), potentially changed to 1 (true) + with computation(FORWARD), interval(...): + if q > constants.QPMIN: + precip_fall = 1 + # end reference Fortran: gfdl_cloud_microphys.F90: subroutine check_column + + with computation(FORWARD), interval(0, 1): + if precip_fall == 0: + precip = 0 + + with computation(FORWARD), interval(1, None): + if precip_fall == 1: + zt = ze - dts * (vt[0, 0, -1] + vt) / 2.0 + + with computation(FORWARD), interval(-1, None): + if precip_fall == 1: + zt[0, 0, 1] = constants.ZS - dts * vt + + with computation(FORWARD), interval(...): + if precip_fall == 1: + if zt[0, 0, 1] >= zt: + zt[0, 0, 1] = zt - constants.DZ_MIN + + +def melting_loop( + # Features required before this can be implemented: + need_2d_temporaries_feature: FloatFieldIJ, + need_double_k_loop_feature: FloatField, +): + + with computation(BACKWARD), interval(...): + # NOTE: code is unfinished, not implemented. do not allow it to run + # this code should never execute anyway. It requries dts > 300, and if this is true + # an error should be triggered in GFDL_1M that prevents execution + disable = True + + if disable == False: + skip_for_now = 1 + """ + # THIS HAS NEVER BEEN TESTED B/C DTS WAS LESS THAN 300 IN THE TEST CASE + if precip_fall == 1 and disable_melt == False: # noqa + # only operate on melted layers + if melting_mask_1 == True and qg1 > constants.QPMIN: # noqa + # initalize exit trigger + stop_melting = False + m: i32 = 0 + while m < k_end and stop_melting == False: # noqa + mplus1: i32 = ( + m + 1 + ) # TODO remove this line only, replace with better solution + # only opterate on previously iterated k-levels + # if melting_mask_2 == True: #noqa + # if zt[0, 0, 1] >= ze.at(K=m): + # stop_melting = ( + # True # if true exit early for ONLY this k level + # ) + # if stop_melting == False: #noqa + # dtime = min(dts, (ze.at(K=m) - ze.at(K=mplus1)) / vtg) + # if ( + # zt < ze.at(K=mplus1) + # and t1.at(K=m) > driver_constants.tice + # ): + # dtime = min(1.0, dtime / tau_g2r) + # sink = min( + # qg1 * dp1 / dp1.at(K=m), + # dtime + # * (t1.at(K=m) - driver_constants.tice) + # / icpk.at(K=m), + # ) + # offset = i32(m - current_k_level) + # hold_t1 = t1.at(K=m) + # t1[0, 0, offset] = hold_t1 - sink * icpk.at(K=m) + # qg1 = qg1 - sink * dp1.at(K=m) / dp1 + # if zt < 0: + # precip_rain = precip_rain + sink * dp1.at( + # K=m + # ) # precip as rain + # else: + # # qr source here will fall next time step + # # (therefore, can evap) + # hold_qr1 = qr1.at(K=m) + # qr1[0, 0, offset] = hold_qr1 + sink + # if stop_melting == False: #noqa + # if qg1 < driver_constants.qpmin: + # stop_melting = True + m = m + 1 + # set current layer as iterated layer + melting_mask_2 = True + """ + + +def update_dm( + dm: FloatField, + dp1: FloatField, + qv1: FloatField, + ql1: FloatField, + qr1: FloatField, + qi1: FloatField, + qs1: FloatField, + qg1: FloatField, + precip_fall: FloatFieldIJ, +): + from __externals__ import do_sedi_w + + with computation(PARALLEL), interval(...): + if precip_fall == 1: + if do_sedi_w == True: # noqa + dm = dp1 * (1.0 + qv1 + ql1 + qr1 + qi1 + qs1 + qg1) + + +def implicit_fall( + q: FloatField, + vt: FloatField, + ze: FloatField, + dp1: FloatField, + m1: FloatField, + m1_sol: FloatField, + precip: FloatFieldIJ, + precip_fall: FloatFieldIJ, +): + from __externals__ import dts, use_ppm + + # begin reference Fortran: gfdl_cloud_microphys.F90: subroutine implicit_fall + # this code computes the time-implicit monotonic scheme + # Fortran author: Shian-Jiann Lin, 2016 + with computation(PARALLEL), interval(...): + if precip_fall == 1: + if use_ppm == False: # noqa + height_diff = ze - ze[0, 0, 1] + dd = dts * vt + q = q * dp1 + + # ----------------------------------------------------------------------- + # sedimentation: non - vectorizable loop + # ----------------------------------------------------------------------- + with computation(FORWARD), interval(0, 1): + if precip_fall == 1: + if use_ppm == False: # noqa + qm = q / (height_diff + dd) + + with computation(FORWARD), interval(1, None): + if precip_fall == 1: + if use_ppm == False: # noqa + qm = (q + dd[0, 0, -1] * qm[0, 0, -1]) / (height_diff + dd) + + # ----------------------------------------------------------------------- + # qm is density at this stage + # ----------------------------------------------------------------------- + with computation(PARALLEL), interval(...): + if precip_fall == 1: + if use_ppm == False: # noqa + qm = qm * height_diff + + # ----------------------------------------------------------------------- + # output mass fluxes: non - vectorizable loop + # ----------------------------------------------------------------------- + with computation(FORWARD), interval(0, 1): + if precip_fall == 1: + if use_ppm == False: # noqa + m1 = q - qm + + with computation(FORWARD), interval(1, None): + if precip_fall == 1: + if use_ppm == False: # noqa + m1 = m1[0, 0, -1] + q - qm + + with computation(FORWARD), interval(-1, None): + if precip_fall == 1: + if use_ppm == False: # noqa + precip = m1 + + # ----------------------------------------------------------------------- + # update: + # ----------------------------------------------------------------------- + with computation(PARALLEL), interval(...): + if precip_fall == 1: + if use_ppm == False: # noqa + q = qm / dp1 + + # end reference Fortran: gfdl_cloud_microphys.F90: subroutine implicit_fall + + with computation(PARALLEL), interval(...): + if precip_fall == 1: + m1_sol = m1_sol + m1 + + +def update_w1( + w1: FloatField, + dm: FloatField, + m1: FloatField, + vt: FloatField, + precip_fall: FloatFieldIJ, +): + from __externals__ import do_sedi_w + + with computation(FORWARD), interval(0, 1): + if precip_fall == 1: + if do_sedi_w: + w1 = (dm * w1 + m1 * vt) / (dm - m1) + + with computation(FORWARD), interval(1, None): + if precip_fall == 1: + if do_sedi_w: + w1 = (dm * w1 - m1[0, 0, -1] * vt[0, 0, -1] + m1 * vt) / ( + dm + m1[0, 0, -1] - m1 + ) + + +def reset( + m1: FloatField, + precip_fall: FloatFieldIJ, +): + # reset masks and temporaries + with computation(FORWARD), interval(0, 1): + precip_fall = 0 + + with computation(PARALLEL), interval(...): + # must be reset to zero everywhere in case there is no precipitate + # in a column and no calculations are performed + m1 = 0.0 + + +@gtscript.function +def prefall_melting( + t: Float, + qv: Float, + ql: Float, + qr: Float, + qg: Float, + qs: Float, + qi: Float, + m1_sol: Float, + is_frozen: bool, + c_air: Float, + c_vap: Float, + d0_vap: Float, + dts: Float, + lv00: Float, + ql_mlt: Float, + tau_imlt: Float, +): + """ + melt cloud ice before fall + + reference Fortran: gfdl_cloud_microphys.F90: subroutine terminal_fall + """ + + fac_imlt = 1.0 - exp(-dts / tau_imlt) + + # ----------------------------------------------------------------------- + # define heat capacity and latent heat coefficient + # ----------------------------------------------------------------------- + + m1_sol = 0.0 + lhl = lv00 + d0_vap * t + lhi = constants.LI00 + constants.DC_ICE * t + q_liq = ql + qr + q_sol = qi + qs + qg + cvm = c_air + qv * c_vap + q_liq * constants.C_LIQ + q_sol * constants.C_ICE + lcpk = lhl / cvm + icpk = lhi / cvm + + # ----------------------------------------------------------------------- + # melting of cloud_ice (before fall) : + # ----------------------------------------------------------------------- + + tc = t - constants.TICE + if is_frozen == False and (qi > constants.QCMIN and tc > 0.0): # noqa + sink = min(qi, fac_imlt * tc / icpk) + if ql_mlt - ql > 0: + ans = ql_mlt - ql + else: + ans = 0 + tmp = min(sink, ans) + ql = ql + tmp + qr = qr + sink - tmp + qi = qi - sink + q_liq = q_liq + sink + q_sol = q_sol - sink + cvm = c_air + qv * c_vap + q_liq * constants.C_LIQ + q_sol * constants.C_ICE + t = t - sink * lhi / cvm + + return t, ql, qr, qi, cvm, lhi, icpk, m1_sol + + +def setup( + t1: FloatField, + qv1: FloatField, + ql1: FloatField, + qr1: FloatField, + qg1: FloatField, + qs1: FloatField, + qi1: FloatField, + dz1: FloatField, + m1_sol: FloatField, + ze: FloatField, + zt: FloatField, + lhi: FloatField, + icpk: FloatField, + cvm: FloatField, + is_frozen: BoolField, +): + """ + calculate terminal fall speed, accounting for + melting of ice, snow, and graupel during fall + + reference Fortran: gfdl_cloud_microphys.F90: subroutine terminal_fall + """ + from __externals__ import ( + dts, + c_air, + c_vap, + d0_vap, + dts, + lv00, + ql_mlt, + tau_imlt, + ) + + # determine frozen levels + # later operations will only be executed if frozen/melted + # initalized to is_frozen = False, True = frozen, False = melted + with computation(PARALLEL), interval(...): + if t1 <= constants.TICE: + is_frozen = True + + # we only want the melting layer closest to the surface + with computation(BACKWARD), interval(0, -1): + if is_frozen[0, 0, 1] == True and is_frozen[0, 0, 0] == False: # noqa + is_frozen = True + + # force surface to "melt" for later calculations + with computation(PARALLEL), interval(-1, None): + is_frozen = False + + with computation(PARALLEL), interval(...): + t1, ql1, qr1, qi1, cvm, lhi, icpk, m1_sol = prefall_melting( + t1, + qv1, + ql1, + qr1, + qg1, + qs1, + qi1, + m1_sol, + is_frozen, + c_air, + c_vap, + d0_vap, + dts, + lv00, + ql_mlt, + tau_imlt, + ) + + # if timestep is too small turn off melting (only calculate at surface) + with computation(PARALLEL), interval(0, -1): + if dts < 300.0: + is_frozen = True + + with computation(PARALLEL), interval(-1, None): + if dts < 300.0: + is_frozen = False + + with computation(BACKWARD), interval(...): + ze = ze[0, 0, 1] - dz1 # dz < 0 + + with computation(FORWARD), interval(0, 1): + zt = ze + + # ----------------------------------------------------------------------- + # update capacity heat and latend heat coefficient + # ----------------------------------------------------------------------- + + with computation(PARALLEL), interval(...): + if is_frozen == False: # noqa + lhi = constants.LI00 + constants.DC_ICE * t1 + icpk = lhi / cvm + + with computation(FORWARD), interval(-1, None): + lhi = constants.LI00 + constants.DC_ICE * t1 + icpk = lhi / cvm diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/terminal_fall/temporaries.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/terminal_fall/temporaries.py new file mode 100644 index 000000000..ff9488bf2 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/terminal_fall/temporaries.py @@ -0,0 +1,22 @@ +import numpy as np +from gt4py.cartesian.gtscript import i32 +from ndsl.dsl.typing import Float +from ndsl.constants import X_DIM, Y_DIM, Z_DIM, Z_INTERFACE_DIM +import pyMoist.GFDL_1M.GFDL_1M_driver.constants as constants +from pyMoist.GFDL_1M.GFDL_1M_driver.config import config +from ndsl.dsl.typing import Int + + +class temporaries: + def __init__(self, quantity_factory): + self.lhi = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.icpk = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.cvm = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.m1 = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + self.dm = quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") + + # used as placeholders to note what features are needed. will be r + self.need_2d_temporaries_feature = quantity_factory.zeros([X_DIM, Y_DIM], "n/a") + self.need_double_k_loop_feature = quantity_factory.zeros( + [X_DIM, Y_DIM, Z_DIM], "n/a" + ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/terminal_fall/terminal_fall.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/terminal_fall/terminal_fall.py new file mode 100644 index 000000000..2c816ab25 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/terminal_fall/terminal_fall.py @@ -0,0 +1,309 @@ +import numpy as np +from gt4py.cartesian.gtscript import i32 + +import pyMoist.GFDL_1M.GFDL_1M_driver.constants as constants +from ndsl import QuantityFactory, StencilFactory, orchestrate +from ndsl.constants import X_DIM, Y_DIM, Z_DIM, Z_INTERFACE_DIM +from ndsl.dsl.typing import Float, FloatField, FloatFieldIJ, Int +from pyMoist.GFDL_1M.GFDL_1M_driver.config import config +from pyMoist.GFDL_1M.GFDL_1M_driver.terminal_fall.temporaries import temporaries +from pyMoist.GFDL_1M.GFDL_1M_driver.terminal_fall.stencils import ( + setup, + check_precip_get_zt, + melting_loop, + update_dm, + implicit_fall, + update_w1, + reset, +) + + +class terminal_fall: + """ + calculate terminal fall speed, accounting for + melting of ice, snow, and graupel during fall + + reference Fortran: gfdl_cloud_microphys.F90: subroutine terminal_fall + """ + + def __init__( + self, + stencil_factory: StencilFactory, + quantity_factory: QuantityFactory, + GFDL_1M_config: config, + config_dependent_constants, + ): + + self.GFDL_1M_config = GFDL_1M_config + + # Initalize temporaries + self.temporaries = temporaries(quantity_factory) + + # Initalize stencils + self._setup = stencil_factory.from_dims_halo( + func=setup, + compute_dims=[X_DIM, Y_DIM, Z_DIM], + externals={ + "dts": config_dependent_constants.DTS, + "tau_imlt": GFDL_1M_config.tau_imlt, + "ql_mlt": GFDL_1M_config.ql_mlt, + "c_air": config_dependent_constants.C_AIR, + "c_vap": config_dependent_constants.C_VAP, + "d0_vap": config_dependent_constants.D0_VAP, + "lv00": config_dependent_constants.LV00, + }, + ) + self._check_precip_get_zt = stencil_factory.from_dims_halo( + func=check_precip_get_zt, + compute_dims=[X_DIM, Y_DIM, Z_DIM], + externals={ + "dts": config_dependent_constants.DTS, + }, + ) + self._melting_loop = stencil_factory.from_dims_halo( + func=melting_loop, + compute_dims=[X_DIM, Y_DIM, Z_DIM], + ) + self._update_dm = stencil_factory.from_dims_halo( + func=update_dm, + compute_dims=[X_DIM, Y_DIM, Z_DIM], + externals={ + "do_sedi_w": GFDL_1M_config.do_sedi_w, + }, + ) + self._implicit_fall = stencil_factory.from_dims_halo( + func=implicit_fall, + compute_dims=[X_DIM, Y_DIM, Z_DIM], + externals={ + "dts": config_dependent_constants.DTS, + "use_ppm": GFDL_1M_config.use_ppm, + }, + ) + self._update_w1 = stencil_factory.from_dims_halo( + func=update_w1, + compute_dims=[X_DIM, Y_DIM, Z_DIM], + externals={ + "do_sedi_w": GFDL_1M_config.do_sedi_w, + }, + ) + self._reset = stencil_factory.from_dims_halo( + func=reset, + compute_dims=[X_DIM, Y_DIM, Z_DIM], + ) + + def __call__( + self, + t1, + qv1, + ql1, + qr1, + qg1, + qs1, + qi1, + dz1, + dp1, + den1, + vtg, + vts, + vti, + precip_rain, + precip_graupel, + precip_snow, + precip_ice, + m1_sol, + w1, + ze, + zt, + is_frozen, + precip_fall, + ): + self._setup( + t1, + qv1, + ql1, + qr1, + qg1, + qs1, + qi1, + dz1, + m1_sol, + ze, + zt, + self.temporaries.lhi, + self.temporaries.icpk, + self.temporaries.cvm, + is_frozen, + ) + + # ----------------------------------------------------------------------- + # ----------------------------------------------------------------------- + # melting of falling cloud ice into rain + # ----------------------------------------------------------------------- + # ----------------------------------------------------------------------- + if self.GFDL_1M_config.vi_fac < 1.0e-5: + precip_ice.view[:] = 0 + else: + self._check_precip_get_zt( + qi1, + vti, + ze, + zt, + precip_ice, + precip_fall, + ) + + # placeholder function. need the listed features to implement + self._melting_loop( + self.temporaries.need_2d_temporaries_feature, + self.temporaries.need_double_k_loop_feature, + ) + + self._update_dm( + self.temporaries.dm, + dp1, + qv1, + ql1, + qr1, + qi1, + qs1, + qg1, + precip_fall, + ) + + self._implicit_fall( + qi1, + vti, + ze, + dp1, + self.temporaries.m1, + m1_sol, + precip_ice, + precip_fall, + ) + + self._update_w1( + w1, + self.temporaries.dm, + self.temporaries.m1, + vti, + precip_fall, + ) + + self._reset( + self.temporaries.m1, + precip_fall, + ) + + # ----------------------------------------------------------------------- + # ----------------------------------------------------------------------- + # melting of falling snow into rain + # ----------------------------------------------------------------------- + # ----------------------------------------------------------------------- + + self._check_precip_get_zt( + qs1, + vts, + ze, + zt, + precip_snow, + precip_fall, + ) + + # placeholder function. need the listed features to implement + self._melting_loop( + self.temporaries.need_2d_temporaries_feature, + self.temporaries.need_double_k_loop_feature, + ) + + self._update_dm( + self.temporaries.dm, + dp1, + qv1, + ql1, + qr1, + qi1, + qs1, + qg1, + precip_fall, + ) + + self._implicit_fall( + qs1, + vts, + ze, + dp1, + self.temporaries.m1, + m1_sol, + precip_snow, + precip_fall, + ) + + self._update_w1( + w1, + self.temporaries.dm, + self.temporaries.m1, + vts, + precip_fall, + ) + + self._reset( + self.temporaries.m1, + precip_fall, + ) + + # ----------------------------------------------------------------------- + # ----------------------------------------------------------------------- + # melting of falling graupel into rain + # ----------------------------------------------------------------------- + # ----------------------------------------------------------------------- + + self._check_precip_get_zt( + qg1, + vtg, + ze, + zt, + precip_graupel, + precip_fall, + ) + + # placeholder function. need the listed features to implement + self._melting_loop( + self.temporaries.need_2d_temporaries_feature, + self.temporaries.need_double_k_loop_feature, + ) + + self._update_dm( + self.temporaries.dm, + dp1, + qv1, + ql1, + qr1, + qi1, + qs1, + qg1, + precip_fall, + ) + + self._implicit_fall( + qg1, + vtg, + ze, + dp1, + self.temporaries.m1, + m1_sol, + precip_graupel, + precip_fall, + ) + + self._update_w1( + w1, + self.temporaries.dm, + self.temporaries.m1, + vtg, + precip_fall, + ) + + self._reset( + self.temporaries.m1, + precip_fall, + ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/warm_rain.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/warm_rain.py index 8bb7c8171..dce5aeb61 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/warm_rain.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/GFDL_1M/GFDL_1M_driver/warm_rain.py @@ -13,11 +13,11 @@ trunc, ) -import pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver_constants as driver_constants +import pyMoist.GFDL_1M.GFDL_1M_driver.constants as constants from ndsl.dsl.typing import Float, FloatField, FloatFieldIJ -GlobalTable_driver_qsat = gtscript.GlobalTable[(Float, (int(driver_constants.LENGTH)))] +GlobalTable_driver_qsat = gtscript.GlobalTable[(Float, (int(constants.LENGTH)))] @gtscript.function @@ -36,7 +36,7 @@ def wqs2( reference Fortran: gfdl_cloud_microphys.F90: function wqs2 """ - tmin = driver_constants.TABLE_ICE - 160.0 + tmin = constants.TABLE_ICE - 160.0 if ta - tmin > 0: ans = ta - tmin @@ -46,7 +46,7 @@ def wqs2( ap1 = min(2621.0, ap1) it = i32(trunc(ap1)) es = table2.A[it - 1] + (ap1 - it) * des2.A[it - 1] - qsat = es / (driver_constants.RVGAS * ta * den) + qsat = es / (constants.RVGAS * ta * den) it = i32( trunc(ap1 - 0.5) ) # check if this rounds or truncates. need truncation here @@ -54,7 +54,7 @@ def wqs2( dqdt = ( 10.0 * (des2.A[it - 1] + (ap1 - it) * (des2.A[it] - des2.A[it - 1])) - / (driver_constants.RVGAS * ta * den) + / (constants.RVGAS * ta * den) ) return qsat, dqdt @@ -104,7 +104,7 @@ def revap_racc( revap = 0.0 - if t1 > driver_constants.T_WFR and qr1 > driver_constants.QPMIN: + if t1 > constants.T_WFR and qr1 > constants.QPMIN: # area and timescale efficiency on revap fac_revp = 1.0 - exp(-half_dt / tau_revp) @@ -118,15 +118,15 @@ def revap_racc( cvm = ( c_air + qv1 * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) lcpk = lhl / cvm tin = t1 - lcpk * ql1 # presence of clouds suppresses the rain evap qpz = qv1 + ql1 qsat, dqsdt = wqs2(tin, den1, table2, des2) - dqh = max(ql1, rh_limited * max(qpz, driver_constants.QCMIN)) + dqh = max(ql1, rh_limited * max(qpz, constants.QCMIN)) dqh = min(dqh, 0.2 * qpz) # new limiter dqv = qsat - qv1 # use this to prevent super - sat the gird box q_minus = qpz - dqh @@ -141,7 +141,7 @@ def revap_racc( # rain evaporation # ----------------------------------------------------------------------- - if dqv > driver_constants.QVMIN and qsat > q_minus: + if dqv > constants.QVMIN and qsat > q_minus: if qsat > q_plus: dq = qsat - qpz else: @@ -166,8 +166,8 @@ def revap_racc( cvm = ( c_air + qv1 * c_vap - + q_liq * driver_constants.C_LIQ - + q_sol * driver_constants.C_ICE + + q_liq * constants.C_LIQ + + q_sol * constants.C_ICE ) t1 = t1 - evap * lhl / cvm revap = evap / half_dt @@ -177,8 +177,8 @@ def revap_racc( # ----------------------------------------------------------------------- if ( - qr1 > driver_constants.QPMIN - and ql1 > driver_constants.QCMIN + qr1 > constants.QPMIN + and ql1 > constants.QCMIN and qsat < q_minus ): sink = half_dt * denfac * cracw * exp(0.95 * log(qr1 * den1)) @@ -269,7 +269,7 @@ def warm_rain( # if it falls anywhere in the column, the entire column becomes true # initalized to 0 (false), potentially changed to 1 (true) with computation(FORWARD), interval(...): - if qr1 > driver_constants.QPMIN: + if qr1 > constants.QPMIN: precip_fall = 1 # end reference Fortran: gfdl_cloud_microphys.F90: subroutine check_column @@ -282,7 +282,7 @@ def warm_rain( with computation(PARALLEL), interval(...): # Use In-Cloud condensates if do_qa == False: # noqa - qadum = max(qa1, driver_constants.QCMIN) + qadum = max(qa1, constants.QCMIN) else: qadum = 1.0 ql1 = ql1 / qadum @@ -292,7 +292,7 @@ def warm_rain( min(1.0, eis / 15.0) ** 2 ) # Estimated inversion strength determine stable regime fac_rc = ( - driver_constants.RC * (rthreshs * fac_rc + rthreshu * (1.0 - fac_rc)) ** 3 + constants.RC * (rthreshs * fac_rc + rthreshu * (1.0 - fac_rc)) ** 3 ) with computation(PARALLEL), interval(...): @@ -301,13 +301,13 @@ def warm_rain( # no subgrid varaibility # ----------------------------------------------------------------------- if qadum > onemsig: - if t1 > driver_constants.T_WFR: + if t1 > constants.T_WFR: qc = fac_rc * ccn / den1 dq = ql1 - qc if dq > 0.0: sink = min( dq, - dts * c_praut * den1 * exp(driver_constants.SO3 * log(ql1)), + dts * c_praut * den1 * exp(constants.SO3 * log(ql1)), ) sink = min(ql0_max, min(ql1, max(0.0, sink))) ql1 = ql1 - sink @@ -358,15 +358,15 @@ def warm_rain( with computation(PARALLEL), interval(...): if irain_f == 0: if z_slope_liq == True: # noqa - dl = max(dl, max(driver_constants.QVMIN, rh_limited * ql1)) + dl = max(dl, max(constants.QVMIN, rh_limited * ql1)) if z_slope_liq == False: # noqa - dl = max(driver_constants.QVMIN, rh_limited * ql1) + dl = max(constants.QVMIN, rh_limited * ql1) # end reference Fortran: gfdl_cloud_microphys.F90: subroutine linear_prof if qadum > onemsig: - if t1 > driver_constants.T_WFR + driver_constants.DT_FR: - dl = min(max(driver_constants.QCMIN, dl), 0.5 * ql1) + if t1 > constants.T_WFR + constants.DT_FR: + dl = min(max(constants.QCMIN, dl), 0.5 * ql1) # -------------------------------------------------------------------- # as in klein's gfdl am2 stratiform scheme (with subgrid variations) # -------------------------------------------------------------------- @@ -386,7 +386,7 @@ def warm_rain( * dts * c_praut * den1 - * exp(driver_constants.SO3 * log(ql1)) + * exp(constants.SO3 * log(ql1)) ) sink = min(ql0_max, min(ql1, max(0.0, sink))) ql1 = ql1 - sink @@ -401,24 +401,24 @@ def warm_rain( # ----------------------------------------------------------------------- if precip_fall == 0: - vtr = driver_constants.VF_MIN + vtr = constants.VF_MIN elif const_vr == True: # noqa vtr = vr_fac # ifs_2016: 4.0 else: qden = qr1 * den1 - if qr1 < driver_constants.THR: - vtr = driver_constants.VR_MIN + if qr1 < constants.THR: + vtr = constants.VR_MIN else: vtr = ( vr_fac - * driver_constants.VCONR - * sqrt(min(10.0, driver_constants.SFCRHO / den1)) - * exp(0.2 * log(qden / driver_constants.NORMR)) + * constants.VCONR + * sqrt(min(10.0, constants.SFCRHO / den1)) + * exp(0.2 * log(qden / constants.NORMR)) ) - vtr = min(vr_max, max(driver_constants.VR_MIN, vtr)) + vtr = min(vr_max, max(constants.VR_MIN, vtr)) with computation(FORWARD), interval(-1, None): - ze[0, 0, 1] = driver_constants.ZS + ze[0, 0, 1] = constants.ZS with computation(BACKWARD), interval(...): ze = ze[0, 0, 1] - dz1 # dz < 0 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/interface/cffi_lib/interface.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/interface/cffi_lib/interface.py index 24bc438e4..0cdbb0286 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/interface/cffi_lib/interface.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/interface/cffi_lib/interface.py @@ -1,3 +1,5 @@ +from distutils.sysconfig import get_config_var + import cffi from mpi4py import MPI diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/interface/wrapper.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/interface/wrapper.py index b9763e1b0..17d8cbe7a 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/interface/wrapper.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/interface/wrapper.py @@ -34,7 +34,7 @@ from ndsl.optional_imports import cupy as cp from pyMoist.aer_activation import AerActivation from pyMoist.GFDL_1M.GFDL_1M import GFDL_1M -from pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver import GFDL_1M_driver +from pyMoist.GFDL_1M.GFDL_1M_driver.driver import driver from pyMoist.interface.flags import gfdl_1m_flags, moist_flags import numpy as np @@ -199,7 +199,7 @@ def __init__( # JIT system for the component of Moist self._aer_activation: Optional[AerActivation] = None self._GFDL_1M_evap: Optional[GFDL_1M] = None - self._GFDL_1M_driver: Optional[GFDL_1M_driver] = None + self._GFDL_1M_driver: Optional[driver] = None # Initalize flags later self.gfdl_1m_flags: Optional[gfdl_1m_flags] = None @@ -213,7 +213,7 @@ def GFDL_1M_driver(self) -> Callable: MPI.COMM_WORLD, self.stencil_config.dace_config, ): - self._GFDL_1M_driver = GFDL_1M_driver( + self._GFDL_1M_driver = driver( self.stencil_factory, self.quantity_factory, self.gfdl_1m_flags.phys_hydrostatic, diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_GFDL_1M_driver.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_GFDL_1M_driver.py index 2c1696d78..d639a551d 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_GFDL_1M_driver.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_GFDL_1M_driver.py @@ -2,7 +2,8 @@ from ndsl.constants import X_DIM, Y_DIM, Z_DIM from ndsl.dsl.typing import Float, Int from ndsl.stencils.testing.translate import TranslateFortranData2Py -from pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver import GFDL_1M_driver +from pyMoist.GFDL_1M.GFDL_1M_driver.driver import driver +from pyMoist.GFDL_1M.GFDL_1M_driver.config import config class TranslateGFDL_1M_driver(TranslateFortranData2Py): @@ -316,13 +317,11 @@ def compute(self, inputs): irain_f = Float(inputs["irain_f"]) mp_print = Float(inputs["mp_print"]) - stencil = GFDL_1M_driver( - self.stencil_factory, - self.quantity_factory, + # Create config class, to be replaced by a proper feature in the pyMoist integration + namelist = config( LPHYS_HYDROSTATIC, LHYDROSTATIC, DT_MOIST, - # Namelist options mp_time, t_min, t_sub, @@ -405,7 +404,14 @@ def compute(self, inputs): mp_print, ) + stencil = driver( + self.stencil_factory, + self.quantity_factory, + namelist, + ) + stencil( + namelist, RAD_QV, RAD_QL, RAD_QR, @@ -452,12 +458,12 @@ def compute(self, inputs): "DQIDTmic": DQIDTmic.view[:], "DQSDTmic": DQSDTmic.view[:], "DQGDTmic": DQGDTmic.view[:], - "PRCP_RAIN": stencil.rain.view[:], - "PRCP_SNOW": stencil.snow.view[:], - "PRCP_ICE": stencil.ice.view[:], - "PRCP_GRAUPEL": stencil.graupel.view[:], - "PFL_LS": stencil.m2_rain.view[:], - "PFI_LS": stencil.m2_sol.view[:], - "REV_LS": stencil.revap.view[:], - "RSU_LS": stencil.isubl.view[:], + "PRCP_RAIN": stencil.outputs.rain.view[:], + "PRCP_SNOW": stencil.outputs.snow.view[:], + "PRCP_ICE": stencil.outputs.ice.view[:], + "PRCP_GRAUPEL": stencil.outputs.graupel.view[:], + "PFL_LS": stencil.outputs.m2_rain.view[:], + "PFI_LS": stencil.outputs.m2_sol.view[:], + "REV_LS": stencil.outputs.revap.view[:], + "RSU_LS": stencil.outputs.isubl.view[:], } diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_GFDL_1M_driver_preloop.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_GFDL_1M_driver_preloop.py index 95f75c3fa..0deb4e2e3 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_GFDL_1M_driver_preloop.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_GFDL_1M_driver_preloop.py @@ -1,8 +1,8 @@ -import pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver_constants as driver_constants +import pyMoist.GFDL_1M.GFDL_1M_driver.constants as constants from ndsl import Namelist, Quantity, StencilFactory, orchestrate from ndsl.constants import X_DIM, Y_DIM, Z_DIM, Z_INTERFACE_DIM from ndsl.stencils.testing.translate import TranslateFortranData2Py -from pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver_core import ( +from pyMoist.GFDL_1M.GFDL_1M_driver.connections import ( fix_negative_values, init_temporaries, ) @@ -82,11 +82,11 @@ def __init__(self, grid, namelist: Namelist, stencil_factory: StencilFactory): # set externals manually to the need for avoid another workaround c_paut = 1 - cpaut = c_paut * 0.104 * driver_constants.GRAV / 1.717e-5 - c_air = driver_constants.CP_AIR - c_vap = driver_constants.CP_VAP - d0_vap = c_vap - driver_constants.C_LIQ - lv00 = driver_constants.HLV0 - d0_vap * driver_constants.T_ICE + cpaut = c_paut * 0.104 * constants.GRAV / 1.717e-5 + c_air = constants.CP_AIR + c_vap = constants.CP_VAP + d0_vap = c_vap - constants.C_LIQ + lv00 = constants.HLV0 - d0_vap * constants.T_ICE # initalize stencils orchestrate(obj=self, config=stencil_factory.config.dace_config) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_GFDL_driver_tables.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_GFDL_driver_tables.py index 95fc4d52b..dbe17a436 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_GFDL_driver_tables.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_GFDL_driver_tables.py @@ -1,6 +1,6 @@ from ndsl import Namelist, StencilFactory from ndsl.stencils.testing.translate import TranslateFortranData2Py -from pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver_tables import get_tables +from pyMoist.GFDL_1M.GFDL_1M_driver.sat_tables import get_tables class TranslateGFDL_driver_tables(TranslateFortranData2Py): diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_fall_speed.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_fall_speed.py index fd4cfcfd7..d94ace8cc 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_fall_speed.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_fall_speed.py @@ -4,7 +4,7 @@ from ndsl.constants import X_DIM, Y_DIM, Z_DIM from ndsl.dsl.typing import Float, FloatField, FloatFieldIJ from ndsl.stencils.testing.translate import TranslateFortranData2Py -from pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver_core import fall_speed_core +from pyMoist.GFDL_1M.GFDL_1M_driver.connections import fall_speed_core class Translatefall_speed(TranslateFortranData2Py): diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_icloud.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_icloud.py index 92737b528..f9a71182c 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_icloud.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_icloud.py @@ -1,11 +1,11 @@ import numpy as np -import pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver_constants as driver_constants +import pyMoist.GFDL_1M.GFDL_1M_driver.constants as constants from ndsl import Namelist, Quantity, StencilFactory, orchestrate from ndsl.constants import X_DIM, Y_DIM, Z_DIM from ndsl.dsl.typing import Float from ndsl.stencils.testing.translate import TranslateFortranData2Py -from pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver_tables import get_tables +from pyMoist.GFDL_1M.GFDL_1M_driver.sat_tables import get_tables from pyMoist.GFDL_1M.GFDL_1M_driver.icloud import icloud @@ -304,15 +304,15 @@ def compute(self, inputs): True # not about to try to serialize this, it should always be true ) if phys_hydrostatic or hydrostatic: - c_air = driver_constants.CP_AIR - c_vap = driver_constants.CP_VAP + c_air = constants.CP_AIR + c_vap = constants.CP_VAP p_nonhydro = False else: - c_air = driver_constants.CV_AIR - c_vap = driver_constants.CV_VAP + c_air = constants.CV_AIR + c_vap = constants.CV_VAP p_nonhydro = True - d0_vap = c_vap - driver_constants.C_LIQ - lv00 = driver_constants.HLV0 - d0_vap * driver_constants.T_ICE + d0_vap = c_vap - constants.C_LIQ + lv00 = constants.HLV0 - d0_vap * constants.T_ICE if hydrostatic: do_sedi_w = False @@ -321,21 +321,21 @@ def compute(self, inputs): # define latent heat coefficient used in wet bulb and bigg mechanism # ----------------------------------------------------------------------- - latv = driver_constants.HLV - lati = driver_constants.HLF + latv = constants.HLV + lati = constants.HLF lats = latv + lati lat2 = lats * lats - lcp = latv / driver_constants.CP_AIR - icp = lati / driver_constants.CP_AIR - tcp = (latv + lati) / driver_constants.CP_AIR + lcp = latv / constants.CP_AIR + icp = lati / constants.CP_AIR + tcp = (latv + lati) / constants.CP_AIR # ----------------------------------------------------------------------- # calculate cloud condensation nuclei (ccn) # the following is based on klein eq. 15 # ----------------------------------------------------------------------- - cpaut = c_paut * 0.104 * driver_constants.GRAV / 1.717e-5 + cpaut = c_paut * 0.104 * constants.GRAV / 1.717e-5 # ----------------------------------------------------------------------- # define conversion scalar / factor for icloud @@ -357,37 +357,37 @@ def compute(self, inputs): # ----------------------------------------------------------------------- cgacs = ( - driver_constants.PISQ - * driver_constants.RNZG - * driver_constants.RNZS - * driver_constants.RHOS + constants.PISQ + * constants.RNZG + * constants.RNZS + * constants.RHOS ) cgacs = cgacs * c_pgacs csacw = ( - driver_constants.PIE - * driver_constants.RNZS + constants.PIE + * constants.RNZS * clin - * driver_constants.GAM325 - / (4.0 * driver_constants.ACT[0] ** 0.8125) + * constants.GAM325 + / (4.0 * constants.ACT[0] ** 0.8125) ) # decreasing csacw to reduce cloud water --- > snow craci = ( - driver_constants.PIE - * driver_constants.RNZR + constants.PIE + * constants.RNZR * alin - * driver_constants.GAM380 - / (4.0 * driver_constants.ACT[1] ** 0.95) + * constants.GAM380 + / (4.0 * constants.ACT[1] ** 0.95) ) csaci = csacw * c_psaci cgacw = ( - driver_constants.PIE - * driver_constants.RNZG - * driver_constants.GAM350 - * driver_constants.GCON - / (4.0 * driver_constants.ACT[5] ** 0.875) + constants.PIE + * constants.RNZG + * constants.GAM350 + * constants.GCON + / (4.0 * constants.ACT[5] ** 0.875) ) cgaci = cgacw * c_pgaci @@ -401,65 +401,65 @@ def compute(self, inputs): cssub[0] = ( 2.0 - * driver_constants.PIE - * driver_constants.VDIFU - * driver_constants.TCOND - * driver_constants.RVGAS - * driver_constants.RNZS + * constants.PIE + * constants.VDIFU + * constants.TCOND + * constants.RVGAS + * constants.RNZS ) cgsub[0] = ( 2.0 - * driver_constants.PIE - * driver_constants.VDIFU - * driver_constants.TCOND - * driver_constants.RVGAS - * driver_constants.RNZG + * constants.PIE + * constants.VDIFU + * constants.TCOND + * constants.RVGAS + * constants.RNZG ) crevp[0] = ( 2.0 - * driver_constants.PIE - * driver_constants.VDIFU - * driver_constants.TCOND - * driver_constants.RVGAS - * driver_constants.RNZR + * constants.PIE + * constants.VDIFU + * constants.TCOND + * constants.RVGAS + * constants.RNZR ) - cssub[1] = 0.78 / np.sqrt(driver_constants.ACT[0]) - cgsub[1] = 0.78 / np.sqrt(driver_constants.ACT[5]) - crevp[1] = 0.78 / np.sqrt(driver_constants.ACT[1]) + cssub[1] = 0.78 / np.sqrt(constants.ACT[0]) + cgsub[1] = 0.78 / np.sqrt(constants.ACT[5]) + crevp[1] = 0.78 / np.sqrt(constants.ACT[1]) cssub[2] = ( 0.31 - * driver_constants.SCM3 - * driver_constants.GAM263 - * np.sqrt(clin / driver_constants.VISK) - / driver_constants.ACT[0] ** 0.65625 + * constants.SCM3 + * constants.GAM263 + * np.sqrt(clin / constants.VISK) + / constants.ACT[0] ** 0.65625 ) cgsub[2] = ( 0.31 - * driver_constants.SCM3 - * driver_constants.GAM275 - * np.sqrt(driver_constants.GCON / driver_constants.VISK) - / driver_constants.ACT[5] ** 0.6875 + * constants.SCM3 + * constants.GAM275 + * np.sqrt(constants.GCON / constants.VISK) + / constants.ACT[5] ** 0.6875 ) crevp[2] = ( 0.31 - * driver_constants.SCM3 - * driver_constants.GAM209 - * np.sqrt(alin / driver_constants.VISK) - / driver_constants.ACT[1] ** 0.725 + * constants.SCM3 + * constants.GAM209 + * np.sqrt(alin / constants.VISK) + / constants.ACT[1] ** 0.725 ) - cssub[3] = driver_constants.TCOND * driver_constants.RVGAS - cssub[4] = driver_constants.HLTS**2 * driver_constants.VDIFU + cssub[3] = constants.TCOND * constants.RVGAS + cssub[4] = constants.HLTS**2 * constants.VDIFU cgsub[3] = cssub[3] crevp[3] = cssub[3] cgsub[4] = cssub[4] - crevp[4] = driver_constants.HLTC**2 * driver_constants.VDIFU + crevp[4] = constants.HLTC**2 * constants.VDIFU cgfr_0 = ( 20.0e2 - * driver_constants.PISQ - * driver_constants.RNZR - * driver_constants.RHOR - / driver_constants.ACT[1] ** 1.75 + * constants.PISQ + * constants.RNZR + * constants.RHOR + / constants.ACT[1] ** 1.75 ) cgfr_1 = 0.66 @@ -486,22 +486,22 @@ def compute(self, inputs): csmlt = np.zeros(5) csmlt[0] = ( 2.0 - * driver_constants.PIE - * driver_constants.TCOND - * driver_constants.RNZS - / driver_constants.HLTF + * constants.PIE + * constants.TCOND + * constants.RNZS + / constants.HLTF ) csmlt[1] = ( 2.0 - * driver_constants.PIE - * driver_constants.VDIFU - * driver_constants.RNZS - * driver_constants.HLTC - / driver_constants.HLTF + * constants.PIE + * constants.VDIFU + * constants.RNZS + * constants.HLTC + / constants.HLTF ) csmlt[2] = cssub[1] csmlt[3] = cssub[2] - csmlt[4] = driver_constants.CH2O / driver_constants.HLTF + csmlt[4] = constants.CH2O / constants.HLTF csmlt_0 = csmlt[0] csmlt_1 = csmlt[1] @@ -514,22 +514,22 @@ def compute(self, inputs): cgmlt = np.zeros(5) cgmlt[0] = ( 2.0 - * driver_constants.PIE - * driver_constants.TCOND - * driver_constants.RNZG - / driver_constants.HLTF + * constants.PIE + * constants.TCOND + * constants.RNZG + / constants.HLTF ) cgmlt[1] = ( 2.0 - * driver_constants.PIE - * driver_constants.VDIFU - * driver_constants.RNZG - * driver_constants.HLTC - / driver_constants.HLTF + * constants.PIE + * constants.VDIFU + * constants.RNZG + * constants.HLTC + / constants.HLTF ) cgmlt[2] = cgsub[1] cgmlt[3] = cgsub[2] - cgmlt[4] = driver_constants.CH2O / driver_constants.HLTF + cgmlt[4] = constants.CH2O / constants.HLTF cgmlt_0 = cgmlt[0] cgmlt_1 = cgmlt[1] @@ -538,7 +538,7 @@ def compute(self, inputs): cgmlt_4 = cgmlt[4] es0 = 6.107799961e2 # ~6.1 mb - ces0 = driver_constants.EPS * es0 + ces0 = constants.EPS * es0 # make temporaries self.TESTVAR_1 = self.quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_terminal_fall.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_terminal_fall.py deleted file mode 100644 index 7d35007a0..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_terminal_fall.py +++ /dev/null @@ -1,231 +0,0 @@ -from ndsl import Namelist, Quantity, StencilFactory, orchestrate -from ndsl.constants import X_DIM, Y_DIM, Z_DIM, Z_INTERFACE_DIM -from ndsl.dsl.typing import Float, Int -from ndsl.stencils.testing.translate import TranslateFortranData2Py -from pyMoist.GFDL_1M.GFDL_1M_driver.terminal_fall import terminal_fall - - -class Translateterminal_fall(TranslateFortranData2Py): - def __init__(self, grid, namelist: Namelist, stencil_factory: StencilFactory): - super().__init__(grid, stencil_factory) - self.stencil_factory = stencil_factory - self.quantity_factory = grid.quantity_factory - self._grid = grid - - # FloatField Inputs - self.in_vars["data_vars"] = { - "t_terminal_fall": {}, - "qv_terminal_fall": {}, - "ql_terminal_fall": {}, - "qr_terminal_fall": {}, - "qg_terminal_fall": {}, - "qs_terminal_fall": {}, - "qi_terminal_fall": {}, - "dz1_terminal_fall": {}, - "dp1_terminal_fall": {}, - "den1_terminal_fall": {}, - "vtg_terminal_fall": {}, - "vts_terminal_fall": {}, - "vti_terminal_fall": {}, - "m1_sol_terminal_fall": {}, - "w1_terminal_fall": {}, - "dts_terminal_fall": {}, - "tau_imlt_terminal_fall": {}, - "ql_mlt_terminal_fall": {}, - "c_air_terminal_fall": {}, - "c_vap_terminal_fall": {}, - "d0_vap_terminal_fall": {}, - "lv00_terminal_fall": {}, - "vi_fac_terminal_fall": {}, - "do_sedi_w_terminal_fall": {}, - "use_ppm_terminal_fall": {}, - "tau_smlt_terminal_fall": {}, - "tau_g2r_terminal_fall": {}, - } - - # FloatField Outputs - self.out_vars = { - "t_terminal_fall": self.grid.compute_dict(), - "qv_terminal_fall": self.grid.compute_dict(), - "ql_terminal_fall": self.grid.compute_dict(), - "qr_terminal_fall": self.grid.compute_dict(), - "qg_terminal_fall": self.grid.compute_dict(), - "qs_terminal_fall": self.grid.compute_dict(), - "qi_terminal_fall": self.grid.compute_dict(), - "dz1_terminal_fall": self.grid.compute_dict(), - "dp1_terminal_fall": self.grid.compute_dict(), - "den1_terminal_fall": self.grid.compute_dict(), - "m1_sol_terminal_fall": self.grid.compute_dict(), - "r1_terminal_fall": self.grid.compute_dict(), - "g1_terminal_fall": self.grid.compute_dict(), - "s1_terminal_fall": self.grid.compute_dict(), - "i1_terminal_fall": self.grid.compute_dict(), - } - - def make_ij_field(self, data) -> Quantity: - qty = self.quantity_factory.empty( - [X_DIM, Y_DIM], - "n/a", - ) - qty.view[:, :] = qty.np.asarray(data[:, :]) - return qty - - def make_ijk_field(self, data) -> Quantity: - qty = self.quantity_factory.empty( - [X_DIM, Y_DIM, Z_DIM], - "n/a", - ) - qty.view[:, :, :] = qty.np.asarray(data[:, :, :]) - return qty - - def compute(self, inputs): - # FloatField Variables - t_terminal_fall = self.make_ijk_field(inputs["t_terminal_fall"]) - t_original = t_terminal_fall - qv_terminal_fall = self.make_ijk_field(inputs["qv_terminal_fall"]) - qv_original = qv_terminal_fall - ql_terminal_fall = self.make_ijk_field(inputs["ql_terminal_fall"]) - ql_original = ql_terminal_fall - qr_terminal_fall = self.make_ijk_field(inputs["qr_terminal_fall"]) - qr_original = qr_terminal_fall - qg_terminal_fall = self.make_ijk_field(inputs["qg_terminal_fall"]) - qg_original = qg_terminal_fall - qs_terminal_fall = self.make_ijk_field(inputs["qs_terminal_fall"]) - qs_original = qs_terminal_fall - qi_terminal_fall = self.make_ijk_field(inputs["qi_terminal_fall"]) - qi_original = qi_terminal_fall - dz1_terminal_fall = self.make_ijk_field(inputs["dz1_terminal_fall"]) - dz1_original = dz1_terminal_fall - dp1_terminal_fall = self.make_ijk_field(inputs["dp1_terminal_fall"]) - dp1_original = dp1_terminal_fall - den1_terminal_fall = self.make_ijk_field(inputs["den1_terminal_fall"]) - den1_original = den1_terminal_fall - vtg_terminal_fall = self.make_ijk_field(inputs["vtg_terminal_fall"]) - vtg_original = vtg_terminal_fall - vts_terminal_fall = self.make_ijk_field(inputs["vts_terminal_fall"]) - vts_original = vts_terminal_fall - vti_terminal_fall = self.make_ijk_field(inputs["vti_terminal_fall"]) - vti_original = vti_terminal_fall - m1_sol_terminal_fall = self.make_ijk_field(inputs["m1_sol_terminal_fall"]) - m1_sol_original = m1_sol_terminal_fall - w1_terminal_fall = self.make_ijk_field(inputs["w1_terminal_fall"]) - w1_original = w1_terminal_fall - - # Float Variables - dts = Float(inputs["dts_terminal_fall"][0]) - tau_imlt = Float(inputs["tau_imlt_terminal_fall"][0]) - ql_mlt = Float(inputs["ql_mlt_terminal_fall"][0]) - c_air = Float(inputs["c_air_terminal_fall"][0]) - c_vap = Float(inputs["c_vap_terminal_fall"][0]) - d0_vap = Float(inputs["d0_vap_terminal_fall"][0]) - lv00 = Float(inputs["lv00_terminal_fall"][0]) - vi_fac = Float(inputs["vi_fac_terminal_fall"][0]) - do_sedi_w = Float(inputs["do_sedi_w_terminal_fall"][0]) - use_ppm = Float(inputs["use_ppm_terminal_fall"][0]) - tau_smlt = Float(inputs["tau_smlt_terminal_fall"][0]) - tau_g2r = Float(inputs["tau_g2r_terminal_fall"][0]) - - # Revert to boolean - if do_sedi_w == 1: - do_sedi_w = True - else: - do_sedi_w = False - if use_ppm == 1: - use_ppm = True - else: - use_ppm = False - - # make masks - self.is_frozen = self.quantity_factory.zeros( - [X_DIM, Y_DIM, Z_DIM], "n/a", dtype=bool - ) - self.precip_fall = self.quantity_factory.zeros([X_DIM, Y_DIM], "n/a") - self.melting_mask_1 = self.quantity_factory.zeros( - [X_DIM, Y_DIM, Z_DIM], "n/a", dtype=bool - ) - self.melting_mask_2 = self.quantity_factory.zeros( - [X_DIM, Y_DIM, Z_DIM], "n/a", dtype=bool - ) - self.current_k_level = self.quantity_factory.zeros( - [X_DIM, Y_DIM, Z_DIM], "n/a", dtype=Int - ) - for k in range(self.current_k_level.view[:].shape[2]): - self.current_k_level.view[:, :, k] = k - - # make temporaries - self.m1 = self.quantity_factory.zeros([X_DIM, Y_DIM, Z_DIM], "n/a") - self.ze = self.quantity_factory.zeros([X_DIM, Y_DIM, Z_INTERFACE_DIM], "n/a") - self.zt = self.quantity_factory.zeros([X_DIM, Y_DIM, Z_INTERFACE_DIM], "n/a") - - # make outputs - self.rain = self.quantity_factory.zeros([X_DIM, Y_DIM], "n/a") - self.snow = self.quantity_factory.zeros([X_DIM, Y_DIM], "n/a") - self.ice = self.quantity_factory.ones([X_DIM, Y_DIM], "n/a") - self.graupel = self.quantity_factory.zeros([X_DIM, Y_DIM], "n/a") - - orchestrate(obj=self, config=self.stencil_factory.config.dace_config) - self._stencil = self.stencil_factory.from_dims_halo( - func=terminal_fall, - compute_dims=[X_DIM, Y_DIM, Z_DIM], - externals={ - "dts": dts, - "tau_imlt": tau_imlt, - "ql_mlt": ql_mlt, - "c_air": c_air, - "c_vap": c_vap, - "d0_vap": d0_vap, - "lv00": lv00, - "vi_fac": vi_fac, - "do_sedi_w": do_sedi_w, - "use_ppm": use_ppm, - "tau_smlt": tau_smlt, - "tau_g2r": tau_g2r, - }, - ) - - self._stencil( - t_terminal_fall, - qv_terminal_fall, - ql_terminal_fall, - qr_terminal_fall, - qg_terminal_fall, - qs_terminal_fall, - qi_terminal_fall, - dz1_terminal_fall, - dp1_terminal_fall, - den1_terminal_fall, - vtg_terminal_fall, - vts_terminal_fall, - vti_terminal_fall, - self.rain, - self.graupel, - self.snow, - self.ice, - m1_sol_terminal_fall, - w1_terminal_fall, - self.ze, - self.zt, - self.is_frozen, - self.precip_fall, - self.melting_mask_1, - self.melting_mask_2, - self.current_k_level, - ) - - return { - "t_terminal_fall": t_terminal_fall.view[:], - "qv_terminal_fall": qv_terminal_fall.view[:], - "ql_terminal_fall": ql_terminal_fall.view[:], - "qr_terminal_fall": qr_terminal_fall.view[:], - "qg_terminal_fall": qg_terminal_fall.view[:], - "qs_terminal_fall": qs_terminal_fall.view[:], - "qi_terminal_fall": qi_terminal_fall.view[:], - "dz1_terminal_fall": dz1_terminal_fall.view[:], - "dp1_terminal_fall": dp1_terminal_fall.view[:], - "den1_terminal_fall": den1_terminal_fall.view[:], - "m1_sol_terminal_fall": m1_sol_terminal_fall.view[:], - "r1_terminal_fall": self.rain.view[:], - "g1_terminal_fall": self.graupel.view[:], - "s1_terminal_fall": self.snow.view[:], - "i1_terminal_fall": self.ice.view[:], - } diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_warm_rain.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_warm_rain.py index 5e29ea331..c2ac0ed6b 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_warm_rain.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/GFDL_1M/GFDL_1M_driver/translate_warm_rain.py @@ -2,7 +2,7 @@ from ndsl.constants import X_DIM, Y_DIM, Z_DIM, Z_INTERFACE_DIM from ndsl.dsl.typing import Float, Int from ndsl.stencils.testing.translate import TranslateFortranData2Py -from pyMoist.GFDL_1M.GFDL_1M_driver.GFDL_1M_driver_tables import get_tables +from pyMoist.GFDL_1M.GFDL_1M_driver.sat_tables import get_tables from pyMoist.GFDL_1M.GFDL_1M_driver.warm_rain import warm_rain diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/__init__.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/__init__.py index eee481c92..cead35d0e 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/__init__.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/__init__.py @@ -7,7 +7,6 @@ TranslateGFDL_driver_tables, ) from .GFDL_1M.GFDL_1M_driver.translate_icloud import Translateicloud -from .GFDL_1M.GFDL_1M_driver.translate_terminal_fall import Translateterminal_fall from .GFDL_1M.GFDL_1M_driver.translate_warm_rain import Translatewarm_rain from .GFDL_1M.translate_GFDL_1M import TranslateGFDL_1M from .translate_aer_activation import TranslateAerActivation diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/scripts/run_tests.sh b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/scripts/run_tests.sh index a9ef7e953..5e50ca437 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/scripts/run_tests.sh +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/scripts/run_tests.sh @@ -15,6 +15,6 @@ python -m pytest -s \ --backend=gt:cpu_kfirst \ --grid=default \ --multimodal_metric \ - --which_modules=icloud \ + --which_modules=terminal_fall \ --which_rank=0 \ ..