Merge branch 'main' into gb/update_nc

.buildkite/longruns/pipeline.yml

@@ -140,6 +140,18 @@ steps:
           slurm_time: 24:00:00
         env:
           JOB_NAME: "longrun_aquaplanet_rhoe_equilmoist_nz63_0M_55km_rs35km_clearsky_tvinsolation"
+
+      - label: ":computer: aquaplanet (ρe_tot) equilmoist clearsky high resolution (nz63) hightop (55km) rayleigh sponge(35e3, 10) Float32 (time-varying insolation) + earth topography"
+        command:
+          - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
+          - julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl --data_dir $$JOB_NAME --out_dir $$JOB_NAME
+          - julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl --nc_dir $$JOB_NAME --fig_dir $$JOB_NAME --case_name aquaplanet
+        artifact_paths: "$$JOB_NAME/*"
+        agents:
+          slurm_ntasks: 64
+          slurm_time: 24:00:00
+        env:
+          JOB_NAME: "longrun_aquaplanet_rhoe_equilmoist_nz63_0M_55km_rs35km_clearsky_tvinsolation_earth"
 
       - label: ":computer: baroclinic wave (ρe_tot) equilmoist high resolution topography (earth)"
         command:

Project.toml

@@ -1,7 +1,7 @@
 name = "ClimaAtmos"
 uuid = "b2c96348-7fb7-4fe0-8da9-78d88439e717"
 authors = ["Climate Modeling Alliance"]
-version = "0.16.2"
+version = "0.16.3"
 
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
@@ -46,13 +46,15 @@ YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
 [compat]
 ArgParse = "1"
 ArtifactWrappers = "0.2"
+Artifacts = "1"
 AtmosphericProfilesLibrary = "0.1"
 CLIMAParameters = "0.7.24"
 ClimaComms = "0.5.5"
 ClimaCore = "0.10.55"
 ClimaTimeSteppers = "0.7.14"
 CloudMicrophysics = "0.15.0"
 Colors = "0.12"
+Dates = "1"
 Dierckx = "0.5"
 DiffEqBase = "6"
 DiffEqCallbacks = "2"
@@ -63,15 +65,21 @@ Insolation = "0.6"
 Interpolations = "0.14"
 IntervalSets = "0.5, 0.6, 0.7"
 Krylov = "0.9"
+LinearAlgebra = "1"
+Logging = "1"
 NCDatasets = "0.11, 0.12, <0.13"
 NVTX = "0.3"
 Pkg = "1.8"
+Printf = "1"
 RRTMGP = "0.9.2"
+Random = "1"
 RootSolvers = "0.2, 0.3, 0.4"
 SciMLBase = "1, 2"
 StaticArrays = "1"
+Statistics = "1"
 StatsBase = "0.33, 0.34"
 SurfaceFluxes = "0.8.0"
+Test = "1"
 Thermodynamics = "0.11"
 YAML = "0.4"
 julia = "1.7"

config/longrun_configs/longrun_aquaplanet_rhoe_equilmoist_nz63_0M_55km_rs35km_clearsky_tvinsolation_earth.yml

@@ -0,0 +1,22 @@
+FLOAT_TYPE: "Float32"
+dt_save_to_disk: "10days"
+dt: "40secs"
+t_end: "200days"
+h_elem: 16
+z_elem: 63
+dz_bottom: 30.0
+dz_top: 3000.0
+z_max: 55000.0
+rayleigh_sponge: true
+topography: "Earth"
+topo_smoothing: true
+kappa_4: 4.0e16
+moist: "equil"
+precip_model: "0M"
+surface_setup: "DefaultMoninObukhov"
+vert_diff: "true"
+rad: "clearsky"
+idealized_insolation: false
+dt_rad: "1hours"
+job_id: "longrun_aquaplanet_rhoe_equilmoist_nz63_0M_55km_rs35km_clearsky_tvinsolation_earth"
+toml: [toml/longrun_aquaplanet_rhoe_equilmoist_nz63_0M_55km_rs35km_clearsky_earth.toml]

