Lowe et al. 2019 Fig S2 (+smoke tests); new optical formulae and prod…

…ucts (optical depth & albedo); LWP product; activation-filtered r_eff product; cleanups (#1256)

Co-authored-by: Sanky <76765634+bhiogade@users.noreply.github.com>
Co-authored-by: Sanket Bhiogade <sbhiogade@agh.edu.pl>
Co-authored-by: Sylwester Arabas <sylwester.arabas@agh.edu.pl>
Co-authored-by: claresinger <clareees@gmail.com>

PySDM/attributes/impl/mapper.py

@@ -51,17 +51,15 @@
     "volume": lambda _, __: Volume,
     "dry volume organic": lambda dynamics, formulae: (
         make_dummy_attribute_factory("dry volume organic")
-        if "Condensation" in dynamics
-        and formulae.surface_tension.__name__ == Constant.__name__
+        if formulae.surface_tension.__name__ == Constant.__name__
         else DryVolumeOrganic
     ),
     "dry volume": lambda dynamics, formulae: (
         DryVolumeDynamic if "AqueousChemistry" in dynamics else DryVolume
     ),
     "dry volume organic fraction": lambda dynamics, formulae: (
         make_dummy_attribute_factory("dry volume organic fraction")
-        if "Condensation" in dynamics
-        and formulae.surface_tension.__name__ == Constant.__name__
+        if formulae.surface_tension.__name__ == Constant.__name__
         else OrganicFraction
     ),
     "kappa times dry volume": lambda _, __: KappaTimesDryVolume,
@@ -124,6 +122,6 @@
 def get_class(name, dynamics, formulae):
     if name not in attributes:
         raise ValueError(
-            f"Unknown attribute name: {name}; valid names: {', '.join(attributes)}"
+            f"Unknown attribute name: {name}; valid names: {', '.join(sorted(attributes))}"
         )
     return attributes[name](dynamics, formulae)

PySDM/formulae.py

@@ -48,12 +48,16 @@ def __init__(  # pylint: disable=too-many-locals
         isotope_equilibrium_fractionation_factors: str = "Null",
         isotope_meteoric_water_line_excess: str = "Null",
         isotope_ratio_evolution: str = "Null",
+        optical_albedo: str = "Null",
+        optical_depth: str = "Null",
         particle_shape_and_density: str = "LiquidSpheres",
         terminal_velocity: str = "GunnKinzer1949",
         handle_all_breakups: bool = False,
     ):
         # initialisation of the fields below is just to silence pylint and to enable code hints
         # in PyCharm and alike, all these fields are later overwritten within this ctor
+        self.optical_albedo = optical_albedo
+        self.optical_depth = optical_depth
         self.condensation_coordinate = condensation_coordinate
         self.saturation_vapour_pressure = saturation_vapour_pressure
         self.hygroscopicity = hygroscopicity

PySDM/physics/__init__.py

@@ -32,6 +32,8 @@
     isotope_meteoric_water_line_excess,
     isotope_ratio_evolution,
     latent_heat,
+    optical_albedo,
+    optical_depth,
     particle_advection,
     particle_shape_and_density,
     saturation_vapour_pressure,

PySDM/physics/constants.py

@@ -34,6 +34,7 @@ def in_unit(value, unit):
 LN_2 = np.log(2)
 ZERO_MASS = 0 * si.kg
 ZERO_VOLUME = 0 * si.m**3
+ONE = 1
 TWO = 2
 THREE = 3
 FOUR = 4

PySDM/physics/constants_defaults.py

@@ -20,6 +20,7 @@
     PPM,
     T0,
     THREE,
+    ONE,
     TWO,
     TWO_THIRDS,
     M,
@@ -317,6 +318,9 @@
 CRAIG_1961_SLOPE_COEFF = 8
 CRAIG_1961_INTERCEPT_COEFF = 10 * PER_MILLE
 
+""" [Bohren 1987](https://doi.org/10.1119/1.15109) """
+asymmetry_g = 0.85  # forward scattering from cloud droplets
+
 
 def compute_derived_values(c: dict):
     c["eps"] = c["Mv"] / c["Md"]
@@ -325,6 +329,6 @@ def compute_derived_values(c: dict):
 
     c["Rd_over_c_pd"] = c["Rd"] / c["c_pd"]
 
-    c["nu_w"] = c["Mv"] / c["rho_w"]
+    c["water_molar_volume"] = c["Mv"] / c["rho_w"]
     c["rho_STP"] = c["p_STP"] / c["Rd"] / c["T_STP"]
     c["H_u"] = c["M"] / c["p_STP"]

PySDM/physics/optical_albedo/__init__.py

@@ -0,0 +1,6 @@
+"""
+alternative formulations of cloud albedo
+"""
+
+from .bohren1987 import Bohren1987
+from .null import Null

PySDM/physics/optical_albedo/bohren1987.py

@@ -0,0 +1,14 @@
+"""
+[Bohren 1987](https://doi.org/10.1119/1.15109) Eq. 14
+"""
+
+
+class Bohren1987:  # pylint: disable=too-few-public-methods
+    def __init__(self, _):
+        pass
+
+    @staticmethod
+    def albedo(const, tau):
+        return ((const.ONE - const.asymmetry_g) * tau) / (
+            const.TWO + (const.ONE - const.asymmetry_g) * tau
+        )

PySDM/physics/optical_albedo/null.py

@@ -0,0 +1,8 @@
+"""
+null (default) class
+"""
+
+
+class Null:  # pylint: disable=too-few-public-methods
+    def __init__(self, _):
+        pass

PySDM/physics/optical_depth/__init__.py

@@ -0,0 +1,6 @@
+"""
+alternative formulations of optical depth
+"""
+
+from .stephens_1978 import Stephens1978
+from .null import Null

PySDM/physics/optical_depth/null.py

@@ -0,0 +1,8 @@
+"""
+null (default) class
+"""
+
+
+class Null:  # pylint: disable=too-few-public-methods
+    def __init__(self, _):
+        pass

PySDM/physics/optical_depth/stephens_1978.py

@@ -0,0 +1,13 @@
+"""
+[Stephens 1978](https://doi.org/10.1175/1520-0469(1978)035%3C2123:RPIEWC%3E2.0.CO;2)
+Eq. 7 for optical depth, where LWP is in g/m^2 and reff is in um.
+"""
+
+
+class Stephens1978:  # pylint: disable=too-few-public-methods
+    def __init__(self, _):
+        pass
+
+    @staticmethod
+    def tau(const, LWP, reff):
+        return (const.ONE_AND_A_HALF * LWP) / (const.rho_w * reff)

PySDM/physics/surface_tension/compressed_film_ruehl.py

@@ -56,7 +56,9 @@ def sigma(const, T, v_wet, v_dry, f_org):  # pylint: disable=too-many-locals
         else:
             # C_bulk is the concentration of the organic in the bulk phase
             # Cb_iso = C_bulk / (1-f_surf)
-            Cb_iso = (f_org * v_dry / const.RUEHL_nu_org) / (v_wet / const.nu_w)
+            Cb_iso = (f_org * v_dry / const.RUEHL_nu_org) / (
+                v_wet / const.water_molar_volume
+            )
 
             # A is the area one molecule of organic occupies at the droplet surface
             # A_iso = A*f_surf (m^2)

PySDM/physics/surface_tension/szyszkowski_langmuir.py

@@ -34,7 +34,9 @@ def sigma(const, T, v_wet, v_dry, f_org):
         else:
             # C_bulk is the concentration of the organic in the bulk phase
             # Cb_iso = C_bulk / (1-f_surf)
-            Cb_iso = (f_org * v_dry / const.RUEHL_nu_org) / (v_wet / const.nu_w)
+            Cb_iso = (f_org * v_dry / const.RUEHL_nu_org) / (
+                v_wet / const.water_molar_volume
+            )
 
             # A is the area that one molecule of organic occupies at the droplet surface
             # A_iso = A*f_surf (m^2)

PySDM/products/__init__.py

@@ -1,7 +1,5 @@
 """
-Simulation output products such as:
-`PySDM.products.size_spectral.particle_size_spectrum.ParticleSizeSpectrum`,
-...
+Simulation output products
 """
 
 from .ambient_thermodynamics import *
@@ -12,3 +10,5 @@
 from .freezing import *
 from .housekeeping import *
 from .size_spectral import *
+from .optical import *
+from .parcel import *

PySDM/products/freezing/__init__.py

@@ -12,5 +12,4 @@
     IceNucleiConcentration,
     SpecificIceNucleiConcentration,
 )
