add vertical fluctuation and nudging

config/default_configs/default_config.yml

@@ -229,6 +229,12 @@ check_precipitation:
 ls_adv:
   help: "Large-scale advection [`nothing` (default), `Bomex`, `LifeCycleTan2018`, `Rico`, `ARM_SGP`, `GATE_III`]"
   value: ~
+vert_fluc:
+  help: "Vertical fluctuation [`nothing` (default), `GCMDriven`]"
+  value: ~
+nudging:
+  help: "Nudging to a mean profile [`nothing` (default), `GCMDriven`]"
+  value: ~
 fps:
   help: "Frames per second for animations"
   value: 5

config/model_configs/prognostic_edmfx_gcmdriven_column.yml

@@ -0,0 +1,40 @@
+job_id: "prognostic_edmfx_gcmdriven_column"
+initial_condition: "Bomex"
+subsidence: "Bomex"
+edmf_coriolis: "Bomex"
+ls_adv: "GCMDriven"
+vert_fluc: "GCMDriven"
+nudging: "GCMDriven"
+surface_setup: "Bomex"
+turbconv: "prognostic_edmfx"
+implicit_diffusion: true
+implicit_sgs_advection: true
+approximate_linear_solve_iters: 2
+max_newton_iters_ode: 3
+edmfx_upwinding: first_order
+edmfx_entr_model: "Generalized"
+edmfx_detr_model: "Generalized"
+edmfx_sgs_mass_flux: true
+edmfx_sgs_diffusive_flux: true
+edmfx_nh_pressure: true
+edmfx_velocity_relaxation: true
+prognostic_tke: true
+moist: "equil"
+config: "column"
+z_max: 3e3
+z_elem: 60
+z_stretch: false
+perturb_initstate: false
+dt: "50secs"
+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, edt, evu]
+    period: 10mins

post_processing/ci_plots.jl

