Change LES default imin index, expose surface area parameter in diag…

…nostic EDMF, set top cell forcing tendency to zero

src/cache/diagnostic_edmf_precomputed_quantities.jl

@@ -104,6 +104,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_bottom_bc!(
     (; ᶠu³⁰, ᶜK⁰) = p.precomputed
 
     thermo_params = CAP.thermodynamics_params(p.params)
+    turbconv_params = CAP.turbconv_params(params)
 
     ρ_int_level = Fields.field_values(Fields.level(Y.c.ρ, 1))
     uₕ_int_level = Fields.field_values(Fields.level(Y.c.uₕ, 1))
@@ -155,7 +156,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_bottom_bc!(
         @. mseʲ_int_level = sgs_scalar_first_interior_bc(
             z_int_level - z_sfc_halflevel,
             ρ_int_level,
-            turbconv_model.a_int,
+            FT(turbconv_params.surface_area),
             h_tot_int_level - K_int_level,
             buoyancy_flux_sfc_halflevel,
             ρ_flux_h_tot_sfc_halflevel,
@@ -166,7 +167,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_bottom_bc!(
         @. q_totʲ_int_level = sgs_scalar_first_interior_bc(
             z_int_level - z_sfc_halflevel,
             ρ_int_level,
-            turbconv_model.a_int,
+            FT(turbconv_params.surface_area),
             q_tot_int_level,
             buoyancy_flux_sfc_halflevel,
             ρ_flux_q_tot_sfc_halflevel,
@@ -191,7 +192,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_bottom_bc!(
             p_int_level,
             Φ_int_level,
         )
-        @. ρaʲ_int_level = ρʲ_int_level * turbconv_model.a_int
+        @. ρaʲ_int_level = ρʲ_int_level * FT(turbconv_params.surface_area)
     end
 
     ρaʲs_int_level = Fields.field_values(Fields.level(ᶜρaʲs, 1))

src/initial_conditions/initial_conditions.jl

@@ -1177,7 +1177,7 @@ end
 function gcm_initial_conditions(external_forcing_file)
     NC.NCDataset(external_forcing_file) do ds
         (  # TODO: Cast to CuVector for GPU compatibility
-            gcm_driven_profile(ds, "thetali_mean")[:, 1],  # 1 is initial time index
+            gcm_driven_profile(ds, "thetali_mean")[:, 1], # 1 is initial time index
             gcm_driven_profile(ds, "u_mean")[:, 1],
             gcm_driven_profile(ds, "v_mean")[:, 1],
             gcm_driven_profile(ds, "qt_mean")[:, 1],

src/prognostic_equations/forcing/external_forcing.jl

@@ -52,14 +52,13 @@ function external_forcing_cache(Y, external_forcing::GCMForcing)
     ᶜls_subsidence = similar(Y.c, FT)
 
     (; external_forcing_file) = external_forcing
-    imin = 100  # TODO: move into `GCMForcing` (and `parsed_args`)
 
     NC.Dataset(external_forcing_file, "r") do ds
         function setvar!(cc_field, varname, colidx, zc_gcm, zc_les)
             parent(cc_field[colidx]) .= interp_vertical_prof(
                 zc_gcm,
                 zc_les,
-                gcm_driven_profile_tmean(ds, varname; imin),  # TODO: time-varying tendencies
+                gcm_driven_profile_tmean(ds, varname),  # TODO: time-varying tendencies
             )
         end
 
@@ -181,6 +180,12 @@ function external_forcing_tendency!(Yₜ, Y, p, t, ::GCMForcing)
         Val{:first_order}(),
     )
     # <-- subsidence
+    # needed to address top boundary condition for forcings. Otherwise upper portion of domain is anomalously warm
+    ρe_tot_top = Fields.level(Yₜ.c.ρe_tot, Spaces.nlevels(axes(Y.c)))
+    @. ρe_tot_top = 0.0
+
+    ρq_tot_top = Fields.level(Yₜ.c.ρq_tot, Spaces.nlevels(axes(Y.c)))
+    @. ρq_tot_top = 0.0
 
     return nothing
 end

src/solver/model_getters.jl

@@ -426,7 +426,7 @@ function get_turbconv_model(FT, parsed_args, turbconv_params)
     elseif turbconv == "diagnostic_edmfx"
         N = parsed_args["updraft_number"]
         TKE = parsed_args["prognostic_tke"]
-        DiagnosticEDMFX{N, TKE}(FT(0.1), turbconv_params.min_area)
+        DiagnosticEDMFX{N, TKE}(turbconv_params.surface_area)
     else
         nothing
     end

src/solver/types.jl

@@ -176,11 +176,10 @@ PrognosticEDMFX{N, TKE}(a_half::FT) where {N, TKE, FT} =
     PrognosticEDMFX{N, TKE, FT}(a_half)
 
 struct DiagnosticEDMFX{N, TKE, FT} <: AbstractEDMF
-    a_int::FT # area fraction of the first interior cell above the surface
     a_half::FT # WARNING: this should never be used outside of divide_by_ρa
 end
-DiagnosticEDMFX{N, TKE}(a_int::FT, a_half::FT) where {N, TKE, FT} =
-    DiagnosticEDMFX{N, TKE, FT}(a_int, a_half)
+DiagnosticEDMFX{N, TKE}(a_half::FT) where {N, TKE, FT} =
+    DiagnosticEDMFX{N, TKE, FT}(a_half)
 
 n_mass_flux_subdomains(::PrognosticEDMFX{N}) where {N} = N
 n_mass_flux_subdomains(::DiagnosticEDMFX{N}) where {N} = N

src/surface_conditions/surface_setups.jl

@@ -252,14 +252,13 @@ end
 function (surface_setup::GCMDriven)(params)
     FT = eltype(params)
     (; external_forcing_file) = surface_setup
-    imin = 100
-    T, lhf, shf = FT.(gcm_surface_conditions(external_forcing_file, imin))
+    T, lhf, shf = FT.(gcm_surface_conditions(external_forcing_file))
     z0 = FT(1e-4)  # zrough
     parameterization = MoninObukhov(; z0, fluxes = HeatFluxes(; lhf, shf))
     return SurfaceState(; parameterization, T)
 end
 
-function gcm_surface_conditions(external_forcing_file, imin)
+function gcm_surface_conditions(external_forcing_file; imin = 793)
     NC.NCDataset(external_forcing_file) do ds
         (
             mean(gcm_driven_timeseries(ds, "surface_temperature")[imin:end]),

src/utils/read_gcm_driven_scm_data.jl

@@ -3,7 +3,7 @@ import NCDatasets as NC
 import StatsBase: mean
 
 """
-    gcm_driven_profile_tmean(ds, varname; imin = 100)
+    gcm_driven_profile_tmean(ds, varname; imin = 793)
 
 Extract time-averaged data (`imin:end`) for `varname` from the "profile" group in the GCM-driven dataset `ds`
 
@@ -14,7 +14,7 @@ Returns a 1D ("z",) `Vector{FT}` of the time-averaged data.
     This method currently assumes the underlying data is `Float64`. 
     If this is not the case, "garbage" data may be returned with no warning.
 """
-function gcm_driven_profile_tmean(ds, varname; imin = 100)
+function gcm_driven_profile_tmean(ds, varname; imin = 793)
     vec(mean(gcm_driven_profile(ds, varname)[:, imin:end]; dims = 2))
 end