add vertical fluctuation and nudging for GCM driven SCM

.buildkite/pipeline.yml

@@ -629,6 +629,14 @@ steps:
         artifact_paths: "prognostic_edmfx_trmm_column_0M/*"
         agents:
           slurm_mem: 20GB
+
+      - label: ":genie: Prognostic EDMFX GCM driven in a column"
+        command: >
+          julia --color=yes --project=examples examples/hybrid/driver.jl
+          --config_file $CONFIG_PATH/prognostic_edmfx_gcmdriven_column.yml
+        artifact_paths: "prognostic_edmfx_gcmdriven_column/output_active/*"
+        agents:
+          slurm_mem: 20GB
 
       - label: ":genie: Prognostic EDMFX Bomex in a box"
         command: >

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: "Bomex"
+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/cache/temporary_quantities.jl

@@ -32,6 +32,7 @@ function temporary_quantities(Y, atmos)
         temp_data_level_3 = Fields.field_values(
             Fields.level(Fields.Field(FT, center_space), 1),
         ), # ρaʲu³ʲ_datah_tot
+        ᶜtemp_C12 = Fields.Field(C12{FT}, center_space), # ᶜuₕ_mean
         ᶜtemp_C3 = Fields.Field(C3{FT}, center_space), # ᶜ∇Φ₃
         ᶜtemp_CT3 = Fields.Field(CT3{FT}, center_space), # ᶜω³, ᶜ∇Φ³
         ᶠtemp_CT3 = Fields.Field(CT3{FT}, face_space), # ᶠuₕ³

src/prognostic_equations/forcing/nudging.jl

@@ -0,0 +1,63 @@
+#####
+##### 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)
+    # TODO: read profiles from LES
+    @. ᶜT_mean = 290
+    @. ᶜq_tot_mean = 0.010
+    @. ᶜu_mean = 5
+    @. ᶜv_mean = 0
+    # TODO: make it a function of z and add timescale to climaparams
+    hr = 3600
+    ᶜτ_scalar = 24hr
+    ᶜτ_wind = 6hr
+    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 = p.scratch.ᶜtemp_C12
+    @. ᶜ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.ρe_tot[colidx] +=
+        Y.c.ρ[colidx] * (
+            TD.cv_m(thermo_params, ᶜts[colidx]) * ᶜdTdt_nudging[colidx] +
+            (
+                cv_v * (TD.air_temperature(thermo_params, ᶜts[colidx]) - T_0) +
+                Lv_0 - R_v * T_0
+            ) * ᶜdqtdt_nudging[colidx]
+        )
+    @. Yₜ.c.ρq_tot[colidx] += Y.c.ρ[colidx] * ᶜdqtdt_nudging[colidx]
+end

src/prognostic_equations/forcing/vertical_fluctuation.jl

@@ -0,0 +1,52 @@
+#####
+##### 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, ::Nothing) = (;)
+function vertical_fluctuation_cache(Y, ::VerticalFluctuation)
+    FT = Spaces.undertype(axes(Y.c))
+    ᶜdTdt_fluc = similar(Y.c, FT)
+    ᶜdqtdt_fluc = similar(Y.c, FT)
+    # TODO: read profiles from LES
+    @. ᶜdTdt_fluc = 0
+    @. ᶜdqtdt_fluc = 0
+    return (; ᶜdTdt_fluc, ᶜdqtdt_fluc)
+end
+
+vertical_fluctuation_tendency!(Yₜ, Y, p, t, colidx, ::Nothing) = nothing
+function vertical_fluctuation_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    colidx,
+    ::VerticalFluctuation,
+)
+    (; params) = p
+    (; ᶜts) = p.precomputed
+    (; ᶜdTdt_fluc, ᶜdqtdt_fluc) = p.vertical_fluctuation
+    thermo_params = CAP.thermodynamics_params(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)
+
+    @. Yₜ.c.ρe_tot[colidx] +=
+        Y.c.ρ[colidx] * (
+            TD.cv_m(thermo_params, ᶜts[colidx]) * ᶜdTdt_fluc[colidx] +
+            (
+                cv_v * (TD.air_temperature(thermo_params, ᶜts[colidx]) - T_0) +
+                Lv_0 - R_v * T_0
+            ) * ᶜdqtdt_fluc[colidx]
+        )
+    @. Yₜ.c.ρq_tot[colidx] += Y.c.ρ[colidx] * ᶜdqtdt_fluc[colidx]
+
+    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

@@ -356,6 +356,26 @@ 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_name = parsed_args["vert_fluc"]
+    @assert vert_fluc_name in (nothing, "GCMDriven")
+    return if isnothing(vert_fluc_name)
+        nothing
+    elseif vert_fluc_name == "GCMDriven"
+        VerticalFluctuation()
+    end
+end
+
+function get_nudging_model(parsed_args)
+    nudging_name = parsed_args["nudging"]
+    @assert nudging_name in (nothing, "GCMDriven")
+    return if isnothing(nudging_name)
+        nothing
+    elseif nudging_name == "GCMDriven"
+        Nudging()
+    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,9 @@ struct LargeScaleAdvection{PT, PQ}
     prof_dTdt::PT # Set large-scale cooling
     prof_dqtdt::PQ # Set large-scale drying
 end
+# maybe need to <: AbstractForcing
+struct VerticalFluctuation end
+struct Nudging end
 
 struct EDMFCoriolis{U, V, FT}
     prof_ug::U
@@ -332,6 +335,8 @@ Base.@kwdef struct AtmosModel{
     S,
     RM,
     LA,
+    VF,
+    NUDGING,
     EC,
     AT,
     TM,
@@ -365,6 +370,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

test/coupler_compatibility.jl

@@ -77,6 +77,8 @@ const T2 = 290
         p.precipitation,
         p.subsidence,
         p.large_scale_advection,
+        p.vertical_fluctuation,
+        p.nudging,
         p.edmf_coriolis,
         p.forcing,
         p.non_orographic_gravity_wave,