Adding the changes used during the upwelling change project.

konrad/convection.py

@@ -263,11 +263,25 @@ 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)
+                #print("diff_neg=",diff_neg)
+
+                return T_con, surfaceT
 
             # to avoid getting stuck in a loop if something weird is going on
             counter += 1
+<<<<<<< HEAD
             if counter == 100:
                 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
 
@@ -421,6 +435,34 @@ 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)
+
+        # Entrain dry air into the convective atmosphere column.
+        T_con = self._entrainment.entrain(T_con, atmosphere)
+
+        T_con_overshoot = T_con
+        # Define overshoot weight T_rad and T_con between Zcontop and Zcontop+10km
+        # Weight for T_con: 1/2^((Z-Zcontop)/1000+7)
+        # Weight for T_rad: 1-1/2^((Z-Zcontop)/1000+7)
+
+        z = atmosphere['z'][0,:]/1000.
+
+        # Combine radiative and convective temperature profiles.
+        if np.any(T_con > T_rad):
+            contop = np.max(np.where(T_con > T_rad))
+            z_contop = (z[contop]+z[contop+1])/2.
+            weight_con = 1./2**(z[contop+1:]-z_contop+7)
+            T_con_overshoot[contop+1:] = (weight_con)*T_con[contop+1:] + (1-weight_con)*T_rad[contop+1:]
+        else:
+            return T_rad
+
+        return T_con_overshoot
 
 class FixedDynamicalHeating(HardAdjustment):
     """Adjustment with a fixed convective (dynamical) heating rate."""

konrad/core.py

