diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index 7974b5f12b..54bf3d9e04 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -96,6 +96,14 @@ steps:
         command: "julia --color=yes --project=examples test/orographic_gravity_wave/ogwd_3d.jl"
         artifact_paths: "orographic_gravity_wave_test_3d/*"
 
+  - group: "Precipitation"
+    steps:
+      - label: ":umbrella: 1-moment precipitation sanity test single column"
+        command: >
+          julia --color=yes --project=examples examples/hybrid/driver.jl
+          --config_file $CONFIG_PATH/single_column_precipitation_test.yml
+        artifact_paths: "single_column_precipitation_test/*"
+
   - group: "Column Examples"
     steps:
 
@@ -240,7 +248,7 @@ steps:
           --job_id sphere_baroclinic_wave_rhoe_equilmoist
           --out_dir sphere_baroclinic_wave_rhoe_equilmoist
         artifact_paths: "sphere_baroclinic_wave_rhoe_equilmoist/*"
-      
+
       - label: ":computer: no lim ARS baroclinic wave (ρe) equilmoist explicit vertdiff"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
@@ -314,7 +322,7 @@ steps:
   - group: "Sphere Examples (Aquaplanet)"
     steps:
 
-      - label: ":computer: aquaplanet (ρe_tot) equilmoist allsky radiation monin_obukhov varying insolation gravity wave (gfdl_restart) high top"
+      - label: ":computer: aquaplanet (ρe_tot) equilmoist allsky radiation monin_obukhov varying insolation gravity wave (gfdl_restart) high top with 1-moment micro"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.yml
@@ -576,7 +584,7 @@ steps:
         artifact_paths: "diagnostic_edmfx_trmm_box/*"
         agents:
           slurm_mem: 20GB
-      
+
       - label: ":genie: Diagnostic EDMFX TRMM stretched grid in a box"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
@@ -619,7 +627,7 @@ steps:
         artifact_paths: "prognostic_edmfx_bomex_fixtke_box/*"
         agents:
           slurm_mem: 20GB
-      
+
       - label: ":genie: Prognostic EDMFX Bomex in a box"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
@@ -627,7 +635,7 @@ steps:
         artifact_paths: "prognostic_edmfx_bomex_box/*"
         agents:
           slurm_mem: 20GB
-      
+
       - label: ":genie: Prognostic EDMFX Bomex stretched grid in a box"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
@@ -643,7 +651,7 @@ steps:
         artifact_paths: "prognostic_edmfx_dycoms_rf01_box/*"
         agents:
           slurm_mem: 20GB
-      
+
       - label: ":genie: Prognostic EDMFX Rico in a column"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
@@ -659,7 +667,7 @@ steps:
         artifact_paths: "prognostic_edmfx_trmm_column/*"
         agents:
           slurm_mem: 20GB
-      
+
       - label: ":genie: Prognostic EDMFX aquaplanet"
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
@@ -750,7 +758,7 @@ steps:
       - label: "GPU: gpu_aquaplanet_dyamond"
         command:
           - mkdir -p gpu_aquaplanet_dyamond
-          - > 
+          - >
             julia --color=yes --project=examples examples/hybrid/driver.jl
             --config_file ${GPU_CONFIG_PATH}gpu_aquaplanet_dyamond.yml
         artifact_paths: "gpu_aquaplanet_dyamond/*"
diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml
index eec06174f0..c678117a67 100644
--- a/config/default_configs/default_config.yml
+++ b/config/default_configs/default_config.yml
@@ -198,6 +198,9 @@ density_upwinding:
 tracer_upwinding:
   help: "Tracer upwinding mode [`none` (default), `first_order` , `third_order`, `boris_book`, `zalesak`]"
   value: none
+precip_upwinding:
+  help: "Precipitation upwinding mode [`none` (default), `first_order` , `third_order`]"
+  value: none
 energy_upwinding:
   help: "Energy upwinding mode [`none` (default), `first_order` , `third_order`, `boris_book`, `zalesak`]"
   value: none
@@ -219,6 +222,9 @@ regression_test:
 check_conservation:
   help: "Check conservation of mass and energy [`false` (default), `true`]"
   value: false
+check_precipitation:
+  help: "Sanity checks for 1-moment precipitation [`false` (default), `true`]"
+  value: false
 ls_adv:
   help: "Large-scale advection [`nothing` (default), `Bomex`, `LifeCycleTan2018`, `Rico`, `ARM_SGP`, `GATE_III`]"
   value: ~
