ft: rcemipii vertical mesh

Author: haakon-e

config/default_configs/default_config.yml

@@ -27,6 +27,9 @@ z_elem:
 z_max:
   help: "Model top height. Default: 30km"
   value: 30000.0
+vertical_mesh:
+  help: "Predefined vertical mesh [`~` (default), `rcemipii`]"
+  value: ~
 vorticity_hyperdiffusion_coefficient:
   help: "Hyperdiffusion coefficient for vorticity (m s-1)"
   value: 0.1857

config/model_configs/rcemipii_box_CRM_1M.yml

@@ -7,30 +7,37 @@ config: box
 rad: allskywithclear
 approximate_linear_solve_iters: 2  # only valid when implicit_diffusion=true
 rayleigh_sponge: true  # set sponge level in toml file
-smagorinsky_lilly: "UV"
+
+# Diffusion
+hyperdiff: "false"
 implicit_diffusion: true
+# smagorinsky_lilly: "UV_W"  # horizontal + vertical
+smagorinsky_lilly: "UV"  # horizontal
+vert_diff: "DecayWithHeightDiffusion"  # vertical
 
 reproducibility_test: true
 
-# Reference diffusion setup
-vert_diff: "DecayWithHeightDiffusion"
-hyperdiff: "false"
-
 # microphysics
+### 0M+Equil option
+# moist: equil
+# precip_model: 0M
 ### !! 1M and 2M !!
 cloud_model: grid_scale   # Ultimately wanted when `moist: nonequil`
 moist: nonequil
 precip_model: 1M
+
 ### spatial discretization ###
-x_max: 96000.0  # 96km
-y_max: 96000.0  # 96km
-# x_max: 6000000.0  # 6000km
-# y_max: 400000.0  # 400km
-z_max: 30000.0  # 30km
 nh_poly: 3  # hor. poly deg -> # quad pts in 1D in a h. elem is Nq = nh + 1
-# z_elem: 43  # Note: Set sponge height (`zd_viscous`) in toml file to cover top ~5 points
-z_elem: 86
-dz_bottom: 30.0
+## RCE_small
+x_max: 96000  # 96km
+y_max: 96000  # 96km
+## RCE_large
+# x_max: 6000000  # 6000km
+# y_max: 400000  # 400km
+## Vertical mesh
+z_max: 33000  # 33km
+z_elem: 74
+vertical_mesh: rcemipii
 
 #~8km grid
 x_elem: 4

src/simulation/grids.jl

@@ -1,4 +1,4 @@
-import ClimaCore: Geometry, Hypsography, Fields, Spaces, Meshes, Grids, CommonGrids
+import ClimaCore: Geometry, Hypsography, Fields, Spaces, Meshes, Grids, CommonGrids, Domains
 using ClimaUtilities: SpaceVaryingInputs.SpaceVaryingInput
 import .AtmosArtifacts as AA
 import ClimaComms
@@ -167,14 +167,27 @@ function BoxGrid(
     topography_damping_factor = 5.0,
     mesh_warp_type::MeshWarpType = LinearWarp(),
     topo_smoothing = false,
+    vertical_mesh = nothing,
 ) where {FT}
     n_quad_points = nh_poly + 1
     stretch =
         z_stretch ? Meshes.HyperbolicTangentStretching{FT}(dz_bottom) : Meshes.Uniform()
     hypsography_fun = hypsography_function_from_topography(
         FT, topography, topography_damping_factor, mesh_warp_type, topo_smoothing,
     )
-    z_mesh = CommonGrids.DefaultZMesh(FT; z_min = 0, z_max, z_elem, stretch)
+    z_mesh = if vertical_mesh == "rcemipii"
+        boundary_layer =
+            FT[0, 37, 112, 194, 288, 395, 520, 667, 843, 1062, 1331, 1664, 2055, 2505]
+        n_bl = length(boundary_layer) - 1
+        free_atmosphere = range(3000, FT(z_max), length = z_elem - n_bl)  # z_elem=74, z_max=33_000m --> 500m spacing
+        CT = Geometry.ZPoint{FT}
+        faces = CT.([boundary_layer; free_atmosphere])
+        z_domain =
+            Domains.IntervalDomain(CT(0), CT(z_max); boundary_names = (:bottom, :top))
+        Meshes.IntervalMesh(z_domain, faces)
+    else
+        CommonGrids.DefaultZMesh(FT; z_min = 0, z_max, z_elem, stretch)
+    end
     grid = CommonGrids.Box3DGrid(
         FT;
         z_elem, x_min = 0, x_max, y_min = 0, y_max, z_min = 0, z_max,

src/solver/type_getters.jl

@@ -764,6 +764,7 @@ function get_grid(parsed_args, params, context)
         z_max = parsed_args["z_max"],
         z_stretch = parsed_args["z_stretch"],
         dz_bottom = parsed_args["dz_bottom"],
+        vertical_mesh = parsed_args["vertical_mesh"],
     )
 
     # Add topography parameters for non-column grids