parameterize sst distribution

config/model_configs/rcemipii_box_diagnostic_edmfx.yml

@@ -13,13 +13,14 @@ edmfx_detr_model: "Generalized"
 edmfx_nh_pressure: true
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
+rayleigh_sponge: true
 moist: equil
 precip_model: 0M
 override_τ_precip: false
 dt: 30secs
 t_end: 3600secs
 dt_save_state_to_disk: 12hours
-toml: [toml/diagnostic_edmfx_0M.toml]
+toml: [toml/rcemipii_diagnostic_edmfx_0M.toml]
 netcdf_interpolation_num_points: [8, 8, 60]
 netcdf_output_at_levels: true
 diagnostics:

config/model_configs/rcemipii_sphere_diagnostic_edmfx.yml

@@ -12,11 +12,12 @@ edmfx_detr_model: "Generalized"
 edmfx_nh_pressure: true
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
+rayleigh_sponge: true
 moist: equil
 precip_model: 0M
 override_τ_precip: false
 dt: 100secs
 t_end: 12hours
 dt_save_state_to_disk: 12hours
-toml: [toml/diagnostic_edmfx_0M.toml]
+toml: [toml/rcemipii_diagnostic_edmfx_0M.toml]
 ode_algo: ARS343

src/parameters/Parameters.jl

@@ -12,8 +12,12 @@ const ACAP = AbstractClimaAtmosParameters
 abstract type AbstractTurbulenceConvectionParameters end
 const ATCP = AbstractTurbulenceConvectionParameters
 
+abstract type AbstractSurfaceTemperatureParameters end
+const ASTP = AbstractSurfaceTemperatureParameters
+
 Base.broadcastable(param_set::ACAP) = tuple(param_set)
 Base.broadcastable(param_set::ATCP) = tuple(param_set)
+Base.broadcastable(param_set::ASTP) = tuple(param_set)
 
 Base.@kwdef struct TurbulenceConvectionParameters{FT} <: ATCP
     surface_area::FT
@@ -44,6 +48,13 @@ Base.@kwdef struct TurbulenceConvectionParameters{FT} <: ATCP
     max_area_limiter_power::FT
 end
 
