diff --git a/NEWS.md b/NEWS.md
index 59ed44c271..e4979d9ab7 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,13 @@ ClimaAtmos.jl Release Notes
 Main
 -------
 
+
+v0.26.0
+-------
+- ![][badge-💥breaking] Add precipitation fluxes to 1M microphysics output.
+  Rename col_integrated_rain (and snow) to surface_rain_flux (and snow)
+  PR [#3084](https://github.com/CliMA/ClimaAtmos.jl/pull/3084).
+
 v0.25.0
 -------
 - ![][badge-💥breaking] Remove reference state from the dycore and the
diff --git a/Project.toml b/Project.toml
index 1bab3bfeb5..e6b0d66690 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ClimaAtmos"
 uuid = "b2c96348-7fb7-4fe0-8da9-78d88439e717"
 authors = ["Climate Modeling Alliance"]
-version = "0.25.0"
+version = "0.26.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
diff --git a/config/model_configs/diagnostic_edmfx_bomex_box.yml b/config/model_configs/diagnostic_edmfx_bomex_box.yml
index 5ae394a8ab..37b294102c 100644
--- a/config/model_configs/diagnostic_edmfx_bomex_box.yml
+++ b/config/model_configs/diagnostic_edmfx_bomex_box.yml
@@ -1,35 +1,35 @@
-initial_condition: "Bomex" 
-subsidence: "Bomex" 
-edmf_coriolis: "Bomex" 
-ls_adv: "Bomex" 
+initial_condition: "Bomex"
+subsidence: "Bomex"
+edmf_coriolis: "Bomex"
+ls_adv: "Bomex"
 surface_setup: "Bomex"
-turbconv: "diagnostic_edmfx" 
+turbconv: "diagnostic_edmfx"
 implicit_diffusion: true
 approximate_linear_solve_iters: 2
 prognostic_tke: true
-edmfx_upwinding: "first_order" 
-edmfx_entr_model: "Generalized" 
-edmfx_detr_model: "Generalized" 
-edmfx_nh_pressure: true 
+edmfx_upwinding: "first_order"
+edmfx_entr_model: "Generalized"
+edmfx_detr_model: "Generalized"
+edmfx_nh_pressure: true
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 moist: "equil"
-config: "box" 
-hyperdiff: "true" 
-x_max: 1e8 
-y_max: 1e8 
-x_elem: 2 
-y_elem: 2 
-z_elem: 60 
-z_max: 3e3 
+config: "box"
+hyperdiff: "true"
+x_max: 1e8
+y_max: 1e8
+x_elem: 2
+y_elem: 2
+z_elem: 60
+z_max: 3e3
 z_stretch: false
-dt: "100secs" 
-t_end: "6hours" 
+dt: "100secs"
+t_end: "6hours"
 dt_save_state_to_disk: "10mins"
 toml: [toml/diagnostic_edmfx.toml]
 netcdf_interpolation_num_points: [8, 8, 60]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, tke, lmix]
     period: 10mins
diff --git a/config/model_configs/diagnostic_edmfx_bomex_stretched_box.yml b/config/model_configs/diagnostic_edmfx_bomex_stretched_box.yml
index a4225b5431..b3fbd64cd4 100644
--- a/config/model_configs/diagnostic_edmfx_bomex_stretched_box.yml
+++ b/config/model_configs/diagnostic_edmfx_bomex_stretched_box.yml
@@ -1,36 +1,36 @@
-initial_condition: "Bomex" 
-subsidence: "Bomex" 
-edmf_coriolis: "Bomex" 
-ls_adv: "Bomex" 
+initial_condition: "Bomex"
+subsidence: "Bomex"
+edmf_coriolis: "Bomex"
+ls_adv: "Bomex"
 surface_setup: "Bomex"
-turbconv: "diagnostic_edmfx" 
+turbconv: "diagnostic_edmfx"
 implicit_diffusion: true
 approximate_linear_solve_iters: 2
 prognostic_tke: true
-edmfx_upwinding: "first_order" 
-edmfx_entr_model: "Generalized" 
-edmfx_detr_model: "Generalized" 
-edmfx_nh_pressure: true 
+edmfx_upwinding: "first_order"
+edmfx_entr_model: "Generalized"
+edmfx_detr_model: "Generalized"
+edmfx_nh_pressure: true
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 moist: "equil"
-config: "box" 
-hyperdiff: "true" 
-x_max: 1e8 
-y_max: 1e8 
-x_elem: 2 
-y_elem: 2 
-z_elem: 30 
-z_max: 3e3 
+config: "box"
+hyperdiff: "true"
+x_max: 1e8
+y_max: 1e8
+x_elem: 2
+y_elem: 2
+z_elem: 30
+z_max: 3e3
 dz_bottom: 50.0
 dz_top: 200.0
-dt: "200secs" 
-t_end: "6hours" 
+dt: "200secs"
+t_end: "6hours"
 dt_save_state_to_disk: "10mins"
 toml: [toml/diagnostic_edmfx.toml]
 netcdf_interpolation_num_points: [8, 8, 30]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, tke, lmix]
     period: 10mins
diff --git a/config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml b/config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml
index 2a30397531..2e954cf50e 100644
--- a/config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml
+++ b/config/model_configs/diagnostic_edmfx_dycoms_rf01_box.yml
@@ -1,35 +1,35 @@
-initial_condition: DYCOMS_RF01 
-subsidence: DYCOMS 
-edmf_coriolis: DYCOMS_RF01 
-rad: DYCOMS_RF01 
+initial_condition: DYCOMS_RF01
+subsidence: DYCOMS
+edmf_coriolis: DYCOMS_RF01
+rad: DYCOMS_RF01
 surface_setup: DYCOMS_RF01
