Account for ERA5/ClimaAtmos topo differences

src/diagnostics/core_diagnostics.jl

@@ -1660,11 +1660,26 @@ function compute_mslp!(out, state, cache, time)
     )
     z_level = Fields.level(Fields.coordinate_field(state.c.ρ).z, 1)
 
-    # compute sea level pressure using the hypsometric equation
+    # Reduce to mean sea level using a standard-lapse hypsometric formulation (ERA-style)
+    # Using constant lapse rate Γ = 6.5 K/km, with virtual temperature
+    # represented via R_m_surf. This more closely matches ERA5's msl computation and
+    # reduces biases over very cold or very warm high-topography regions.
+    FT = Spaces.undertype(Fields.axes(state.c.ρ))
+    Γ = FT(6.5e-3) # K m^-1
+    ϵ = FT(1e-6)   # small floor to avoid non-positive base
+
+    # bracket = max(ϵ, 1 + Γ z / T(z))
+    bracket = similar(t_level)
+    @. bracket = max(ϵ, 1 + Γ * z_level / t_level)
+
+    # exponent = g / (R_m Γ)
+    exponent = similar(R_m_surf)
+    @. exponent = g / (R_m_surf * Γ)
+
     if isnothing(out)
-        return @. p_level * exp(g * z_level / (R_m_surf * t_level))
+        return p_level .* (bracket .^ exponent)
     else
-        @. out = p_level * exp(g * z_level / (R_m_surf * t_level))
+        out .= p_level .* (bracket .^ exponent)
     end
 end
 

src/initial_conditions/initial_conditions.jl

@@ -371,6 +371,136 @@ function overwrite_initial_conditions!(
     )
 end
 
