Merge pull request #3416 from CliMA/as/update-configs

Update configs (default topo smoothing) + add docs

config/default_configs/default_config.yml

@@ -51,7 +51,7 @@ topo_smoothing:
   value: false
 topography_damping_factor:
   help: "Factor by which smallest resolved length-scale is to be damped"
-  value: 512
+  value: 5
 # ODE
 use_newton_rtol:
   help: "Whether to check if the current iteration of Newton's method has an error within a relative tolerance, instead of always taking the maximum number of iterations (only for ClimaTimeSteppers.jl)"

docs/make.jl

@@ -33,6 +33,7 @@ makedocs(;
         "Diagnostic EDMF Equations" => "diagnostic_edmf_equations.md",
         "Gravity Wave Drag Parameterizations" => "gravity_wave.md",
         "Ocean Surface Albedo Parameterization" => "surface_albedo.md",
+        "Topography Representation" => "topography.md",
         "Tracers" => "tracers.md",
         "Radiative Equilibrium" => "radiative_equilibrium.md",
         "Restarts and checkpoints" => "restarts.md",

docs/src/topography.md

@@ -0,0 +1,12 @@
+### Topography in ClimaAtmos
+
+Dataset source : https://www.ncei.noaa.gov/products/etopo-global-relief-model
+ClimaArtifact  : https://github.com/CliMA/ClimaArtifacts/tree/main/earth_orography 
+
+We currently use the ClimaUtilities `SpaceVaryingInput` tool to regrid (using linear interpolation) the ETOPO2022 ice-surface elevation dataset (see ClimaArtifacts) onto the required spectral element horizontal grid. The file `examples/topography_spectra.jl` provides simple tools to generate such regridded fields (and their spectra) on user-defined horizontal spaces. For existing ClimaAtmos simulation data, users may access the `orog.nc` dataset from the default diagnostic outputs to visualize or examine this data. As an example, we include plots of the generated topography on a cubed sphere with 16 elements and 64 elements per panel edge, and compare instances of `unsmoothed` and `smoothed` datasets to visualize the effect of our topography smoothing methods. 
+
+- Elevation data (elems per panel = 16)
+![](assets/smoothing_16elem.png)
+
+- Elevation data (elems per panel = 64)
+![](assets/smoothing_64elem.png)

examples/topography_spectra.jl

@@ -44,7 +44,6 @@ function generate_spaces(;
     Δh_scale = Spaces.node_horizontal_length_scale(h_space)
     @assert h_space isa CC.Spaces.SpectralElementSpace2D
     coords = CC.Fields.coordinate_field(h_space)
-    target_field = CC.Fields.zeros(h_space)
     elev_from_file = SpaceVaryingInputs.SpaceVaryingInput(
         AA.earth_orography_file_path(; context = comms_ctx),
         "z",
@@ -129,7 +128,7 @@ function generate_all_spectra(; h_elem = 16)
         xgridvisible = true,
         ygridvisible = false,
     )
-    for ii in (1, 2, 4, 8, 16, 32, 64, 128, 256, 512)
+    for ii in collect(range(1, 10, length = 10))
         n_attenuation = ii
         test_var = generate_spaces(; h_elem, n_attenuation)
         Δh_scale = Spaces.node_horizontal_length_scale(Spaces.axes(test_var))
@@ -142,3 +141,84 @@ function generate_all_spectra(; h_elem = 16)
     CairoMakie.Legend(fig[1, 2], ax1)
     return fig
 end
+
+"""
+    compare_elevation(;h_elem, n_attenuation)
+Computes the difference between the (regridded) elevation profile
+on a spectral element grid with `h_elem` without any smoothing against
+the corresponding regridded and smoothed surface elevation field. 
+Returns regridded arrays of target topography, and differences
+on cubed-sphere with `h_elem` elements per panel.
+(This field can then be plotted with `ClimaCorePlots`)
+"""
+function compare_elevation(;
+    h_elem = 16,
+    n_attenuation = 5,
+    planet_radius = 6.378e6,
+)
+    FT = Float32
+    space1 = generate_spaces(; h_elem, n_attenuation)
+    space2 = generate_spaces(; h_elem, n_attenuation = 1)
+    Δh_scale = Spaces.node_horizontal_length_scale(Spaces.axes(space2))
+    raw_space_diff = space2 .- space1
+    # Linear Interpolation for visualization
+    npts = Int(round(2π * planet_radius / Δh_scale))
+    npts = ifelse(rem(npts, 2) == 0, npts, npts + 1)
+    longpts = range(FT(-180), FT(180.0), Int(npts))
+    latpts = range(FT(-90), FT(90), Int(npts // 2))
+    hcoords =
+        [Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts]
+    regridded_space_diff = remap_to_array(raw_space_diff, hcoords)
+    regridded_tgt_space = remap_to_array(space1, hcoords)
+    return (
+        raw_space_diff,
+        regridded_space_diff,
+        longpts,
+        latpts,
+        regridded_tgt_space,
+    )
+end
+
+(hspace_diff, regrid_space_diff, lon, lat, regridded_target_space) =
+    compare_elevation(; h_elem = 16, n_attenuation = 5);
+
+"""
+    generate_fig_Δelevation(;h_elem, n_attenuation)
+Given elevation data regridded onto a cubed-sphere grid with `h_elem²` panels per
+face, generate a figure with the following plots: 
+(1) Difference in regridded elevation data (raw cubed-sphere data - smoothed cubed-sphere data)
+(2) Surface elevation as seen by the model (smoothed cubed-sphere data)
+The "regridding" step implied here consists of interpolating (linear interoplation by default)the ETOPO2022 
+raw dataset (via ClimaArtifacts) onto the required ClimaAtmos grid.
+"""
+function generate_fig_Δelevation(; h_elem = 16, n_attenuation = 5)
+    fig = Figure(; size = (2000, 1000))
+    ax1 = Axis(
+        fig[1, 1],
+        xlabel = "lon",
+        ylabel = "lat",
+        title = "Difference in regridded surface elevation (raw - smoothed)",
+    )
+    ax2 = Axis(
+        fig[2, 1],
+        xlabel = "lon",
+        ylabel = "lat",
+        title = "Regridded surface elevation (smoothed)",
+    )
+    color_levels = 25
+    cmap = cgrad(:curl, color_levels; categorical = true)
+    raw_Δelevation, regridded_Δelevation, lon, lat, regridded_tgt =
+        compare_elevation(; h_elem, n_attenuation)
+    figdata1 = CairoMakie.contourf!(
+        ax1,
+        lon,
+        lat,
+        regridded_Δelevation;
+        colormap = cmap,
+    )
+    figdata2 =
+        CairoMakie.contourf!(ax2, lon, lat, regridded_tgt; colormap = cmap)
+    Colorbar(fig[1, 2], figdata1, label = "Elevation diff [m]")
+    Colorbar(fig[2, 2], figdata2, label = "Regridded elevation [m]")
+    return fig
+end