-turbconv: diagnostic_edmfx 
+turbconv: diagnostic_edmfx
 implicit_diffusion: true
 approximate_linear_solve_iters: 2
 prognostic_tke: true
-edmfx_upwinding: first_order 
-edmfx_entr_model: "Generalized" 
-edmfx_detr_model: "Generalized" 
-edmfx_nh_pressure: true 
+edmfx_upwinding: first_order
+edmfx_entr_model: "Generalized"
+edmfx_detr_model: "Generalized"
+edmfx_nh_pressure: true
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 moist: equil
-config: box 
-hyperdiff: "true" 
-x_max: 1e8 
-y_max: 1e8 
-x_elem: 2 
-y_elem: 2 
-z_elem: 30 
-z_max: 1500 
+config: box
+hyperdiff: "true"
+x_max: 1e8
+y_max: 1e8
+x_elem: 2
+y_elem: 2
+z_elem: 30
+z_max: 1500
 z_stretch: false
 dt: 100secs
-t_end: 4hours 
+t_end: 4hours
 dt_save_state_to_disk: 10mins
 toml: [toml/diagnostic_edmfx.toml]
 netcdf_interpolation_num_points: [8, 8, 30]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, tke, lmix]
     period: 10mins
diff --git a/config/model_configs/diagnostic_edmfx_dycoms_rf01_explicit_box.yml b/config/model_configs/diagnostic_edmfx_dycoms_rf01_explicit_box.yml
index 700b9445dc..3754110fd1 100644
--- a/config/model_configs/diagnostic_edmfx_dycoms_rf01_explicit_box.yml
+++ b/config/model_configs/diagnostic_edmfx_dycoms_rf01_explicit_box.yml
@@ -1,33 +1,33 @@
-initial_condition: DYCOMS_RF01 
-subsidence: DYCOMS 
-edmf_coriolis: DYCOMS_RF01 
-rad: DYCOMS_RF01 
+initial_condition: DYCOMS_RF01
+subsidence: DYCOMS
+edmf_coriolis: DYCOMS_RF01
+rad: DYCOMS_RF01
 surface_setup: DYCOMS_RF01
-turbconv: diagnostic_edmfx 
+turbconv: diagnostic_edmfx
 prognostic_tke: true
-edmfx_upwinding: first_order 
-edmfx_entr_model: "Generalized" 
-edmfx_detr_model: "Generalized" 
-edmfx_nh_pressure: true 
+edmfx_upwinding: first_order
+edmfx_entr_model: "Generalized"
+edmfx_detr_model: "Generalized"
+edmfx_nh_pressure: true
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 moist: equil
-config: box 
-hyperdiff: "true" 
-x_max: 1e8 
-y_max: 1e8 
-x_elem: 2 
-y_elem: 2 
-z_elem: 30 
-z_max: 1500 
+config: box
+hyperdiff: "true"
+x_max: 1e8
+y_max: 1e8
+x_elem: 2
+y_elem: 2
+z_elem: 30
+z_max: 1500
 z_stretch: false
-dt: 20secs 
-t_end: 4hours 
+dt: 20secs
+t_end: 4hours
 dt_save_state_to_disk: 10mins
 toml: [toml/diagnostic_edmfx.toml]
 netcdf_interpolation_num_points: [8, 8, 30]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, tke, lmix]
     period: 10mins
diff --git a/config/model_configs/diagnostic_edmfx_gabls_box.yml b/config/model_configs/diagnostic_edmfx_gabls_box.yml
index 1b2be79983..8ccb60b1ba 100644
--- a/config/model_configs/diagnostic_edmfx_gabls_box.yml
+++ b/config/model_configs/diagnostic_edmfx_gabls_box.yml
@@ -6,9 +6,9 @@ implicit_diffusion: true
 approximate_linear_solve_iters: 2
 prognostic_tke: true
 edmfx_upwinding: first_order
-edmfx_entr_model: "Generalized" 
-edmfx_detr_model: "Generalized" 
-edmfx_nh_pressure: true 
+edmfx_entr_model: "Generalized"
+edmfx_detr_model: "Generalized"
+edmfx_nh_pressure: true
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
 moist: equil
@@ -27,7 +27,7 @@ dt_save_state_to_disk: 10mins
 toml: [toml/diagnostic_edmfx.toml]
 netcdf_interpolation_num_points: [8, 8, 8]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, tke, lmix]
     period: 10mins
diff --git a/config/model_configs/diagnostic_edmfx_rico_box.yml b/config/model_configs/diagnostic_edmfx_rico_box.yml
index 1e50ac5486..5ea89cb68c 100644
--- a/config/model_configs/diagnostic_edmfx_rico_box.yml
+++ b/config/model_configs/diagnostic_edmfx_rico_box.yml
@@ -30,7 +30,7 @@ dt_save_state_to_disk: 10mins
 toml: [toml/diagnostic_edmfx.toml]
 netcdf_interpolation_num_points: [8, 8, 100]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, tke, lmix]
     period: 10mins
diff --git a/config/model_configs/diagnostic_edmfx_test_box.yml b/config/model_configs/diagnostic_edmfx_test_box.yml
index 18eaf61205..f4879a20f5 100644
--- a/config/model_configs/diagnostic_edmfx_test_box.yml
+++ b/config/model_configs/diagnostic_edmfx_test_box.yml
@@ -14,7 +14,7 @@ moist: "equil"
 toml: [toml/diagnostic_edmfx.toml]
 netcdf_interpolation_num_points: [8, 8, 45]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, tke, lmix]
     period: 10mins
