0-moment micro with diag EDMFX

CliMA · Aug 28, 2023 · e610a5b · e610a5b
1 parent 1ecbf6d
commit e610a5b
Show file tree

Hide file tree

Showing 9 changed files with 235 additions and 70 deletions.
diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
@@ -194,7 +194,7 @@ steps:
 
       - label: ":computer: hydrostatic balance (ρe) float64"
         command: >
-          julia --color=yes --project=examples examples/hybrid/driver.jl 
+          julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_hydrostatic_balance_rhoe_ft64.yml
 
           julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
@@ -449,7 +449,7 @@ steps:
 
       - label: ":computer: held suarez (ρe) equilmoist topography (dcmip)"
         command: >
-          julia --color=yes --project=examples examples/hybrid/driver.jl 
+          julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_held_suarez_rhoe_equilmoist_topography_dcmip.yml
 
           julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
@@ -464,8 +464,8 @@ steps:
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_ssp_baroclinic_wave_rhoe_equilmoist_dcmip200.yml
 
-          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl 
-          --data_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_dcmip200 
+          julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
+          --data_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_dcmip200
           --out_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_dcmip200
 
           julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl
@@ -478,14 +478,14 @@ steps:
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth.yml
-          
+
           julia --color=yes --project=examples post_processing/remap/remap_pipeline.jl
-          --data_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth 
+          --data_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth
           --out_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth
-          
-          julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl 
-          --nc_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth 
-          --fig_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth 
+
+          julia --color=yes --project=examples post_processing/plot/plot_pipeline.jl
+          --nc_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth
+          --fig_dir sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth
           --case_name elevation_spectrum
         artifact_paths: "sphere_ssp_baroclinic_wave_rhoe_equilmoist_earth/*"
 
@@ -616,7 +616,7 @@ steps:
         artifact_paths: "diagnostic_edmfx_test_box/*"
         agents:
           slurm_mem: 20GB
-      
+
       - label: ":genie: Diagnostic EDMFX GABLS in a box"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
@@ -632,7 +632,7 @@ steps:
         artifact_paths: "diagnostic_edmfx_bomex_box/*"
         agents:
           slurm_mem: 20GB
-      
+
       - label: ":genie: Diagnostic EDMFX Bomex stretched grid in a box"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
@@ -656,7 +656,7 @@ steps:
         artifact_paths: "diagnostic_edmfx_rico_box/*"
         agents:
           slurm_mem: 20GB
-      
+
       - label: ":genie: Diagnostic EDMFX TRMM in a box"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
@@ -744,7 +744,7 @@ steps:
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/edmf_nieuwstadt_jfnk_imex.yml
-          
+
           julia --color=yes --project=examples regression_tests/test_mse.jl
           --job_id edmf_nieuwstadt_jfnk_imex --out_dir edmf_nieuwstadt_jfnk_imex
         artifact_paths: "edmf_nieuwstadt_jfnk_imex/*"
@@ -863,7 +863,7 @@ steps:
         key: "gpu_baroclinic_wave_rhoe"
         command: >
           julia --project -e 'using CUDA; CUDA.versioninfo()'
-          
+
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $PERF_CONFIG_PATH/gpu_baroclinic_wave_rhoe.yml
         artifact_paths: "gpu_implicit_barowave_ref/*"
@@ -884,7 +884,7 @@ steps:
         key: "gpu_held_suarez_rhoe_hightop"
         command: >
           julia --project -e 'using CUDA; CUDA.versioninfo()'
-          
+
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $PERF_CONFIG_PATH/gpu_held_suarez_rhoe_hightop.yml
         artifact_paths: "gpu_held_suarez_rhoe_hightop/*"

diff --git a/config/model_configs/diagnostic_edmfx_rico_box.yml b/config/model_configs/diagnostic_edmfx_rico_box.yml
@@ -1,27 +1,27 @@
 
 job_id: diagnostic_edmfx_rico_box
-initial_condition: Rico 
-subsidence: Rico 
-ls_adv: Rico 
+initial_condition: Rico
+subsidence: Rico
+ls_adv: Rico
 surface_setup: Rico
-turbconv: diagnostic_edmfx 
+turbconv: diagnostic_edmfx
 prognostic_tke: true
-edmfx_upwinding: first_order 
-edmfx_entr_detr: true 
-edmfx_nh_pressure: true 
+edmfx_upwinding: first_order
+edmfx_entr_detr: true
+edmfx_nh_pressure: true
 edmfx_sgs_flux: true
 moist: equil
-precip_model: "0M" 
-config: box 
+precip_model: "0M"
+config: box
 hyperdiff: "true"
 kappa_4: 1e12
-x_max: 1e5 
-y_max: 1e5 
-x_elem: 2 
-y_elem: 2 
-z_elem: 100 
-z_max: 4e3 
+x_max: 1e5
+y_max: 1e5
+x_elem: 2
+y_elem: 2
+z_elem: 100
+z_max: 4e3
 z_stretch: false
