Merge pull request #541 from SpeedyWeather/mk/lsc_prev

Large-scale condensation at previous time step

docs/src/examples_3D.md

@@ -212,6 +212,72 @@ And the comparison looks like
 ![Aquaplanet, no deep convection](aquaplanet_nodeepconvection.png)
 ![Aquaplanet, no convection](aquaplanet_noconvection.png)
 
+## Large-scale vs convective precipitation
+
+```@example precipitation
+using SpeedyWeather
+
+# components
+spectral_grid = SpectralGrid(trunc=31, nlev=8)
+large_scale_condensation = ImplicitCondensation(spectral_grid)
+convection = SimplifiedBettsMiller(spectral_grid)
+
+# create model, initialize, run
+model = PrimitiveWetModel(; spectral_grid, large_scale_condensation, convection)
+simulation = initialize!(model)
+model.feedback.verbose = false # hide
+run!(simulation, period=Day(10))
+nothing # hide
+```
+
+We run the default `PrimitiveWetModel` with `ImplicitCondensation` as large-scale condensation
+(see [Implicit large-scale condensation](@ref)) and the `SimplifiedBettsMiller`
+for convection (see [Simplified Betts-Miller](@ref BettsMiller)). These schemes
+have some additional parameters, we leave them as default for now, but you could
+do `ImplicitCondensation(spectral_grid, relative_humidity_threshold = 0.8)` to
+let it rain at 80% instead of 100% relative humidity. We now want to analyse
+the precipitation that comes from these parameterizations
+
+```@example precipitation
+using CairoMakie
+
+(; precip_large_scale, precip_convection) = simulation.diagnostic_variables.surface
+m2mm = 1000     # convert from [m] to [mm]
+heatmap(m2mm*precip_large_scale, title="Large-scale precipiation [mm]: Accumulated over 10 days", colormap=:dense)
+save("large-scale_precipitation_acc.png", ans) # hide
+nothing # hide
+```
+![Large-scale precipitation](large-scale_precipitation_acc.png)
+
+Precipitation (both large-scale and convective) are written into the
+`simulation.diagnostic_variables.surface` which, however, accumulate all precipitation
+until netCDF output is written. As we didn't specify `output=true` here, they accumulated
+precipitation over the last 10 days of the simulation which includes the initial
+adjustment of humidity in the tropics (the large band of precipitation here).
+So let us reset these accumulators and integrate for another 6 hours to get the
+precipitation only in the period.
+
+```@example precipitation
+# reset accumulators and simulate 6 hours
+simulation.diagnostic_variables.surface.precip_large_scale .= 0
+simulation.diagnostic_variables.surface.precip_convection .= 0
+run!(simulation, period=Hour(6))
+
+# visualise, precip_* arrays are flat copies, no need to read them out again!
+m2mm_hr = (1000*Hour(1)/Hour(6))    # convert from [m] to [mm/hr]
+heatmap(m2mm_hr*precip_large_scale, title="Large-scale precipiation [mm/hr]", colormap=:dense)
+save("large-scale_precipitation.png", ans) # hide
+heatmap(m2mm_hr*precip_convection, title="Convective precipiation [mm/hr]", colormap=:dense)
+save("convective_precipitation.png", ans) # hide
+nothing # hide
+```
+![Large-scale precipitation](large-scale_precipitation.png)
+![Convective precipitation](convective_precipitation.png)
+
+As the precipitation fields are accumulated meters over the integration period
+(in the case of no output) we divide by 6 hours to get a precipitation rate ``[m/s]``
+but then multiply with 1 hour and 1000 to get the typical precipitation unit of ``[mm/hr]``.
+
 ## References
 
 [^JW06]: Jablonowski, C. and Williamson, D.L. (2006), A baroclinic instability test case for atmospheric model dynamical cores. Q.J.R. Meteorol. Soc., 132: 2943-2975. DOI:[10.1256/qj.06.12](https://doi.org/10.1256/qj.06.12)

src/dynamics/diagnostic_variables.jl

@@ -35,7 +35,8 @@ Base.@kwdef struct GridVariables{NF<:AbstractFloat, Grid<:AbstractGrid{NF}}
     temp_grid       ::Grid = zeros(Grid, nlat_half)  # absolute temperature [K]
     temp_grid_prev  ::Grid = zeros(Grid, nlat_half)  # absolute temperature of previous time step [K]
     temp_virt_grid  ::Grid = zeros(Grid, nlat_half)  # virtual tempereature [K]  
-    humid_grid      ::Grid = zeros(Grid, nlat_half)  # specific_humidity
+    humid_grid      ::Grid = zeros(Grid, nlat_half)  # specific_humidity [kg/kg]
+    humid_grid_prev ::Grid = zeros(Grid, nlat_half)  # specific_humidity at previous time step
     u_grid          ::Grid = zeros(Grid, nlat_half)  # zonal velocity *coslat [m/s]
     v_grid          ::Grid = zeros(Grid, nlat_half)  # meridional velocity *coslat [m/s]
     u_grid_prev     ::Grid = zeros(Grid, nlat_half)  # zonal velocity *coslat of previous time step [m/s]

src/dynamics/tendencies.jl

@@ -766,12 +766,14 @@ function SpeedyTransforms.gridded!(
     diagn::DiagnosticVariables,
     progn::PrognosticVariables,
     lf::Int,
-    model::ModelSetup)
+    model::ModelSetup;
+    kwargs...
+)
 
     # all variables on layers
     for (progn_layer, diagn_layer) in zip(progn.layers, diagn.layers)
         progn_layer_lf = progn_layer.timesteps[lf]