diff --git a/config/model_configs/diagnostic_edmfx_trmm_box.yml b/config/model_configs/diagnostic_edmfx_trmm_box.yml
index e2a519e512..7629585317 100644
--- a/config/model_configs/diagnostic_edmfx_trmm_box.yml
+++ b/config/model_configs/diagnostic_edmfx_trmm_box.yml
@@ -30,7 +30,7 @@ FLOAT_TYPE: "Float64"
 toml: [toml/diagnostic_edmfx.toml]
 netcdf_interpolation_num_points: [8, 8, 82]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, tke, lmix]
     period: 10mins
diff --git a/config/model_configs/diagnostic_edmfx_trmm_box_0M.yml b/config/model_configs/diagnostic_edmfx_trmm_box_0M.yml
index 2afe16ddd1..7307e4a4a0 100644
--- a/config/model_configs/diagnostic_edmfx_trmm_box_0M.yml
+++ b/config/model_configs/diagnostic_edmfx_trmm_box_0M.yml
@@ -32,7 +32,7 @@ toml: [toml/diagnostic_edmfx_0M.toml]
 netcdf_output_at_levels: true
 netcdf_interpolation_num_points: [8, 8, 82]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, tke, lmix]
     period: 10mins
diff --git a/config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml b/config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml
index e9e28398c7..035c052e8f 100644
--- a/config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml
+++ b/config/model_configs/diagnostic_edmfx_trmm_stretched_box.yml
@@ -33,7 +33,7 @@ FLOAT_TYPE: "Float64"
 toml: [toml/diagnostic_edmfx.toml]
 netcdf_interpolation_num_points: [8, 8, 30]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, tke, lmix]
     period: 10mins
diff --git a/config/model_configs/gpu_prognostic_edmfx_aquaplanet.yml b/config/model_configs/gpu_prognostic_edmfx_aquaplanet.yml
index 52eebd2d0b..4ac20f61d4 100644
--- a/config/model_configs/gpu_prognostic_edmfx_aquaplanet.yml
+++ b/config/model_configs/gpu_prognostic_edmfx_aquaplanet.yml
@@ -21,7 +21,7 @@ dt_save_to_disk: 600secs
 toml: [toml/prognostic_edmfx.toml]
 output_default_diagnostics: false
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hfes, hur, hus, cl, clw, cli, evspsbl, rsd, rsu, rld, rlu]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hfes, hur, hus, cl, clw, cli, evspsbl, pr, rsd, rsu, rld, rlu]
     reduction_time: average
     period: 1hours
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, tke, lmix]
diff --git a/config/model_configs/prognostic_edmfx_adv_test_column.yml b/config/model_configs/prognostic_edmfx_adv_test_column.yml
index a53d91d1a5..9fb8af18ee 100644
--- a/config/model_configs/prognostic_edmfx_adv_test_column.yml
+++ b/config/model_configs/prognostic_edmfx_adv_test_column.yml
@@ -18,7 +18,7 @@ FLOAT_TYPE: "Float64"
 toml: [toml/prognostic_edmfx_advection.toml]
 netcdf_interpolation_num_points: [2, 2, 63]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke, lmix]
     period: 10mins
diff --git a/config/model_configs/prognostic_edmfx_aquaplanet.yml b/config/model_configs/prognostic_edmfx_aquaplanet.yml
index f9895b9a0a..0c5d0975fe 100644
--- a/config/model_configs/prognostic_edmfx_aquaplanet.yml
+++ b/config/model_configs/prognostic_edmfx_aquaplanet.yml
@@ -1,23 +1,23 @@
-surface_setup: DefaultMoninObukhov 
+surface_setup: DefaultMoninObukhov
 rad: clearsky
-turbconv: prognostic_edmfx 
+turbconv: prognostic_edmfx
 prognostic_tke: true
-edmfx_upwinding: first_order 
-edmfx_entr_model: "Generalized" 
-edmfx_detr_model: "Generalized" 
-edmfx_nh_pressure: true 
+edmfx_upwinding: first_order
+edmfx_entr_model: "Generalized"
+edmfx_detr_model: "Generalized"
+edmfx_nh_pressure: true
 edmfx_filter: true
 edmfx_sgs_mass_flux: true
 edmfx_sgs_diffusive_flux: true
-moist: equil 
+moist: equil
 precip_model: 0M
-dt: 10secs 
-t_end: 1hours 
+dt: 10secs
+t_end: 1hours
 dt_save_state_to_disk: 600secs
 toml: [toml/prognostic_edmfx.toml]
 output_default_diagnostics: false
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hfes, hur, hus, cl, clw, cli, evspsbl, rsd, rsu, rld, rlu]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hfes, hur, hus, cl, clw, cli, evspsbl, pr, rsd, rsu, rld, rlu]
     reduction_time: average
     period: 1hours
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, tke, lmix]
diff --git a/config/model_configs/prognostic_edmfx_bomex_box.yml b/config/model_configs/prognostic_edmfx_bomex_box.yml
index 673b5d036c..4b7c339b17 100644
--- a/config/model_configs/prognostic_edmfx_bomex_box.yml
+++ b/config/model_configs/prognostic_edmfx_bomex_box.yml
@@ -33,7 +33,7 @@ dt_save_state_to_disk: "30mins"
 toml: [toml/prognostic_edmfx_bomex.toml]
 netcdf_interpolation_num_points: [8, 8, 60]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
     period: 10mins
