separate entrainment and detrainment and add options

.buildkite/pipeline.yml

@@ -632,7 +632,6 @@ steps:
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/diagnostic_edmfx_aquaplanet_tke.yml
         artifact_paths: "diagnostic_edmfx_aquaplanet_tke/*"
-        soft_fail: true
         agents:
           slurm_mem: 20GB
 

config/default_configs/default_edmf_config.yml

@@ -23,9 +23,12 @@ edmfx_nh_pressure:
 turbconv_case:
   help: "The case run by Turbulence convection scheme [`Bomex` (default), `Bomex`, `DYCOMS_RF01`, `TRMM_LBA`, `GABLS`]"
   value: ~
-edmfx_entr_detr:
-  help: "If set to true, it switches on EDMFX entrainment/detrainment closure.  [`true`, `false` (default)]"
-  value: false
+edmfx_entr_model:
+  help: "EDMFX entrainment closure.  [`nothing` (default), `PiGroups`, `ConstantCoefficent`, `ConstantTimescale`]"
+  value: ~
+edmfx_detr_model:
+  help: "EDMFX detrainment closure.  [`nothing` (default), `PiGroups`, `ConstantCoefficent`]"
+  value: ~
 edmfx_upwinding:
   help: "EDMFX upwinding mode [`none` (default), `first_order` , `third_order`, `boris_book`, `zalesak`]"
   value: none

config/longrun_configs/longrun_aquaplanet_amip.yml

@@ -13,7 +13,8 @@ surface_setup: "DefaultMoninObukhov"
 turbconv: "diagnostic_edmfx" 
 prognostic_tke: true
 edmfx_upwinding: "first_order" 
-edmfx_entr_detr: true 
+edmfx_entr_model: "ConstantTimescale" 
+edmfx_detr_model: "ConstantCoefficient" 
 edmfx_nh_pressure: true 
 edmfx_sgs_flux: true
 rayleigh_sponge: true

config/model_configs/diagnostic_edmfx_aquaplanet_tke.yml

@@ -5,11 +5,12 @@ rad: gray
 turbconv: diagnostic_edmfx 
 prognostic_tke: true
 edmfx_upwinding: first_order 
-edmfx_entr_detr: true 
+edmfx_entr_model: "ConstantTimescale" 
+edmfx_detr_model: "ConstantCoefficient" 
 edmfx_nh_pressure: true 
 edmfx_sgs_flux: true
 moist: equil 
 precip_model: 0M
-dt: 50secs 
-t_end: 3hours 
+dt: 10secs 
+t_end: 1hours 
 dt_save_to_disk: 600secs

config/model_configs/diagnostic_edmfx_bomex_box.yml

@@ -8,7 +8,8 @@ surface_setup: "Bomex"
 turbconv: "diagnostic_edmfx" 
 prognostic_tke: true
 edmfx_upwinding: "first_order" 
-edmfx_entr_detr: true 
+edmfx_entr_model: "ConstantTimescale" 
+edmfx_detr_model: "ConstantCoefficient" 
 edmfx_nh_pressure: true 
 edmfx_sgs_flux: true
 moist: "equil"

config/model_configs/diagnostic_edmfx_bomex_stretched_box.yml

@@ -8,7 +8,8 @@ surface_setup: "Bomex"
 turbconv: "diagnostic_edmfx" 
 prognostic_tke: true
 edmfx_upwinding: "first_order" 
-edmfx_entr_detr: true 
+edmfx_entr_model: "ConstantTimescale" 
+edmfx_detr_model: "ConstantCoefficient" 
 edmfx_nh_pressure: true 
 edmfx_sgs_flux: true
 moist: "equil"

config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml

@@ -7,7 +7,8 @@ surface_setup: DYCOMS_RF01
 turbconv: diagnostic_edmfx 
 prognostic_tke: true
 edmfx_upwinding: first_order 
-edmfx_entr_detr: true 
+edmfx_entr_model: "ConstantTimescale" 
+edmfx_detr_model: "ConstantCoefficient" 
 edmfx_nh_pressure: true 
 edmfx_sgs_flux: true
 moist: equil