-from .ice_water_content import IceWaterContent, SpecificIceWaterContent
 from .total_unfrozen_immersed_surface_area import TotalUnfrozenImmersedSurfaceArea

PySDM/products/optical/__init__.py

@@ -0,0 +1,4 @@
+""" cloud optical properties """
+
+from .cloud_optical_depth import CloudOpticalDepth
+from .cloud_albedo import CloudAlbedo

PySDM/products/optical/cloud_albedo.py

@@ -0,0 +1,13 @@
+"""
+cloud albedo
+"""
+
+from PySDM.products.impl.product import Product
+
+
+class CloudAlbedo(Product):
+    def __init__(self, *, unit="dimensionless", name=None):
+        super().__init__(name=name, unit=unit)
+
+    def _impl(self, **kwargs):
+        return self.formulae.optical_albedo.albedo(kwargs["optical_depth"])

PySDM/products/optical/cloud_optical_depth.py

@@ -0,0 +1,16 @@
+"""
+cloud optical depth
+"""
+
+from PySDM.products.impl.product import Product
+
+
+class CloudOpticalDepth(Product):
+    def __init__(self, *, unit="dimensionless", name=None):
+        super().__init__(name=name, unit=unit)
+
+    def _impl(self, **kwargs):
+        return self.formulae.optical_depth.tau(
+            kwargs["liquid_water_path"],
+            kwargs["effective_radius"],
+        )