diff --git a/config/model_configs/prognostic_edmfx_bomex_column.yml b/config/model_configs/prognostic_edmfx_bomex_column.yml
index f9114908a6..4ad5c00bbf 100644
--- a/config/model_configs/prognostic_edmfx_bomex_column.yml
+++ b/config/model_configs/prognostic_edmfx_bomex_column.yml
@@ -27,7 +27,7 @@ dt_save_to_disk: "10mins"
 toml: [toml/prognostic_edmfx_bomex.toml]
 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]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
     period: 10mins
diff --git a/config/model_configs/prognostic_edmfx_bomex_fixtke_column.yml b/config/model_configs/prognostic_edmfx_bomex_fixtke_column.yml
index 5bdecdabca..4d61f75c5d 100644
--- a/config/model_configs/prognostic_edmfx_bomex_fixtke_column.yml
+++ b/config/model_configs/prognostic_edmfx_bomex_fixtke_column.yml
@@ -26,7 +26,7 @@ dt_save_state_to_disk: "10mins"
 toml: [toml/prognostic_edmfx_bomex.toml]
 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]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
     period: 10mins
diff --git a/config/model_configs/prognostic_edmfx_bomex_implicit_column.yml b/config/model_configs/prognostic_edmfx_bomex_implicit_column.yml
index dbfab47ba7..1517ea6d81 100644
--- a/config/model_configs/prognostic_edmfx_bomex_implicit_column.yml
+++ b/config/model_configs/prognostic_edmfx_bomex_implicit_column.yml
@@ -28,7 +28,7 @@ dt_save_to_disk: "10mins"
 toml: [toml/prognostic_edmfx_bomex.toml]
 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]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
     period: 10mins
diff --git a/config/model_configs/prognostic_edmfx_bomex_pigroup_column.yml b/config/model_configs/prognostic_edmfx_bomex_pigroup_column.yml
index 681af9e256..3cc3f60f16 100644
--- a/config/model_configs/prognostic_edmfx_bomex_pigroup_column.yml
+++ b/config/model_configs/prognostic_edmfx_bomex_pigroup_column.yml
@@ -25,7 +25,7 @@ dt_save_to_disk: "10mins"
 toml: [toml/prognostic_edmfx_bomex_pigroup.toml]
 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]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
     period: 10mins
diff --git a/config/model_configs/prognostic_edmfx_bomex_stretched_column.yml b/config/model_configs/prognostic_edmfx_bomex_stretched_column.yml
index 840611a78e..47a5fa9b96 100644
--- a/config/model_configs/prognostic_edmfx_bomex_stretched_column.yml
+++ b/config/model_configs/prognostic_edmfx_bomex_stretched_column.yml
@@ -31,7 +31,7 @@ dt_save_state_to_disk: "10mins"
 toml: [toml/prognostic_edmfx_bomex.toml]
 netcdf_interpolation_num_points: [2, 2, 30]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
     period: 10mins
diff --git a/config/model_configs/prognostic_edmfx_dycoms_rf01_column.yml b/config/model_configs/prognostic_edmfx_dycoms_rf01_column.yml
index 718e7df823..fd68f3cf5b 100644
--- a/config/model_configs/prognostic_edmfx_dycoms_rf01_column.yml
+++ b/config/model_configs/prognostic_edmfx_dycoms_rf01_column.yml
@@ -30,7 +30,7 @@ dt_save_state_to_disk: 10mins
 toml: [toml/prognostic_edmfx.toml]
 netcdf_interpolation_num_points: [2, 2, 30]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
     period: 10mins
diff --git a/config/model_configs/prognostic_edmfx_gabls_column.yml b/config/model_configs/prognostic_edmfx_gabls_column.yml
index 4aa7374c2e..cdac772292 100644
--- a/config/model_configs/prognostic_edmfx_gabls_column.yml
+++ b/config/model_configs/prognostic_edmfx_gabls_column.yml
@@ -28,7 +28,7 @@ perturb_initstate: false
 toml: [toml/prognostic_edmfx.toml]
 netcdf_interpolation_num_points: [2, 2, 8]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
     period: 10mins
diff --git a/config/model_configs/prognostic_edmfx_gcmdriven_column.yml b/config/model_configs/prognostic_edmfx_gcmdriven_column.yml
index d3aa32a2a2..d5a51d2c26 100644
--- a/config/model_configs/prognostic_edmfx_gcmdriven_column.yml
+++ b/config/model_configs/prognostic_edmfx_gcmdriven_column.yml
@@ -27,7 +27,7 @@ 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]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
     period: 10mins
diff --git a/config/model_configs/prognostic_edmfx_rico_column.yml b/config/model_configs/prognostic_edmfx_rico_column.yml
index b2e3ceb7d4..5c13a7a3f1 100644
--- a/config/model_configs/prognostic_edmfx_rico_column.yml
+++ b/config/model_configs/prognostic_edmfx_rico_column.yml
@@ -29,7 +29,7 @@ dt_save_state_to_disk: "10mins"
 toml: [toml/prognostic_edmfx.toml]
 netcdf_interpolation_num_points: [8, 8, 100]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
     period: 10mins
diff --git a/config/model_configs/prognostic_edmfx_simpleplume_column.yml b/config/model_configs/prognostic_edmfx_simpleplume_column.yml
index 1acc5a0ef5..f0cc8f2f1a 100644
--- a/config/model_configs/prognostic_edmfx_simpleplume_column.yml
+++ b/config/model_configs/prognostic_edmfx_simpleplume_column.yml
@@ -26,7 +26,7 @@ dt_save_to_disk: "10mins"
 toml: [toml/prognostic_edmfx_simpleplume.toml]
 netcdf_interpolation_num_points: [2, 2, 80]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
     period: 10mins