config/model_configs/diagnostic_edmfx_gabls_box.yml

@@ -5,7 +5,8 @@ surface_setup: GABLS
 turbconv: diagnostic_edmfx
 prognostic_tke: true
 edmfx_upwinding: first_order
-edmfx_entr_detr: true
+edmfx_entr_model: "ConstantTimescale" 
+edmfx_detr_model: "ConstantCoefficient" 
 edmfx_nh_pressure: true 
 edmfx_sgs_flux: true
 moist: equil

config/model_configs/diagnostic_edmfx_rico_box.yml

@@ -7,7 +7,8 @@ surface_setup: Rico
 turbconv: diagnostic_edmfx
 prognostic_tke: true
 edmfx_upwinding: first_order
-edmfx_entr_detr: true
+edmfx_entr_model: "ConstantTimescale" 
+edmfx_detr_model: "ConstantCoefficient" 
 edmfx_nh_pressure: true
 edmfx_sgs_flux: true
 moist: equil

config/model_configs/diagnostic_edmfx_trmm_box.yml

@@ -6,7 +6,8 @@ surface_setup: TRMM_LBA
 turbconv: diagnostic_edmfx
 prognostic_tke: true
 edmfx_upwinding: first_order
-edmfx_entr_detr: true
+edmfx_entr_model: "ConstantTimescale" 
+edmfx_detr_model: "ConstantCoefficient" 
 edmfx_nh_pressure: true
 edmfx_sgs_flux: true
 moist: equil

config/model_configs/edmfx_bomex_box_v1.yml

@@ -8,7 +8,8 @@ dt: "1secs"
 edmfx_upwinding: first_order
 surface_setup: "Bomex"
 perturb_initstate: false
-edmfx_entr_detr: true
+edmfx_entr_model: "ConstantCoefficient" 
+edmfx_detr_model: "ConstantCoefficient" 
 t_end: "10800secs"
 y_max: 100000.0
 y_elem: 2

config/model_configs/edmfx_bomex_box_v2.yml

@@ -8,7 +8,8 @@ dt: "1secs"
 edmfx_upwinding: first_order
 surface_setup: "Bomex"
 perturb_initstate: false
-edmfx_entr_detr: true
+edmfx_entr_model: "ConstantCoefficient" 
+edmfx_detr_model: "ConstantCoefficient" 
 t_end: "10800secs"
 y_max: 100000.0
 y_elem: 2

config/model_configs/sphere_baroclinic_wave_rhoe_equilmoist_edmfx.yml

@@ -1,7 +1,8 @@
 dt_save_to_disk: "30secs"
 initial_condition: "MoistBaroclinicWaveWithEDMF"
 dt: "1secs"
-edmfx_entr_detr: true
+edmfx_entr_model: "ConstantCoefficient" 
+edmfx_detr_model: "ConstantCoefficient" 
 t_end: "6mins"
 turbconv: "edmfx"
 ode_algo: "SSP33ShuOsher"

src/ClimaAtmos.jl

@@ -80,6 +80,7 @@ include(joinpath("prognostic_equations", "hyperdiffusion.jl"))
 include(joinpath("prognostic_equations", "edmf_coriolis.jl"))
 include(joinpath("prognostic_equations", "buoyancy_gradients.jl"))
 include(joinpath("prognostic_equations", "edmfx_closures.jl"))
+include(joinpath("prognostic_equations", "edmfx_entr_detr.jl"))
 include(joinpath("prognostic_equations", "edmfx_tke.jl"))
 include(joinpath("prognostic_equations", "edmfx_sgs_flux.jl"))
 include(joinpath("prognostic_equations", "edmfx_boundary_condition.jl"))

src/cache/diagnostic_edmf_precomputed_quantities.jl

@@ -140,7 +140,8 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
         ᶜq_totʲs,
         ᶜtsʲs,
         ᶜρʲs,
-        ᶜentr_detrʲs,
+        ᶜentrʲs,
+        ᶜdetrʲs,
         ᶜnh_pressureʲs,
         ᶜS_q_totʲs,
         ᶜS_e_totʲs_helper,