-        gridded!(diagn_layer, progn_layer_lf, model)
+        gridded!(diagn_layer, progn_layer_lf, model; kwargs...)
     end
 
     # surface only for ShallowWaterModel or PrimitiveEquation
@@ -788,9 +790,12 @@ $(TYPEDSIGNATURES)
 Propagate the spectral state of the prognostic variables `progn` to the
 diagnostic variables in `diagn` for the barotropic vorticity model.
 Updates grid vorticity, spectral stream function and spectral and grid velocities u, v."""
-function SpeedyTransforms.gridded!( diagn::DiagnosticVariablesLayer,   
-                                    progn::PrognosticVariablesLayer,
-                                    model::Barotropic)
+function SpeedyTransforms.gridded!( 
+    diagn::DiagnosticVariablesLayer,   
+    progn::PrognosticVariablesLayer,
+    model::Barotropic;
+    kwargs...
+)
 
     (; vor_grid, u_grid, v_grid ) = diagn.grid_variables
     (; vor ) = progn                    # relative vorticity
@@ -818,10 +823,12 @@ Propagate the spectral state of the prognostic variables `progn` to the
 diagnostic variables in `diagn` for the shallow water model. Updates grid vorticity,
 grid divergence, grid interface displacement (`pres_grid`) and the velocities
 u, v."""
-function SpeedyTransforms.gridded!( diagn::DiagnosticVariablesLayer,
-                                    progn::PrognosticVariablesLayer,
-                                    model::ShallowWater,                # everything that's constant
-                                    )
+function SpeedyTransforms.gridded!( 
+    diagn::DiagnosticVariablesLayer,
+    progn::PrognosticVariablesLayer,
+    model::ShallowWater;
+    kwargs...
+)
 
     (; vor_grid, div_grid, u_grid, v_grid ) = diagn.grid_variables
     (; vor, div) = progn