diff --git a/config/model_configs/single_column_precipitation_test.yml b/config/model_configs/single_column_precipitation_test.yml
new file mode 100644
index 0000000000..8caac6f508
--- /dev/null
+++ b/config/model_configs/single_column_precipitation_test.yml
@@ -0,0 +1,32 @@
+config: "column"
+initial_condition: "PrecipitatingColumn"
+surface_setup: "DefaultExchangeCoefficients"
+z_elem: 200
+z_max: 10000.0
+z_stretch: false
+dt: "10secs"
+t_end: "1500secs"
+dt_save_to_disk: "500secs"
+moist: "nonequil"
+precip_model: "1M"
+precip_upwinding: "first_order"
+hyperdiff: "false"
+regression_test: false
+check_precipitation: true
+job_id: "single_column_precipitation_test"
+toml: [toml/single_column_precipitation_test.toml]
+diagnostics:
+  - short_name: hus
+    period: 500secs
+  - short_name: clw
+    period: 500secs
+  - short_name: cli
+    period: 500secs
+  - short_name: husra
+    period: 500secs
+  - short_name: hussn
+    period: 500secs
+  - short_name: ta
+    period: 500secs
+  - short_name: wa
+    period: 500secs
diff --git a/config/model_configs/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.yml b/config/model_configs/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.yml
index 66d727c5ea..d9a0f27897 100644
--- a/config/model_configs/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.yml
+++ b/config/model_configs/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.yml
@@ -11,7 +11,8 @@ dz_bottom: 300.0
 vert_diff: "true"
 idealized_insolation: false
 z_max: 45000.0
-precip_model: "0M"
+precip_model: "1M"
+precip_upwinding: "first_order"
 job_id: "sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res"
 moist: "equil"
 toml: [toml/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.toml]
diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl
index 85f42466a1..1c1f5e52ef 100644
--- a/examples/hybrid/driver.jl
+++ b/examples/hybrid/driver.jl
@@ -36,6 +36,7 @@ include(joinpath(pkgdir(CA), "post_processing", "post_processing_funcs.jl"))
 include(
     joinpath(pkgdir(CA), "post_processing", "define_tc_quicklook_profiles.jl"),
 )
+include(joinpath(pkgdir(CA), "post_processing", "plot_single_column_precip.jl"))
 
 ref_job_id = config.parsed_args["reference_job_id"]
 reference_job_id = isnothing(ref_job_id) ? simulation.job_id : ref_job_id
@@ -123,3 +124,55 @@ if config.parsed_args["check_conservation"]
           (p.net_energy_flux_sfc[][] - p.net_energy_flux_toa[][]) ≈
           sum(sol.u[end].c.ρe_tot) rtol = 100 * eps(FT)
 end
