Merge #961

961: Compartmentalize grid construction into standalone function. r=ilopezgp a=ilopezgp This PR creates a standalone function for the construction of the SCM grid from the namelist. This allows users to construct the grid without having to run the model. Co-authored-by: Ignacio <[email protected]>
CliMA · May 1, 2022 · af1fa3c · af1fa3c
2 parents b739536 + e5ea345
commit af1fa3c
Showing 1 changed file with 50 additions and 45 deletions.
diff --git a/driver/main.jl b/driver/main.jl
@@ -72,51 +72,7 @@ function Simulation1d(namelist)
     cfl_limit = namelist["time_stepping"]["cfl_limit"]
     dt_min = namelist["time_stepping"]["dt_min"]
 
-    truncated_gcm_mesh = TC.parse_namelist(namelist, "grid", "stretch", "flag"; default = false)
-
-    if Cases.get_case(namelist) == Cases.LES_driven_SCM()
-        Δz = get(namelist["grid"], "dz", nothing)
-        nz = get(namelist["grid"], "nz", nothing)
-        @assert isnothing(Δz) ⊻ isnothing(nz) string(
-            "LES_driven_SCM supports nz or Δz, not both.",
-            "The domain height is enforced to be the same as in LES.",
-        )
-
-        les_filename = namelist["meta"]["lesfile"]
-        zmax = NC.Dataset(les_filename, "r") do data
-            Array(TC.get_nc_data(data, "zc"))[end]
-        end
-        nz = isnothing(nz) ? Int(zmax ÷ Δz) : Int(nz)
-        Δz = isnothing(Δz) ? FT(zmax ÷ nz) : FT(Δz)
-    else
-        Δz = FT(namelist["grid"]["dz"])
-        nz = namelist["grid"]["nz"]
-    end
-
-    z₀, z₁ = FT(0), FT(nz * Δz)
-    if truncated_gcm_mesh
-        nzₛ = namelist["grid"]["stretch"]["nz"]
-        Δzₛ_surf = FT(namelist["grid"]["stretch"]["dz_surf"])
-        Δzₛ_top = FT(namelist["grid"]["stretch"]["dz_toa"])
-        zₛ_toa = FT(namelist["grid"]["stretch"]["z_toa"])
-        stretch = CC.Meshes.GeneralizedExponentialStretching(Δzₛ_surf, Δzₛ_top)
-        domain = CC.Domains.IntervalDomain(
-            CC.Geometry.ZPoint{FT}(z₀),
-            CC.Geometry.ZPoint{FT}(zₛ_toa),
-            boundary_tags = (:bottom, :top),
-        )
-        gcm_mesh = CC.Meshes.IntervalMesh(domain, stretch; nelems = nzₛ)
-        mesh = TC.TCMeshFromGCMMesh(gcm_mesh; z_max = z₁)
-    else
-        CC.Meshes.Uniform()
-        domain = CC.Domains.IntervalDomain(
-            CC.Geometry.ZPoint{FT}(z₀),
-            CC.Geometry.ZPoint{FT}(z₁),
-            boundary_tags = (:bottom, :top),
-        )
-        mesh = CC.Meshes.IntervalMesh(domain, nelems = nz)
-    end
-    grid = TC.Grid(mesh)
+    grid = construct_grid(namelist; FT = FT)
 
     Stats = skip_io ? nothing : NetCDFIO_Stats(namelist, grid)
 
@@ -284,6 +240,55 @@ function initialize(sim::Simulation1d)
     return
 end
 
+function construct_grid(namelist; FT = Float64)
+
+    truncated_gcm_mesh = TC.parse_namelist(namelist, "grid", "stretch", "flag"; default = false)
+
+    if Cases.get_case(namelist) == Cases.LES_driven_SCM()
+        Δz = get(namelist["grid"], "dz", nothing)
+        nz = get(namelist["grid"], "nz", nothing)
+        @assert isnothing(Δz) ⊻ isnothing(nz) string(
+            "LES_driven_SCM supports nz or Δz, not both.",
+            "The domain height is enforced to be the same as in LES.",
+        )
+
+        les_filename = namelist["meta"]["lesfile"]
+        zmax = NC.Dataset(les_filename, "r") do data
+            Array(TC.get_nc_data(data, "zc"))[end]
+        end
+        nz = isnothing(nz) ? Int(zmax ÷ Δz) : Int(nz)
+        Δz = isnothing(Δz) ? FT(zmax ÷ nz) : FT(Δz)
+    else
+        Δz = FT(namelist["grid"]["dz"])
+        nz = namelist["grid"]["nz"]
+    end
+
+    z₀, z₁ = FT(0), FT(nz * Δz)
+    if truncated_gcm_mesh
+        nzₛ = namelist["grid"]["stretch"]["nz"]
+        Δzₛ_surf = FT(namelist["grid"]["stretch"]["dz_surf"])
+        Δzₛ_top = FT(namelist["grid"]["stretch"]["dz_toa"])
+        zₛ_toa = FT(namelist["grid"]["stretch"]["z_toa"])
+        stretch = CC.Meshes.GeneralizedExponentialStretching(Δzₛ_surf, Δzₛ_top)
+        domain = CC.Domains.IntervalDomain(
+            CC.Geometry.ZPoint{FT}(z₀),
+            CC.Geometry.ZPoint{FT}(zₛ_toa),
+            boundary_tags = (:bottom, :top),
+        )
+        gcm_mesh = CC.Meshes.IntervalMesh(domain, stretch; nelems = nzₛ)
+        mesh = TC.TCMeshFromGCMMesh(gcm_mesh; z_max = z₁)
+    else
+        CC.Meshes.Uniform()
+        domain = CC.Domains.IntervalDomain(
+            CC.Geometry.ZPoint{FT}(z₀),
+            CC.Geometry.ZPoint{FT}(z₁),
+            boundary_tags = (:bottom, :top),
+        )
+        mesh = CC.Meshes.IntervalMesh(domain, nelems = nz)
+    end
+    return TC.Grid(mesh)
+end
+
 include("callbacks.jl")
 
 function solve_args(sim::Simulation1d)