Account for ERA5/ClimaAtmos topo differences #4098
Conversation
…re. Account for lapse rate in mslp sea level reduction (diagnostics).
costachris
force-pushed
the
cc/improve_wxmodel_pressure_init_diag
branch
from
6f20da6 to
1ce4418
Compare
| @assert all(map(x -> x in (keys(ncin)), in_dims)) "Source file $source_file is missing subset of the required dimensions: $in_dims" | ||
|
|
||
| # assert ncin has required variables | ||
| req_vars = ["u", "v", "w", "t", "q", "skt", "sp"] |
There was a problem hiding this comment.
If it breaks without surface_geopotential catch it with the @assert here
| # Choose attributes; for z_sfc, set clean altitude attributes | ||
| var_attrib = if dst_name == "z_sfc" | ||
| Dict( | ||
| "standard_name" => "surface_altitude", |
There was a problem hiding this comment.
I assume the data you're writing is geopotential based on the ERA5 name? If that is the case the units are gravity * height (m^2/s^2) and would be good to at least rename the long_name and the units. I would also use zg_sfc even though ERA5 uses z for geopotential because then it's very clear.
There was a problem hiding this comment.
Ah ok I see what is happening now, you're correct. Not sure I have an immediate suggestion other that it is a little confusing because there are two if statements in this loop to catch z_sfc since it needs a new var_attrib and needs to be divided by gravity. Possibly separating it out could be cleaner / easier to read if a little longer.
| # Convert geopotential to meters if necessary | ||
| if dst_name == "z_sfc" | ||
| data2d .= data2d ./ grav | ||
| end |
There was a problem hiding this comment.
| end | |
| end | |
| # Duplicate 2D surface field across all target vertical levels |
I think this comment would be helpful here
| 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)) |
There was a problem hiding this comment.
Do we need a coalesce? Are there instances when the data is NaN and our model will still run? If not then we could consider asserting not NaN
| if !has_z_sfc | ||
| @warn "Skipping topographic correction because variable `$var_name` is missing from $(file_path)." | ||
| return false | ||
| end |
There was a problem hiding this comment.
I think this haskey check is maybe better done in overwrite_initial_conditions!? Perhaps also passing the variable name z_sfc as a function argument would make it more generalizable
| z_surface_arr = Array(Fields.field2array(ᶠz_surface)) | ||
| z_model_arr = Array(Fields.field2array(ᶠz_model_surface)) | ||
| Δz_arr = Array(Fields.field2array(ᶠΔz)) |
There was a problem hiding this comment.
Shouldn't field2array already return an array - i.e., does Array() actually do anything here?
| 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), | ||
| ) |
There was a problem hiding this comment.
Repeated code. either delete the first bit or simplify the second
| 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) |
There was a problem hiding this comment.
This is written confusingly. Could you either comment to explain why we need so many steps (e.g. why is floor and finite checks needed?) Why is 10 chosen as a floor?
| p_sfc_arr = Array(Fields.field2array(p_sfc)) | ||
| p_sfc_finite = filter(isfinite, p_sfc_arr) |
There was a problem hiding this comment.
Again is Array needed after field2array? Is isfinite needed if its surface pressure - shouldn't it be finite? I also think one line and a clearer comment would help
| long_name = "Mean Sea Level Pressure", | ||
| standard_name = "mean_sea_level_pressure", | ||
| units = "Pa", | ||
| comments = "Mean sea level pressure computed from the hypsometric equation", |
There was a problem hiding this comment.
Maybe update the comment to note that this formulation is similar to ERA5 and uses a lapse-rate temperature dependent hypsometric calculation instead of just the standard one
| 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 * Γ) |
There was a problem hiding this comment.
This will allocate new arrays which is slow. Maybe @dennisYatunin could recommend a fix?
| 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) |
There was a problem hiding this comment.
Could you write out the whole equation here for clarity? Then it matters less how you break up the computation for performance. e.g., "p_msl = p(z1)(1+\Gamma z1 / t)^(g/(Rm*\Gamma) is the ERA5 standard calculation for MSLP...".
Julians42
left a comment
Julians42
left a comment
There was a problem hiding this comment.
Left some comments to improve clarity and then one performance optimization we should do in the diagnostics for Dennis. Glad this is looking so much better!
Adjust for ERA5/ClimaAtmos topo differences when initializing pressure. Account for lapse rate in mslp sea level reduction (diagnostics).
The colorbars vary, but note each improvement significantly reduces the RMSE/Bias.
Current 😭
Topography correction in surface pressure init 👀
Topography correction + account for lapse rate in MSLP diagnostic 🥳