diff --git a/config/model_configs/prognostic_edmfx_trmm_column.yml b/config/model_configs/prognostic_edmfx_trmm_column.yml
index 59383b8de6..c953ea5847 100644
--- a/config/model_configs/prognostic_edmfx_trmm_column.yml
+++ b/config/model_configs/prognostic_edmfx_trmm_column.yml
@@ -29,7 +29,7 @@ FLOAT_TYPE: "Float64"
 toml: [toml/prognostic_edmfx.toml]
 netcdf_interpolation_num_points: [8, 8, 82]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
     period: 10mins
diff --git a/config/model_configs/prognostic_edmfx_trmm_column_0M.yml b/config/model_configs/prognostic_edmfx_trmm_column_0M.yml
index 384b2df555..9f78337ab3 100644
--- a/config/model_configs/prognostic_edmfx_trmm_column_0M.yml
+++ b/config/model_configs/prognostic_edmfx_trmm_column_0M.yml
@@ -29,7 +29,7 @@ FLOAT_TYPE: "Float64"
 toml: [toml/prognostic_edmfx.toml]
 netcdf_interpolation_num_points: [8, 8, 82]
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
     period: 10mins
   - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke, lmix]
     period: 10mins
diff --git a/config/model_configs/single_column_precipitation_test.yml b/config/model_configs/single_column_precipitation_test.yml
index f105d951d0..28e2d0a5d7 100644
--- a/config/model_configs/single_column_precipitation_test.yml
+++ b/config/model_configs/single_column_precipitation_test.yml
@@ -18,3 +18,5 @@ toml: [toml/single_column_precipitation_test.toml]
 diagnostics:
   - short_name: [hus, clw, cli, husra, hussn, ta, wa]
     period: 500secs
+  - short_name: [pr, ]
+    period: 10secs
diff --git a/post_processing/ci_plots.jl b/post_processing/ci_plots.jl
index 6dc05468c5..ddc069924b 100644
--- a/post_processing/ci_plots.jl
+++ b/post_processing/ci_plots.jl
@@ -42,7 +42,7 @@ import CairoMakie
 import CairoMakie.Makie
 import ClimaAnalysis
 import ClimaAnalysis: Visualize as viz
-import ClimaAnalysis: SimDir, slice
+import ClimaAnalysis: SimDir, slice, read_var
 import ClimaAnalysis.Utils: kwargs as ca_kwargs
 
 import ClimaCoreSpectra: power_spectrum_2d
@@ -515,7 +515,6 @@ function make_plots(
     ::Val{:single_column_precipitation_test},
     output_paths::Vector{<:AbstractString},
 )
-
     simdirs = SimDir.(output_paths)
 
     # TODO: Move this plotting code into the same framework as the other ones