@@ -318,7 +319,8 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
             ᶜtsʲ = ᶜtsʲs.:($j)
             ᶜρʲ = ᶜρʲs.:($j)
             ᶜq_totʲ = ᶜq_totʲs.:($j)
-            ᶜentr_detrʲ = ᶜentr_detrʲs.:($j)
+            ᶜentrʲ = ᶜentrʲs.:($j)
+            ᶜdetrʲ = ᶜdetrʲs.:($j)
             ᶜnh_pressureʲ = ᶜnh_pressureʲs.:($j)
             ᶜS_q_totʲ = ᶜS_q_totʲs.:($j)
             ᶜS_e_totʲ_helper = ᶜS_e_totʲs_helper.:($j)
@@ -341,8 +343,8 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
                 Fields.field_values(Fields.level(ᶜq_totʲ, i - 1))
             ρʲ_prev_level = Fields.field_values(Fields.level(ᶜρʲ, i - 1))
             tsʲ_prev_level = Fields.field_values(Fields.level(ᶜtsʲ, i - 1))
-            entr_detrʲ_prev_level =
-                Fields.field_values(Fields.level(ᶜentr_detrʲ, i - 1))
+            entrʲ_prev_level = Fields.field_values(Fields.level(ᶜentrʲ, i - 1))
+            detrʲ_prev_level = Fields.field_values(Fields.level(ᶜdetrʲ, i - 1))
             nh_pressureʲ_prev_level =
                 Fields.field_values(Fields.level(ᶜnh_pressureʲ, i - 1))
             scale_height =
@@ -352,9 +354,8 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
             S_e_totʲ_helper_prev_level =
                 Fields.field_values(Fields.level(ᶜS_e_totʲ_helper, i - 1))
 
-            @. entr_detrʲ_prev_level = pi_groups_entr_detr(
+            @. entrʲ_prev_level = entrainment(
                 params,
-                p.atmos.edmfx_entr_detr,
                 z_prev_level,
                 z_sfc_halflevel,
                 p_prev_level,
@@ -374,6 +375,31 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
                 TD.relative_humidity(thermo_params, ts_prev_level),
                 ᶜphysical_buoyancy(params, ρ_ref_prev_level, ρ_prev_level),
                 dt,
+                p.atmos.edmfx_entr_model,
+            )
+
+            @. detrʲ_prev_level = detrainment(
+                params,
+                z_prev_level,
+                z_sfc_halflevel,
+                p_prev_level,
+                ρ_prev_level,
+                buoyancy_flux_sfc_halflevel,
+                ρaʲ_prev_level / ρʲ_prev_level,
+                get_physical_w(
+                    u³ʲ_prev_halflevel,
+                    local_geometry_prev_halflevel,
+                ),
+                TD.relative_humidity(thermo_params, tsʲ_prev_level),
+                ᶜphysical_buoyancy(params, ρ_ref_prev_level, ρʲ_prev_level),
+                get_physical_w(
+                    u³⁰_prev_halflevel,
+                    local_geometry_prev_halflevel,
+                ),
+                TD.relative_humidity(thermo_params, ts_prev_level),
+                ᶜphysical_buoyancy(params, ρ_ref_prev_level, ρ_prev_level),
+                dt,
+                p.atmos.edmfx_detr_model,
             )
 
             # TODO: use updraft top instead of scale height
@@ -414,10 +440,7 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
                 (1 / local_geometry_halflevel.J) * (
                     local_geometry_prev_level.J *
                     ρaʲ_prev_level *
-                    (
-                        entr_detrʲ_prev_level.entr -
-                        entr_detrʲ_prev_level.detr + S_q_totʲ_prev_level
-                    )
+                    (entrʲ_prev_level - detrʲ_prev_level + S_q_totʲ_prev_level)
                 )
 
             # Using constant exponents in broadcasts allocate, so we use
@@ -457,8 +480,8 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
                     local_geometry_prev_level.J *
                     2 *
                     (
-                        entr_detrʲ_prev_level.entr * u³⁰_data_prev_halflevel -
-                        entr_detrʲ_prev_level.entr * u³ʲ_data_prev_halflevel
+                        entrʲ_prev_level * u³⁰_data_prev_halflevel -
+                        entrʲ_prev_level * u³ʲ_data_prev_halflevel
                     )
                 )
 