+
+if config.parsed_args["check_precipitation"]
+
+    # plot results of the single column precipitation test
+    plot_single_column_precip(simulation.output_dir, reference_job_id)
+
+    # run some simple tests based on the output
+    FT = Spaces.undertype(axes(sol.u[end].c.ρ))
+    Yₜ = similar(sol.u[end])
+
+    Yₜ_ρ = similar(Yₜ.c.ρq_rai)
+    Yₜ_ρqₚ = similar(Yₜ.c.ρq_rai)
+    Yₜ_ρqₜ = similar(Yₜ.c.ρq_rai)
+
+    CA.remaining_tendency!(Yₜ, sol.u[end], sol.prob.p, sol.t[end])
+
+    @. Yₜ_ρqₚ = -Yₜ.c.ρq_rai - Yₜ.c.ρq_sno
+    @. Yₜ_ρqₜ = Yₜ.c.ρq_tot
+    @. Yₜ_ρ = Yₜ.c.ρ
+
+    Fields.bycolumn(axes(sol.u[end].c.ρ)) do colidx
+
+        # no nans
+        @assert !any(isnan, Yₜ.c.ρ[colidx])
+        @assert !any(isnan, Yₜ.c.ρq_tot[colidx])
+        @assert !any(isnan, Yₜ.c.ρe_tot[colidx])
+        @assert !any(isnan, Yₜ.c.ρq_rai[colidx])
+        @assert !any(isnan, Yₜ.c.ρq_sno[colidx])
+        @assert !any(isnan, sol.prob.p.precipitation.ᶜwᵣ[colidx])
+        @assert !any(isnan, sol.prob.p.precipitation.ᶜwₛ[colidx])
+
+        # treminal velocity is positive
+        @test minimum(sol.prob.p.precipitation.ᶜwᵣ[colidx]) >= FT(0)
+        @test minimum(sol.prob.p.precipitation.ᶜwₛ[colidx]) >= FT(0)
+
+        # checking for water budget conservation
+        # in the presence of precipitation sinks
+        # (This test only works without surface flux of q_tot)
+        @test all(
+            ClimaCore.isapprox(
+                Yₜ_ρqₜ[colidx],
+                Yₜ_ρqₚ[colidx],
+                rtol = 1e2 * eps(FT),
+            ),
+        )
+
+        # mass budget consistency
+        @test all(
+            ClimaCore.isapprox(Yₜ_ρ[colidx], Yₜ_ρqₜ[colidx], rtol = eps(FT)),
+        )
+    end
+end
diff --git a/post_processing/plot_single_column_precip.jl b/post_processing/plot_single_column_precip.jl
new file mode 100644
index 0000000000..41b8d75d94
--- /dev/null
+++ b/post_processing/plot_single_column_precip.jl
@@ -0,0 +1,110 @@
+import ClimaAtmos as CA
+import ClimaCore as CO
+import CairoMakie as MK
+
+function plot_single_column_precip(output_dir, job_id)
+
+    fig = MK.Figure(resolution = (1200, 600))
+    ax1 = MK.Axis(fig[1, 1], ylabel = "z [km]", xlabel = "q_tot [g/kg]")
+    ax4 = MK.Axis(fig[2, 1], ylabel = "z [km]", xlabel = "T [K]")
+    ax2 = MK.Axis(fig[1, 2], xlabel = "q_liq [g/kg]")
+    ax3 = MK.Axis(fig[1, 3], xlabel = "q_ice [g/kg]")
+    ax5 = MK.Axis(fig[2, 2], xlabel = "q_rai [g/kg]")
+    ax6 = MK.Axis(fig[2, 3], xlabel = "q_sno [g/kg]")
+
+    path = joinpath(pkgdir(CA), output_dir)
+
+    col = Dict(
+        "0" => :navy,
+        "500" => :blue2,
+        "1000" => :royalblue,
+        "1500" => :skyblue1,
+    )
+
+    for time in ["0", "500", "1000", "1500"]
+
+        fqₜ = CO.InputOutput.HDF5Reader(
+            joinpath(path, "hus_inst_" * time * ".0.h5"),
+        )
+        fqₗ = CO.InputOutput.HDF5Reader(
+            joinpath(path, "clw_inst_" * time * ".0.h5"),
+        )
+        fqᵢ = CO.InputOutput.HDF5Reader(
+            joinpath(path, "cli_inst_" * time * ".0.h5"),
+        )
+        fqᵣ = CO.InputOutput.HDF5Reader(
+            joinpath(path, "husra_inst_" * time * ".0.h5"),
+        )
+        fqₛ = CO.InputOutput.HDF5Reader(
+            joinpath(path, "hussn_inst_" * time * ".0.h5"),
+        )
+        fTₐ = CO.InputOutput.HDF5Reader(
+            joinpath(path, "ta_inst_" * time * ".0.h5"),
+        )
+        fwₐ = CO.InputOutput.HDF5Reader(
+            joinpath(path, "wa_inst_" * time * ".0.h5"),
+        )
+
+        qₜ = CO.InputOutput.read_field(fqₜ, "hus_inst")
+        qₗ = CO.InputOutput.read_field(fqₗ, "clw_inst")
+        qᵢ = CO.InputOutput.read_field(fqᵢ, "cli_inst")
+        qᵣ = CO.InputOutput.read_field(fqᵣ, "husra_inst")
+        qₛ = CO.InputOutput.read_field(fqₛ, "hussn_inst")
+        Tₐ = CO.InputOutput.read_field(fTₐ, "ta_inst")
+        wₐ = CO.InputOutput.read_field(fwₐ, "wa_inst")
+
+        qₜ_col = CO.Fields.column(qₜ, 1, 1, 1)
+        qₗ_col = CO.Fields.column(qₗ, 1, 1, 1)
+        qᵢ_col = CO.Fields.column(qᵢ, 1, 1, 1)
+        qᵣ_col = CO.Fields.column(qᵣ, 1, 1, 1)
+        qₛ_col = CO.Fields.column(qₛ, 1, 1, 1)
+        Tₐ_col = CO.Fields.column(Tₐ, 1, 1, 1)
+        wₐ_col = CO.Fields.column(wₐ, 1, 1, 1)
+        z = CO.Fields.coordinate_field(qₜ_col).z
+
+        MK.lines!(
+            ax1,
+            vec(parent(qₜ_col)) .* 1e3,
+            vec(parent(z)) ./ 1e3,
+            color = col[time],
+        )
+        MK.lines!(
+            ax2,
+            vec(parent(qₗ_col)) .* 1e3,
+            vec(parent(z)) ./ 1e3,
+            color = col[time],
+        )
+        MK.lines!(
+            ax3,
+            vec(parent(qᵢ_col)) .* 1e3,
+            vec(parent(z)) ./ 1e3,
+            color = col[time],
+        )
+        MK.lines!(
+            ax4,
+            vec(parent(Tₐ_col)),
+            vec(parent(z)) ./ 1e3,
+            color = col[time],
+        )
+        MK.lines!(
+            ax5,
+            vec(parent(qᵣ_col)) .* 1e3,
+            vec(parent(z)) ./ 1e3,
+            color = col[time],
+        )
+        MK.lines!(
+            ax6,
+            vec(parent(qₛ_col)) .* 1e3,
+            vec(parent(z)) ./ 1e3,
+            color = col[time],
+        )
+
+        for fid in [fqₜ, fqₗ, fqᵢ, fqᵣ, fqₛ, fTₐ, fwₐ]
+            close(fid)
+        end
+    end
+
+    @info("Saving precipitation plots to " * path)
+
+    MK.save(joinpath(path, "single_column_precipitation.png"), fig)
+end
diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
index 59253291c1..ddb73bde86 100644
--- a/src/cache/precomputed_quantities.jl
+++ b/src/cache/precomputed_quantities.jl
@@ -227,7 +227,7 @@ function thermo_state(
             eltype(thermo_params)(0.003),
         )
     get_ts(ρ::Real, ::Nothing, ::Nothing, e_int::Real, ::Nothing, q_pt) =