@@ -564,6 +563,14 @@ function make_plots(
         end
     end
 
+    # surface_precipitation
+    surface_precip = read_var(simdir.variable_paths["pr"]["inst"]["10s"])
+    viz.line_plot1D!(
+        fig,
+        slice(surface_precip, x = 0.0, y = 0.0);
+        p_loc = [3, 1:3],
+    )
+
     file_path = joinpath(output_paths[1], "summary.pdf")
     CairoMakie.save(file_path, fig)
 end
@@ -822,7 +829,7 @@ function make_plots(
     simdirs = SimDir.(output_paths)
 
     short_names_3D, reduction = ["ua", "ta", "hus"], "average"
-    short_names_2D = ["hfes", "evspsbl"]
+    short_names_2D = ["hfes", "evspsbl", "pr"]
     vars_3D = map_comparison(simdirs, short_names_3D) do simdir, short_name
         get(simdir; short_name, reduction) |> ClimaAnalysis.average_lon
     end
@@ -861,6 +868,7 @@ function make_plots(
         "rlus",
         "hfes",
         "evspsbl",
+        "pr",
     ]
     available_periods = ClimaAnalysis.available_periods(
         simdirs[1];
@@ -935,6 +943,7 @@ function make_plots(
         "hfes",
         "evspsbl",
         "ts",
+        "pr",
     ]
     vars_3D = map_comparison(simdirs, short_names_3D) do simdir, short_name
         get(simdir; short_name, reduction, period) |> ClimaAnalysis.average_lon
@@ -984,6 +993,7 @@ function make_plots(::AquaplanetPlots, output_paths::Vector{<:AbstractString})
         "rlus",
         "hfes",
         "evspsbl",
+        "pr",
     ]
     available_periods = ClimaAnalysis.available_periods(
         simdirs[1];
@@ -1041,6 +1051,7 @@ function make_plots(::Aquaplanet1MPlots, output_paths::Vector{<:AbstractString})
         "rlus",
         "hfes",
         "evspsbl",
+        "pr",
     ]
     available_periods = ClimaAnalysis.available_periods(
         simdirs[1];
diff --git a/src/diagnostics/Diagnostics.jl b/src/diagnostics/Diagnostics.jl
index 04d01941b0..3f2094286e 100644
--- a/src/diagnostics/Diagnostics.jl
+++ b/src/diagnostics/Diagnostics.jl
@@ -23,6 +23,7 @@ import ..EquilMoistModel
 import ..NonEquilMoistModel
 
 # precip_model
+import ..NoPrecipitation
 import ..Microphysics0Moment
 import ..Microphysics1Moment
 
diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl
index 5ddd5cb6cc..520b6ac942 100644
--- a/src/diagnostics/core_diagnostics.jl
+++ b/src/diagnostics/core_diagnostics.jl
@@ -599,21 +599,30 @@ add_diagnostic_variable!(
 
 ###
 # Precipitation (2d)
-# TODO: add precipitation flux for the 1-moment microphysics
 ###
 compute_pr!(out, state, cache, time) =
     compute_pr!(out, state, cache, time, cache.atmos.precip_model)
 compute_pr!(_, _, _, _, precip_model::T) where {T} =
     error_diagnostic_variable("pr", precip_model)
 
-function compute_pr!(out, state, cache, time, precip_model::Microphysics0Moment)
+function compute_pr!(
+    out,
+    state,
+    cache,
+    time,
+    precip_model::Union{
+        NoPrecipitation,
+        Microphysics0Moment,
+        Microphysics1Moment,
+    },
+)
     if isnothing(out)
-        return cache.precipitation.col_integrated_rain .+
-               cache.precipitation.col_integrated_snow
+        return cache.precipitation.surface_rain_flux .+
+               cache.precipitation.surface_snow_flux
     else
         out .=
-            cache.precipitation.col_integrated_rain .+
-            cache.precipitation.col_integrated_snow
+            cache.precipitation.surface_rain_flux .+
+            cache.precipitation.surface_snow_flux
     end
 end
 
diff --git a/src/diagnostics/default_diagnostics.jl b/src/diagnostics/default_diagnostics.jl
index 592a5faa32..b8b6064761 100644
--- a/src/diagnostics/default_diagnostics.jl
+++ b/src/diagnostics/default_diagnostics.jl
@@ -164,8 +164,8 @@ function default_diagnostics(
     t_end;
     output_writer,
 ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}}
-    moist_diagnostics = ["hur", "hus", "cl", "clw", "cli", "hussfc", "evspsbl"]
-
+    moist_diagnostics =
+        ["hur", "hus", "cl", "clw", "cli", "hussfc", "evspsbl", "pr"]
     average_func = frequency_averages(t_start, t_end)
     return [average_func(moist_diagnostics...; output_writer, t_start)...]
 end
diff --git a/src/parameterized_tendencies/microphysics/precipitation.jl b/src/parameterized_tendencies/microphysics/precipitation.jl
index a2b9db6273..d2557f586b 100644
--- a/src/parameterized_tendencies/microphysics/precipitation.jl
+++ b/src/parameterized_tendencies/microphysics/precipitation.jl
@@ -17,7 +17,13 @@ precipitation_cache(Y, atmos::AtmosModel) =
 ##### No Precipitation
 #####
 
-precipitation_cache(Y, precip_model::NoPrecipitation) = (;)
+function precipitation_cache(Y, precip_model::NoPrecipitation)
+    FT = Spaces.undertype(axes(Y.c))
+    return (;
+        surface_rain_flux = zeros(axes(Fields.level(Y.f, half))),
+        surface_snow_flux = zeros(axes(Fields.level(Y.f, half))),
+    )
+end
 precipitation_tendency!(Yₜ, Y, p, t, ::NoPrecipitation, _) = nothing
 
 #####
@@ -31,12 +37,8 @@ function precipitation_cache(Y, precip_model::Microphysics0Moment)
         ᶜS_ρe_tot = similar(Y.c, FT),
         ᶜ3d_rain = similar(Y.c, FT),
         ᶜ3d_snow = similar(Y.c, FT),
-        col_integrated_rain = zeros(
-            axes(Fields.level(Geometry.WVector.(Y.f.u₃), half)),
-        ),
-        col_integrated_snow = zeros(
-            axes(Fields.level(Geometry.WVector.(Y.f.u₃), half)),
-        ),
+        surface_rain_flux = zeros(axes(Fields.level(Y.f, half))),
+        surface_snow_flux = zeros(axes(Fields.level(Y.f, half))),
     )
 end
 
@@ -126,6 +128,32 @@ function compute_precipitation_cache!(
     end
 end
 
+function compute_precipitation_surface_fluxes!(
+    Y,
+    p,
+    precip_model::Microphysics0Moment,
+)
+    (; ᶜT) = p.core
+    (; ᶜts) = p.precomputed  # assume ᶜts has been updated
+    (; ᶜ3d_rain, ᶜ3d_snow, ᶜS_ρq_tot, ᶜS_ρe_tot) = p.precipitation
+    (; surface_rain_flux, surface_snow_flux) = p.precipitation
+    (; col_integrated_precip_energy_tendency,) = p.conservation_check
+
+    # update total column energy source for surface energy balance
+    Operators.column_integral_definite!(
+        col_integrated_precip_energy_tendency,
+        ᶜS_ρe_tot,
+    )
+    # update surface precipitation fluxes in cache for coupler's use
+    thermo_params = CAP.thermodynamics_params(p.params)
+    T_freeze = TD.Parameters.T_freeze(thermo_params)
+    @. ᶜT = TD.air_temperature(thermo_params, ᶜts)
+    @. ᶜ3d_rain = ifelse(ᶜT >= T_freeze, ᶜS_ρq_tot, 0)
+    @. ᶜ3d_snow = ifelse(ᶜT < T_freeze, ᶜS_ρq_tot, 0)
+    Operators.column_integral_definite!(surface_rain_flux, ᶜ3d_rain)
+    Operators.column_integral_definite!(surface_snow_flux, ᶜ3d_snow)
+end
+
 function precipitation_tendency!(
     Yₜ,
     Y,
@@ -134,36 +162,19 @@ function precipitation_tendency!(
     precip_model::Microphysics0Moment,
     _,
 )
-    (; ᶜT, ᶜΦ) = p.core
-    (; ᶜts) = p.precomputed  # assume ᶜts has been updated
-    (; params) = p
     (; turbconv_model) = p.atmos
-    (; ᶜ3d_rain, ᶜ3d_snow, ᶜS_ρq_tot, ᶜS_ρe_tot) = p.precipitation
-    (; col_integrated_rain, col_integrated_snow) = p.precipitation
-    (; col_integrated_precip_energy_tendency,) = p.conservation_check
-    thermo_params = CAP.thermodynamics_params(params)
+    (; ᶜS_ρq_tot, ᶜS_ρe_tot) = p.precipitation
 
     # Compute the ρq_tot and ρe_tot precipitation source terms
     compute_precipitation_cache!(Y, p, precip_model, turbconv_model)
+    # Compute surface precipitation flux
+    compute_precipitation_surface_fluxes!(Y, p, precip_model)
 
     # Add the source terms to the tendencies
     @. Yₜ.c.ρq_tot += ᶜS_ρq_tot
     @. Yₜ.c.ρ += ᶜS_ρq_tot
     @. Yₜ.c.ρe_tot += ᶜS_ρe_tot
 
-    # update total column energy source for surface energy balance
-    Operators.column_integral_definite!(
-        col_integrated_precip_energy_tendency,
-        ᶜS_ρe_tot,
-    )
-    # update precip in cache for coupler's use
-    # 3d rain and snow
-    T_freeze = TD.Parameters.T_freeze(thermo_params)
-    @. ᶜT = TD.air_temperature(thermo_params, ᶜts)
-    @. ᶜ3d_rain = ifelse(ᶜT >= T_freeze, ᶜS_ρq_tot, 0)
-    @. ᶜ3d_snow = ifelse(ᶜT < T_freeze, ᶜS_ρq_tot, 0)
-    Operators.column_integral_definite!(col_integrated_rain, ᶜ3d_rain)
-    Operators.column_integral_definite!(col_integrated_snow, ᶜ3d_snow)
     return nothing
 end
 
@@ -178,6 +189,8 @@ function precipitation_cache(Y, precip_model::Microphysics1Moment)
         ᶜSqᵣᵖ = similar(Y.c, FT),
         ᶜSqₛᵖ = similar(Y.c, FT),
         ᶜSeₜᵖ = similar(Y.c, FT),
+        surface_rain_flux = zeros(axes(Fields.level(Y.f, half))),
+        surface_snow_flux = zeros(axes(Fields.level(Y.f, half))),
     )
 end
 
@@ -235,7 +248,6 @@ function compute_precipitation_cache!(Y, p, ::Microphysics1Moment, _)
     # first term of eq 36 from Raymond 2013
     compute_precipitation_heating!(ᶜSeₜᵖ, ᶜwᵣ, ᶜwₛ, ᶜu, ᶜqᵣ, ᶜqₛ, ᶜts, ᶜ∇T, thp)
 end
-
 function compute_precipitation_cache!(
     Y,
     p,
@@ -283,6 +295,34 @@ function compute_precipitation_cache!(
     compute_precipitation_heating!(ᶜSeₜᵖ, ᶜwᵣ, ᶜwₛ, ᶜu, ᶜqᵣ, ᶜqₛ, ᶜts, ᶜ∇T, thp)
 end
 
+function compute_precipitation_surface_fluxes!(
+    Y,
+    p,
+    precip_model::Microphysics1Moment,
+)
+    (; surface_rain_flux, surface_snow_flux) = p.precipitation
+    (; col_integrated_precip_energy_tendency,) = p.conservation_check
+    (; ᶜqᵣ, ᶜqₛ, ᶜwᵣ, ᶜwₛ) = p.precomputed
+
+    (; ᶠtemp_scalar) = p.scratch
+    slg = Fields.level(Fields.local_geometry_field(ᶠtemp_scalar), Fields.half)
+
+    # Constant extrapolation: - put values from bottom cell center to bottom cell face
+    ˢρ = Fields.Field(Fields.field_values(Fields.level(Y.c.ρ, 1)), axes(slg))
+    # For density this is equivalent with ᶠwinterp(ᶜJ, Y.c.ρ) and therefore
+    # consistent with the way we do vertical advection
+    ˢqᵣ = Fields.Field(Fields.field_values(Fields.level(ᶜqᵣ, 1)), axes(slg))
+    ˢqₛ = Fields.Field(Fields.field_values(Fields.level(ᶜqₛ, 1)), axes(slg))
+    ˢwᵣ = Fields.Field(Fields.field_values(Fields.level(ᶜwᵣ, 1)), axes(slg))
+    ˢwₛ = Fields.Field(Fields.field_values(Fields.level(ᶜwₛ, 1)), axes(slg))
+
+    # Project the flux to CT3 vector and convert to physical units.
+    @. surface_rain_flux =
+        -projected_vector_data(CT3, ˢρ * ˢqᵣ * Geometry.WVector(ˢwᵣ), slg)
+    @. surface_snow_flux =
+        -projected_vector_data(CT3, ˢρ * ˢqₛ * Geometry.WVector(ˢwₛ), slg)
+end
+
 function precipitation_tendency!(
     Yₜ,
     Y,
@@ -292,10 +332,12 @@ function precipitation_tendency!(
     _,
 )
     (; turbconv_model) = p.atmos
-    compute_precipitation_cache!(Y, p, precip_model, turbconv_model)
-
     (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ) = p.precipitation
 
+    # Populate the cache and precipitation surface fluxes
+    compute_precipitation_cache!(Y, p, precip_model, turbconv_model)
+    compute_precipitation_surface_fluxes!(Y, p, precip_model)
+
     # Update grid mean tendencies
     @. Yₜ.c.ρ += Y.c.ρ * ᶜSqₜᵖ
     @. Yₜ.c.ρq_tot += Y.c.ρ * ᶜSqₜᵖ
@@ -322,7 +364,9 @@ function precipitation_tendency!(
 
     (; ᶜρaʲs) = p.precomputed
 
+    # Populate the cache and precipitation surface fluxes
     compute_precipitation_cache!(Y, p, precip_model, turbconv_model)
+    compute_precipitation_surface_fluxes!(Y, p, precip_model)
 
     # Update from environment precipitation sources
     # and the grid mean precipitation sinks
@@ -357,7 +401,9 @@ function precipitation_tendency!(
     # Grid mean precipitation sinks
     (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ) = p.precipitation
 
+    # Populate the cache and precipitation surface fluxes
     compute_precipitation_cache!(Y, p, precip_model, turbconv_model)
+    compute_precipitation_surface_fluxes!(Y, p, precip_model)
 
     # Update from environment precipitation sources
     # and the grid mean precipitation sinks
diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl
index 8675b16ee0..44dcbf6ba5 100644
--- a/src/prognostic_equations/advection.jl
+++ b/src/prognostic_equations/advection.jl
@@ -47,6 +47,7 @@ end
 NVTX.@annotate function horizontal_tracer_advection_tendency!(Yₜ, Y, p, t)
     n = n_mass_flux_subdomains(p.atmos.turbconv_model)
     (; ᶜu) = p.precomputed
+
     if p.atmos.turbconv_model isa PrognosticEDMFX
         (; ᶜuʲs) = p.precomputed
     end
@@ -62,7 +63,6 @@ NVTX.@annotate function horizontal_tracer_advection_tendency!(Yₜ, Y, p, t)
                 Y.c.sgsʲs.:($$j).q_tot * wdivₕ(ᶜuʲs.:($$j))
         end
     end
-
     return nothing
 end
 
diff --git a/src/prognostic_equations/implicit/implicit_tendency.jl b/src/prognostic_equations/implicit/implicit_tendency.jl
index 3dd83984f2..2fa603db2a 100644
--- a/src/prognostic_equations/implicit/implicit_tendency.jl
+++ b/src/prognostic_equations/implicit/implicit_tendency.jl
@@ -151,17 +151,12 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t)
         # is done with other passive tracers in the explicit tendency.
         # Here we add the advection with precipitation terminal velocity
         # using downward biasing and free outflow bottom boundary condition
-
-        ᶠlg = Fields.local_geometry_field(Y.f)
+        (; ᶜqᵣ, ᶜqₛ, ᶜwᵣ, ᶜwₛ) = p.precomputed
         @. Yₜ.c.ρq_rai -= ᶜprecipdivᵥ(
-            CT3(unit_basis_vector_data(CT3, ᶠlg)) *
-            ᶠwinterp(ᶜJ, Y.c.ρ) *
-            ᶠright_bias(-p.precomputed.ᶜwᵣ * ᶜspecific.q_rai),
+            ᶠwinterp(ᶜJ, Y.c.ρ) * ᶠright_bias(Geometry.WVector(-(ᶜwᵣ)) * ᶜqᵣ),
         )
         @. Yₜ.c.ρq_sno -= ᶜprecipdivᵥ(
-            CT3(unit_basis_vector_data(CT3, ᶠlg)) *
-            ᶠwinterp(ᶜJ, Y.c.ρ) *
-            ᶠright_bias(-p.precomputed.ᶜwₛ * ᶜspecific.q_sno),
+            ᶠwinterp(ᶜJ, Y.c.ρ) * ᶠright_bias(Geometry.WVector(-(ᶜwₛ)) * ᶜqₛ),
         )
     end
 
diff --git a/src/prognostic_equations/surface_temp.jl b/src/prognostic_equations/surface_temp.jl
index 39ae52d726..2f8ad703cd 100644
--- a/src/prognostic_equations/surface_temp.jl
+++ b/src/prognostic_equations/surface_temp.jl
@@ -51,8 +51,8 @@ function surface_temp_tendency!(Yₜ, Y, p, t, slab::PrognosticSurfaceTemperatur
         sfc_turb_w_flux = p.precomputed.sfc_conditions.ρ_flux_q_tot
 
         # precipitation
-        P_liq = p.precipitation.col_integrated_rain
-        P_snow = p.precipitation.col_integrated_snow
+        P_liq = p.precipitation.surface_rain_flux
+        P_snow = p.precipitation.surface_snow_flux
 
         @. Yₜ.sfc.water -=
             P_liq +
diff --git a/test/parameterized_tendencies/microphysics/precipitation.jl b/test/parameterized_tendencies/microphysics/precipitation.jl
index 4f061aa52b..9adb233d08 100644
--- a/test/parameterized_tendencies/microphysics/precipitation.jl
+++ b/test/parameterized_tendencies/microphysics/precipitation.jl
@@ -38,8 +38,8 @@ include("../../test_helpers.jl")
         :ᶜS_ρq_tot,
         :ᶜ3d_rain,
         :ᶜ3d_snow,
-        :col_integrated_rain,
-        :col_integrated_snow,
+        :surface_rain_flux,
+        :surface_snow_flux,
     )
     for var_name in test_varnames
         @test var_name ∈ propertynames(precip_cache)
@@ -72,7 +72,8 @@ include("../../test_helpers.jl")
     (; turbconv_model) = p.atmos
     precip_cache = CA.precipitation_cache(Y, precip_model)
     ᶜYₜ = Y .* FT(0)
-    test_varnames = (:ᶜSqₜᵖ, :ᶜSqᵣᵖ, :ᶜSqₛᵖ, :ᶜSeₜᵖ)
+    test_varnames =
+        (:ᶜSqₜᵖ, :ᶜSqᵣᵖ, :ᶜSqₛᵖ, :ᶜSeₜᵖ, :surface_rain_flux, :surface_snow_flux)
     for var_name in test_varnames
         @test var_name ∈ propertynames(precip_cache)
     end