@@ -850,21 +857,27 @@ Propagate the spectral state of the prognostic variables `progn` to the
 diagnostic variables in `diagn` for primitive equation models. Updates grid vorticity,
 grid divergence, grid temperature, pressure (`pres_grid`) and the velocities
 u, v."""
-function SpeedyTransforms.gridded!( diagn::DiagnosticVariablesLayer,
-                                    progn::PrognosticVariablesLayer,
-                                    model::PrimitiveEquation)           # everything that's constant
+function SpeedyTransforms.gridded!( 
+    diagn::DiagnosticVariablesLayer,
+    progn::PrognosticVariablesLayer,
+    model::PrimitiveEquation;
+    initialize=false
+)
 
     (; vor_grid, div_grid, u_grid, v_grid ) = diagn.grid_variables
     (; temp_grid, humid_grid ) = diagn.grid_variables
-    (; temp_grid_prev, u_grid_prev, v_grid_prev) = diagn.grid_variables
+    (; temp_grid_prev, humid_grid_prev, u_grid_prev, v_grid_prev) = diagn.grid_variables
     (; vor, div, temp, humid) = progn
     U = diagn.dynamics_variables.a      # reuse work arrays for velocities spectral
     V = diagn.dynamics_variables.b      # U = u*coslat, V=v*coslat
 
-    # retain previous time step for vertical advection
-    @. temp_grid_prev = temp_grid
-    @. u_grid_prev = u_grid
-    @. v_grid_prev = v_grid
+    # retain previous time step for vertical advection and some parameterizations
+    if initialize == false              # only store prev after initial step
+        @. temp_grid_prev = temp_grid   # this is temperature anomaly wrt to implicit reference profile!
+        @. humid_grid_prev = humid_grid
+        @. u_grid_prev = u_grid
+        @. v_grid_prev = v_grid
+    end
 
     S = model.spectral_transform
     wet_core = model isa PrimitiveWet
@@ -874,23 +887,31 @@ function SpeedyTransforms.gridded!( diagn::DiagnosticVariablesLayer,
     # V = v*coslat =  coslat*∂ϕ/∂lat + ∂Ψ/dlon
     UV_from_vordiv!(U, V, vor, div, S)
 
-    gridded!(vor_grid, vor, S)                # get vorticity on grid from spectral vor
-    gridded!(div_grid, div, S)                # get divergence on grid from spectral div
-    gridded!(temp_grid, temp, S)              # (absolute) temperature
+    gridded!(vor_grid, vor, S)                  # get vorticity on grid from spectral vor
+    gridded!(div_grid, div, S)                  # get divergence on grid from spectral div
+    gridded!(temp_grid, temp, S)                # (absolute) temperature
 
-    if wet_core                             # specific humidity (wet core only)
+    if wet_core                                 # specific humidity (wet core only)
         gridded!(humid_grid, humid, S)        
         hole_filling!(humid_grid, model.hole_filling, model)  # remove negative humidity
     end
 
     # include humidity effect into temp for everything stability-related
     temperature_average!(diagn, temp, S)
-    virtual_temperature!(diagn, temp, model)  # temp = virt temp for dry core
+    virtual_temperature!(diagn, temp, model)    # temp = virt temp for dry core
 
     # transform from U, V in spectral to u, v on grid (U, V = u, v*coslat)
     gridded!(u_grid, U, S, unscale_coslat=true)
     gridded!(v_grid, V, S, unscale_coslat=true)
 
+    if initialize   # at initial step store prev <- current
+        # subtract the reference temperature profile as temp_grid is too after every time step
+        @. temp_grid_prev = temp_grid - model.implicit.temp_profile[diagn.k]
+        @. humid_grid_prev = humid_grid
+        @. u_grid_prev = u_grid
+        @. v_grid_prev = v_grid
+    end
+
     return nothing
 end
 
@@ -903,7 +924,6 @@ function temperature_average!(
     temp::LowerTriangularMatrix,
     S::SpectralTransform,
 )
-
     # average from l=m=0 harmonic divided by norm of the sphere
     diagn.temp_average[] = real(temp[1])/S.norm_sphere
 end 

src/dynamics/time_integration.jl

@@ -392,7 +392,7 @@ function time_stepping!(
     # propagate spectral state to grid variables for initial condition output
     (; output, feedback) = model
     lf = 1                                  # use first leapfrog index
-    gridded!(diagn, progn, lf, model)
+    gridded!(diagn, progn, lf, model, initialize=true)
     initialize!(progn.particles, progn, diagn, model.particle_advection)
     initialize!(output, feedback, time_stepping, clock, diagn, model)
     initialize!(model.callbacks, progn, diagn, model)

src/physics/column_variables.jl

@@ -12,7 +12,8 @@ function get_column!(
         model.planet,
         model.orography,
         model.land_sea_mask,
-        model.albedo)
+        model.albedo,
+        model.implicit)
 end
 
 """
@@ -30,9 +31,11 @@ function get_column!(
     orography::AbstractOrography,
     land_sea_mask::AbstractLandSeaMask,
     albedo::AbstractAlbedo,
+    implicit::AbstractImplicit,
 )
 
     (; σ_levels_full, ln_σ_levels_full) = geometry
+    (; temp_profile) = implicit     # reference temperature on this layer
 
     @boundscheck C.nlev == D.nlev || throw(BoundsError)
 
@@ -57,6 +60,12 @@ function get_column!(
         C.temp[k] = layer.grid_variables.temp_grid[ij]
         C.temp_virt[k] = layer.grid_variables.temp_virt_grid[ij]    # actually diagnostic
         C.humid[k] = layer.grid_variables.humid_grid[ij] 
+
+        # and at previous time step, add temp reference profile back in as temp_grid_prev is anomaly
+        C.u_prev[k] = layer.grid_variables.u_grid_prev[ij]
+        C.v_prev[k] = layer.grid_variables.v_grid_prev[ij]
+        C.temp_prev[k] = layer.grid_variables.temp_grid_prev[ij] + temp_profile[k]
+        C.humid_prev[k] = layer.grid_variables.humid_grid_prev[ij] 
     end
 
     # TODO skin = surface approximation for now

src/physics/define_column.jl

@@ -21,10 +21,16 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV
     orography::NF = 0                       # orography height [m]
 
     # PROGNOSTIC VARIABLES
-    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)         # absolute temperature [K]
-    const humid::Vector{NF} = zeros(NF, nlev)        # specific humidity [kg/kg]
+    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)        # absolute temperature [K]
+    const humid::Vector{NF} = zeros(NF, nlev)       # specific humidity [kg/kg]
+
+    # PROGNOSTIC VARIABLES at previous time step
+    const u_prev::Vector{NF} = zeros(NF, nlev)      # zonal velocity [m/s]
+    const v_prev::Vector{NF} = zeros(NF, nlev)      # meridional velocity [m/s]
+    const temp_prev::Vector{NF} = zeros(NF, nlev)   # absolute temperature [K]
+    const humid_prev::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(Pa)]