-dt: 10secs 
-t_end: 6hours 
-dt_save_to_disk: 10mins
+dt: 5secs
+t_end: 6hours
+dt_save_to_disk: 10mins
diff --git a/config/model_configs/diagnostic_edmfx_trmm_box.yml b/config/model_configs/diagnostic_edmfx_trmm_box.yml
@@ -1,28 +1,28 @@
 
 job_id: diagnostic_edmfx_trmm_box
-initial_condition: TRMM_LBA 
-rad: TRMM_LBA 
+initial_condition: TRMM_LBA
+rad: TRMM_LBA
 surface_setup: TRMM_LBA
-turbconv: diagnostic_edmfx 
+turbconv: diagnostic_edmfx
 prognostic_tke: true
-edmfx_upwinding: first_order 
-edmfx_entr_detr: true 
-edmfx_nh_pressure: true 
+edmfx_upwinding: first_order
+edmfx_entr_detr: true
+edmfx_nh_pressure: true
 edmfx_sgs_flux: true
 moist: equil
 apply_limiter: false
 precip_model: "0M"
-config: box 
+config: box
 hyperdiff: "true"
 kappa_4: 1e12
-x_max: 1e5 
-y_max: 1e5 
-x_elem: 2 
-y_elem: 2 
-z_elem: 82 
-z_max: 16400 
+x_max: 1e5
+y_max: 1e5
+x_elem: 2
+y_elem: 2
+z_elem: 82
+z_max: 16400
 z_stretch: false
-dt: 10secs 
-t_end: 6hours 
-dt_save_to_disk: 10mins 
-FLOAT_TYPE: "Float64"
+dt: 10secs
+t_end: 6hours
+dt_save_to_disk: 10mins
+FLOAT_TYPE: "Float64"
diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl
@@ -87,6 +87,13 @@ include(joinpath("prognostic_equations", "cloud_fraction.jl"))
 include(
     joinpath("parameterized_tendencies", "microphysics", "precipitation.jl"),
 )
+include(
+    joinpath(
+        "parameterized_tendencies",
+        "microphysics",
+        "microphysics_wrappers.jl",
+    ),
+)
 include(
     joinpath("prognostic_equations", "vertical_diffusion_boundary_layer.jl"),
 )