+Base.@kwdef struct SurfaceTemperatureParameters{FT} <: ASTP
+    SST_mean::FT
+    SST_delta::FT
+    SST_wavelength::FT
+    SST_wavelength_latitude::FT
+end
+
 Base.@kwdef struct ClimaAtmosParameters{
     FT,
     TP,
@@ -54,6 +65,7 @@ Base.@kwdef struct ClimaAtmosParameters{
     WP,
     SFP,
     TCP,
+    STP,
 } <: ACAP
     thermodynamics_params::TP
     rrtmgp_params::RP
@@ -63,6 +75,7 @@ Base.@kwdef struct ClimaAtmosParameters{
     water_params::WP
     surface_fluxes_params::SFP
     turbconv_params::TCP
+    surface_temp_params::STP
     Omega::FT
     f_plane_coriolis_frequency::FT
     planet_radius::FT

src/parameters/create_parameters.jl

@@ -41,6 +41,18 @@ function TurbulenceConvectionParameters(toml_dict::CP.AbstractTOMLDict)
     CAP.TurbulenceConvectionParameters{FT}(; parameters...)
 end
 
+function SurfaceTemperatureParameters(toml_dict::CP.AbstractTOMLDict)
+    name_map = (;
+        :SST_mean => :SST_mean,
+        :SST_delta => :SST_delta,
+        :SST_wavelength => :SST_wavelength,
+        :SST_wavelength_latitude => :SST_wavelength_latitude,
+    )
+    parameters = CP.get_parameter_values(toml_dict, name_map, "ClimaAtmos")
+    FT = CP.float_type(toml_dict)
+    CAP.SurfaceTemperatureParameters{FT}(; parameters...)
+end
+
 function create_parameter_set(config::AtmosConfig)
     (; toml_dict, parsed_args) = config
     FT = CP.float_type(toml_dict)
@@ -64,6 +76,9 @@ function create_parameter_set(config::AtmosConfig)
         SF.Parameters.SurfaceFluxesParameters(toml_dict, UF.BusingerParams)
     SFP = typeof(surface_fluxes_params)
 
+    surface_temp_params = SurfaceTemperatureParameters(toml_dict)
+    STP = typeof(surface_temp_params)
+
     moisture_model = parsed_args["moist"]
     microphysics_cloud_params = if moisture_model == "nonequil"
         (;
@@ -128,7 +143,7 @@ function create_parameter_set(config::AtmosConfig)
         :optics_lookup_temperature_max => :optics_lookup_temperature_max,
     )
     parameters = CP.get_parameter_values(toml_dict, name_map, "ClimaAtmos")
-    return CAP.ClimaAtmosParameters{FT, TP, RP, IP, MPC, MPP, WP, SFP, TCP}(;
+    return CAP.ClimaAtmosParameters{FT, TP, RP, IP, MPC, MPP, WP, SFP, TCP, STP}(;
         parameters...,
         thermodynamics_params,
         rrtmgp_params,
@@ -138,5 +153,6 @@ function create_parameter_set(config::AtmosConfig)
         water_params,
         surface_fluxes_params,
         turbconv_params,
+        surface_temp_params,
     )
 end

src/surface_conditions/surface_conditions.jl

@@ -18,7 +18,8 @@ function update_surface_conditions!(Y, p, t)
     (; ᶜts, ᶜu, sfc_conditions) = p.precomputed
     (; params, sfc_setup, atmos) = p
     thermo_params = CAP.thermodynamics_params(params)
-    surface_params = CAP.surface_fluxes_params(params)
+    surface_fluxes_params = CAP.surface_fluxes_params(params)
+    surface_temp_params = CAP.surface_temp_params(params)
     int_ts_values = Fields.field_values(Fields.level(ᶜts, 1))
     int_u_values = Fields.field_values(Fields.level(ᶜu, 1))
     int_z_values =
@@ -36,7 +37,8 @@ function update_surface_conditions!(Y, p, t)
         projected_vector_data(CT2, int_u_values, int_local_geometry_values),
         int_z_values,
         thermo_params,
-        surface_params,
+        surface_fluxes_params,
+        surface_temp_params,
         atmos,
         sfc_temp_var,
         t,
@@ -136,7 +138,8 @@ ifelsenothing(x::Nothing, default) = default
         interior_v,
         interior_z,
         thermo_params,
-        surface_params,
+        surface_fluxes_params,
+        surface_temp_params,
         atmos,
         sfc_prognostic_temp,
         t,
@@ -153,7 +156,8 @@ function surface_state_to_conditions(
     interior_v,
     interior_z,
     thermo_params,
-    surface_params,
+    surface_fluxes_params,
+    surface_temp_params,
     atmos,
     sfc_prognostic_temp,
     t,
@@ -169,7 +173,11 @@ function surface_state_to_conditions(
 
     T = if isnothing(sfc_prognostic_temp)
         if isnothing(surf_state.T)
-            surface_temperature(atmos.sfc_temperature, coordinates)
+            surface_temperature(
+                atmos.sfc_temperature,
+                coordinates,
+                surface_temp_params,
+            )
         else
             surf_state.T
         end
@@ -277,7 +285,7 @@ function surface_state_to_conditions(
             end
             if isnothing(surf_state.gustiness)
                 buoyancy_flux = SF.compute_buoyancy_flux(
-                    surface_params,
+                    surface_fluxes_params,
                     shf,
                     lhf,
                     interior_ts,
@@ -320,7 +328,7 @@ function surface_state_to_conditions(
     end
 
     return atmos_surface_conditions(
-        SF.surface_conditions(surface_params, inputs),
+        SF.surface_conditions(surface_fluxes_params, inputs),
         ts,
         surface_local_geometry,
         atmos,
@@ -332,28 +340,33 @@ end
 function surface_temperature(
     ::RCEMIPIISST,
     coordinates::Union{Geometry.LatLongZPoint, Geometry.LatLongPoint},
+    surface_temp_params,
 )
     (; lat) = coordinates
+    (; SST_mean, SST_delta, SST_wavelength_latitude) = surface_temp_params
     FT = eltype(lat)
-    T = FT(300) + FT(1.25) / 2 * cosd(360 * lat / 54)
+    T = SST_mean + SST_delta / 2 * cosd(360 * lat / SST_wavelength_latitude)
     return T
 end
 
-#Plane SST distribution from Wing et al. (2023) https://gmd.copernicus.org/preprints/gmd-2023-235/
+#Box SST distribution from Wing et al. (2023) https://gmd.copernicus.org/preprints/gmd-2023-235/
 function surface_temperature(
     ::RCEMIPIISST,
     coordinates::Union{Geometry.XZPoint, Geometry.XYZPoint},
+    surface_temp_params,
 )
     (; x) = coordinates
+    (; SST_mean, SST_delta, SST_wavelength) = surface_temp_params
     FT = eltype(x)
-    T = FT(300) + FT(1.25) / 2 * cos(2 * FT(pi) * x / 6000)
+    T = SST_mean - SST_delta / 2 * cos(2 * FT(pi) * x / SST_wavelength)
     return T
 end
 
 #For non-RCEMIPII box models with prescribed surface temp, assume that the latitude is 0.
 function surface_temperature(
     ::Union{ZonallySymmetricSST, ZonallyAsymmetricSST},
     coordinates::Union{Geometry.XZPoint, Geometry.XYZPoint},
+    surface_temp_params,
 )
     (; x) = coordinates
     FT = eltype(x)
@@ -363,6 +376,7 @@ end
 function surface_temperature(
     ::ZonallySymmetricSST,
     coordinates::Geometry.LatLongZPoint,
+    surface_temp_params,
 )
     (; lat, z) = coordinates
     FT = eltype(lat)
@@ -373,6 +387,7 @@ end
 function surface_temperature(
     ::ZonallyAsymmetricSST,
     coordinates::Geometry.LatLongZPoint,
+    surface_temp_params,
 )
     (; lat, long, z) = coordinates
     FT = eltype(lat)

toml/rcemipii_diagnostic_edmfx_0M.toml

@@ -0,0 +1,33 @@
+[precipitation_timescale]
+value = 600
+
+[entr_inv_tau]
+value = 0.002
+
+[entr_coeff]
+value = 0
+
+[detr_inv_tau]
+value = 0
+
+[detr_vertdiv_coeff]
+value = 0.6
+
+[detr_buoy_coeff]
+value = 0.12
+
+[min_area_limiter_scale]
+value = 0
+
+[max_area_limiter_scale]
+value = 0
+
+[angular_velocity_planet_rotation]
+value = 0
+
+[ocean_surface_albedo]
+value = 0.07
+
+[mean_sea_level_pressure]
+value = 101480
+