PySDM_examples.Rozanski_and_Sonntag_1982.multibox

 1"""
 2definition of a multi-box environment (parcel iteratively adjusting to a stationary state)
 3plus helper equation for isotopic ratio evolution
 4"""
 5
 6# pylint: disable=invalid-name
 7
 8import numpy as np
 9
10from PySDM.environments import Parcel
11
12
13def Rv_prim(*, Rl, Nl, Rv, Nv, dNl, Rr, K, a):
14    """eq. (2) in [Rozanski and Sonntag 1982](https://doi.org/10.3402/tellusa.v34i2.10795)"""
15    return (Rl * Nl + Rv * Nv + dNl * Rr * K) / ((Nl + dNl * K) * a + Nv)
16
17
18class MultiBox(Parcel):
19    """iterative parcel model in which each new iteration operates with ambient isotopic profile
20    resultant from the previous iteration, leading to a stationary state
21    """
22
23    def __init__(
24        self,
25        *,
26        isotopes,
27        delta_nl,
28        rain_isotope_ratios,
29        nt,
30        autoconversion_mixrat_threshold,
31        isotope_exchange_factor,
32        **kwargs,
33    ):
34        super().__init__(
35            variables=[
36                f"{ratio}_{isotope}" for ratio in ("Rv", "Rr") for isotope in isotopes
37            ],
38            **kwargs,
39        )
40        self.isotopes = isotopes
41        self.delta_nl = delta_nl
42        self.rain_isotope_ratios = rain_isotope_ratios
43        self.nt = nt
44        self.autoconversion_mixrat_threshold = autoconversion_mixrat_threshold
45        self.isotope_exchange_factor = isotope_exchange_factor
46
47    def advance_parcel_vars(self):
48        """explicit Euler integration of isotope-ratio time derivative"""
49        assert self.delta_liquid_water_mixing_ratio >= 0
50        self._recalculate_temperature_pressure_relative_humidity(self._tmp)
51
52        alpha_old = {}
53        dRv__dt = {}
54        for isotope in self.isotopes:
55            alpha_fun = getattr(
56                self.particulator.formulae.isotope_equilibrium_fractionation_factors,
57                f"alpha_l_{isotope}",
58            )
59            alpha_old[isotope] = alpha_fun(self["T"][0])
60            alpha_new = alpha_fun(self._tmp["T"][0])
61
62            dRv__dt[isotope] = self[f"Rv_{isotope}"][
63                0
64            ] * self.particulator.formulae.isotope_ratio_evolution.d_Rv_over_Rv(
65                alpha=alpha_old[isotope],
66                d_alpha=(alpha_new - alpha_old[isotope]) / self.dt,
67                n_vapour=self["water_vapour_mixing_ratio"][0],
68                d_n_vapour=-self.delta_liquid_water_mixing_ratio / self.dt,
69                n_liquid=self.autoconversion_mixrat_threshold,  # TODO #1207
70            )
71
72        super().advance_parcel_vars()
73
74        for isotope in self.isotopes:
75            self.particulator.backend.explicit_euler(
76                self._tmp[f"Rv_{isotope}"], self.particulator.dt, dRv__dt[isotope]
77            )
78            level = self.particulator.n_steps
79            if self.delta_nl is not None:
80                self._tmp[f"Rv_{isotope}"][:] = Rv_prim(
81                    Rl=alpha_old[isotope] * self._tmp[f"Rv_{isotope}"][0],
82                    Nl=self.autoconversion_mixrat_threshold,
83                    Rv=self._tmp[f"Rv_{isotope}"][0],
84                    Nv=self["water_vapour_mixing_ratio"][0],
85                    dNl=np.sum(self.delta_nl[level:]),
86                    Rr=self.rain_isotope_ratios[isotope][
87                        min(
88                            level + 2, self.nt  # TODO #1207: this warrants a comment...
89                        )
90                    ],
91                    K=self.isotope_exchange_factor,
92                    a=alpha_old[isotope],
93                )
94            self._tmp[f"Rr_{isotope}"][:] = (
95                alpha_old[isotope] * self._tmp[f"Rv_{isotope}"][0]
96            )