wip

config/longrun_configs/longrun_aquaplanet_rhoe_equil_55km_nz63_clearsky_progedmf_0M.yml

@@ -0,0 +1,36 @@
+h_elem: 16
+z_max: 55000.0
+z_elem: 63
+dz_bottom: 30.0
+dz_top: 3000.0
+moist: "equil"
+precip_model: "0M"
+override_τ_precip: false
+rad: "clearsky"
+dt_rad: "6hours"
+surface_setup: "DefaultMoninObukhov"
+turbconv: "prognostic_edmfx"
+implicit_diffusion: true
+approximate_linear_solve_iters: 2
+max_newton_iters_ode: 3
+prognostic_tke: true
+edmfx_upwinding: "first_order"
+edmfx_entr_model: "Generalized"
+edmfx_detr_model: "Generalized"
+edmfx_nh_pressure: true
+edmfx_velocity_relaxation: true
+edmfx_sgs_mass_flux: false
+edmfx_sgs_diffusive_flux: true
+rayleigh_sponge: true
+dt_save_state_to_disk: "3hours"
+dt: "50secs"
+t_end: "2days"
+job_id: "longrun_aquaplanet_rhoe_equil_55km_nz63_clearsky_progedmf_0M"
+toml: [toml/longrun_aquaplanet_rhoe_equil_55km_nz63_clearsky_progedmf_0M.toml]
+diagnostics:
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+    period: 3hours
+  - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
+    period: 3hours
+  - short_name: [entr, detr, lmix, bgrad, strain, edt, evu]
+    period: 3hours

config/model_configs/prognostic_edmfx_bomex_column.yml

@@ -5,25 +5,26 @@ edmf_coriolis: "Bomex"
 ls_adv: "Bomex"
 surface_setup: "Bomex"
 turbconv: "prognostic_edmfx"
-implicit_diffusion: true
-implicit_sgs_advection: true
+implicit_diffusion: false
+implicit_sgs_advection: false
 approximate_linear_solve_iters: 2
 max_newton_iters_ode: 3
 edmfx_upwinding: first_order
-edmfx_entr_model: "Generalized"
-edmfx_detr_model: "Generalized"
+edmfx_entr_model: "PiGroups"
+edmfx_detr_model: "PiGroups"
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 edmfx_nh_pressure: true
 edmfx_velocity_relaxation: true
 prognostic_tke: true
+cloud_model: "grid_scale"
 moist: "equil"
 config: "column"
 z_max: 3e3
 z_elem: 60
 z_stretch: false
 perturb_initstate: false
-dt: "50secs"
+dt: "1secs"
 t_end: "6hours"
 dt_save_to_disk: "10mins"
 toml: [toml/prognostic_edmfx_bomex.toml]
@@ -34,5 +35,5 @@ diagnostics:
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
     period: 10mins
-  - short_name: [entr, detr, lmix, bgrad, strain, edt, evu]
+  - short_name: [entr, detr, lmix, bgrad, strain]
     period: 10mins

config/model_configs/prognostic_edmfx_bomex_explicit_column.yml

@@ -0,0 +1,39 @@
+job_id: "prognostic_edmfx_bomex_column"
+initial_condition: "Bomex"
+subsidence: "Bomex"
+edmf_coriolis: "Bomex"
+ls_adv: "Bomex"
+surface_setup: "Bomex"
+turbconv: "prognostic_edmfx"
+implicit_diffusion: false
+implicit_sgs_advection: false
+approximate_linear_solve_iters: 2
+max_newton_iters_ode: 3
+edmfx_upwinding: first_order
+edmfx_entr_model: "PiGroups"
+edmfx_detr_model: "PiGroups"
+edmfx_sgs_mass_flux: true
+edmfx_sgs_diffusive_flux: true
+edmfx_nh_pressure: true
+edmfx_velocity_relaxation: true
+prognostic_tke: true
+cloud_model: "grid_scale"
+moist: "equil"
+config: "column"
+z_max: 3e3
+z_elem: 60
+z_stretch: false
+perturb_initstate: false
+dt: "1secs"
+t_end: "6hours"
+dt_save_to_disk: "10mins"
+toml: [toml/prognostic_edmfx_bomex.toml]
+netcdf_output_at_levels: true
+netcdf_interpolation_num_points: [2, 2, 60]
+diagnostics:
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+    period: 10mins
+  - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
+    period: 10mins
+  - short_name: [entr, detr, lmix, bgrad, strain]
+    period: 10mins

post_processing/ci_plots.jl