src/physics/large_scale_condensation.jl

@@ -11,11 +11,14 @@ export ImplicitCondensation
 Large scale condensation as with implicit precipitation.
 $(TYPEDFIELDS)"""
 Base.@kwdef struct ImplicitCondensation{NF<:AbstractFloat} <: AbstractCondensation
-    "Flux limiter for latent heat release [K] per timestep"
-    max_heating::NF = 0.2
+    "Relative humidity threshold [1 = 100%] to trigger condensation"
+    relative_humidity_threshold::NF = 1
+
+    "Flux limiter for latent heat release [W/m²] per timestep"
+    max_heating::NF = 1
 
-    "Time scale in multiples of time step Δt"
-    time_scale::NF = 9
+    "Time scale in multiples of time step Δt, the larger the less immediate"
+    time_scale::NF = 3
 end
 
 ImplicitCondensation(SG::SpectralGrid; kwargs...) = ImplicitCondensation{SG.NF}(; kwargs...)
@@ -58,18 +61,20 @@ relative humidity. Calculates the tendencies for specific humidity
 and temperature from latent heat release and integrates the
 large-scale precipitation vertically for output."""
 function large_scale_condensation!(
-    column::ColumnVariables{NF},
+    column::ColumnVariables,
     scheme::ImplicitCondensation,
     clausius_clapeyron::AbstractClausiusClapeyron,
     geometry::Geometry,
     planet::AbstractPlanet,
     atmosphere::AbstractAtmosphere,
     time_stepping::AbstractTimeStepper,
-) where NF
+)
 
-    (; temp, humid, pres) = column           # prognostic vars: specific humidity, pressure
-    (; temp_tend, humid_tend) = column       # tendencies to write into
-    (; sat_humid) = column                   # intermediate variable, calculated in thermodynamics!
+    (; pres) = column                       # prognostic vars: pressure
+    temp = column.temp_prev                 # but use temp, humid from
+    humid = column.humid_prev               # from previous time step for numerical stability
+    (; temp_tend, humid_tend) = column      # tendencies to write into
+    (; sat_humid) = column                  # intermediate variable, calculated in thermodynamics!
 
     # precompute scaling constant for precipitation output
     pₛ = pres[end]                          # surface pressure
@@ -80,10 +85,10 @@ function large_scale_condensation!(
     (; Lᵥ, cₚ, Lᵥ_Rᵥ) = clausius_clapeyron
     Lᵥ_cₚ = Lᵥ/cₚ                           # latent heat of vaporization over heat capacity
     max_heating = scheme.max_heating/Δt_sec
-    time_scale = scheme.time_scale
+    (; time_scale, relative_humidity_threshold) = scheme
 
     @inbounds for k in eachindex(column)
-        if humid[k] > sat_humid[k]
+        if humid[k] > sat_humid[k]*relative_humidity_threshold
 
             # tendency for Implicit humid = sat_humid, divide by leapfrog time step below
             δq = sat_humid[k] - humid[k]

src/physics/tendencies.jl

@@ -36,6 +36,9 @@ function parameterization_tendencies!(
     end
 end
 
+"""$(TYPEDSIGNATURES)
+Calls for `column` one physics parameterization after another
+and convert fluxes to tendencies."""
 function parameterization_tendencies!(
     column::ColumnVariables,
     model::PrimitiveEquation

src/physics/thermodynamics.jl

@@ -148,7 +148,11 @@ function saturation_humidity!(
     column::ColumnVariables,
     clausius_clapeyron::AbstractClausiusClapeyron,
 )
-    (; sat_humid, pres, temp) = column
+    (; sat_humid, pres) = column
+
+    # use previous time step for temperature for stability of large-scale condensation
+    # TODO also use previous pressure, but sat_humid is only weakly dependent on it, skip for now
+    temp = column.temp_prev
 
     for k in eachlayer(column)
         sat_humid[k] = saturation_humidity(temp[k], pres[k], clausius_clapeyron)