config/model_configs/prognostic_edmfx_bomex_box.yml

@@ -7,7 +7,7 @@ surface_setup: "Bomex"
 turbconv: "prognostic_edmfx"
 edmfx_upwinding: first_order
 edmfx_entr_model: "Generalized"
-edmfx_detr_model: "Generalized"
+edmfx_detr_model: "ConstantArea"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
@@ -27,7 +27,7 @@ perturb_initstate: false
 dt: "5secs"
 t_end: "6hours"
 dt_save_to_disk: "10mins"
-toml: [toml/prognostic_edmfx_box.toml]
+toml: [toml/prognostic_edmfx_bomex_box.toml]
 diagnostics:
   - short_name: ts
     period: 10mins

config/model_configs/prognostic_edmfx_bomex_column.yml

@@ -7,7 +7,7 @@ surface_setup: "Bomex"
 turbconv: "prognostic_edmfx"
 edmfx_upwinding: first_order
 edmfx_entr_model: "Generalized"
-edmfx_detr_model: "Generalized"
+edmfx_detr_model: "ConstantArea"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
@@ -22,4 +22,4 @@ perturb_initstate: false
 dt: "5secs"
 t_end: "6hours"
 dt_save_to_disk: "10mins"
-toml: [toml/prognostic_edmfx_box.toml]
+toml: [toml/prognostic_edmfx_bomex_box.toml]

config/model_configs/prognostic_edmfx_bomex_tke_box.yml

@@ -7,7 +7,7 @@ surface_setup: "Bomex"
 turbconv: "prognostic_edmfx"
 edmfx_upwinding: first_order
 edmfx_entr_model: "Generalized"
-edmfx_detr_model: "Generalized"
+edmfx_detr_model: "ConstantArea"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
@@ -27,4 +27,4 @@ perturb_initstate: false
 dt: "5secs"
 t_end: "6hours"
 dt_save_to_disk: "10mins"
-toml: [toml/prognostic_edmfx_box.toml]
+toml: [toml/prognostic_edmfx_bomex_box.toml]

config/model_configs/prognostic_edmfx_bomex_tke_column.yml

@@ -7,7 +7,7 @@ surface_setup: "Bomex"
 turbconv: "prognostic_edmfx"
 edmfx_upwinding: first_order
 edmfx_entr_model: "Generalized"
-edmfx_detr_model: "Generalized"
+edmfx_detr_model: "ConstantArea"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
@@ -22,4 +22,4 @@ perturb_initstate: false
 dt: "5secs"
 t_end: "6hours"
 dt_save_to_disk: "10mins"
-toml: [toml/prognostic_edmfx_box.toml]
+toml: [toml/prognostic_edmfx_bomex_box.toml]

src/cache/diagnostic_edmf_precomputed_quantities.jl

@@ -392,6 +392,8 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
                 ),
                 TD.relative_humidity(thermo_params, ts_prev_level),
                 FT(0),
+                FT(0), # ᶜentr, not implemented
+                FT(0), # ᶜvert_div, not implemented
                 dt,
                 p.atmos.edmfx_detr_model,
             )

src/cache/precomputed_quantities.jl

@@ -44,6 +44,8 @@ function precomputed_quantities(Y, atmos)
         !(atmos.moisture_model isa DryModel) &&
         atmos.energy_form isa TotalEnergy
     ) || !(atmos.turbconv_model isa PrognosticEDMFX)
+    @assert !(atmos.edmfx_detr_model isa ConstantAreaDetrainment) ||
+            !(atmos.turbconv_model isa DiagnosticEDMFX)
     TST = thermo_state_type(atmos.moisture_model, FT)
     SCT = SurfaceConditions.surface_conditions_type(atmos, FT)
     n = n_mass_flux_subdomains(atmos.turbconv_model)
@@ -344,8 +346,8 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
     end
 
     if turbconv_model isa PrognosticEDMFX