@@ -1040,26 +1040,26 @@ function make_plots(
 
     short_names = [
         "ua",
-        "wa",
+        "va",
+        "tke",
+        "lmix",
+        "bgrad",
+        "strain",
         "thetaa",
         "thetaaup",
         "ta",
         "taup",
         "ha",
         "haup",
+        "wa",
         "waup",
-        "tke",
-        "lmix",
         "arup",
         "hus",
         "husup",
         "hur",
         "hurup",
-        "cl",
         "clw",
         "clwup",
-        "cli",
-        "cliup",
         precip_names...,
     ]
     reduction = "inst"
@@ -1081,6 +1081,13 @@ function make_plots(
         ([slice(v, time = LAST_SNAP) for v in group]...,) for
         group in var_groups_zt
     ]
+    # var_groups_z = [
+    #     (
+    #         [
+    #             ClimaAnalysis.window(v, "time", left = 14400, right = 21600) |> ClimaAnalysis.average_time for v in group
+    #         ]...,
+    #     ) for group in var_groups_zt
+    # ]
 
     tmp_file = make_plots_generic(
         output_paths,

src/cache/precomputed_quantities.jl

@@ -38,7 +38,7 @@ function precomputed_quantities(Y, atmos)
             !(atmos.turbconv_model isa DiagnosticEDMFX)
     @assert !(atmos.moisture_model isa DryModel) ||
             !(atmos.turbconv_model isa PrognosticEDMFX)
-    @assert isnothing(atmos.turbconv_model) || isnothing(atmos.vert_diff)
+    # @assert isnothing(atmos.turbconv_model) || isnothing(atmos.vert_diff)
     TST = thermo_state_type(atmos.moisture_model, FT)
     SCT = SurfaceConditions.surface_conditions_type(atmos, FT)
     cspace = axes(Y.c)

src/cache/prognostic_edmf_precomputed_quantities.jl

@@ -182,7 +182,7 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
         ᶜK_h,
         ρatke_flux,
     ) = p.precomputed
-    (; ᶜuʲs, ᶜtsʲs, ᶠu³ʲs, ᶜρʲs, ᶜentrʲs, ᶜdetrʲs) = p.precomputed
+    (; ᶜuʲs, ᶜtsʲs, ᶠu³ʲs, ᶜρʲs, ᶜentrʲs, ᶜdetrʲs, ᶜtke⁰) = p.precomputed
     (; ustar, obukhov_length, buoyancy_flux) = p.precomputed.sfc_conditions
 
     ᶜz = Fields.coordinate_field(Y.c).z
@@ -200,14 +200,14 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             z_sfc,
             ᶜp,
             Y.c.ρ,
-            buoyancy_flux,
             draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)),
             get_physical_w(ᶜuʲs.:($$j), ᶜlg),
             TD.relative_humidity(thermo_params, ᶜtsʲs.:($$j)),
             ᶜphysical_buoyancy(params, Y.c.ρ, ᶜρʲs.:($$j)),
             get_physical_w(ᶜu, ᶜlg),
             TD.relative_humidity(thermo_params, ᶜts⁰),
             FT(0),
+            max(ᶜtke⁰, 0),
             p.atmos.edmfx_entr_model,
         )
         @. ᶜentrʲs.:($$j) = limit_entrainment(
@@ -224,7 +224,6 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             z_sfc,
             ᶜp,
             Y.c.ρ,
-            buoyancy_flux,
             Y.c.sgsʲs.:($$j).ρa,
             draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)),
             get_physical_w(ᶜuʲs.:($$j), ᶜlg),
@@ -236,6 +235,7 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
             ᶜentrʲs.:($$j),
             ᶜvert_div,
             ᶜmassflux_vert_div,
+            max(ᶜtke⁰, 0),
             p.atmos.edmfx_detr_model,
         )
         @. ᶜdetrʲs.:($$j) = limit_detrainment(
@@ -303,7 +303,7 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t)
 
     turbconv_params = CAP.turbconv_params(params)
     c_m = CAP.tke_ed_coeff(turbconv_params)
-    @. ᶜK_u = c_m * ᶜmixing_length * sqrt(max(ᶜtke⁰, 0))
+    @. ᶜK_u = max(c_m * ᶜmixing_length * sqrt(max(ᶜtke⁰, 0)), 0)
     @. ᶜK_h = ᶜK_u / ᶜprandtl_nvec
 
     ρatke_flux_values = Fields.field_values(ρatke_flux)

src/initial_conditions/initial_conditions.jl