-        TD.PhaseNonEquil_ρeq(thermo_params, ρ, e_int, q_pt)
+        TD.PhaseNonEquil(thermo_params, e_int, ρ, q_pt)
     get_ts(::Nothing, p::Real, θ::Real, ::Nothing, ::Nothing, ::Nothing) =
         TD.PhaseDry_pθ(thermo_params, p, θ)
     get_ts(::Nothing, p::Real, θ::Real, ::Nothing, q_tot::Real, ::Nothing) =
diff --git a/src/diagnostics/Diagnostics.jl b/src/diagnostics/Diagnostics.jl
index b9ea4868ac..65b7e78dd7 100644
--- a/src/diagnostics/Diagnostics.jl
+++ b/src/diagnostics/Diagnostics.jl
@@ -21,6 +21,7 @@ import ..NonEquilMoistModel
 
 # precip_model
 import ..Microphysics0Moment
+import ..Microphysics1Moment
 
 # radiation
 import ClimaAtmos.RRTMGPInterface as RRTMGPI
diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl
index 24a21b79cf..969a2d9ed8 100644
--- a/src/diagnostics/core_diagnostics.jl
+++ b/src/diagnostics/core_diagnostics.jl
@@ -523,7 +523,9 @@ add_diagnostic_variable!(
 )
 
 ###
-# Precipitation (2d) - TODO: change to kg m^-2 s^-1
+# Precipitation (2d)
+# TODO: change to kg m^-2 s^-1
+# 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)
@@ -549,3 +551,68 @@ add_diagnostic_variable!(
     comments = "Total precipitation including rain and snow",
     compute! = compute_pr!,
 )
