Final modifications for the paper. This is the working version for th…

…e paper.

konrad/atmosphere.py

@@ -377,14 +377,14 @@ def get_cold_point_plev(self, interpolate=False):
             # Return the single coldest point on the actual pressure grid.
             return self["plev"][self.get_cold_point_index()]
 
-    def get_thermal_tropopause_index(self):
+    def get_thermal_and_ozone_tropopause_index(self):
         """Return the model level index of the thermal tropopause.
         
         Returns:
             int: Model level index at the thermal tropopause as defined by the WMO.
         """
         # Calculates the lapse-rate in Kelvin per kilometre
-        lapse_rate = np.gradient(self["T"][0, :], self["z"][0, :] / 1000)
+        lapse_rate_K_per_km = np.gradient(self["T"][0, :], self["z"][0, :] / 1000)
 
         # Gets the indices where the lapse rate decreases at less than two Kelvin per
         # kilometre
@@ -401,23 +401,30 @@ def get_thermal_tropopause_index(self):
 
             index = mask[i]
 
-            limit_test = np.where(self["z"][0, :]-self["z"][0, index] >= 2)[0][0]
-            region = lapse_rate[index:limit_test + 1]
+            limit_test = np.where(self["z"][0, :] - self["z"][0, index] >= 2)[0][0]
+            region = lapse_rate_K_per_km[index: limit_test + 1]
 
             if np.all(region >= -2):
                 thermal_tropopause = index
             i += 1
-
-        return thermal_tropopause
-
-    def get_thermal_tropopause_plev(self):
-        """Return the thermal tropopause pressure.
 
-        Returns:
-            float: Pressure at the thermal tropopause [Pa].
-        """
+        ozone_tropopause = np.where(
+            self["z"][0, :]- self["z"][0, thermal_tropopause] > -1
+        )[0][0]
+        ozone_tropopause_sub = np.where(
+            self["z"][0, :]- self["z"][0, ozone_tropopause] > -2
+        )[0][0]
+        ozone_tropopause_super = np.where(
+            self["z"][0, :]- self["z"][0, thermal_tropopause] > 2
+        )[0][0]
 
-        return self['plev'][self.get_thermal_tropopause_index()]
+
+        return (
+            thermal_tropopause,
+            ozone_tropopause_sub,
+            ozone_tropopause,
+            ozone_tropopause_super,
+        )
 
     def get_triple_point_index(self):
         """Return the model level index at the triple point.

konrad/convection.py

@@ -263,7 +263,7 @@ def convective_adjustment(self, atmosphere, lapse, surface, timestep=0.1):
             else:
                 diff_neg = diff
                 surfaceTneg = surfaceT
-                
+
             if  np.abs(surfaceTpos - surfaceTneg) < 1e-7:
                 #print("counter=",counter)
                 #print("diff_pos=",diff_pos)
@@ -273,15 +273,8 @@ def convective_adjustment(self, atmosphere, lapse, surface, timestep=0.1):
 
             # to avoid getting stuck in a loop if something weird is going on
             counter += 1
-<<<<<<< HEAD
-            if counter == 100:
+            if counter == 10000:
                 raise ValueError("No energy conserving convective profile can be found")
-=======
-            if counter == 3000:
-                raise ValueError(
-                    "No energy conserving convective profile can be found"
-                )
->>>>>>> b978d81 (Adding the changes used during the upwelling change project.)
 
         return T_con, surfaceT
 
@@ -435,14 +428,9 @@ def update_convective_top_height(self, z, lim=0.2):
         self.create_variable("convective_top_height", np.array([contop_z]))
         return
 
-class HardAdjustment_Overshoot(HardAdjustment):
-
-    def convective_profile(self, atmosphere, surfaceT, lapse, **kwargs):
-        # The convective adjustment is only applied to the atmospheric profile,
-        # if it causes heating somewhere
-        T_rad = atmosphere["T"][-1]
-        T_con = self.get_moist_adiabat(atmosphere, surfaceT, lapse)
 
+class FixedDynamicalHeating(HardAdjustment):
+    """Adjustment with a fixed convective (dynamical) heating rate."""
         # Entrain dry air into the convective atmosphere column.
         T_con = self._entrainment.entrain(T_con, atmosphere)
 
@@ -464,9 +452,6 @@ def convective_profile(self, atmosphere, surfaceT, lapse, **kwargs):
 
         return T_con_overshoot
 
-class FixedDynamicalHeating(HardAdjustment):
-    """Adjustment with a fixed convective (dynamical) heating rate."""
-
     def __init__(self, heating=None, *args, **kwargs):
         """
         Parameters:

