AbstractSurfacePerturbation (#631)

SpeedyWeather · Dec 6, 2024 · aa279a0 · aa279a0
1 parent bdac940
commit aa279a0
Show file tree

Hide file tree

Showing 6 changed files with 164 additions and 19 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Unreleased
 
+- AbstractSurfacePerturbation introduced [#631](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/631)
 - AbstractStochasticPhysics and SPPT implemented [#629](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/629)
 - Introduced seperate extended tests that are not run on every commit [#628](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/628)
 - One-band longwave radiation [#624](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/624)

diff --git a/docs/src/stochastic_physics.md b/docs/src/stochastic_physics.md
@@ -3,9 +3,9 @@
 Stochastic physics introduces stochasticity into the parameterizations of physical process.
 There is conceptually several classes of how this can be done
 
-    - Stochastic perturbations of the tendencies
-    - Stochastic parameter perturbations
-    - Stochastic perturbations to the inputs of parameterizations
+- Stochastic perturbations of the tendencies
+- Stochastic parameter perturbations
+- Stochastic perturbations to the inputs of parameterizations
 
 all of these use random numbers created from some random processes
 (white noise or with some autocorrelation in space or time or sampled from some other distribution)

diff --git a/src/physics/convection.jl b/src/physics/convection.jl
@@ -1,19 +1,51 @@
+abstract type AbstractSurfacePerturbation <: AbstractParameterization end
+
+# subtypes don't require an initialize! method defined by default
+initialize!(::AbstractSurfacePerturbation, ::PrimitiveEquation) = nothing
+
+export NoSurfacePerturbation
+
+"""Returns the surface temperature and humidity without
+perturbation from the lowermost layer."""
+struct NoSurfacePerturbation <: AbstractSurfacePerturbation end
+
+"""$(TYPEDSIGNATURES)
+Returns the surface temperature and humidity without
+perturbation from the lowermost layer. Used as a functor."""
+function (SP::NoSurfacePerturbation)(
+    column::ColumnVariables,
+    model::PrimitiveEquation,
+)
+    (; temp, humid, nlayers) = column
+    return temp[nlayers], humid[nlayers]
+end
+
+# only return temperature for primitive dry model
+function (SP::NoSurfacePerturbation)(
+    column::ColumnVariables,
+    model::PrimitiveDry,
+)
+    (; temp, nlayers) = column
+    return temp[nlayers]
+end
+
 abstract type AbstractConvection <: AbstractParameterization end
 
 export NoConvection
+
+"""Dummy type to disable convection."""
 struct NoConvection <: AbstractConvection end
 NoConvection(::SpectralGrid) = NoConvection()
 initialize!(::NoConvection, ::PrimitiveEquation) = nothing
 convection!(::ColumnVariables, ::NoConvection, ::PrimitiveEquation) = nothing
 
 export SimplifiedBettsMiller
 
-"""
-The simplified Betts-Miller convection scheme from Frierson, 2007,
+"""The simplified Betts-Miller convection scheme from Frierson, 2007,
 https://doi.org/10.1175/JAS3935.1. This implements the qref-formulation
 in their paper. Fields and options are
 $(TYPEDFIELDS)"""
-@kwdef struct SimplifiedBettsMiller{NF} <: AbstractConvection
+@kwdef struct SimplifiedBettsMiller{NF, SP} <: AbstractConvection
     "number of vertical layers"
     nlayers::Int
 
@@ -22,10 +54,23 @@ $(TYPEDFIELDS)"""
 
     "[OPTION] Relative humidity for reference profile"
     relative_humidity::NF = 0.7
+
+    "[OPTION] Surface perturbation of temp, humid to calculate the moist pseudo adiabat"
+    surface_temp_humid::SP  # don't specify default here but in generator below
+end
+
+# generator function 
+function SimplifiedBettsMiller(
+    SG::SpectralGrid;
+    surface_temp_humid = NoSurfacePerturbation(),
+    kwargs...)
+    # type inference happens here as @kwdef only defines no/all types specified
+    return SimplifiedBettsMiller{SG.NF, typeof(surface_temp_humid)}(
+        nlayers=SG.nlayers; surface_temp_humid, kwargs...)
 end
 
-SimplifiedBettsMiller(SG::SpectralGrid; kwargs...) = SimplifiedBettsMiller{SG.NF}(nlayers=SG.nlayers; kwargs...)
-initialize!(::SimplifiedBettsMiller, ::PrimitiveEquation) = nothing
+# nothing to initialize but maybe the surface temp/humid functor?
+initialize!(SBM::SimplifiedBettsMiller, model::PrimitiveEquation) = initialize!(SBM.surface_temp_humid, model)
 
 # function barrier for all AbstractConvection
 function convection!(
@@ -42,7 +87,7 @@ function convection!(
     model::PrimitiveWet,
 )
     convection!(column, scheme, model.clausius_clapeyron,
-                    model.geometry, model.planet, model.atmosphere, model.time_stepping)
+                    model.geometry, model.planet, model.atmosphere, model.time_stepping, model)
 end
 
 """
@@ -60,6 +105,7 @@ function convection!(
     planet::AbstractPlanet,
     atmosphere::AbstractAtmosphere,
     time_stepping::AbstractTimeStepper,
+    model::PrimitiveWet,
 ) where NF
 
     σ = geometry.σ_levels_full
@@ -74,8 +120,8 @@ function convection!(
     humid_ref_profile = column.b    # specific humidity [kg/kg] profile to adjust to
 
     # CONVECTIVE CRITERIA AND FIRST GUESS RELAXATION
-    temp_parcel = temp[nlayers]
-    humid_parcel = humid[nlayers]
+    # force conversion to NF here
+    temp_parcel, humid_parcel = convert.(NF, SBM.surface_temp_humid(column, model))
     level_zero_buoyancy = pseudo_adiabat!(temp_ref_profile,
                                             temp_parcel, humid_parcel,
                                             temp_virt, geopot, pₛ, σ,
@@ -248,24 +294,37 @@ The simplified Betts-Miller convection scheme from Frierson, 2007,
 https://doi.org/10.1175/JAS3935.1 but with humidity set to zero.
 Fields and options are
 $(TYPEDFIELDS)"""
-@kwdef struct DryBettsMiller{NF} <: AbstractConvection
+@kwdef struct DryBettsMiller{NF, SP} <: AbstractConvection
     "number of vertical layers/levels"
     nlayers::Int
 
     "[OPTION] Relaxation time for profile adjustment"
     time_scale::Second = Hour(4)
+
+    "[OPTION] Surface perturbation of temp to calculate the dry adiabat"
+    surface_temp::SP    # don't set default here but in generator below
 end
 
-DryBettsMiller(SG::SpectralGrid; kwargs...) = DryBettsMiller{SG.NF}(nlayers=SG.nlayers; kwargs...)
-initialize!(::DryBettsMiller, ::PrimitiveEquation) = nothing
+# generator function
+function DryBettsMiller(
+    SG::SpectralGrid;
+    surface_temp = NoSurfacePerturbation(),
+    kwargs...)
+    # infer type here as @kwdef only defines constructors for no/all types specified
+    return DryBettsMiller{SG.NF, typeof(surface_temp)}(nlayers=SG.nlayers; surface_temp, kwargs...)
+end
+
+# nothing to initialize but maybe the surface temp functor?
+initialize!(DBM::DryBettsMiller, model::PrimitiveEquation) = initialize!(DBM.surface_temp, model)
 
 # function barrier to unpack model
 function convection!(
     column::ColumnVariables,
     scheme::DryBettsMiller,
     model::PrimitiveEquation,
 )
-    convection!(column, scheme, model.geometry, model.atmosphere)
+    # TODO check whether unpacking makes a performance difference? (it shouldn't?)
+    convection!(column, scheme, model.geometry, model.atmosphere, model)
 end
 
 """
@@ -281,6 +340,7 @@ function convection!(
     DBM::DryBettsMiller,
     geometry::Geometry,
     atmosphere::AbstractAtmosphere,
+    model::PrimitiveEquation,
 ) where NF
 
     σ = geometry.σ_levels_full
@@ -292,8 +352,12 @@ function convection!(
     temp_ref_profile = column.a     # temperature [K] reference profile to adjust to
 
     # CONVECTIVE CRITERIA AND FIRST GUESS RELAXATION
+    # force conversion to NF here, add "," to ignore additional returns like humidity
+    temp_parcel, = convert.(NF, DBM.surface_temp(column, model))
     level_zero_buoyancy = dry_adiabat!(temp_ref_profile,
-                                            temp, σ,
+                                            temp, 
+                                            temp_parcel,
+                                            σ,
                                             atmosphere)
 
     local PT::NF = 0        # precipitation due to coolinga
@@ -331,6 +395,7 @@ set to NaN instead and should be skipped in the relaxation."""
 function dry_adiabat!(
     temp_ref_profile::AbstractVector,
     temp_environment::AbstractVector,
+    temp_parcel::Real,
     σ::AbstractVector,
     atmosphere::AbstractAtmosphere,
 )
@@ -342,7 +407,6 @@ function dry_adiabat!(
         length(σ) == length(temp_environment) || throw(BoundsError)
 
     nlayers = length(temp_ref_profile)      # number of vertical levels
-    temp_parcel = temp_environment[nlayers] # parcel is at lowermost layer temperature
     temp_ref_profile .= NaN                 # reset profile from any previous calculation
     temp_ref_profile[nlayers] = temp_parcel    # start profile with parcel temperature
 

diff --git a/src/physics/longwave_radiation.jl b/src/physics/longwave_radiation.jl
@@ -16,7 +16,7 @@ end
 export UniformCooling
 
 """
-Uniform cooling following Pauluis and Garner, 2006. JAS. https://doi.org/10.1175/JAS3705.1
+Uniform cooling following Paulius and Garner, 2006. JAS. https://doi.org/10.1175/JAS3705.1
 imposing a default temperature tendency of -1.5K/day (=1K/16hours for a `time_scale` of 16 hours)
 on every level except for the stratosphere (diagnosed as `temp < temp_min`) where a
 relaxation term with `time_scale_stratosphere` towards `temp_stratosphere` is applied.
@@ -199,9 +199,12 @@ function longwave_radiation!(
         column.flux_temp_upward[k] += U         # accumulate that flux
     end
 
+    # store outgoing longwave radiation (OLR) for diagnostics
+    column.outgoing_longwave_radiation = U
+
     # DOWNWARD flux D
     local D::NF = 0                             # top boundary condition of longwave flux
-                                                # non need to accumulate 0 at top downward flux
+                                                # no need to accumulate 0 at top downward flux
     @inbounds for k in 1:nlayers
         # Frierson et al. 2006, equation 7
         D += dτ[k]*(B[k] - D)

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -48,6 +48,7 @@ include("optical_depth.jl")
 include("longwave_radiation.jl")
 # include("shortwave_radiation.jl")
 include("random_process.jl")
+include("stochastic_physics.jl")
 
 # INITIALIZATION AND INTEGRATION
 include("run_speedy.jl")

diff --git a/test/stochastic_physics.jl b/test/stochastic_physics.jl
@@ -0,0 +1,76 @@
+@testset "Surface perturbation wet/dry model" begin
+    @kwdef struct MySurfacePerturbation <: SpeedyWeather.AbstractSurfacePerturbation
+        noise_scale_temp::Float32 = 1            # standard deviation of perturbation, in Kelvin
+    end
+
+    # implement the actual perturbation as a functor 
+    function (SP::MySurfacePerturbation)(
+        column::ColumnVariables,
+        model::PrimitiveEquation,
+    )
+        (; temp, humid, nlayers) = column
+        T, q = temp[nlayers], humid[nlayers]    # unperturbed lowermost layer variables
+
+        # perturb temperature additive only
+        r = column.random_value
+        T += SP.noise_scale_temp*r
+
+        return T, q
+    end
+
+    spectral_grid = SpectralGrid()
+    my_surface_perturbation = MySurfacePerturbation(noise_scale_temp=1)
+
+    # construct convection with that perturbation
+    convection = SimplifiedBettsMiller(spectral_grid, surface_temp_humid=my_surface_perturbation)
+
+    # don't forget the random process otherwise your `column.random_value` is always 0!
+    random_process = SpectralAR1Process(spectral_grid, seed=1)
+
+    # construct model
+    model = PrimitiveWetModel(spectral_grid; random_process, convection)
+    simulation = initialize!(model)
+
+    run!(simulation, period=Day(5))
+
+    k = spectral_grid.nlayers
+    surface_humid = simulation.diagnostic_variables.grid.humid_grid[:, k]
+
+    # reconstruct model with different seed
+    random_process = SpectralAR1Process(spectral_grid, seed=2)
+    model = PrimitiveWetModel(spectral_grid; random_process, convection)
+    # reinitialize, rerun check whether it's different!
+    simulation = initialize!(model)     # choose a different random seed for random_process
+    run!(simulation, period=Day(5))
+    surface_humid2 = simulation.diagnostic_variables.grid.humid_grid[:, k]
+
+    for ij in eachindex(surface_humid, surface_humid2)
+        if surface_humid[ij] != 0
+            @test surface_humid[ij] != surface_humid2[ij]
+        end
+    end
+
+    # for the dry model too
+    convection = DryBettsMiller(spectral_grid, surface_temp=my_surface_perturbation)
+    model = PrimitiveDryModel(spectral_grid; random_process, convection)
+    simulation = initialize!(model)
+    run!(simulation, period=Day(5))
+
+    k = spectral_grid.nlayers
+    surface_temp = simulation.diagnostic_variables.grid.temp_grid[:, k]
+
+    # reinitialize, rerun check whether it's different!
+    simulation = initialize!(model)     # choose a different random seed for random_process
+    run!(simulation, period=Day(5))
+    surface_temp2 = simulation.diagnostic_variables.grid.temp_grid[:, k]
+
+    # count the number of grid points that are different in temperature
+    N = length(surface_temp)
+    local n::Int64 = 0
+    for ij in eachindex(surface_temp, surface_temp2)
+        n += surface_temp[ij] != surface_temp2[ij]
+    end
+
+    # 75% of grid points must be different
+    @test n/N > 3/4
+end