+
+###
+# Precipitation (3d)
+###
+compute_husra!(out, state, cache, time) =
+    compute_husra!(out, state, cache, time, cache.atmos.precip_model)
+compute_husra!(_, _, _, _, model::T) where {T} =
+    error_diagnostic_variable("husra", model)
+
+function compute_husra!(
+    out,
+    state,
+    cache,
+    time,
+    precip_model::Microphysics1Moment,
+)
+    if isnothing(out)
+        return state.c.ρq_rai ./ state.c.ρ
+    else
+        out .= state.c.ρq_rai ./ state.c.ρ
+    end
+end
+
+add_diagnostic_variable!(
+    short_name = "husra",
+    long_name = "Mass Fraction of Rain",
+    standard_name = "mass_fraction_of_rain_in_air",
+    units = "kg kg^-1",
+    comments = """
+    This is calculated as the mass of rain water in the grid cell divided by
+    the mass of air (dry air + water vapor + cloud condensate) in the grid cells.
+    """,
+    compute! = compute_husra!,
+)
+
+compute_hussn!(out, state, cache, time) =
+    compute_hussn!(out, state, cache, time, cache.atmos.precip_model)
+compute_hussn!(_, _, _, _, model::T) where {T} =
+    error_diagnostic_variable("hussn", model)
+
+function compute_hussn!(
+    out,
+    state,
+    cache,
+    time,
+    precip_model::Microphysics1Moment,
+)
+    if isnothing(out)
+        return state.c.ρq_sno ./ state.c.ρ
+    else
+        out .= state.c.ρq_sno ./ state.c.ρ
+    end
+end
+
+add_diagnostic_variable!(
+    short_name = "hussn",
+    long_name = "Mass Fraction of Snow",
+    standard_name = "mass_fraction_of_snow_in_air",
+    units = "kg kg^-1",
+    comments = """
+    This is calculated as the mass of snow in the grid cell divided by
+    the mass of air (dry air + water vapor + cloud condensate) in the grid cells.
+    """,
+    compute! = compute_hussn!,
+)
diff --git a/src/initial_conditions/atmos_state.jl b/src/initial_conditions/atmos_state.jl
index 98ebbcf5b9..c4250362bd 100644
--- a/src/initial_conditions/atmos_state.jl
+++ b/src/initial_conditions/atmos_state.jl
@@ -102,8 +102,10 @@ moisture_variables(ls, ::NonEquilMoistModel) = (;
 # TODO: Remove perf_mode. Currently, adding tracers hurts performance.
 precip_variables(ls, ::NoPrecipitation, ::PerfStandard) = (;)
 precip_variables(ls, ::Microphysics0Moment, ::PerfStandard) = (;)
-precip_variables(ls, ::Microphysics1Moment, ::PerfStandard) =
-    (; ρq_rai = zero(eltype(ls)), ρq_sno = zero(eltype(ls)))
+precip_variables(ls, ::Microphysics1Moment, ::PerfStandard) = (;
+    ρq_rai = ls.ρ * ls.precip_state.q_rai,
+    ρq_sno = ls.ρ * ls.precip_state.q_sno,
+)
 precip_variables(ls, _, ::PerfExperimental) =
     (; ρq_rai = zero(eltype(ls)), ρq_sno = zero(eltype(ls)))
 
diff --git a/src/initial_conditions/initial_conditions.jl b/src/initial_conditions/initial_conditions.jl
index 45a6b3fbd8..5d650800f1 100644
--- a/src/initial_conditions/initial_conditions.jl
+++ b/src/initial_conditions/initial_conditions.jl
@@ -946,3 +946,46 @@ function (initial_condition::TRMM_LBA)(params)
     end
     return local_state
 end
+
+"""
+    PrecipitatingColumn
+
+A 1-dimensional precipitating column test
+"""
+struct PrecipitatingColumn <: InitialCondition end
+
+prescribed_prof(::Type{FT}, z_mid, z_max, val) where {FT} =
+    z -> z < z_max ? FT(val) * exp(-(z - FT(z_mid))^2 / 2 / FT(1e3)^2) : FT(0)
+
+function (initial_condition::PrecipitatingColumn)(params)
+    FT = eltype(params)
+    thermo_params = CAP.thermodynamics_params(params)
+    p_0 = FT(101300.0)
+    qᵣ = prescribed_prof(FT, 2000, 5000, 1e-6)
+    qₛ = prescribed_prof(FT, 5000, 8000, 2e-6)
+    qₗ = prescribed_prof(FT, 4000, 5500, 2e-5)
+    qᵢ = prescribed_prof(FT, 6000, 9000, 1e-5)
+    θ = APL.Rico_θ_liq_ice(FT)
+    q_tot = APL.Rico_q_tot(FT)
+    u = prescribed_prof(FT, 0, Inf, 0)
+    v = prescribed_prof(FT, 0, Inf, 0)
+    p = hydrostatic_pressure_profile(; thermo_params, p_0, θ, q_tot)
+    function local_state(local_geometry)
+        (; z) = local_geometry.coordinates
+        ts = TD.PhaseNonEquil_pθq(
+            thermo_params,
+            p(z),
+            θ(z),
+            TD.PhasePartition(q_tot(z), qₗ(z), qᵢ(z)),
+        )
+        return LocalState(;
+            params,
+            geometry = local_geometry,
+            thermo_state = ts,
+            velocity = Geometry.UVVector(u(z), v(z)),
+            turbconv_state = nothing,
+            precip_state = PrecipState1M(; q_rai = qᵣ(z), q_sno = qₛ(z)),
+        )
+    end
+    return local_state
+end
diff --git a/src/initial_conditions/local_state.jl b/src/initial_conditions/local_state.jl
index 163e3297a2..bc40924086 100644
--- a/src/initial_conditions/local_state.jl
+++ b/src/initial_conditions/local_state.jl
@@ -7,7 +7,15 @@ an `AtmosModel`.
 abstract type TurbconvState{FT} end
 
 """
-    LocalState(; params, geometry, thermo_state, velocity, turbconv_state)
+    PrecipState{FT}
+
+A collection of values that are required to initialize the `precip_model` of
+an `AtmosModel`.
+"""
+abstract type PrecipState{FT} end
+
+"""
+    LocalState(; params, geometry, thermo_state, velocity, turbconv_state, precip_state)
 
 A generic representation of all the data required to initialize an `AtmosModel`
 at some point in the domain. If `velocity` or `turbconv_state` are omitted, they
@@ -20,6 +28,7 @@ struct LocalState{
     TS <: TD.ThermodynamicState{FT},
     V <: Geometry.LocalVector{FT},
     TC <: TurbconvState{FT},
+    PS <: PrecipState{FT},
     TP,
 }
     params::P
@@ -27,6 +36,7 @@ struct LocalState{
     thermo_state::TS
     velocity::V
     turbconv_state::TC
+    precip_state::PS
 
     # commonly used values that can be inferred from the values above
     thermo_params::TP
@@ -39,6 +49,7 @@ function LocalState(;
     thermo_state,
     velocity = nothing,
     turbconv_state = nothing,
+    precip_state = nothing,
 )
     FT = eltype(params)
     return LocalState(
@@ -47,6 +58,7 @@ function LocalState(;
         thermo_state,
         isnothing(velocity) ? Geometry.UVVector(FT(0), FT(0)) : velocity,
         isnothing(turbconv_state) ? NoTurbconvState{FT}() : turbconv_state,
+        isnothing(precip_state) ? NoPrecipState{FT}() : precip_state,
         CAP.thermodynamics_params(params),
         TD.air_density(CAP.thermodynamics_params(params), thermo_state),
     )
@@ -78,3 +90,26 @@ struct EDMFState{FT} <: TurbconvState{FT}
 end
 EDMFState(; tke, draft_area = 0, velocity = Geometry.WVector(0)) =
     EDMFState{typeof(tke)}(tke, draft_area, velocity)
+
+"""
+    NoPrecipState{FT}()
+
+Indicates that no initial conditions are available for the `precip_model`. Any
+values required by the `precip_model` are set to 0.
+"""
+struct NoPrecipState{FT} <: PrecipState{FT} end
+
+@inline Base.getproperty(::NoPrecipState{FT}, ::Symbol) where {FT} = FT(0)
+
+"""
+    PrecipState1M(; q_rai, q_sno)
+
+Stores the values of `ρq_rai` and `ρq_sno` for the `precip_model`.
+If no values are provided, they are set to zero.
+"""
+struct PrecipState1M{FT} <: PrecipState{FT}
+    q_rai::FT
+    q_sno::FT
+end
+PrecipState1M(; q_rai = 0, q_sno = 0) =
+    PrecipState1M{typeof(q_rai)}(q_rai, q_sno)
diff --git a/src/parameterized_tendencies/microphysics/precipitation.jl b/src/parameterized_tendencies/microphysics/precipitation.jl
index a38162f3cd..d64271f665 100644
--- a/src/parameterized_tendencies/microphysics/precipitation.jl
+++ b/src/parameterized_tendencies/microphysics/precipitation.jl
@@ -20,7 +20,7 @@ precipitation_cache(Y, precip_model::NoPrecipitation) = (;)
 precipitation_tendency!(Yₜ, Y, p, t, colidx, ::NoPrecipitation) = nothing
 
 #####
-##### 0-Moment without sgs scheme
+##### 0-Moment without sgs scheme or with diagnostic edmf
 #####
 
 function precipitation_cache(Y, precip_model::Microphysics0Moment)
@@ -144,3 +144,61 @@ function precipitation_tendency!(
     end
     return nothing
 end
+
+#####
+##### 1-Moment without sgs scheme
+#####
+
+function precipitation_cache(Y, precip_model::Microphysics1Moment)
+    FT = Spaces.undertype(axes(Y.c))
+    return (;
+        ᶜSqₜᵖ = similar(Y.c, FT),
+        ᶜSqᵣᵖ = similar(Y.c, FT),
+        ᶜSqₛᵖ = similar(Y.c, FT),
+        ᶜSeₜᵖ = similar(Y.c, FT),
+        ᶜwᵣ = similar(Y.c, FT),
+        ᶜwₛ = similar(Y.c, FT),
+    )
+end
+
+function compute_precipitation_cache!(Y, p, colidx, ::Microphysics1Moment, _)
+    FT = Spaces.undertype(axes(Y.c))
+
+    (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ, ᶜwᵣ, ᶜwₛ) = p.precipitation
+
+    # TODO in the next PR
+    @. ᶜwᵣ[colidx] = FT(0)
+    @. ᶜwₛ[colidx] = FT(0)
+    @. ᶜSqₜᵖ[colidx] = FT(0)
+    @. ᶜSeₜᵖ[colidx] = FT(0)
+    @. ᶜSqᵣᵖ[colidx] = FT(0)
+    @. ᶜSqₛᵖ[colidx] = FT(0)
+end
+
+function precipitation_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    colidx,
+    precip_model::Microphysics1Moment,
+)
+    FT = Spaces.undertype(axes(Y.c))
+    (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ) = p.precipitation
+
+    compute_precipitation_cache!(
+        Y,
+        p,
+        colidx,
+        precip_model,
+        p.atmos.turbconv_model,
+    )
+
+    @. Yₜ.c.ρ[colidx] += Y.c.ρ[colidx] * ᶜSqₜᵖ[colidx]
+    @. Yₜ.c.ρq_tot[colidx] += Y.c.ρ[colidx] * ᶜSqₜᵖ[colidx]
+    @. Yₜ.c.ρe_tot[colidx] += Y.c.ρ[colidx] * ᶜSeₜᵖ[colidx]
+    @. Yₜ.c.ρq_rai[colidx] += Y.c.ρ[colidx] * ᶜSqᵣᵖ[colidx]
+    @. Yₜ.c.ρq_sno[colidx] += Y.c.ρ[colidx] * ᶜSqₛᵖ[colidx]
+
+    return nothing
+end
diff --git a/src/parameters/create_parameters.jl b/src/parameters/create_parameters.jl
index 65c3e55a3e..fb438de6f6 100644
--- a/src/parameters/create_parameters.jl
+++ b/src/parameters/create_parameters.jl
@@ -46,6 +46,16 @@ function create_parameter_set(config::AtmosConfig)
         nothing
     elseif precip_model == "0M"
         CM.Parameters.Parameters0M(FT, toml_dict)
+    elseif precip_model == "1M"
+        (;
+            cl = CM.Parameters.CloudLiquid(FT, toml_dict),
+            ci = CM.Parameters.CloudIce(FT, toml_dict),
+            pr = CM.Parameters.Rain(FT, toml_dict),
+            ps = CM.Parameters.Snow(FT, toml_dict),
+            ce = CM.Parameters.CollisionEff(FT, toml_dict),
+            tv = CM.Parameters.Blk1MVelType(FT, toml_dict),
+            aps = CM.Parameters.AirProperties(FT, toml_dict),
+        )
     else
         error("Invalid precip_model $(precip_model)")
     end
diff --git a/src/prognostic_equations/implicit/implicit_solver.jl b/src/prognostic_equations/implicit/implicit_solver.jl
index 3ef3ac671f..e8e3eb1332 100644
--- a/src/prognostic_equations/implicit/implicit_solver.jl
+++ b/src/prognostic_equations/implicit/implicit_solver.jl
@@ -241,7 +241,8 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx)
     (; matrix, enthalpy_flag) = A
     (; ᶜspecific, ᶠu³, ᶜK, ᶜp, ∂ᶜK_∂ᶠu₃) = p
     (; ᶜΦ, ᶠgradᵥ_ᶜΦ, ᶜρ_ref, ᶜp_ref, params) = p
-    (; energy_upwinding, tracer_upwinding, density_upwinding) = p.atmos.numerics
+    (; energy_upwinding, density_upwinding) = p.atmos.numerics
+    (; tracer_upwinding, precip_upwinding) = p.atmos.numerics
 
     FT = Spaces.undertype(axes(Y.c))
     one_ATC3 = CT3(FT(1))'
@@ -320,6 +321,8 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx)
         (@name(c.ρq_ice), @name(ᶜspecific.q_ice), tracer_upwinding),
         (@name(c.ρq_rai), @name(ᶜspecific.q_rai), tracer_upwinding),
         (@name(c.ρq_sno), @name(ᶜspecific.q_sno), tracer_upwinding),
+        (@name(c.ρq_rai), @name(ᶜspecific.q_rai), precip_upwinding),
+        (@name(c.ρq_sno), @name(ᶜspecific.q_sno), precip_upwinding),
     )
     available_tracer_info =
         MatrixFields.unrolled_filter(tracer_info) do (ρχ_name, _, _)
diff --git a/src/prognostic_equations/implicit/implicit_tendency.jl b/src/prognostic_equations/implicit/implicit_tendency.jl
index 1eb6b34bcf..5a4db2cef3 100644
--- a/src/prognostic_equations/implicit/implicit_tendency.jl
+++ b/src/prognostic_equations/implicit/implicit_tendency.jl
@@ -23,28 +23,30 @@ end
 # the implicit tendency function. Since dt >= dtγ, we can safely use dt for now.
 # TODO: Can we rewrite ᶠfct_boris_book and ᶠfct_zalesak so that their broadcast
 # expressions are less convoluted?
-vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:none}) =
-    @. ᶜρχₜ += -(ᶜadvdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠu³ * ᶠinterp(ᶜχ)))
-vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:first_order}) =
-    @. ᶜρχₜ += -(ᶜadvdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ)))
-vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:third_order}) =
-    @. ᶜρχₜ += -(ᶜadvdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind3(ᶠu³, ᶜχ)))
-vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:boris_book}) =
-    @. ᶜρχₜ += -(ᶜadvdivᵥ(
+vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding) =
+    vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding, ᶜadvdivᵥ)
+vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:none}, ᶜdivᵥ) =
+    @. ᶜρχₜ += -(ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠu³ * ᶠinterp(ᶜχ)))
+vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:first_order}, ᶜdivᵥ) =
+    @. ᶜρχₜ += -(ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ)))
+vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:third_order}, ᶜdivᵥ) =
+    @. ᶜρχₜ += -(ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind3(ᶠu³, ᶜχ)))
+vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:boris_book}, ᶜdivᵥ) =
+    @. ᶜρχₜ += -(ᶜdivᵥ(
         ᶠwinterp(ᶜJ, ᶜρ) * (
             ᶠupwind1(ᶠu³, ᶜχ) + ᶠfct_boris_book(
                 ᶠupwind3(ᶠu³, ᶜχ) - ᶠupwind1(ᶠu³, ᶜχ),
-                ᶜχ / dt - ᶜadvdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ)) / ᶜρ,
+                ᶜχ / dt - ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ)) / ᶜρ,
             )
         ),
     ))
-vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:zalesak}) =
-    @. ᶜρχₜ += -(ᶜadvdivᵥ(
+vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:zalesak}, ᶜdivᵥ) =
+    @. ᶜρχₜ += -(ᶜdivᵥ(
         ᶠwinterp(ᶜJ, ᶜρ) * (
             ᶠupwind1(ᶠu³, ᶜχ) + ᶠfct_zalesak(
                 ᶠupwind3(ᶠu³, ᶜχ) - ᶠupwind1(ᶠu³, ᶜχ),
                 ᶜχ / dt,
-                ᶜχ / dt - ᶜadvdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ)) / ᶜρ,
+                ᶜχ / dt - ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ)) / ᶜρ,
             )
         ),
     ))
@@ -57,8 +59,9 @@ vertical_advection!(ᶜρχₜ, ᶠu³, ᶜχ, ::Val{:third_order}) =
     @. ᶜρχₜ -= ᶜadvdivᵥ(ᶠupwind3(ᶠu³, ᶜχ)) - ᶜχ * ᶜadvdivᵥ(ᶠu³)
 
 function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx)
-    (; energy_upwinding, tracer_upwinding, density_upwinding) = p.atmos.numerics
-    (; turbconv_model, rayleigh_sponge) = p.atmos
+    (; energy_upwinding, density_upwinding) = p.atmos.numerics
+    (; tracer_upwinding, precip_upwinding) = p.atmos.numerics
+    (; turbconv_model, rayleigh_sponge, precip_model) = p.atmos
     (; dt) = p.simulation
     n = n_mass_flux_subdomains(turbconv_model)
     ᶜJ = Fields.local_geometry_field(Y.c).J
@@ -102,6 +105,56 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx)
         )
     end
 
+    if precip_model isa Microphysics1Moment
+        # Advection of precipitation with the mean flow
+        # is done with other tracers above.
+        # Here we add the advection with precipitation terminal velocity
+        # using first order upwind and free outflow bottom boundary condition
+
+        ᶠu³ₚ = p.scratch.ᶠtemp_CT3
+        ᶜqₚ = p.scratch.ᶜtemp_scalar
+        lgf = Fields.local_geometry_field(Y.f)
+        FT = Spaces.undertype(axes(Y.c))
+
+        @. ᶠu³ₚ[colidx] =
+            FT(-1) *
+            ᶠinterp(p.precipitation.ᶜwᵣ[colidx]) *
+            CT3(unit_basis_vector_data(CT3, lgf[colidx]))
+        @. ᶜqₚ[colidx] = Y.c.ρq_rai[colidx] / Y.c.ρ[colidx]
+
+        ᶜdivᵥ_ρqₚ = Operators.DivergenceF2C(
+            top = Operators.SetValue(C3(FT(0))),
+            bottom = Operators.Extrapolate(),
+        )
+
+        vertical_transport!(
+            Yₜ.c.ρq_rai[colidx],
+            ᶜJ[colidx],
+            Y.c.ρ[colidx],
+            ᶠu³ₚ[colidx],
+            ᶜqₚ[colidx],
+            dt,
+            precip_upwinding,
+            ᶜdivᵥ_ρqₚ,
+        )
+
+        @. ᶠu³ₚ[colidx] =
+            FT(-1) *
+            ᶠinterp(p.precipitation.ᶜwₛ[colidx]) *
+            CT3(unit_basis_vector_data(CT3, lgf[colidx]))
+        @. ᶜqₚ[colidx] = Y.c.ρq_sno[colidx] / Y.c.ρ[colidx]
+        vertical_transport!(
+            Yₜ.c.ρq_sno[colidx],
+            ᶜJ[colidx],
+            Y.c.ρ[colidx],
+            ᶠu³ₚ[colidx],
+            ᶜqₚ[colidx],
+            dt,
+            precip_upwinding,
+            ᶜdivᵥ_ρqₚ,
+        )
+    end
+
     @. Yₜ.f.u₃[colidx] +=
         -(
             ᶠgradᵥ(ᶜp[colidx] - ᶜp_ref[colidx]) +
diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl
index 362f1928a0..49c66c8cb1 100644
--- a/src/solver/type_getters.jl
+++ b/src/solver/type_getters.jl
@@ -86,6 +86,7 @@ function get_numerics(parsed_args)
 
     energy_upwinding = Val(Symbol(parsed_args["energy_upwinding"]))
     tracer_upwinding = Val(Symbol(parsed_args["tracer_upwinding"]))
+    precip_upwinding = Val(Symbol(parsed_args["precip_upwinding"]))
     density_upwinding = Val(Symbol(parsed_args["density_upwinding"]))
     edmfx_upwinding = Val(Symbol(parsed_args["edmfx_upwinding"]))
     edmfx_sgsflux_upwinding =
@@ -98,6 +99,7 @@ function get_numerics(parsed_args)
     numerics = AtmosNumerics(;
         energy_upwinding,
         tracer_upwinding,
+        precip_upwinding,
         density_upwinding,
         edmfx_upwinding,
         edmfx_sgsflux_upwinding,
@@ -323,6 +325,7 @@ function get_initial_condition(parsed_args)
         "DryDensityCurrentProfile",
         "RisingThermalBubbleProfile",
         "ScharProfile",
+        "PrecipitatingColumn",
     ]
         return getproperty(ICs, Symbol(parsed_args["initial_condition"]))()
     else
diff --git a/src/solver/types.jl b/src/solver/types.jl
index ab882bcba8..442388434f 100644
--- a/src/solver/types.jl
+++ b/src/solver/types.jl
@@ -265,6 +265,7 @@ Base.broadcastable(x::AbstractPerformanceMode) = tuple(x)
 Base.@kwdef struct AtmosNumerics{
     EN_UP,
     TR_UP,
+    PR_UP,
     DE_UP,
     ED_UP,
     ED_SG_UP,
@@ -275,6 +276,7 @@ Base.@kwdef struct AtmosNumerics{
     """Enable specific upwinding schemes for specific equations"""
     energy_upwinding::EN_UP
     tracer_upwinding::TR_UP
+    precip_upwinding::PR_UP
     density_upwinding::DE_UP
     edmfx_upwinding::ED_UP
     edmfx_sgsflux_upwinding::ED_SG_UP
@@ -469,7 +471,7 @@ AtmosConfig(::Nothing; comms_ctx = nothing) = AtmosConfig(Dict(); comms_ctx)
 
 """
     AtmosConfig(config::Dict; comms_ctx = nothing)
-Constructs the AtmosConfig from the Dict passed in. This Dict overrides all of 
+Constructs the AtmosConfig from the Dict passed in. This Dict overrides all of
 the default configurations set in `default_config_dict()`.
 """
 function AtmosConfig(config::Dict; comms_ctx = nothing)
diff --git a/toml/single_column_precipitation_test.toml b/toml/single_column_precipitation_test.toml
new file mode 100644
index 0000000000..e39b3da960
--- /dev/null
+++ b/toml/single_column_precipitation_test.toml
@@ -0,0 +1,4 @@
+[C_H]
+alias = "C_H"
+value = 0.0
+type = "float"