Precipitation rate for coupling (#669)

* precipitation rate in m/s

* Precipitation rate for convection

* update changelog

* set convection rate back to zero for every column

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Unreleased
 
+- Precipitation rate from large-scale condensation/convection for coupling [#669](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/669)
 - Fixed a bug, so that sum now works for LowerTraingularArrays [#668](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/668)
 - Require only AbstractVectors in the interpolation [#665](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/665)
 - Intitialize! simulation restructured [#666](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/666)

src/dynamics/diagnostic_variables.jl

@@ -266,10 +266,10 @@ $(TYPEDFIELDS)"""
     "Accumulated large-scale precipitation [m]"
     precip_convection::GridVariable2D = zeros(GridVariable2D, nlat_half)
 
-    "Rate of large-scale precipitation [m], over time step"
+    "Rate of large-scale precipitation [m/s], instantaneous"
     precip_rate_large_scale::GridVariable2D = zeros(GridVariable2D, nlat_half)
 
-    "Rate of large-scale precipitation [m], over time step"
+    "Rate of large-scale precipitation [m/s], instantaneous"
     precip_rate_convection::GridVariable2D = zeros(GridVariable2D, nlat_half)
 
     "Cloud top [m]"

src/physics/column_variables.jl

@@ -136,9 +136,9 @@ function write_column_tendencies!(
     diagn.physics.precip_large_scale[ij] += column.precip_large_scale
     diagn.physics.precip_convection[ij] += column.precip_convection
 
-    # precipitation rate [m] over this time step (i.e. overwrite, do not accumulate)
-    diagn.physics.precip_rate_large_scale[ij] = column.precip_large_scale
-    diagn.physics.precip_rate_convection[ij] = column.precip_convection
+    # precipitation rate [m/s], instantaneous (i.e. overwrite, do not accumulate)
+    diagn.physics.precip_rate_large_scale[ij] = column.precip_rate_large_scale
+    diagn.physics.precip_rate_convection[ij] = column.precip_rate_convection
 
     # Cloud top in height [m] from geopotential height divided by gravity, 0 for no clouds
     diagn.physics.cloud_top[ij] = column.cloud_top == nlayers+1 ? 0 : column.geopot[column.cloud_top]
@@ -186,8 +186,10 @@ function reset_column!(column::ColumnVariables{NF}) where NF
 
     # Convection and precipitation
     column.cloud_top = column.nlayers+1
-    column.precip_convection = 0
+    column.precip_convection = 0            # set back to zero to accumulate in the vertical
     column.precip_large_scale = 0
+    column.precip_rate_convection = 0       # instantaneously overwritten, but convection may escape early
+    column.precip_rate_large_scale = 0
 
     # radiation
     column.outgoing_longwave_radiation = 0

src/physics/convection.jl

@@ -203,7 +203,9 @@ function convection!(
     (; Δt_sec) = time_stepping
     pₛΔt_gρ = (pₛ * Δt_sec / gravity / water_density) * deep_convection # enfore no precip for shallow conv 
     column.precip_convection *= pₛΔt_gρ                                 # convert to [m] of rain during Δt
+    column.precip_rate_convection = column.precip_convection / Δt_sec   # rate: convert to [m/s] of rain
     column.cloud_top = min(column.cloud_top, level_zero_buoyancy)       # clouds reach to top of convection
+    return nothing
 end
 
 """

src/physics/define_column.jl

@@ -86,6 +86,8 @@ $(TYPEDFIELDS)"""
     cloud_top::Int = nlayers+1              # layer index k of top-most layer with clouds
     precip_convection::NF = 0               # Precipitation due to convection [m]
     precip_large_scale::NF = 0              # precipitation due to large-scale condensation [m]
+    precip_rate_convection::NF = 0          # Precipitation rate due to convection [m/s]
+    precip_rate_large_scale::NF = 0         # precipitation rate due to large-scale condensation [m/s]
 
     # RADIATION
     cos_zenith::NF = 0                      # cosine of solar zenith angle [1]

src/physics/large_scale_condensation.jl

@@ -75,7 +75,7 @@ function large_scale_condensation!(
     pₛ = pres[end]                          # surface pressure
     (; Δt_sec) = time_stepping
     Δσ = geometry.σ_levels_thick
-    pₛΔt_gρ = pₛ*Δt_sec/(planet.gravity * atmosphere.water_density)
+    pₛΔt_gρ = (pₛ * Δt_sec)/(planet.gravity * atmosphere.water_density)
 
     (; Lᵥ, cₚ, Lᵥ_Rᵥ) = clausius_clapeyron
     Lᵥ_cₚ = Lᵥ/cₚ                           # latent heat of vaporization over heat capacity
@@ -101,12 +101,16 @@ function large_scale_condensation!(
 
             # 2. Precipitation due to large-scale condensation [kg/m²/s] /ρ for [m/s]
             # += for vertical integral
-            ΔpₖΔt_gρ = Δσ[k] * pₛΔt_gρ                      # Formula 4 *Δt for [m] of rain during Δt
-            column.precip_large_scale -= ΔpₖΔt_gρ * δq      # Formula 25, unit [m]
+            precip = Δσ[k] * pₛΔt_gρ * -δq          # precipitation [m] on layer k, Formula 4
+            column.precip_large_scale += precip     # integrate vertically, Formula 25, unit [m]
 
             # only accumulate into humid_tend now to allow humid_tend != 0 before this scheme is called
             humid_tend[k] += δq
             temp_tend[k] += δT
         end
     end
+
+    # convert to rain rate [m/s]
+    column.precip_rate_large_scale = column.precip_large_scale / Δt_sec
+    return nothing
 end