add perturbed ERA5 forcing (#1291)

Author: kmdeck

experiments/long_runs/low_res_snowy_land.jl

@@ -62,11 +62,12 @@ function setup_model(FT, start_date, stop_date, Δt, domain, earth_param_set)
         context,
         lowres = true,
     )
-    atmos, radiation = ClimaLand.prescribed_forcing_era5(
+    atmos, radiation = ClimaLand.prescribed_perturbed_forcing_era5(
         era5_ncdata_path,
         surface_space,
         start_date,
         earth_param_set,
+        5.0,
         FT;
         max_wind_speed = 25.0,
         time_interpolation_method = LinearInterpolation(PeriodicCalendar()),

src/shared_utilities/drivers.jl

@@ -37,6 +37,7 @@ export AbstractAtmosphericDrivers,
     make_update_drivers,
     prescribed_lai_era5,
     prescribed_forcing_era5,
+    prescribed_perturbed_forcing_era5,
     prescribed_analytic_forcing,
     default_zenith_angle
 
@@ -1031,16 +1032,53 @@ function specific_humidity_from_dewpoint(
         sim_FT(T_dew_air) - _T_freeze,
         sim_FT(T_air) - _T_freeze,
     )
-
-    # In the future, we will use Thermodynamics functions for all of the below. Note that the type of T_air must match the type of the parameters in thermo_params.
-    esat::sim_FT = Thermodynamics.saturation_vapor_pressure(
+    q = Thermodynamics.q_vap_from_RH_liquid(
         thermo_params,
+        sim_FT(P_air),
         sim_FT(T_air),
-        Thermodynamics.Liquid(),
+        rh,
     )
+    return q
+end
+
+
+"""
+    perturbed_specific_humidity_from_dewpoint(dewpoint_temperature, temperature, air_pressure, earth_param_set, ΔT)
 
-    e = rh * esat
-    q = data_FT(0.622 * e / (P_air - 0.378 * e))
+Estimates the perturbed specific humidity given the true dewpoint temperature, true temperature of the air
+in Kelvin, and true air pressure in Pa, along with a perturbation to the temperature ΔT, and the ClimaLand
+`earth_param_set`. We first compute relative humidity, and then using RH and the perturbed air temperature,
+equal to  true temperature+ΔT, compute a perturbed `q`.
+
+Relative humidity is computed using the Magnus formula.
+
+For more information on the Magnus formula, see e.g.
+Lawrence, Mark G. (1 February 2005).
+"The Relationship between Relative Humidity and the Dewpoint Temperature in Moist Air:
+A Simple Conversion and Applications".
+Bulletin of the American Meteorological Society. 86 (2): 225–234.
+"""
+function perturbed_specific_humidity_from_dewpoint(
+    T_dew_air::data_FT,
+    T_air::data_FT,
+    P_air::data_FT,
+    earth_param_set,
+    ΔT,
+) where {data_FT <: Real}
+    thermo_params = LP.thermodynamic_parameters(earth_param_set)
+    _T_freeze = LP.T_freeze(earth_param_set)
+    sim_FT = typeof(_T_freeze)
+    # Obtain the relative humidity. This function requires temperatures in Celsius
+    rh::sim_FT = rh_from_dewpoint(
+        sim_FT(T_dew_air) - _T_freeze,
+        sim_FT(T_air) - _T_freeze,
+    )
+    q = Thermodynamics.q_vap_from_RH_liquid(
+        thermo_params,
+        sim_FT(P_air),
+        sim_FT(T_air + ΔT),
+        rh,
+    )
     return q
 end
 
@@ -1049,6 +1087,8 @@ end
 
 Returns the relative humidity given the dewpoint temperature in Celsius and the
 air temperature in Celsius, using the Magnus formula.
+Since this formula is not restricted to the [0,1] range,
+we enforce it.
 
 For more information on the Magnus formula, see e.g.
 Lawrence, Mark G. (1 February 2005).
@@ -1061,7 +1101,7 @@ function rh_from_dewpoint(Td_C::FT, T_C::FT) where {FT <: Real}
     b = FT(17.625) # unitless
     γ = Td_C * b / (c + Td_C)
     rh = exp(γ - b * T_C / (c + T_C))
-    return rh
+    return max(FT(0), min(rh, FT(1.0)))
 end
 
 """
@@ -1499,27 +1539,28 @@ function prescribed_forcing_era5(
     # Pass a list of files in all cases
     era5_ncdata_path isa String && (era5_ncdata_path = [era5_ncdata_path])
 
-    # Precip is provide as a mass flux; convert to volume flux of liquid water with ρ =1000 kg/m^3
+    _ρ_liq = LP.ρ_cloud_liq(earth_param_set)
+    # Precip is provided as a mass flux; convert to volume flux of liquid water with ρ = 1000 kg/m^3
     precip = TimeVaryingInput(
         [era5_ncdata_path, era5_ncdata_path],
         ["mtpr", "msr"],
         surface_space;
         start_date,
         regridder_type,
         regridder_kwargs = (; interpolation_method),
-        file_reader_kwargs = (; preprocess_func = (data) -> -data / 1000,),
+        file_reader_kwargs = (; preprocess_func = (data) -> -data / _ρ_liq,),
         method = time_interpolation_method,
         compose_function = (mtpr, msr) -> min.(mtpr .- msr, Float32(0)),
     )
-    # Precip is provide as a mass flux; convert to volume flux of liquid water with ρ =1000 kg/m^3
+    # Precip is provided as a mass flux; convert to volume flux of liquid water with ρ = 1000 kg/m^3
     snow_precip = TimeVaryingInput(
         era5_ncdata_path,
         "msr",
         surface_space;
         start_date,
         regridder_type,
         regridder_kwargs = (; interpolation_method),
-        file_reader_kwargs = (; preprocess_func = (data) -> -data / 1000,),
+        file_reader_kwargs = (; preprocess_func = (data) -> -data / _ρ_liq,),
         method = time_interpolation_method,
     )
     if max_wind_speed isa Nothing
@@ -1642,6 +1683,194 @@ function prescribed_forcing_era5(
     return (; atmos, radiation)
 end
 
+"""
+     prescribed_perturbed_forcing_era5(era5_ncdata_path,
+                             surface_space,
+                             start_date,
+                             earth_param_set,
+                             ΔT,
+                             FT;
+                             gustiness=1,
+                             max_wind_speed = nothing,
+                             c_co2 = TimeVaryingInput((t) -> 4.2e-4),
+                             time_interpolation_method = LinearInterpolation(PeriodicCalendar()),
+                             regridder_type = :InterpolationsRegridder,
+                             interpolation_method = Interpolations.Constant(),)
+
+A helper function which constructs the `PrescribedAtmosphere` and `PrescribedRadiativeFluxes`
+from a file path pointing to the ERA5 data in a netcdf file, the surface_space, the start date,
+and the earth_param_set, applying a change in the instantaneous temperature at each point
+of ΔT, while keeping the relative humidity fixed.
+
+Please see the documentation for `prescribed_forcing_era5` for more information.
+"""
+function prescribed_perturbed_forcing_era5(
+    era5_ncdata_path,
+    surface_space,
+    start_date,
+    earth_param_set,
+    ΔT,
+    FT;
+    gustiness = 1,
+    max_wind_speed = nothing,
+    c_co2 = TimeVaryingInput((t) -> 4.2e-4),
+    time_interpolation_method = LinearInterpolation(PeriodicCalendar()),
+    regridder_type = :InterpolationsRegridder,
+    interpolation_method = Interpolations.Constant(),
+)
+    # Pass a list of files in all cases
+    era5_ncdata_path isa String && (era5_ncdata_path = [era5_ncdata_path])
+    _ρ_liq = LP.ρ_cloud_liq(earth_param_set)
+    # Precip is provided as a mass flux; convert to volume flux of liquid water with ρ = 1000 kg/m^3
+    precip = TimeVaryingInput(
+        [era5_ncdata_path, era5_ncdata_path],
+        ["mtpr", "msr"],
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        file_reader_kwargs = (; preprocess_func = (data) -> -data / _ρ_liq,),
+        method = time_interpolation_method,
+        compose_function = (mtpr, msr) -> min.(mtpr .- msr, Float32(0)),
+    )
+    # Precip is provided as a mass flux; convert to volume flux of liquid water with ρ = 1000 kg/m^3
+    snow_precip = TimeVaryingInput(
+        era5_ncdata_path,
+        "msr",
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        file_reader_kwargs = (; preprocess_func = (data) -> -data / _ρ_liq,),
+        method = time_interpolation_method,
+    )
+    if max_wind_speed isa Nothing
+        wind_function = (u, v) -> sqrt.(u .^ 2 .+ v .^ 2)
+    else
+        wind_function = (u, v) -> min.(sqrt.(u .^ 2 .+ v .^ 2), max_wind_speed)
+    end
+
+    u_atmos = TimeVaryingInput(
+        [era5_ncdata_path, era5_ncdata_path],
+        ["u10", "v10"],
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        compose_function = wind_function,
+        method = time_interpolation_method,
+    )
+    specific_humidity(Td, T, P; params = earth_param_set, ΔT = ΔT) =
+        ClimaLand.perturbed_specific_humidity_from_dewpoint.(
+            Td,
+            T,
+            P,
+            params,
+            ΔT,
+        )
+    q_atmos = TimeVaryingInput(
+        [era5_ncdata_path, era5_ncdata_path, era5_ncdata_path],
+        ["d2m", "t2m", "sp"],
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        compose_function = specific_humidity,
+        method = time_interpolation_method,
+    )
+    P_atmos = TimeVaryingInput(
+        era5_ncdata_path,
+        "sp",
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        method = time_interpolation_method,
+    )
+    perturb_temp(T; ΔT = ΔT) = T .+ ΔT
+    T_atmos = TimeVaryingInput(
+        era5_ncdata_path,
+        "t2m",
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        file_reader_kwargs = (; preprocess_func = perturb_temp),
+        method = time_interpolation_method,
+    )
+    h_atmos = FT(10)
+
+    atmos = PrescribedAtmosphere(
+        precip,
+        snow_precip,
+        T_atmos,
+        u_atmos,
+        q_atmos,
+        P_atmos,
+        start_date,
+        h_atmos,
+        earth_param_set;
+        gustiness = FT(gustiness),
+        c_co2 = c_co2,
+    )
+
+    SW_d = TimeVaryingInput(
+        era5_ncdata_path,
+        "msdwswrf",
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        method = time_interpolation_method,
+    )
+    function compute_diffuse_fraction(total, direct)
+        diff = max(total - direct, Float32(0))
+        return min(diff / (total + eps(Float32)), Float32(1))
+    end
+    function compute_diffuse_fraction_broadcasted(total, direct)
+        return @. compute_diffuse_fraction(total, direct)
+    end
+
+    frac_diff = TimeVaryingInput(
+        [era5_ncdata_path, era5_ncdata_path],
+        ["msdwswrf", "msdrswrf"],
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        method = time_interpolation_method,
+        compose_function = compute_diffuse_fraction_broadcasted,
+    )
+    LW_d = TimeVaryingInput(
+        era5_ncdata_path,
+        "msdwlwrf",
+        surface_space;
+        start_date,
+        regridder_type,
+        regridder_kwargs = (; interpolation_method),
+        method = time_interpolation_method,
+    )
+    zenith_angle =
+        (t, s) -> default_zenith_angle(
+            t,
+            s;
+            latitude = ClimaCore.Fields.coordinate_field(surface_space).lat,
+            longitude = ClimaCore.Fields.coordinate_field(surface_space).long,
+            insol_params = earth_param_set.insol_params,
+        )
+
+    radiation = PrescribedRadiativeFluxes(
+        FT,
+        SW_d,
+        LW_d,
+        start_date;
+        θs = zenith_angle,
+        earth_param_set = earth_param_set,
+        frac_diff = frac_diff,
+    )
+    return (; atmos, radiation)
+end
+
 
 """
      prescribed_lai_era5(era5_lai_ncdata_path,