+"""
+    correct_surface_pressure_for_topography!(
+        p_sfc,
+        file_path,
+        face_space,
+        Y,
+        ᶜT,
+        ᶜq_tot,
+        thermo_params,
+        regridder_kwargs,
+    )
+
+Adjusts the surface pressure field `p_sfc` to account for mismatches between
+ERA5 (file) surface altitude and the model orography when specifying pressure.
+
+    Δz = z_model_surface - z_sfc
+
+and applies a hydrostatic correction at the surface using the local moist gas
+constant and temperature at the surface:
+
+    p_sfc .= p_sfc .* exp.(-Δz * g ./ (R_m_sfc .* T_sfc))
+
+where:
+- `g` is gravitational acceleration from `thermo_params`
+- `R_m_sfc` is the moist-air gas constant evaluated from `ᶜq_tot` at the surface
+- `T_sfc` is the air temperature from `ᶜT` at the surface
+
+A small floor is applied to the denominator `R_m_sfc .* T_sfc` to avoid
+numerical blow-ups, and non-finite `Δz` entries leave `p_sfc` unchanged.
+
+Returns `true` if the correction is applied. If `z_sfc` is missing in the
+provided file, the function logs a warning and returns `false` without
+modifying `p_sfc`.
+
+Arguments
+- `p_sfc`: face field of surface pressure to be corrected (modified in-place)
+- `file_path`: path to the ERA5-derived initialization NetCDF file
+- `face_space`: face space of the model grid (for reading/regridding)
+- `Y`: prognostic state, used to obtain model surface height
+- `ᶜT`: center field of temperature
+- `ᶜq_tot`: center field of total specific humidity
+- `thermo_params`: thermodynamics parameter set
+- `regridder_kwargs`: keyword arguments forwarded to the regridder
+
+Notes
+- Expects `z_sfc` units in meters; if the variable is not present, the
+  correction is skipped.
+"""
+function correct_surface_pressure_for_topography!(
+    p_sfc,
+    file_path,
+    face_space,
+    Y,
+    ᶜT,
+    ᶜq_tot,
+    thermo_params,
+    regridder_kwargs,
+)
+    var_name = "z_sfc"
+    ᶠz_surface = nothing
+    # Explicitly check if `z_sfc` exists in the input file; skip if not found
+    has_z_sfc = NC.NCDataset(file_path) do ds
+        haskey(ds, var_name)
+    end
+    if !has_z_sfc
+        @warn "Skipping topographic correction because variable `$var_name` is missing from $(file_path)."
+        return false
+    end
+    ᶠz_surface = Fields.level(
+        SpaceVaryingInputs.SpaceVaryingInput(
+            file_path,
+            var_name,
+            face_space,
+            regridder_kwargs = regridder_kwargs,
+        ),
+        Fields.half,
+    )
+
+    if ᶠz_surface === nothing
+        return false
+    end
+
+    FT = eltype(thermo_params)
+    grav = thermo_params.grav
+
+    # Log initial surface pressure statistics
+    p_sfc_arr = Array(Fields.field2array(p_sfc))
+    p_sfc_finite = filter(isfinite, p_sfc_arr)
+
+    ᶠz_model_surface = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
+    ᶠΔz = zeros(face_space)
+    @. ᶠΔz = ᶠz_model_surface - ᶠz_surface
+
+    # Log surface height statistics
+    z_surface_arr = Array(Fields.field2array(ᶠz_surface))
+    z_model_arr = Array(Fields.field2array(ᶠz_model_surface))
+    Δz_arr = Array(Fields.field2array(ᶠΔz))
+
+    ᶜphase_partition = TD.PhasePartition.(ᶜq_tot)
+    ᶜR_m = TD.gas_constant_air.(thermo_params, ᶜphase_partition)
+    ᶠR_m = ᶠinterp.(ᶜR_m)
+    ᶠR_m_sfc = Fields.level(ᶠR_m, Fields.half)
+
+    ᶠT = ᶠinterp.(ᶜT)
+    ᶠT_sfc = Fields.level(ᶠT, Fields.half)
+
+    denom = copy(ᶠR_m_sfc)
+    @. denom = ᶠR_m_sfc * ᶠT_sfc
+    denom_floor = FT(10)
+    @. denom = ifelse(isfinite(denom) && denom > denom_floor, denom, denom_floor)
+
+    correction_factor = similar(ᶠΔz)
+    @. correction_factor = ifelse(
+        isfinite(ᶠΔz),
+        exp(FT(-1) * ᶠΔz * grav / denom),
+        one(FT),
+    )
+
+    @. p_sfc =
+        p_sfc *
+        ifelse(
+            isfinite(ᶠΔz),
+            exp(FT(-1) * ᶠΔz * grav / denom),
+            one(FT),
+        )
+
+    @info "Adjusted surface pressure to account for ERA5/model surface-height differences."
+    return true
+end
+
 # WeatherModel function using the shared implementation
 function overwrite_initial_conditions!(
     initial_condition::WeatherModel,
@@ -421,6 +551,16 @@ function overwrite_initial_conditions!(
         center_space,
         regridder_kwargs = regridder_kwargs,
     )
+    correct_surface_pressure_for_topography!(
+        p_sfc,
+        file_path,
+        face_space,
+        Y,
+        ᶜT,
+        ᶜq_tot,
+        thermo_params,
+        regridder_kwargs,
+    )
 
     # With the known temperature (ᶜT) and moisture (ᶜq_tot) profile,
     # recompute the pressure levels assuming hydrostatic balance is maintained.

src/utils/weather_model.jl

@@ -159,12 +159,28 @@ function to_z_levels(source_file, target_file, target_levels, FT)
 
     # Write 2D surface variables - extend to all levels (TODO: accept 2D variables in atmos)
     # Simply repeat the surface values for all levels
-    surf_map = Dict("skt" => "skt", "sp" => "p")
+    surf_map = Dict("skt" => "skt", "sp" => "p", "surface_geopotential" => "z_sfc")
     for (src_name, dst_name) in surf_map
-        var_obj =
-            defVar(ncout, dst_name, FT, ("lon", "lat", "z"), attrib = ncin[src_name].attrib)
+        # Choose attributes; for z_sfc, set clean altitude attributes
+        var_attrib = if dst_name == "z_sfc"
+            Dict(
+                "standard_name" => "surface_altitude",
+                "long_name" => "surface altitude derived from ERA5",
+                "units" => "m",
+                "source_variable" => src_name,
+            )
+        else
+            ncin[src_name].attrib
+        end
+        var_obj = defVar(ncout, dst_name, FT, ("lon", "lat", "z"), attrib = var_attrib)
+        # Read first time slice and coalesce; follow same convention as sp (use [:, :, 1])
+        data2d = FT.(coalesce.(ncin[src_name][:, :, 1], NaN))
+        # Convert geopotential to meters if necessary
+        if dst_name == "z_sfc"
+            data2d .= data2d ./ grav
+        end
         for k in 1:length(target_levels)
-            var_obj[:, :, k] = FT.(ncin[src_name][:, :, 1])
+            var_obj[:, :, k] = data2d
         end
     end