isotopes: Rozanski & Sonntag 1982 example (#1209)

Author: slayoo

PySDM/dynamics/isotopic_fractionation.py

@@ -35,4 +35,4 @@ def register(self, builder):
             builder.request_attribute(f"moles_{isotope}")
 
     def __call__(self):
-        self.particulator.isotopic_fractionation()
+        self.particulator.isotopic_fractionation(self.isotopes)

PySDM/environments/impl/moist.py

@@ -12,6 +12,9 @@ def __init__(self, dt, mesh, variables, mixed_phase=False):
         variables += ["water_vapour_mixing_ratio", "thd", "T", "p", "RH"]
         if mixed_phase:
             variables += ["a_w_ice"]
+        all_vars_unique = len(variables) == len(set(variables))
+        assert all_vars_unique
+
         self.particulator = None
         self.dt = dt
         self.mesh = mesh
@@ -57,11 +60,7 @@ def get_predicted(self, index):
             )
         return self._values["predicted"][index]
 
-    def sync(self):
-        target = self._tmp
-        target["water_vapour_mixing_ratio"].ravel(self.get_water_vapour_mixing_ratio())
-        target["thd"].ravel(self.get_thd())
-
+    def _recalculate_temperature_pressure_relative_humidity(self, target):
         self.particulator.backend.temperature_pressure_rh(
             rhod=target["rhod"],
             thd=target["thd"],
@@ -70,6 +69,14 @@ def sync(self):
             p=target["p"],
             RH=target["RH"],
         )
+
+    def sync(self):
+        target = self._tmp
+        target["water_vapour_mixing_ratio"].ravel(self.get_water_vapour_mixing_ratio())
+        target["thd"].ravel(self.get_thd())
+
+        self._recalculate_temperature_pressure_relative_humidity(target)
+
         if "a_w_ice" in self.variables:
             self.particulator.backend.a_w_ice(
                 T=target["T"],

PySDM/environments/parcel.py

@@ -2,6 +2,8 @@
 Zero-dimensional adiabatic parcel framework
 """
 
+from typing import List, Optional
+
 import numpy as np
 
 from PySDM.environments.impl.moist import Moist
@@ -24,10 +26,10 @@ def __init__(
         w: [float, callable],
         z0: float = 0,
         mixed_phase=False,
+        variables: Optional[List[str]] = None,
     ):
-        super().__init__(
-            dt, Mesh.mesh_0d(), ["rhod", "z", "t"], mixed_phase=mixed_phase
-        )
+        variables = (variables or []) + ["rhod", "z", "t"]
+        super().__init__(dt, Mesh.mesh_0d(), variables, mixed_phase=mixed_phase)
 
         self.p0 = p0
         self.initial_water_vapour_mixing_ratio = initial_water_vapour_mixing_ratio
@@ -38,7 +40,7 @@ def __init__(
         self.w = w if callable(w) else lambda _: w
 
         self.formulae = None
-        self.delta_liquid_water_mixing_ratio = None
+        self.delta_liquid_water_mixing_ratio = np.nan
         self.params = None
 
     @property
@@ -116,6 +118,7 @@ def advance_parcel_vars(self):
             - self.delta_liquid_water_mixing_ratio / 2
         )
 
+        # derivate evaluated at p_old, T_old, mixrat_mid, w_mid
         drho_dz = self.formulae.hydrostatics.drho_dz(
             g=self.formulae.constants.g_std,
             p=p,

PySDM/physics/isotope_ratio_evolution/merlivat_and_jouzel_1979.py

@@ -1,5 +1,7 @@
 """
 see derivation of eq. (12) in [Merlivat and Jouzel 1979](https://doi.org/10.1029/JC084iC08p05029)
+(for constant alpha leads to eq. (13) in
+[Gedzelman & Arnold 1994](https://doi.org/10.1029/93jd03518))
 """
 
 

examples/PySDM_examples/Rozanski_and_Sonntag_1982/__init__.py

@@ -0,0 +1,6 @@
+"""
+based on Rozanski and Sonntag 1982 (Tellus)
+https://doi.org/10.3402/tellusa.v34i2.10795
+"""
+
+# pylint: disable=invalid-name