@@ -513,8 +536,8 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
                     local_geometry_prev_level.J *
                     ρaʲ_prev_level *
                     (
-                        entr_detrʲ_prev_level.entr * h_tot⁰_prev_level -
-                        entr_detrʲ_prev_level.detr * h_totʲ_prev_level +
+                        entrʲ_prev_level * h_tot⁰_prev_level -
+                        detrʲ_prev_level * h_totʲ_prev_level +
                         S_q_totʲ_prev_level * S_e_totʲ_helper_prev_level
                     )
                 )
@@ -541,8 +564,8 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
                     local_geometry_prev_level.J *
                     ρaʲ_prev_level *
                     (
-                        entr_detrʲ_prev_level.entr * q_tot_prev_level -
-                        entr_detrʲ_prev_level.detr * q_totʲ_prev_level +
+                        entrʲ_prev_level * q_tot_prev_level -
+                        detrʲ_prev_level * q_totʲ_prev_level +
                         S_q_totʲ_prev_level
                     )
                 )
@@ -604,15 +627,18 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
         ᶠu³ʲ = ᶠu³ʲs.:($j)
         ᶜuʲ = ᶜuʲs.:($j)
         ᶠu³ʲ = ᶠu³ʲs.:($j)
-        ᶜentr_detrʲ = ᶜentr_detrʲs.:($j)
+        ᶜentrʲ = ᶜentrʲs.:($j)
+        ᶜdetrʲ = ᶜdetrʲs.:($j)
         ᶜS_q_totʲ = ᶜS_q_totʲs.:($j)
         ᶜS_e_totʲ_helper = ᶜS_e_totʲs_helper.:($j)
 
         u³ʲ_halflevel = Fields.field_values(Fields.level(ᶠu³ʲ, i_top + half))
         @. u³ʲ_halflevel = CT3(0)
 
-        entr_detrʲ_level = Fields.field_values(Fields.level(ᶜentr_detrʲ, i_top))
-        fill!(entr_detrʲ_level, RecursiveApply.rzero(eltype(entr_detrʲ_level)))
+        entrʲ_level = Fields.field_values(Fields.level(ᶜentrʲ, i_top))
+        detrʲ_level = Fields.field_values(Fields.level(ᶜdetrʲ, i_top))
+        fill!(entrʲ_level, RecursiveApply.rzero(eltype(entrʲ_level)))
+        fill!(detrʲ_level, RecursiveApply.rzero(eltype(detrʲ_level)))
         @. ᶜuʲ = C123(Y.c.uₕ) + ᶜinterp(C123(ᶠu³ʲ))
 
         S_q_totʲ_level = Fields.field_values(Fields.level(ᶜS_q_totʲ, i_top))
@@ -671,7 +697,7 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
     # TKE equation, where using ᶜu⁰ results in allocation
     for j in 1:n
         @. ᶜtke_exch +=
-            ᶜρaʲs.:($$j) * ᶜentr_detrʲs.:($$j).detr / Y.c.ρ *
+            ᶜρaʲs.:($$j) * ᶜdetrʲs.:($$j) / Y.c.ρ *
             (1 / 2 * norm_sqr(ᶜinterp(ᶠu³⁰) - ᶜinterp(ᶠu³ʲs.:($$j))) - ᶜtke⁰)
     end
 

src/cache/edmf_precomputed_quantities.jl

@@ -25,7 +25,17 @@ function set_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
     (; ᶜspecific, ᶜp, ᶜΦ, ᶜh_tot, ᶜρ_ref) = p
     (; ᶜspecific⁰, ᶜρa⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰, ᶜh_tot⁰) = p
     (; ᶜmixing_length, ᶜlinear_buoygrad, ᶜstrain_rate_norm, ᶜK_u, ᶜK_h) = p