@@ -310,7 +310,7 @@ def is_converged(self):
             logger.info(d_txt.format(self.newDN, 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}')
         # If the equilibrium is larger than the threshold count, it declares
         #     convergence
         return self.counteq > self.post_count
@@ -405,6 +405,7 @@ def run(self):
 
             # Applies cooling due to stratospheric upwelling
             self.upwelling.cool(
+                surface=self.surface,
                 atmosphere=self.atmosphere,
                 convection=self.convection,
                 timestep=self.timestep_days,

konrad/ozone.py

@@ -43,6 +43,8 @@
 
 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",
@@ -358,7 +360,6 @@ def __init__(self, w=0, norm_level=None,
         """
         super().__init__(w=w, is_coupled_upwelling=is_coupled_upwelling)
         self.norm_level = norm_level
-        self._new_norm_level = None
         self._a = None
         self._profile_coupling = profile_coupling
 
@@ -375,18 +376,44 @@ def get_param_interpolations(self, atmosphere):
             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 / self.norm_level, a, bounds_error = False,
+            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._new_norm_level))
+            alist.append(self._a[param_num - 1](p / self.norm_level))
         return alist
 
     def __call__(self, atmosphere, timestep, upwelling, **kwargs):
@@ -400,14 +427,11 @@ def __call__(self, atmosphere, timestep, upwelling, **kwargs):
                 "It can be found here:\n"
                 "https://gitlab.com/simple-stratospheric-chemistry/simotrostra"
             )
-
-        if self.norm_level is None:
-            self.norm_level = self.get_norm_level(atmosphere)
 
         if self._a is None:
             self._a = self.get_param_interpolations(atmosphere)
 
-        self._new_norm_level = self.get_norm_level(atmosphere) 
+        self.norm_level = self.get_norm_level(atmosphere) 
 
         T = atmosphere['T'][0, :]
         p = atmosphere['plev']  # [Pa]

konrad/upwelling.py

@@ -79,10 +79,10 @@ class Upwelling(Component, metaclass=abc.ABCMeta):
     """Base class to define abstract methods for all upwelling handlers."""
 
     @abc.abstractmethod
-    def cool(self, atmosphere, convection, timestep):
-        """Cool the atmosphere according to an upwelling.
-
+    def cool(self, surface, atmosphere, convection, timestep):
+        """ Cool the atmosphere according to an upwelling.
         Parameters:
+            surface (konrad.surface.Surface): Surface model.
             atmosphere (konrad.atmosphere.Atmosphere): Atmosphere model.
             convection (konrad.convection): Convection model.
             timestep (float): Timestep width [day].
@@ -110,11 +110,12 @@ def __init__(self, w=0.2, lowest_level=None):
         self._w = w * meters_per_day  # in m/day
         self._lowest_level = lowest_level
 
-    def cool(self, atmosphere, convection, timestep):
+    def cool(self, surface, atmosphere, convection, timestep):
         """Apply cooling above the convective top (level where the net
         radiative heating becomes small).
 
         Parameters:
+            surface (konrad.surface.Surface): Surface model.
             atmosphere (konrad.atmosphere.Atmosphere): Atmosphere model.
             convection (konrad.convection): Convection model.
             timestep (float): Timestep width [day].
@@ -147,6 +148,88 @@ def cool(self, atmosphere, convection, timestep):
 
         self["cooling_rates"] = (("time", "plev"), -Q.reshape(1, -1))
 
+class Coupled_StratosphericUpwelling(Upwelling):
+    """Include an evolving upwelling in a logistic fashion"""
+    def __init__(self, w=0.2, lowest_level=None, w_max=0.25, rate=1):
+        self._w = w * meters_per_day
+        self._lowest_level = lowest_level
+        self._w_o = w
+        self._w_f = w_max
+        self._rate = rate
+        self._Tso = None
+
+    def _update_upwelling(self, surface):
+        """Updates the upwelling speed.
+        
+        Parameters:
+            surface (konrad.surface.Surface): Surface model.
+        
+        Notes:
+            
+            This updating model considers the relationship between warmer surface tem-
+            perature and increasing stratospheric upwelling speed. The model links
+            warming and upwelling speed through a logistic model, where there is an
+            initial pre-industrial upwelling and a final forcing-dependent upwelling.
+            The rate at which the upwelling increases is prescribed.
+            
+            I use the solutions of the equation
+            
+            .. math:: \\frac{d\, w}{d(\Delta T)} = r w \\left(1-\\frac{w}{w_{f}}\\right)
+            
+            with the constraint :math:w(0)=w_{o} and :math:r is the rate of growth.
+        """
+        new_w = np.exp(- self._rate * (surface["temperature"][-1] - self._Tso))
+        new_w *= ((self._w_f / self._w_o) - 1)
+        new_w += 1
+        new_w = 1 / new_w
+        new_w *= self._w_f
+        return new_w
+
+    def cool(self, surface, atmosphere, convection, timestep):
+        """Apply cooling above the convective top (level where the net
+        radiative heating becomes small).
+        
+        Parameters:
+            surface (konrad.surface.Surface): Surface model.
+            atmosphere (konrad.atmosphere.Atmosphere): Atmosphere model.
+            convection (konrad.convection): Convection model.
+            timestep (float): Timestep width [day].
+        """
+        # Get the first surface temperature to calculate the surface temperature change
+        if self._Tso is None:
+            self._Tso = surface["temperature"][-1]
+
+        # Calculate the new upwelling speed
+        self._w = self._update_upwelling(surface) * meters_per_day
+
+        # Ensure that the variable is initialized even if the upwelling is not
+        # applied. A proper initaliziation is needed to write netCDF output.
+        self.create_variable(
+            'cooling_rates',
+            dims=('time', 'plev'),
+            data=np.zeros_like(atmosphere["plev"]).reshape(1, -1),
+        )
+
+        T = atmosphere['T'][0, :]
+        z = atmosphere['z'][0, :]
+        Cp = atmosphere.get_heat_capacity()
+
+        if self._lowest_level is not None:
+            above_level_index = self._lowest_level
+        else:
+            above_level_index = convection.get('convective_top_index')[0]
+            if np.isnan(above_level_index):
+                # if convection hasn't been applied and a lowest level for the
+                # upwelling has not been specified, upwelling is not applied
+                return
+        above_level_index = int(np.round(above_level_index))
+
+        Q = cooling_rates(T, z, self._w, Cp, above_level_index)
+
+        atmosphere['T'][0, :] += Q * timestep
+
+        self['cooling_rates'] = (('time', 'plev'), -Q.reshape(1, -1))
+
 
 class SpecifiedCooling(Upwelling):
     """Include an upwelling with specified cooling"""