konrad/core.py

@@ -283,7 +283,7 @@ def is_converged(self):
         self.newDDN = np.abs(self.newDN) - np.abs(self.oldDN)
 
         # Checks whether the difference is below the threshold
-        test1 = np.abs(self.newDN) <= self.delta
+        test1 = (np.abs(self.newDN) / self.timestep_days) <= self.delta
         # Checks whether the second order difference is below the threshold
         test2 = (np.abs(self.newDDN) / self.timestep_days ** 2) <= self.delta2
 
@@ -307,10 +307,12 @@ def is_converged(self):
             d_txt = "Days within equilibrium conditions: {0:3.2f}"
             logger.debug(d_txt.format(self.counteq.total_seconds() / (60 * 60 * 24)))
             d_txt = "Delta N (TOA): {0:2.2e} (Threshold: {1:2.2e})"
-            logger.info(d_txt.format(self.newDN, self.delta))
+            logger.debug(d_txt.format(self.newDN / self.timestep_days, self.delta))
             d_txt = "Delta (Delta N (TOA)): {0:2.2e} (Threshold: {1:2.2e})"
-            logger.info(d_txt.format(self.newDDN, self.delta2))
-            logger.info(f'Upwelling speed: {self.upwelling._w}')
+            logger.debug(
+                d_txt.format(self.newDDN / self.timestep_days ** 2, self.delta2)
+            )
+
         # If the equilibrium is larger than the threshold count, it declares
         #     convergence
         return self.counteq > self.post_count

konrad/entrainment.py

@@ -84,6 +84,7 @@ def entrain(self, T_con_adiabat, atmosphere):
 
         RH_hlev = interp1d(np.log(p), RH, fill_value="extrapolate")(np.log(phlev[:-1]))
 
+        entr = self.entr
         deltaT = np.ones_like(p) * 0.0
         k_cb = np.max(np.where(p >= 96000.0))
 

konrad/ozone.py

@@ -43,8 +43,6 @@
 
 from konrad import constants
 from konrad.component import Component