PySDM/products/parcel/__init__.py

@@ -0,0 +1,3 @@
+""" products specific to the parcel environment """
+
+from .cloud_water_path import ParcelLiquidWaterPath

PySDM/products/parcel/cloud_water_path.py

@@ -0,0 +1,55 @@
+"""
+cloud water path integrated over parcel displacement taking into account changes
+in parcel volume along the way
+"""
+
+from PySDM.environments.parcel import Parcel
+
+from PySDM.products.impl.activation_filtered_product import _ActivationFilteredProduct
+from PySDM.products.impl.moment_product import MomentProduct
+
+
+class ParcelLiquidWaterPath(MomentProduct, _ActivationFilteredProduct):
+    def __init__(
+        self,
+        count_unactivated: bool,
+        count_activated: bool,
+        name=None,
+        unit="kg/m^2",
+    ):
+        MomentProduct.__init__(self, unit=unit, name=name)
+        _ActivationFilteredProduct.__init__(
+            self, count_activated=count_activated, count_unactivated=count_unactivated
+        )
+        self.previous = {"z": 0.0, "cwc": 0.0}
+        self.cwp = 0.0
+
+    def register(self, builder):
+        if not isinstance(builder.particulator.environment, Parcel):
+            raise NotImplementedError()
+        _ActivationFilteredProduct.register(self, builder)
+        MomentProduct.register(self, builder)
+        self.particulator.observers.append(self)
+
+    def notify(self):
+        _ActivationFilteredProduct.impl(self, attr="water mass", rank=1)
+        avg_mass = self.buffer.copy()
+
+        _ActivationFilteredProduct.impl(self, attr="water mass", rank=0)
+        tot_numb = self.buffer.copy()
+
+        self._download_to_buffer(self.particulator.environment["z"])
+        current_z = self.buffer.copy()
+
+        cwc = avg_mass * tot_numb / self.particulator.mesh.dv
+        dz = current_z - self.previous["z"]
+        cwc_mean = (cwc + self.previous["cwc"]) / 2
+
+        if self.previous["cwc"] > 0:
+            self.cwp += cwc_mean * dz
+
+        self.previous["z"] = current_z
+        self.previous["cwc"] = cwc
+
+    def _impl(self, **kwargs):
+        return self.cwp

PySDM/products/size_spectral/__init__.py

@@ -7,6 +7,7 @@
     ZerothMoment,
 )
 from .effective_radius import EffectiveRadius
+from .effective_radius_activated import ActivatedEffectiveRadius
 from .mean_radius import MeanRadius
 from .mean_radius_activated import ActivatedMeanRadius
 from .mean_volume_radius import MeanVolumeRadius
@@ -34,3 +35,11 @@
 from .total_particle_concentration import TotalParticleConcentration
 from .total_particle_specific_concentration import TotalParticleSpecificConcentration
 from .water_mixing_ratio import WaterMixingRatio
+from .cloud_water_content import (
+    CloudWaterContent,
+    SpecificCloudWaterContent,
+    LiquidWaterContent,
+    SpecificLiquidWaterContent,
+    IceWaterContent,
+    SpecificIceWaterContent,
+)