PR #473: from firedrakeproject/fix_SWSaturationAdjustment

Remove L_v from SWSaturationAdjustment

examples/shallow_water/moist_thermal_williamson5.py

@@ -32,9 +32,8 @@
 NP = -20*epsilon
 mu1 = 0.05
 mu2 = 0.98
-L_v = 10
 q0 = 135  # chosen to give an initial max vapour of approx 0.02
-beta2 = 1
+beta2 = 10
 qprecip = 10e-4
 gamma_r = 10e-3
 # topography parameters
@@ -94,10 +93,10 @@ def sat_func(x_in):
 def gamma_v(x_in):
     h = x_in.split()[1]
     b = x_in.split()[2]
-    return (1 + L_v*(20*q0/(g*h + g*tpexpr) * exp(20*(1 - b/g))))**(-1)
+    return (1 + beta2*(20*q0/(g*h + g*tpexpr) * exp(20*(1 - b/g))))**(-1)
 
 
-SWSaturationAdjustment(eqns, sat_func, L_v, time_varying_saturation=True,
+SWSaturationAdjustment(eqns, sat_func, time_varying_saturation=True,
                        parameters=parameters, thermal_feedback=True,
                        beta2=beta2, gamma_v=gamma_v,
                        time_varying_gamma_v=True)

gusto/physics.py

@@ -938,7 +938,7 @@ class SWSaturationAdjustment(PhysicsParametrisation):
 
     """
 
-    def __init__(self, equation, saturation_curve, L_v=None,
+    def __init__(self, equation, saturation_curve,
                  time_varying_saturation=False, vapour_name='water_vapour',
                  cloud_name='cloud_water', convective_feedback=False,
                  beta1=None, thermal_feedback=False, beta2=None, gamma_v=1,
@@ -956,11 +956,6 @@ def __init__(self, equation, saturation_curve, L_v=None,
                 prognostic field.
             time_varying_saturation (bool, optional): set this to True if the
                 saturation curve is changing in time. Defaults to False.
-            L_v (float, optional): The air expansion factor multiplied by the
-                latent heat due to phase change divided by the specific heat
-                capacity. For the atmosphere we take L_v to be 10, following A.2
-                in Zerroukat and Allen (2015). Defaults to None but must be
-                specified if using thermal feedback.
             vapour_name (str, optional): name of the water vapour variable.
                 Defaults to 'water_vapour'.
             cloud_name (str, optional): name of the cloud variable. Defaults to
@@ -975,7 +970,8 @@ def __init__(self, equation, saturation_curve, L_v=None,
                 affect the buoyancy equation. Defaults to False.
             beta2 (float, optional): Condensation proportionality constant
                 for thermal feedback. Defaults to None, but must be specified
-                if thermal_feedback is True.
+                if thermal_feedback is True. This is equivalent to the L_v
+                parameter in Zerroukat and Allen (2015).
             gamma_v (ufl expression or :class: `function`): The proportion of
                 moist species that is converted when a conversion between
                 vapour and cloud is taking place. Defaults to one, in which
@@ -1015,7 +1011,6 @@ def __init__(self, equation, saturation_curve, L_v=None,
         if self.thermal_feedback:
             assert "b" in equation.field_names, "Buoyancy field must exist for thermal feedback"
             assert beta2 is not None, "If thermal feedback is used, beta2 parameter must be specified"
-            assert L_v is not None, "If thermal feedback is used, L_v parameter must be specified"
 
         # Obtain function spaces and functions
         W = equation.function_space
@@ -1086,7 +1081,7 @@ def __init__(self, equation, saturation_curve, L_v=None,
         if convective_feedback:
             factors.append(self.gamma_v*beta1)
         if thermal_feedback:
-            factors.append(parameters.g*L_v*self.gamma_v*beta2)
+            factors.append(parameters.g*self.gamma_v*beta2)
 
         # Add terms to equations and make interpolators
         self.source = [Function(Vc) for factor in factors]

integration-tests/physics/test_sw_saturation_adjustment.py

@@ -30,8 +30,7 @@ def run_sw_cond_evap(dirname, process):
     theta_c = pi
     lamda_c = pi/2
     rc = R/4
-    L_v = 10
-    beta2 = 1
+    beta2 = 10
 
     # Domain
     mesh = IcosahedralSphereMesh(radius=R, refinement_level=3, degree=2)
@@ -61,7 +60,7 @@ def run_sw_cond_evap(dirname, process):
             diagnostic_fields=[Sum('water_vapour', 'cloud_water')])
 
     # Physics schemes
-    physics_schemes = [(SWSaturationAdjustment(eqns, sat, L_v=L_v,
+    physics_schemes = [(SWSaturationAdjustment(eqns, sat,
                                                parameters=parameters,
                                                thermal_feedback=True,
                                                beta2=beta2),
@@ -89,7 +88,7 @@ def run_sw_cond_evap(dirname, process):
         v_true = Function(v0.function_space()).interpolate(sat*(0.96+0.005*pert))
         c_true = Function(c0.function_space()).interpolate(Constant(0.0))
         # gain buoyancy
-        factor = parameters.g*L_v
+        factor = parameters.g*beta2
         sat_adj_expr = (v0 - sat) / dt
         sat_adj_expr = conditional(sat_adj_expr < 0,
                                    max_value(sat_adj_expr, -c0 / dt),
@@ -104,7 +103,7 @@ def run_sw_cond_evap(dirname, process):
         v_true = Function(v0.function_space()).interpolate(Constant(sat))
         c_true = Function(c0.function_space()).interpolate(v0 - sat)
         # lose buoyancy
-        factor = parameters.g*L_v
+        factor = parameters.g*beta2
         sat_adj_expr = (v0 - sat) / dt
         sat_adj_expr = conditional(sat_adj_expr < 0,
                                    max_value(sat_adj_expr, -c0 / dt),