-        set_prognostic_edmf_precomputed_quantities_environment!(Y, p, ᶠuₕ³, t)
         set_prognostic_edmf_precomputed_quantities_draft_and_bc!(Y, p, ᶠuₕ³, t)
+        set_prognostic_edmf_precomputed_quantities_environment!(Y, p, ᶠuₕ³, t)
         set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
     end
 

src/cache/prognostic_edmf_precomputed_quantities.jl

@@ -148,6 +148,7 @@ function set_prognostic_edmf_precomputed_quantities_draft_and_bc!(Y, p, ᶠuₕ
             Fields.field_values(Fields.level(Y.c.sgsʲs.:($j).ρa, 1))
 
         turbconv_params = CAP.turbconv_params(params)
+        @. sgsʲs_ρ_int_val = TD.air_density(thermo_params, ᶜtsʲ_int_val)
         @. sgsʲs_ρa_int_val =
             $(FT(turbconv_params.surface_area)) *
             TD.air_density(thermo_params, ᶜtsʲ_int_val)
@@ -191,6 +192,7 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
     ᶜdz = Fields.Δz_field(axes(Y.c))
     ᶜlg = Fields.local_geometry_field(Y.c)
 
+    ᶜvert_div = p.scratch.ᶜtemp_scalar
     for j in 1:n
         @. ᶜentrʲs.:($$j) = entrainment(
             params,
@@ -209,6 +211,7 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             dt,
             p.atmos.edmfx_entr_model,
         )
+        @. ᶜvert_div = ᶜdivᵥ(ᶠinterp(ᶜρʲs.:($$j)) * ᶠu³ʲs.:($$j)) / ᶜρʲs.:($$j)
         @. ᶜdetrʲs.:($$j) = detrainment(
             params,
             ᶜz,
@@ -223,6 +226,8 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             get_physical_w(ᶜu, ᶜlg),
             TD.relative_humidity(thermo_params, ᶜts⁰),
             FT(0),
+            ᶜentrʲs.:($$j),
+            ᶜvert_div,
             dt,
             p.atmos.edmfx_detr_model,
         )

src/prognostic_equations/edmfx_entr_detr.jl

@@ -184,6 +184,8 @@ function detrainment(
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
     dt::FT,
+    ᶜentr::FT,
+    ᶜvert_div::FT,
     ::NoDetrainment,
 ) where {FT}
     return FT(0)
@@ -203,6 +205,8 @@ function detrainment(
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
+    ᶜentr::FT,
+    ᶜvert_div::FT,
     dt::FT,
     ::PiGroupsDetrainment,
 ) where {FT}
@@ -255,6 +259,8 @@ function detrainment(
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
+    ᶜentr::FT,
+    ᶜvert_div::FT,
     dt::FT,
     ::GeneralizedDetrainment,
 ) where {FT}
@@ -294,6 +300,8 @@ function detrainment(
     ᶜw⁰::FT,
     ᶜRH⁰::FT,
     ᶜbuoy⁰::FT,
+    ᶜentr::FT,
+    ᶜvert_div::FT,
     dt::FT,
     ::GeneralizedHarmonicsDetrainment,
 ) where {FT}
@@ -319,6 +327,29 @@ function detrainment(
     return detr * FT(2) * hm_limiter(ᶜaʲ)
 end
 
+function detrainment(
+    params,
+    ᶜz::FT,
+    z_sfc::FT,
+    ᶜp::FT,
+    ᶜρ::FT,
+    buoy_flux_surface::FT,
+    ᶜaʲ::FT,
+    ᶜwʲ::FT,
+    ᶜRHʲ::FT,
+    ᶜbuoyʲ::FT,
+    ᶜw⁰::FT,
+    ᶜRH⁰::FT,
+    ᶜbuoy⁰::FT,
+    ᶜentr::FT,
+    ᶜvert_div::FT,
+    dt::FT,
+    ::ConstantAreaDetrainment,
+) where {FT}
+    detr = max(min(ᶜentr - ᶜvert_div, 1 / dt), 0)
+    return detr
+end
+
 edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model) = nothing
 
 function edmfx_entr_detr_tendency!(

src/prognostic_equations/edmfx_tke.jl

@@ -22,6 +22,7 @@ function edmfx_tke_tendency!(
     (; ᶠu³⁰, ᶜstrain_rate_norm, ᶜlinear_buoygrad, ᶜtke⁰, ᶜmixing_length) =
         p.precomputed
     (; ᶜK_u, ᶜK_h, ρatke_flux) = p.precomputed
+    (; dt) = p.simulation
     ᶜρa⁰ = turbconv_model isa PrognosticEDMFX ? p.precomputed.ᶜρa⁰ : Y.c.ρ
     ᶠgradᵥ = Operators.GradientC2F()
 
@@ -55,8 +56,10 @@ function edmfx_tke_tendency!(
         # pressure work
         # dissipation
         @. Yₜ.c.sgs⁰.ρatke[colidx] -=
-            ᶜρa⁰[colidx] * c_d * max(ᶜtke⁰[colidx], 0)^(FT(3) / 2) /
-            ᶜmixing_length[colidx]
+            ᶜρa⁰[colidx] * c_d * max(ᶜtke⁰[colidx], 0)^(FT(3) / 2) / max(
+                ᶜmixing_length[colidx],
+                c_d * dt * sqrt(max(ᶜtke⁰[colidx], eps(FT))),
+            )
     end
 
     return nothing

src/solver/model_getters.jl

@@ -376,6 +376,8 @@ function get_detrainment_model(parsed_args)
         GeneralizedDetrainment()
     elseif detr_model == "GeneralizedHarmonics"
         GeneralizedHarmonicsDetrainment()
+    elseif detr_model == "ConstantArea"
+        ConstantAreaDetrainment()
     else
         error("Invalid detr_model $(detr_model)")
     end

src/solver/types.jl

@@ -188,6 +188,7 @@ struct NoDetrainment <: AbstractDetrainmentModel end
 struct PiGroupsDetrainment <: AbstractDetrainmentModel end
 struct GeneralizedDetrainment <: AbstractDetrainmentModel end
 struct GeneralizedHarmonicsDetrainment <: AbstractDetrainmentModel end
+struct ConstantAreaDetrainment <: AbstractDetrainmentModel end
 
 abstract type AbstractQuadratureType end
 struct LogNormalQuad <: AbstractQuadratureType end

toml/longrun_aquaplanet_rhoe_equilmoist_nz63_0M_55km_rs35km_clearsky_earth.toml

@@ -0,0 +1,10 @@
+[alpha_rayleigh_w]
+alias = "alpha_rayleigh_w"
+value = 10.0
+type = "float"
+
+[zd_rayleigh]
+alias = "zd_rayleigh"
+value = 35000.0
+type = "float"
+

toml/prognostic_edmfx_bomex_box.toml

@@ -0,0 +1,49 @@
+[EDMF_surface_area]
+alias = "surface_area"
+value = 0.1
+type = "float"
+
+[EDMF_min_area]
+alias = "min_area"
+value = 1.0e-5
+type = "float"
+
+[EDMF_max_area]
+alias = "max_area"
+value = 0.6
+type = "float"
+
+[entr_tau]
+alias = "entr_tau"
+value = 0.002
+type = "float"
+
+[entr_coeff]
+alias = "entr_coeff"
+value = 0
+type = "float"
+
+[min_area_limiter_scale]
+alias = "min_area_limiter_scale"
+value = 0
+type = "float"
+
+[detr_tau]
+alias = "detr_tau"
+value = 0
+type = "float"
+
+[detr_coeff]
+alias = "detr_coeff"
+value = 0
+type = "float"
+
+[detr_buoy_coeff]
+alias = "detr_buoy_coeff"
+value = 0
+type = "float"
+
+[max_area_limiter_scale]
+alias = "max_area_limiter_scale"
+value = 0
+type = "float"