precip units in [m] to [mm] for output

src/SpeedyWeather.jl

@@ -168,14 +168,14 @@ include("physics/thermodynamics.jl")
 include("physics/boundary_layer.jl")
 include("physics/temperature_relaxation.jl")
 include("physics/vertical_diffusion.jl")
+include("physics/large_scale_condensation.jl")
 include("physics/pretty_printing.jl")
 
 # MODELS
 include("dynamics/models.jl")
 
 # # PHYSICS
 # include("physics/convection.jl")
-include("physics/large_scale_condensation.jl")
 # include("physics/longwave_radiation.jl")
 # include("physics/shortwave_radiation.jl")
 

src/output/output.jl

@@ -442,9 +442,9 @@ function write_netcdf_variables!(   output::OutputWriter,
     precip_large_scale .= 0
     precip_convection .= 0
 
-    # convert from [m/s] to [mm/Δt], where Δt is the output time step
-    precip_cond *= 1000*output.output_dt_sec
-    precip_conv *= 1000*output.output_dt_sec
+    # convert from [m] to [mm] within output time step (e.g. 6hours)
+    precip_cond *= 1000
+    precip_conv *= 1000
 
     if Model <: PrimitiveEquation
         @. pres = exp(pres)/100     # convert from log(pₛ) to surface pressure pₛ [hPa]

src/physics/column_variables.jl

@@ -1,6 +1,5 @@
 """
-    get_column!(C,D,ij,G)
-
+$(TYPEDSIGNATURES)
 Update `C::ColumnVariables` by copying the prognostic variables from `D::DiagnosticVariables`
 at gridpoint index `ij`. Provide `G::Geometry` for coordinate information."""
 function get_column!(   C::ColumnVariables,
@@ -31,6 +30,7 @@ function get_column!(   C::ColumnVariables,
         C.humid[k] = layer.grid_variables.humid_grid[ij]
 
         # as well as geopotential (not actually prognostic though)
+        # TODO geopot on the grid is currently not computed in dynamics
         C.geopot[k] = layer.grid_variables.geopot_grid[ij]
     end
 end
@@ -47,8 +47,7 @@ function get_column!(   C::ColumnVariables,
 end
 
 """
-    write_column_tendencies!(D,C,ij,G)
-
+$(TYPEDSIGNATURES)
 Write the parametrization tendencies from `C::ColumnVariables` into the horizontal fields
 of tendencies stored in `D::DiagnosticVariables` at gridpoint index `ij`."""
 function write_column_tendencies!(  D::DiagnosticVariables,
@@ -72,27 +71,26 @@ function write_column_tendencies!(  D::DiagnosticVariables,
 end
 
 """
-    reset_column!(column::ColumnVariables)
-
+$(TYPEDSIGNATURES)
 Set the accumulators (tendencies but also vertical sums and similar) back to zero
 for `column` to be reused at other grid points."""
 function reset_column!(column::ColumnVariables{NF}) where NF
 
     # set tendencies to 0 for += accumulation
-    fill!(column.u_tend,0)
-    fill!(column.v_tend,0)
-    fill!(column.temp_tend,0)
-    fill!(column.humid_tend,0)
+    column.u_tend .= 0
+    column.v_tend .= 0
+    column.temp_tend .= 0
+    column.humid_tend .= 0
 
     # set fluxes to 0 for += accumulation
-    fill!(column.flux_u_upward,0)
-    fill!(column.flux_u_downward,0)
-    fill!(column.flux_v_upward,0)
-    fill!(column.flux_v_downward,0)
-    fill!(column.flux_humid_upward,0)
-    fill!(column.flux_humid_downward,0)
-    fill!(column.flux_temp_upward,0)
-    fill!(column.flux_temp_downward,0)
+    column.flux_u_upward .= 0
+    column.flux_u_downward .= 0
+    column.flux_v_upward .= 0
+    column.flux_v_downward .= 0
+    column.flux_humid_upward .= 0
+    column.flux_humid_downward .= 0
+    column.flux_temp_upward .= 0
+    column.flux_temp_downward .= 0
 
     # # Convection
     # column.cloud_top = column.nlev+1

src/physics/define_column.jl

@@ -17,14 +17,14 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV
     latd::NF = 0                            # latitude, needed for shortwave radiation
 
     # PROGNOSTIC VARIABLES
-    const u::Vector{NF} = zeros(NF,nlev)            # zonal velocity
-    const v::Vector{NF} = zeros(NF,nlev)            # meridional velocity
-    const temp::Vector{NF} = zeros(NF,nlev)         # temperature
-    const humid::Vector{NF} = zeros(NF,nlev)        # specific humidity
+    const u::Vector{NF} = zeros(NF,nlev)            # zonal velocity [m/s]
+    const v::Vector{NF} = zeros(NF,nlev)            # meridional velocity [m/s]
+    const temp::Vector{NF} = zeros(NF,nlev)         # temperature [K]
+    const humid::Vector{NF} = zeros(NF,nlev)        # specific humidity [kg/kg]
 
     # (log) pressure per layer, surface is prognostic, last element here, but precompute other layers too
-    const ln_pres::Vector{NF} = zeros(NF,nlev+1)    # logarithm of pressure [log(hPa)]
-    const pres::Vector{NF} = zeros(NF,nlev+1)       # pressure [hPa]
+    const ln_pres::Vector{NF} = zeros(NF,nlev+1)    # logarithm of pressure [log(Pa)]
+    const pres::Vector{NF} = zeros(NF,nlev+1)       # pressure [Pa]
 
     # TENDENCIES to accumulate the parametrizations into
     const u_tend::Vector{NF} = zeros(NF,nlev)                   # zonal velocity [m]
@@ -49,8 +49,8 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV
     const flux_humid_downward::Vector{NF} = zeros(NF,nlev+1)
 
     # THERMODYNAMICS
-    const sat_humid::Vector{NF} = zeros(NF,nlev)                # Saturation specific humidity
-    const sat_vap_pres::Vector{NF} = zeros(NF,nlev)             # Saturation vapour pressure
+    const sat_humid::Vector{NF} = zeros(NF,nlev)                # Saturation specific humidity [kg/kg]
+    const sat_vap_pres::Vector{NF} = zeros(NF,nlev)             # Saturation vapour pressure [Pa]
     const dry_static_energy::Vector{NF} = zeros(NF,nlev)        # Dry static energy
     const moist_static_energy::Vector{NF} = zeros(NF,nlev)      # Moist static energy
 
@@ -62,18 +62,18 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV
     sat_moist_static_energy_half::Vector{NF} = zeros(NF,nlev)  # Saturation moist static energy interpolated to half-levels
 
     # Convection
-    conditional_instability::Bool = false                    # Whether a conditional instability exists in this column (condition 1)
-    activate_convection::Bool = false                        # Whether convection should be activated in this column (condition 2)
-    cloud_top::Int = nlev+1                                  # Top-of-convection layer
-    excess_humidity::NF = 0                                  # Excess humidity due to convection
-    cloud_base_mass_flux::NF = 0                             # Mass flux at the top of the PBL
-    precip_convection::NF = 0                                # Precipitation due to convection
+    conditional_instability::Bool = false                       # in this column? (condition 1)
+    activate_convection::Bool = false                           # in this column? (condition 2)
+    cloud_top::Int = nlev+1                                     # Top-of-convection layer
+    excess_humidity::NF = 0                                     # due to convection
+    cloud_base_mass_flux::NF = 0                                # Mass flux at the top of the PBL
+    precip_convection::NF = 0                                   # Precipitation due to convection [m]
     net_flux_humid::Vector{NF} = zeros(NF,nlev)              # Net fluxes of moisture in this column
     net_flux_dry_static_energy::Vector{NF} = zeros(NF,nlev)  # Net fluxes of dry static energy in this column
     entrainment_profile::Vector{NF} = zeros(NF,nlev)         # Entrainment coefficients
 
     # Large-scale condensation
-    precip_large_scale::NF = 0  # Precipitation due to large-scale condensation
+    precip_large_scale::NF = 0                                  # precipitation [m]
 
     # Longwave radiation
     ## New vars in radlw_down!

src/physics/large_scale_condensation.jl

@@ -71,7 +71,7 @@ function large_scale_condensation!(
     model::PrimitiveWet,
 )
     large_scale_condensation!(column,model.large_scale_condensation,
-        model.geometry,model.constants,model.atmosphere)
+        model.geometry,model.constants,model.atmosphere,model.time_stepping)
 end
 
 
@@ -81,6 +81,7 @@ function large_scale_condensation!(
     geometry::Geometry,
     constants::DynamicsConstants,
     atmosphere::AbstractAtmosphere,
+    time_stepping::TimeStepper,
     ) where NF
 
     (;relative_threshold,humid_tend_max) = scheme
@@ -94,6 +95,9 @@ function large_scale_condensation!(
     pₛ = pres[end]                         # surface pressure
 
     (;gravity, water_density) = constants
+    (;Δt_sec) = time_stepping
+    Δt_gρ = Δt_sec/gravity/water_density
+
     (;σ_levels_thick) = geometry
     latent_heat = convert(NF, atmosphere.latent_heat_condensation/atmosphere.cₚ)
 
@@ -106,16 +110,16 @@ function large_scale_condensation!(
         if humid[k] > humid_threshold
             # accumulate in tendencies (nothing is added if humidity not above threshold)
             humid_tend_k = -(humid[k] - humid_threshold) / time_scale                   # Formula 22
-            temp_tend[k] += -latent_heat * min(humid_tend_k, humid_tend_max[k]*pres[k])   # Formula 23
+            temp_tend[k] += -latent_heat * min(humid_tend_k, humid_tend_max[k]*pres[k]) # Formula 23
 
             # If there is large-scale condensation at a level higher (i.e. smaller k) than
             # the cloud-top previously diagnosed due to convection, then increase the cloud-top
-            column.cloud_top = min(column.cloud_top, k)             # Page 7 (last sentence)
+            column.cloud_top = min(column.cloud_top, k)         # Page 7 (last sentence)
 
             # 2. Precipitation due to large-scale condensation [kg/m²/s] /ρ for [m/s]
             # += for vertical integral
-            Δpₖ_g = pₛ*σ_levels_thick[k]/gravity/water_density      # Formula 4
-            column.precip_large_scale += -Δpₖ_g * humid_tend_k       # Formula 25
+            Δpₖ_g = pₛ*σ_levels_thick[k]*Δt_gρ                  # Formula 4 *Δt for [m] of rain during Δt
+            column.precip_large_scale += -Δpₖ_g * humid_tend_k  # Formula 25, unit [m]
 
             # only write into humid_tend now to allow humid_tend != 0 before this scheme is called
             humid_tend[k] += humid_tend_k