Add support for interpolation of topography

[skip ci][ci skip]

config/default_configs/default_config.yml

@@ -249,3 +249,6 @@ log_params:
 output_default_diagnostics:
   help: "Output the default diagnostics associated to the selected atmospheric model"
   value: false
+netcdf_interpolate_topography:
+  help: "Interpolate diagnostics in the NetCDF files in such a way that `z` is the elevation from the sea level (as opposed to the elevation from the surface)"
+  value: false

src/diagnostics/netcdf_writer.jl

@@ -414,6 +414,9 @@ struct NetCDFWriter{T, TS}
 
     # NetCDF files that are currently open. Only the root process uses this field.
     open_files::Dict{String, NCDatasets.NCDataset}
+
+    # Whether to treat z as altitude over the surface or the sea level
+    interpolate_topography::Bool
 end
 
 """
@@ -458,9 +461,9 @@ function NetCDFWriter(;
     # common paradigm for all the other dimensions. Moreover, with topography, we need need
     # to perform an interpolation for the surface.
 
-    if hypsography isa Spaces.Flat
-        # We don't have information about the type and space, so we set interpolated_surface
-        # to nothing, and deal with it later.
+    # We have to deal with the surface only if we are not interpolating the topography and
+    # if our topography is non-trivial
+    if hypsography isa Spaces.Flat || interpolate_topography
         interpolated_surface = nothing
     elseif hypsography isa Hypsography.LinearAdaption
         horizontal_space = axes(hypsography.surface)
@@ -472,7 +475,11 @@ function NetCDFWriter(;
         )
         vcoords = []
         remapper = Remapper(hcoords, vcoords, horizontal_space)
-        interpolated_surface = interpolate(remapper, hypsography.surface)
+        interpolated_surface = interpolate(
+            remapper,
+            hypsography.surface,
+            interpolate_topography = false,
+        )
     else
         error("Cannot process hysography $hypsography")
     end
@@ -483,6 +490,7 @@ function NetCDFWriter(;
         compression_level,
         interpolated_surface,
         Dict(),
+        interpolate_topography,
     )
 end
 
@@ -532,7 +540,8 @@ function write_field!(writer::NetCDFWriter, field, diagnostic, integrator)
 
     # Now we can interpolate onto the target points
     # There's an MPI call in here (to aggregate the results)
-    interpolated_field = interpolate(remapper, field)
+    interpolated_field =
+        interpolate(remapper, field; writer.interpolate_topography)
 
     # Only the root process has to write
     ClimaComms.iamroot(ClimaComms.context(field)) || return

src/solver/type_getters.jl

@@ -632,7 +632,10 @@ function get_diagnostics(parsed_args, atmos_model, hypsography)
     )
 
     hdf5_writer = CAD.HDF5Writer()
-    netcdf_writer = CAD.NetCDFWriter(; hypsography)
+    netcdf_writer = CAD.NetCDFWriter(;
+        hypsography,
+        interpolate_topography = parsed_args["netcdf_interpolate_topography"],
+    )
     writers = (hdf5_writer, netcdf_writer)
 
     # The default writer is HDF5