@@ -945,6 +945,7 @@ EDMFBoxPlots = Union{
     Val{:prognostic_edmfx_gabls_column},
     Val{:prognostic_edmfx_bomex_fixtke_column},
     Val{:prognostic_edmfx_bomex_column},
+    Val{:prognostic_edmfx_gcmdriven_column},
     Val{:prognostic_edmfx_bomex_stretched_column},
     Val{:prognostic_edmfx_bomex_explicit_column},
     Val{:prognostic_edmfx_dycoms_rf01_column},

src/ClimaAtmos.jl

@@ -50,6 +50,8 @@ include(joinpath("prognostic_equations", "implicit", "implicit_solver.jl"))
 
 include(joinpath("prognostic_equations", "remaining_tendency.jl"))
 include(joinpath("prognostic_equations", "forcing", "large_scale_advection.jl")) # TODO: should this be in tendencies/?
+include(joinpath("prognostic_equations", "forcing", "vertical_fluctuation.jl"))
+include(joinpath("prognostic_equations", "forcing", "nudging.jl"))
 include(joinpath("prognostic_equations", "forcing", "subsidence.jl"))
 
 include(joinpath("prognostic_equations", "surface_temp.jl"))

src/cache/cache.jl

@@ -19,6 +19,8 @@ struct AtmosCache{
     PR,
     SUB,
     LSAD,
+    VERTFLUC,
+    NUDGING,
     EDMFCOR,
     FOR,
     NONGW,
@@ -81,6 +83,8 @@ struct AtmosCache{
     precipitation::PR
     subsidence::SUB
     large_scale_advection::LSAD
+    vertical_fluctuation::VERTFLUC
+    nudging::NUDGING
     edmf_coriolis::EDMFCOR
     forcing::FOR
     non_orographic_gravity_wave::NONGW
@@ -225,6 +229,8 @@ function build_cache(Y, atmos, params, surface_setup, sim_info, aerosol_names)
     precipitation = precipitation_cache(Y, atmos)
     subsidence = subsidence_cache(Y, atmos)
     large_scale_advection = large_scale_advection_cache(Y, atmos)
+    vertical_fluctuation = vertical_fluctuation_cache(Y, atmos)
+    nudging = nudging_cache(Y, atmos)
     edmf_coriolis = edmf_coriolis_cache(Y, atmos)
     forcing = forcing_cache(Y, atmos)
     non_orographic_gravity_wave = non_orographic_gravity_wave_cache(Y, atmos)
@@ -253,6 +259,8 @@ function build_cache(Y, atmos, params, surface_setup, sim_info, aerosol_names)
         precipitation,
         subsidence,
         large_scale_advection,
+        vertical_fluctuation,
+        nudging,
         edmf_coriolis,
         forcing,
         non_orographic_gravity_wave,

src/prognostic_equations/forcing/nudging.jl

@@ -0,0 +1,64 @@
+#####
+##### Nudging/relaxation tendency for GCM-driven SCM
+#####
+
+import ClimaCore.Fields as Fields
+
+nudging_cache(Y, atmos::AtmosModel) =
+    nudging_cache(Y, atmos.nudging)
+
+nudging_cache(Y, ::Nothing) = (;)
+nudging_tendency!(Yₜ, Y, p, t, colidx, ::Nothing) = nothing
+
+function nudging_cache(Y, ::Nudging)
+    FT = Spaces.undertype(axes(Y.c))
+    ᶜT_mean = similar(Y.c, FT)
+    ᶜq_tot_mean = similar(Y.c, FT)
+    ᶜu_mean = similar(Y.c, FT)
+    ᶜv_mean = similar(Y.c, FT)
+    @. ᶜT_mean = 290
+    @. ᶜq_tot_mean = 0.010
+    @. ᶜu_mean = 5
+    @. ᶜv_mean = 0
+    hr = 3600
+    ᶜτ_scalar = 6hr
+    ᶜτ_wind = 1hr
+    return (;
+        ᶜτ_scalar,
+        ᶜτ_wind,
+        ᶜT_mean,
+        ᶜq_tot_mean,
+        ᶜu_mean,
+        ᶜv_mean,
+    )
+end
+
+function nudging_tendency!(Yₜ, Y, p, t, colidx, ::Nudging)
+
+    thermo_params = CAP.thermodynamics_params(p.params)
+    T_0 = TD.Parameters.T_0(thermo_params)
+    Lv_0 = TD.Parameters.LH_v0(thermo_params)
+    cv_v = TD.Parameters.cv_v(thermo_params)
+    R_v = TD.Parameters.R_v(thermo_params)
+
+    (; ᶜu_mean, ᶜv_mean, ᶜτ_wind, ᶜτ_scalar, ᶜT_mean, ᶜq_tot_mean) = p.nudging
+    (; ᶜspecific, ᶜts) = p.precomputed
+    ᶜlg = Fields.local_geometry_field(Y.c)
+    ᶜuₕ_mean_colidx = C12(Geometry.UVVector(ᶜu_mean[colidx], ᶜv_mean[colidx]), ᶜlg[colidx])
+    @. Yₜ.c.uₕ[colidx] -= (Y.c.uₕ[colidx] - ᶜuₕ_mean_colidx) / ᶜτ_wind
+
+    dTdt_nudging = p.scratch.ᶜtemp_scalar
+    dqtdt_nudging = p.scratch.ᶜtemp_scalar_2
+    @. dTdt_nudging[colidx] = - (TD.air_temperature(thermo_params, ᶜts[colidx]) - ᶜT_mean[colidx]) / ᶜτ_scalar
+    @. dqtdt_nudging[colidx] = - (ᶜspecific.q_tot[colidx] - ᶜq_tot_mean[colidx]) / ᶜτ_scalar
+
+    @. Yₜ.c.ρq_tot[colidx] += Y.c.ρ[colidx] * dqtdt_nudging[colidx]
+
+    @. Yₜ.c.ρe_tot[colidx] +=
+        Y.c.ρ[colidx] * (
+            TD.cp_m(thermo_params, ᶜts) * ᶜdTdt_nudging[colidx] + (
+                cv_v * (TD.air_temperature(thermo_params, ᶜts[colidx]) - T_0) + Lv_0 -
+                R_v * T_0
+            ) * ᶜdqtdt_nudging[colidx]
+        )
+end

src/prognostic_equations/forcing/vertical_fluctuation.jl

@@ -0,0 +1,55 @@
+#####
+##### Vertical fluctuation (for single column experiments)
+#####
+
+import Thermodynamics as TD
+import ClimaCore.Spaces as Spaces
+import ClimaCore.Fields as Fields
+
+vertical_fluctuation_cache(Y, atmos::AtmosModel) =
+    vertical_fluctuation_cache(Y, atmos.vert_fluc)
+
+vertical_fluctuation_cache(Y, vert_fluc::Nothing) = (;)
+function vertical_fluctuation_cache(Y, vert_fluc::VerticalFluctuation)
+    FT = Spaces.undertype(axes(Y.c))
+    return (; ᶜdqtdt_fluc = similar(Y.c, FT), ᶜdTdt_fluc = similar(Y.c, FT))
+end
+
+vertical_fluctuation_tendency!(Yₜ, Y, p, t, colidx, ::Nothing) = nothing
+function vertical_fluctuation_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    colidx,
+    vert_fluc::VerticalFluctuation,
+)
+    (; prof_dTdt, prof_dqtdt) = vert_fluc
+
+    thermo_params = CAP.thermodynamics_params(p.params)
+    ᶜts = p.precomputed.ᶜts[colidx]
+    ᶜdqtdt_fluc = p.vertical_fluctuation.ᶜdqtdt_fluc[colidx]
+    ᶜdTdt_fluc = p.vertical_fluctuation.ᶜdTdt_fluc[colidx]
+    z = Fields.coordinate_field(axes(ᶜdqtdt_fluc)).z
+
+    @. ᶜdTdt_fluc = prof_dTdt(thermo_params, ᶜts, t, z)
+    @. ᶜdqtdt_fluc = prof_dqtdt(thermo_params, ᶜts, t, z)
+
+    T_0 = TD.Parameters.T_0(thermo_params)
+    Lv_0 = TD.Parameters.LH_v0(thermo_params)
+    cv_v = TD.Parameters.cv_v(thermo_params)
+    R_v = TD.Parameters.R_v(thermo_params)
+
+    @. Yₜ.c.ρq_tot[colidx] += Y.c.ρ[colidx] * ᶜdqtdt_fluc
+    # TODO: should `hv` be a thermo function?
+    #     (hv = cv_v * (ᶜT - T_0) + Lv_0 - R_v * T_0)
+    @. Yₜ.c.ρe_tot[colidx] +=
+        Y.c.ρ[colidx] * (
+            TD.cp_m(thermo_params, ᶜts) * ᶜdTdt_fluc +
+            (
+                cv_v * (TD.air_temperature(thermo_params, ᶜts) - T_0) + Lv_0 -
+                R_v * T_0
+            ) * ᶜdqtdt_fluc
+        )
+    return nothing
+end

src/prognostic_equations/remaining_tendency.jl

@@ -33,6 +33,8 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
         subsidence_tendency!(Yₜ, Y, p, t, colidx, p.atmos.subsidence)
         edmf_coriolis_tendency!(Yₜ, Y, p, t, colidx, p.atmos.edmf_coriolis)
         large_scale_advection_tendency!(Yₜ, Y, p, t, colidx, p.atmos.ls_adv)
+        vertical_fluctuation_tendency!(Yₜ, Y, p, t, colidx, p.atmos.vert_fluc)
+        nudging_tendency!(Yₜ, Y, p, t, colidx, p.atmos.nudging)
 
         if p.atmos.sgs_adv_mode == Explicit()
             edmfx_sgs_vertical_advection_tendency!(

src/solver/model_getters.jl

@@ -333,18 +333,22 @@ function get_large_scale_advection_model(parsed_args, ::Type{FT}) where {FT}
         (APL.ARM_SGP_dTdt(FT), APL.ARM_SGP_dqtdt(FT))
     elseif ls_adv == "GATE_III"
         (APL.GATE_III_dTdt(FT), APL.GATE_III_dqtdt(FT))
+    elseif ls_adv == "GCMDriven"
+        (APL.Bomex_dTdt(FT), APL.Bomex_dqtdt(FT))
     else
         error("Uncaught case")
     end
     # See https://clima.github.io/AtmosphericProfilesLibrary.jl/dev/
     # for which functions accept which arguments.
-    prof_dqtdt = if ls_adv in ("Bomex", "LifeCycleTan2018", "Rico", "GATE_III")
-        (thermo_params, ᶜts, t, z) -> prof_dqtdt₀(z)
-    elseif ls_adv == "ARM_SGP"
-        (thermo_params, ᶜts, t, z) ->
-            prof_dqtdt₀(TD.exner(thermo_params, ᶜts), t, z)
-    end
-    prof_dTdt = if ls_adv in ("Bomex", "LifeCycleTan2018", "Rico")
+    prof_dqtdt =
+        if ls_adv in
+           ("Bomex", "LifeCycleTan2018", "Rico", "GATE_III", "GCMDriven")
+            (thermo_params, ᶜts, t, z) -> prof_dqtdt₀(z)
+        elseif ls_adv == "ARM_SGP"
+            (thermo_params, ᶜts, t, z) ->
+                prof_dqtdt₀(TD.exner(thermo_params, ᶜts), t, z)
+        end
+    prof_dTdt = if ls_adv in ("Bomex", "LifeCycleTan2018", "Rico", "GCMDriven")
         (thermo_params, ᶜts, t, z) ->
             prof_dTdt₀(TD.exner(thermo_params, ᶜts), z)
     elseif ls_adv == "ARM_SGP"
@@ -356,6 +360,39 @@ function get_large_scale_advection_model(parsed_args, ::Type{FT}) where {FT}
     return LargeScaleAdvection(prof_dTdt, prof_dqtdt)
 end
 
+function get_vertical_fluctuation_model(parsed_args, ::Type{FT}) where {FT}
+    vert_fluc = parsed_args["vert_fluc"]
+    vert_fluc == nothing && return nothing
+
+    (prof_dTdt₀, prof_dqtdt₀) = if vert_fluc == "GCMDriven"
+        (APL.Bomex_dTdt(FT), APL.Bomex_dqtdt(FT))
+    else
+        error("Uncaught case")
+    end
+    # See https://clima.github.io/AtmosphericProfilesLibrary.jl/dev/
+    # for which functions accept which arguments.
+    prof_dqtdt = if vert_fluc == "GCMDriven"
+        (thermo_params, ᶜts, t, z) -> prof_dqtdt₀(z)
+    end
+    prof_dTdt = if vert_fluc == "GCMDriven"
+        (thermo_params, ᶜts, t, z) ->
+            prof_dTdt₀(TD.exner(thermo_params, ᶜts), z)
+    end
+
+    return VerticalFluctuation(prof_dTdt, prof_dqtdt)
+end
+
+function get_nudging_model(parsed_args)
+    nudging_name = parsed_args["nudging"]
+    return if isnothing(nudging_name)
+        nothing
+    elseif nudging_name == "GCMDriven"
+        Nudging()
+    else
+        error("Uncaught nudging model `$nudging_name`.")
+    end
+end
+
 function get_edmf_coriolis(parsed_args, ::Type{FT}) where {FT}
     edmf_coriolis = parsed_args["edmf_coriolis"]
     edmf_coriolis == nothing && return nothing

src/solver/type_getters.jl

@@ -60,6 +60,8 @@ function get_atmos(config::AtmosConfig, params)
         radiation_mode,
         subsidence = get_subsidence_model(parsed_args, radiation_mode, FT),
         ls_adv = get_large_scale_advection_model(parsed_args, FT),
+        vert_fluc = get_vertical_fluctuation_model(parsed_args, FT),
+        nudging = get_nudging_model(parsed_args),
         edmf_coriolis = get_edmf_coriolis(parsed_args, FT),
         advection_test,
         tendency_model = get_tendency_model(parsed_args),

src/solver/types.jl

@@ -114,6 +114,11 @@ struct LargeScaleAdvection{PT, PQ}
     prof_dTdt::PT # Set large-scale cooling
     prof_dqtdt::PQ # Set large-scale drying
 end
+struct VerticalFluctuation{PT, PQ}
+    prof_dTdt::PT # Set large-scale cooling
+    prof_dqtdt::PQ # Set large-scale drying
+end
+struct Nudging end  # maybe need to <: AbstractForcing
 
 struct EDMFCoriolis{U, V, FT}
     prof_ug::U
@@ -332,6 +337,8 @@ Base.@kwdef struct AtmosModel{
     S,
     RM,
     LA,
+    VF,
+    NUDGING,
     EC,
     AT,
     TM,
@@ -365,6 +372,8 @@ Base.@kwdef struct AtmosModel{
     subsidence::S = nothing
     radiation_mode::RM = nothing
     ls_adv::LA = nothing
+    vert_fluc::VF = nothing
+    nudging::NUDGING = nothing
     edmf_coriolis::EC = nothing
     advection_test::AT = nothing
     tendency_model::TM = nothing