diff --git a/src/TurbulenceConvection_deprecated/closures/microphysics.jl b/src/TurbulenceConvection_deprecated/closures/microphysics.jl
@@ -54,7 +54,9 @@ function precipitation_formation(
             qi_tendency += S_qt * (1 - λ)
             θ_liq_ice_tendency -=
                 S_qt / Π_m / c_pm * (L_v0 * λ + L_s0 * (1 - λ))
-            e_tot_tendency += (λ * I_l + (1 - λ) * I_i + Φ) * S_qt
+            e_tot_tendency +=
+                S_qt *
+                e_tot_0M_precipitation_sources_helper(thermo_params, ts, Φ)
         end
 
         if precip_model isa Microphysics1Moment

diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl
@@ -83,7 +83,7 @@ end
 Updates the precomputed quantities stored in `p` for diagnostic edmfx.
 """
 function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
-    (; moisture_model, turbconv_model) = p.atmos
+    (; moisture_model, turbconv_model, precip_model) = p.atmos
     FT = eltype(Y)
     n = n_mass_flux_subdomains(turbconv_model)
     ᶜz = Fields.coordinate_field(Y.c).z
@@ -105,10 +105,14 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
         ᶜρʲs,
         ᶜentr_detrʲs,
         ᶜnh_pressureʲs,
+        ᶜS_q_totʲs,
+        ᶜS_e_totʲs_helper,
     ) = p
-    (; ᶠu³⁰, ᶜu⁰, ᶜtke⁰, ᶜlinear_buoygrad, ᶜshear², ᶜmixing_length) = p
+    (; ᶠu³⁰, ᶜu⁰, ᶜtke⁰, ᶜlinear_buoygrad, ᶜshear², ᶜmixing_length, ᶜS_q_tot⁰) =
+        p
     (; ᶜK_h, ᶜK_u, ρatke_flux) = p
     thermo_params = CAP.thermodynamics_params(params)
+    microphys_params = CAP.microphysics_params(params)
     ᶜlg = Fields.local_geometry_field(Y.c)
 
     @. ᶜtke⁰ = Y.c.sgs⁰.ρatke / Y.c.ρ
@@ -237,6 +241,7 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
         end_index = fieldcount(eltype(∂x∂ξ_level)) # This will be 4 in 2D and 9 in 3D.
         ∂x³∂ξ³_level = ∂x∂ξ_level.:($end_index)
 
+        Φ_prev_level = Fields.field_values(Fields.level(ᶜΦ, i - 1))
         ρ_ref_prev_level = Fields.field_values(Fields.level(ᶜρ_ref, i - 1))
         ∇Φ³_prev_level = Fields.field_values(Fields.level(ᶜ∇Φ³, i - 1))
         ∇Φ³_data_prev_level = ∇Φ³_prev_level.components.data.:1
@@ -267,6 +272,8 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
             ᶜq_totʲ = ᶜq_totʲs.:($j)
             ᶜentr_detrʲ = ᶜentr_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)
 
             ρaʲ_level = Fields.field_values(Fields.level(ᶜρaʲ, i))
             u³ʲ_halflevel = Fields.field_values(Fields.level(ᶠu³ʲ, i - half))
@@ -292,6 +299,10 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
                 Fields.field_values(Fields.level(ᶜnh_pressureʲ, i - 1))
             scale_height =
                 CAP.R_d(params) * CAP.T_surf_ref(params) / CAP.grav(params)
+            S_q_totʲ_prev_level =
+                Fields.field_values(Fields.level(ᶜS_q_totʲ, i - 1))
+            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(
                 params,
@@ -331,6 +342,18 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
             nh_pressureʲ_data_prev_level =
                 nh_pressureʲ_prev_level.components.data.:1
 
+            # Updraft q_tot sources from precipitation formation
+            # To be applied in updraft continuity, moisture and energy
+            # for updrafts and grid mean
+            @. S_q_totʲ_prev_level = q_tot_precipitation_sources(
+                precip_model,
+                thermo_params,
+                microphys_params,
+                dt,
+                q_totʲ_prev_level,
+                tsʲ_prev_level,
+            )
+
             @. ρaʲu³ʲ_data =
                 (1 / local_geometry_halflevel.J) * (
                     local_geometry_prev_halflevel.J *
@@ -342,7 +365,10 @@ 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)
+                    (
+                        entr_detrʲ_prev_level.entr -
+                        entr_detrʲ_prev_level.detr + S_q_totʲ_prev_level
+                    )
                 )
 
             # Using constant exponents in broadcasts allocate, so we use
@@ -420,6 +446,12 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
                 ρaʲu³ʲ_data / sqrt(max(0, u³ʲ_datau³ʲ_data)),
             )
 
+            @. S_e_totʲ_helper_prev_level =
+                e_tot_0M_precipitation_sources_helper(
+                    thermo_params,
+                    tsʲ_prev_level,
+                    Φ_prev_level,
+                )
             @. ρaʲu³ʲ_datah_tot =
                 (1 / local_geometry_halflevel.J) * (
                     local_geometry_prev_halflevel.J *
@@ -433,7 +465,8 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
                     ρaʲ_prev_level *
                     (
                         entr_detrʲ_prev_level.entr * h_tot_prev_level -
-                        entr_detrʲ_prev_level.detr * h_totʲ_prev_level
+                        entr_detrʲ_prev_level.detr * h_totʲ_prev_level +
+                        S_q_totʲ_prev_level * S_e_totʲ_helper_prev_level
                     )
                 )
             @. h_totʲ_level = ifelse(
@@ -460,7 +493,8 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
                     ρaʲ_prev_level *
                     (
                         entr_detrʲ_prev_level.entr * q_tot_prev_level -
-                        entr_detrʲ_prev_level.detr * q_totʲ_prev_level
+                        entr_detrʲ_prev_level.detr * q_totʲ_prev_level +
+                        S_q_totʲ_prev_level
                     )
                 )
             @. q_totʲ_level = ifelse(
@@ -507,16 +541,27 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
     i_top = Spaces.nlevels(axes(Y.c))
     u³⁰_halflevel = Fields.field_values(Fields.level(ᶠu³⁰, i_top + half))
     @. u³⁰_halflevel = CT3(0)
+
     for j in 1:n
         ᶠu³ʲ = ᶠu³ʲs.:($j)
         ᶜuʲ = ᶜuʲs.:($j)
         ᶠu³ʲ = ᶠu³ʲs.:($j)
         ᶜentr_detrʲ = ᶜentr_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)))
         @. ᶜuʲ = C123(Y.c.uₕ) + ᶜinterp(C123(ᶠu³ʲ))
+
+        S_q_totʲ_level = Fields.field_values(Fields.level(ᶜS_q_totʲ, i_top))
+        S_e_totʲ_helper_level =
+            Fields.field_values(Fields.level(ᶜS_e_totʲ_helper, i_top))
+        @. S_q_totʲ_level = FT(0)
+        @. S_e_totʲ_helper_level = FT(0)
     end
 
     @. ᶜu⁰ = C123(Y.c.uₕ) + ᶜinterp(C123(ᶠu³⁰))
@@ -609,5 +654,15 @@ function set_diagnostic_edmf_precomputed_quantities!(Y, p, t)
         sfc_local_geometry_values,
     )
 
+    # Environment precipitation sources (to be applied to grid mean)
+    @. S_q_tot⁰ = q_tot_precipitation_sources(
+        precip_model,
+        thermo_params,
+        microphys_params,
+        dt,
+        q_tot,
+        ᶜts,
+    )
+
     return nothing
 end
diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
@@ -106,6 +106,9 @@ function precomputed_quantities(Y, atmos)
                 NTuple{n, NamedTuple{(:entr, :detr), NTuple{2, 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),
+            ᶜS_e_totʲs_helper = similar(Y.c, NTuple{n, FT}),
             ᶠu³⁰ = similar(Y.f, CT3{FT}),
             ᶜu⁰ = similar(Y.c, C123{FT}),
             ᶜtke⁰ = similar(Y.c, FT),