-from konrad.constants import earth_standard_gravity as g
-from typhon.physics import density
 
 __all__ = [
     "Ozone",
@@ -352,108 +350,52 @@ def __init__(self, w=0, norm_level=None,
                  is_coupled_upwelling=False, profile_coupling="cold_point"):
         """
         Parameters:
-            w (ndarray / int / float): upwelling velocity [mm / s]
-            is_coupled_upwelling (bool): whether to couple to upwelling model
-            profile_coupling (str): Mechanism to determine the normalisation level:
-                    * "cold_point": Coupled to the cold-point tropopause
-                    * "thermal_tropopause": Coupled to the WMO thermal tropopause
+            p (ndarray): pressure levels [Pa]
+            T (ndarray): temperature values [K]
+        Returns:
+            ndarray: density of molecules [number of molecules / m3]
         """
-        super().__init__(w=w, is_coupled_upwelling=is_coupled_upwelling)
-        self.norm_level = norm_level
-        self._a = None
-        self._profile_coupling = profile_coupling
-
-    def get_norm_level(self, atmosphere):
-        if self._profile_coupling == 'cold_point':
-            return atmosphere.get_cold_point_plev(interpolate=True)
-        elif self._profile_coupling == 'thermal_tropopause':
-            return atmosphere.get_thermal_tropopause_plev()
-        else:
-            raise ValueError("Invalid coupling mechanism.")
-
-    def get_param_interpolations(self, atmosphere):
-        cariolle_data = Dataset(
-            os.path.join(os.path.dirname(__file__),
-                         'data/Cariolle_data.nc'))
-        p_data = cariolle_data['p'][:]
-        T_data = cariolle_data['A5'][:]
-
-        # Calculate normalisation level from the Cariolle data
-        if self._profile_coupling == "cold_point":
-            cp_idx = np.where(T_data == np.min(T_data[p_data > 1000]))[0][0]
-            p_data /= p_data[cp_idx]
-        elif self._profile_coupling == "thermal_tropopause":
-            z_data = np.cumsum(-np.gradient(np.flip(p_data)) /
-                 (density(np.flip(p_data), np.flip(Tp_data)) * g)
-            )
-            lapse_rate_km = np.gradient(T_data, z_data / 1000)
-            mask = np.where(lapse_rate_km > -2)[0]
-            tt_idx = False
-            i=0
-            while not tt_idx:
-                index = mask[i]
-                limit_test = np.where(z_data - z_data[index] >= 2)[0][0]
-                region = lapse_rate_km[index:limit_test + 1]
-                if np.all(region >= -2):
-                    tt_idx = index
-                i += 1
-            p_data /= p_data[tt_idx]
-        else:
-            raise ValueError("Invalid coupling mechanism.")
-
-        # Normalise the Cariolle pressure profile and create the interpolators
-        alist = []
-        for param_num in range(1, 8):
-            a = cariolle_data[f'A{param_num}'][:]
-            alist.append(interp1d(p_data, a, bounds_error = False,
-                                  fill_value=(a[0],a[-1]))
-            )
-        return alist
-
-    def get_params(self,p):
-        alist = []
-        for param_num in range(1, 8):
-            alist.append(self._a[param_num - 1](p / self.norm_level))
-        return alist
+        return (constants.avogadro * p) / (constants.molar_gas_constant_dry_air * T)
 
-    def __call__(self, atmosphere, timestep, upwelling, **kwargs):
-
-        try:
-            from simotrostra.utils import overhead_molecules
-        except ModuleNotFoundError:
-            raise ModuleNotFoundError(
-                "No module named 'simotrostra'. You need to install it in order"
-                " to use this ozone class.\n"
-                "It can be found here:\n"
-                "https://gitlab.com/simple-stratospheric-chemistry/simotrostra"
-            )
+    def overhead_molecules(self, gas, p, z, T):
+        """
+        Parameters:
+            gas (ndarray): gas concentration [ppv] corresponding to levels p
+            p (ndarray): pressure levels [Pa]
+            z (ndarray): height values [m]
+            T (ndarray): temperature values [K]
+        Returns:
+            ndarray: number of molecules per m2 in overhead column
+        """
+        # molecules / m2 of air in each layer
+        molecules_air = self.density_of_molecules(p, T) * np.gradient(z)
+        # overhead column in molecules / m2
+        col = np.flipud(np.cumsum(np.flipud(gas * molecules_air)))
 
-        if self._a is None:
-            self._a = self.get_param_interpolations(atmosphere)
-
-        self.norm_level = self.get_norm_level(atmosphere) 
+        return col
 
-        T = atmosphere['T'][0, :]
-        p = atmosphere['plev']  # [Pa]
-        o3 = atmosphere['O3'][0, :]  # moles of ozone / moles of air
-        z = atmosphere['z'][0, :]  # m
+    def __call__(self, atmosphere, timestep, upwelling, **kwargs):
+        T = atmosphere["T"][0, :]
+        p = atmosphere["plev"]  # [Pa]
+        o3 = atmosphere["O3"][0, :]  # moles of ozone / moles of air
+        z = atmosphere["z"][0, :]  # m
 
-        o3col = overhead_molecules(o3, p, z, T
-                                   ) * 10 ** -4  # in molecules / cm2
+        o3col = self.overhead_molecules(o3, p, z, T) * 10 ** -4  # in molecules / cm2
 
         A1, A2, A3, A4, A5, A6, A7 = self.get_params(p)
         # A7 is in molecules / cm2
         # tendency of ozone volume mixing ratio per second
-        do3dt = A1 + A2*(o3 - A3) + A4*(T - A5) + A6*(o3col - A7)
+        do3dt = A1 + A2 * (o3 - A3) + A4 * (T - A5) + A6 * (o3col - A7)
 
         # transport term
         transport_ox = self.ozone_transport(o3, z, upwelling)
 
-        atmosphere['O3'] = (
-            ('time', 'plev'),
-            (o3 + (do3dt * 24 * 60**2 + transport_ox) * timestep).reshape(1, -1)
+        atmosphere["O3"] = (
+            ("time", "plev"),
+            (o3 + (do3dt * 24 * 60 ** 2 + transport_ox) * timestep).reshape(1, -1),
         )
 
+
 class Simotrostra(Cariolle):
     """Wrapper for Ed Charlesworth's simple chemistry scheme."""