-    (; ᶜspecificʲs, ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜtsʲs, ᶜρʲs, ᶜh_totʲs, ᶜentr_detrʲs) = p
+    (;
+        ᶜspecificʲs,
+        ᶜuʲs,
+        ᶠu³ʲs,
+        ᶜKʲs,
+        ᶜtsʲs,
+        ᶜρʲs,
+        ᶜh_totʲs,
+        ᶜentrʲs,
+        ᶜdetrʲs,
+    ) = p
     (; ustar, obukhov_length, buoyancy_flux) = p.sfc_conditions
 
     @. ᶜspecific⁰ = specific_full_sgs⁰(Y.c, turbconv_model)
@@ -150,9 +160,8 @@ function set_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
     ᶜlg = Fields.local_geometry_field(Y.c)
 
     for j in 1:n
-        @. ᶜentr_detrʲs.:($$j) = pi_groups_entr_detr(
+        @. ᶜentrʲs.:($$j) = entrainment(
             params,
-            p.atmos.edmfx_entr_detr,
             ᶜz,
             z_sfc,
             ᶜp,
@@ -166,6 +175,24 @@ function set_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
             TD.relative_humidity(thermo_params, ᶜts⁰),
             ᶜphysical_buoyancy(params, ᶜρ_ref, ᶜρ⁰),
             dt,
+            p.atmos.edmfx_entr_model,
+        )
+        @. ᶜdetrʲs.:($$j) = detrainment(
+            params,
+            ᶜz,
+            z_sfc,
+            ᶜp,
+            Y.c.ρ,
+            buoyancy_flux,
+            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, ᶜρ_ref, ᶜρʲs.:($$j)),
+            get_physical_w(ᶜu⁰, ᶜlg),
+            TD.relative_humidity(thermo_params, ᶜts⁰),
+            ᶜphysical_buoyancy(params, ᶜρ_ref, ᶜρ⁰),
+            dt,
+            p.atmos.edmfx_detr_model,
         )
     end
 
@@ -204,7 +231,7 @@ function set_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
     @. ᶜtke_exch = 0
     for j in 1:n
         @. ᶜtke_exch +=
-            Y.c.sgsʲs.:($$j).ρa * ᶜentr_detrʲs.:($$j).detr / ᶜρa⁰ * (
+            Y.c.sgsʲs.:($$j).ρa * ᶜdetrʲs.:($$j) / ᶜρa⁰ * (
                 1 / 2 *
                 (
                     get_physical_w(ᶜuʲs.:($$j), ᶜlg) - get_physical_w(ᶜu⁰, ᶜlg)

src/cache/precomputed_quantities.jl

@@ -78,10 +78,8 @@ function precomputed_quantities(Y, atmos)
             ᶜKʲs = similar(Y.c, NTuple{n, FT}),
             ᶜtsʲs = similar(Y.c, NTuple{n, TST}),
             ᶜρʲs = similar(Y.c, NTuple{n, FT}),
-            ᶜentr_detrʲs = similar(
-                Y.c,
-                NTuple{n, NamedTuple{(:entr, :detr), NTuple{2, FT}}},
-            ),
+            ᶜentrʲs = similar(Y.c, NTuple{n, FT}),
+            ᶜdetrʲs = similar(Y.c, NTuple{n, FT}),
             (
                 atmos.energy_form isa TotalEnergy ?
                 (;
@@ -101,10 +99,8 @@ function precomputed_quantities(Y, atmos)
             ᶜρʲs = similar(Y.c, NTuple{n, FT}),
             ᶜh_totʲs = similar(Y.c, NTuple{n, FT}),
             ᶜq_totʲs = similar(Y.c, NTuple{n, FT}),
-            ᶜentr_detrʲs = similar(
-                Y.c,
-                NTuple{n, NamedTuple{(:entr, :detr), NTuple{2, FT}}},
-            ),
+            ᶜentrʲs = similar(Y.c, NTuple{n, FT}),
+            ᶜdetrʲs = similar(Y.c, NTuple{n, FT}),
             ᶜnh_pressureʲs = similar(Y.c, NTuple{n, CT3{FT}}),
             ᶜS_q_totʲs = similar(Y.c, NTuple{n, FT}),
             ᶜS_q_tot⁰ = similar(Y.c, FT),