@@ -774,6 +774,13 @@ end
 The `InitialCondition` described in [Brown2002](@cite), but with a
 hydrostatically balanced pressure profile.
 """
+
+bomex_tke(::Type{FT}) where {FT} = z -> if (z <= 2500.0)
+    FT(1) - z / 3000
+else
+    FT(0)
+end
+
 Base.@kwdef struct ARM_SGP <: InitialCondition
     prognostic_tke::Bool = false
 end
@@ -811,7 +818,7 @@ for IC in (:Soares, :Bomex, :LifeCycleTan2018, :ARM_SGP)
                 ),
                 velocity = Geometry.UVector(u(z)),
                 turbconv_state = EDMFState(;
-                    tke = prognostic_tke ? FT(0) : tke(z),
+                    tke = prognostic_tke ? bomex_tke(FT)(z) : tke(z),
                 ),
             )
         end

src/prognostic_equations/edmfx_closures.jl

@@ -110,8 +110,37 @@ function edmfx_nh_pressure_tendency!(
     ᶠlg = Fields.local_geometry_field(Y.f)
 
     scale_height = CAP.R_d(params) * CAP.T_surf_ref(params) / CAP.grav(params)
+    FT = eltype(Y)
+    ᶜz = Fields.coordinate_field(Y.c).z
+    z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
+    z_sfc_data = Fields.field_values(z_sfc[colidx])
+    turbconv_params = CAP.turbconv_params(params)
+    a_min = CAP.min_area(turbconv_params)
+    ᶜupdraft_top = p.scratch.ᶜtemp_scalar
 
     for j in 1:n
+        # # look for updraft top
+        # ᶜupdraft_top_data = Fields.field_values(ᶜupdraft_top[colidx])
+        # @. ᶜupdraft_top_data = FT(0)
+        # for level in 1:Spaces.nlevels(axes(ᶜz))
+        #     ρaʲ_lev = Fields.field_values(
+        #         Spaces.level(Y.c.sgsʲs.:($j).ρa[colidx], level),
+        #     )
+        #     ρʲ_lev =
+        #         Fields.field_values(Spaces.level(ᶜρʲs.:($j)[colidx], level))
+        #     ᶜz_lev = Fields.field_values(Spaces.level(ᶜz[colidx], level))
+        #     @. ᶜupdraft_top_data = ifelse(
+        #         draft_area(ρaʲ_lev, ρʲ_lev) > a_min,
+        #         ᶜz_lev,
+        #         ᶜupdraft_top_data,
+        #     )
+        # end
+        # @. ᶜupdraft_top_data = ᶜupdraft_top_data - z_sfc_data
+
+        # # There's only one updraft_top per column, so it's
+        # # safe to use at cell centers and cell faces, correct?:
+        # ᶠupdraft_top = Fields.Field(ᶜupdraft_top_data, axes(ᶠu₃⁰[colidx]))
+
         @. ᶠnh_pressure₃ʲs.:($$j)[colidx] = ᶠupdraft_nh_pressure(
             params,
             p.atmos.edmfx_nh_pressure,
@@ -193,6 +222,7 @@ function mixing_length(
 
     # compute the maximum mixing length at height z
     l_z = ᶜz - z_sfc
+    l_z = 1000000
 
     # compute the l_W - the wall constraint mixing length
     # which imposes an upper limit on the size of eddies near the surface
@@ -242,14 +272,20 @@ function mixing_length(
     )
 
     # add limiters
+    # l = SA.SVector(
+    #     l_N > l_z ? l_z : l_N,
+    #     l_TKE > l_z ? l_z : l_TKE,
+    #     l_W > l_z ? l_z : l_W,
+    # )
     l = SA.SVector(
-        l_N > l_z ? l_z : l_N,
-        l_TKE > l_z ? l_z : l_TKE,
-        l_W > l_z ? l_z : l_W,
+        (l_N < eps(FT) || l_N > l_z) ? l_z : l_N,
+        (l_TKE < eps(FT) || l_TKE > l_z) ? l_z : l_TKE,
+        (l_W < eps(FT) || l_W > l_z) ? l_z : l_W,
     )
     # get soft minimum
     l_smin = lamb_smooth_minimum(l, smin_ub, smin_rm)
     l_limited = max(l_smag, min(l_smin, l_z))
+    #l_limited = min(200, l_z)
 
     return l_limited
 end
@@ -305,11 +341,25 @@ function edmfx_velocity_relaxation_tendency!(
 
     n = n_mass_flux_subdomains(turbconv_model)
     (; dt) = p
+    (; ᶜK, ᶜh_tot, ᶜspecific) = p.precomputed
+    ᶜu₃ʲ = p.scratch.ᶜtemp_C3
 
     if p.atmos.edmfx_velocity_relaxation
         for j in 1:n
             @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] -=
                 C3(min(Y.f.sgsʲs.:($$j).u₃[colidx].components.data.:1, 0)) / dt
+            # @. Yₜ.c.sgsʲs.:($$j).ρa[colidx] -= min(Y.c.sgsʲs.:($$j).ρa[colidx], 0) / dt
+            # @. ᶜu₃ʲ[colidx] = ᶜinterp(Y.f.sgsʲs.:($$j).u₃[colidx])
+            # @. Yₜ.c.sgsʲs.:($$j).mse[colidx] -= ifelse(
+            #     ᶜu₃ʲ[colidx].components.data.:1 < 0,
+            #     (Y.c.sgsʲs.:($$j).mse[colidx] - (ᶜh_tot[colidx] - ᶜK[colidx])) / dt,
+            #     0,
+            # )
+            # @. Yₜ.c.sgsʲs.:($$j).q_tot[colidx] -= ifelse(
+            #     ᶜu₃ʲ[colidx].components.data.:1 < 0,
+            #     (Y.c.sgsʲs.:($$j).q_tot[colidx] - ᶜspecific.q_tot[colidx]) / dt,
+            #     0,
+            # )
         end
     end
 end