diff --git a/src/SpeedyWeather.jl b/src/SpeedyWeather.jl
index 178809ba1..a43de646b 100644
--- a/src/SpeedyWeather.jl
+++ b/src/SpeedyWeather.jl
@@ -196,8 +196,8 @@ include("physics/large_scale_condensation.jl")
 include("physics/surface_fluxes.jl")
 include("physics/convection.jl")
 include("physics/zenith.jl")
-# include("physics/shortwave_radiation.jl")
-# include("physics/longwave_radiation.jl")
+include("physics/shortwave_radiation.jl")
+include("physics/longwave_radiation.jl")
 include("physics/pretty_printing.jl")
 
 # OCEAN AND LAND
diff --git a/src/dynamics/initial_conditions.jl b/src/dynamics/initial_conditions.jl
index 7b1e06a9c..47ce7cd53 100644
--- a/src/dynamics/initial_conditions.jl
+++ b/src/dynamics/initial_conditions.jl
@@ -455,7 +455,7 @@ function initialize_humidity!(  progn::PrognosticVariables,
     temp_grid = gridded(progn.layers[end].timesteps[1].temp,model.spectral_transform)
     humid_surf_grid = zero(pres_surf_grid)
     for ij in eachgridpoint(humid_surf_grid)
-        q_sat = saturation_humidity(temp_grid[ij],exp(pres_surf_grid[ij]),model.clausis_clapeyron)
+        q_sat = saturation_humidity(temp_grid[ij],exp(pres_surf_grid[ij]),model.clausius_clapeyron)
         humid_surf_grid[ij] = relhumid_ref*q_sat
     end
 
diff --git a/src/dynamics/models.jl b/src/dynamics/models.jl
index 42e904240..54546b7f8 100644
--- a/src/dynamics/models.jl
+++ b/src/dynamics/models.jl
@@ -163,13 +163,15 @@ Base.@kwdef mutable struct PrimitiveDryModel{NF<:AbstractFloat, D<:AbstractDevic
 
     # PHYSICS/PARAMETERIZATIONS
     physics::Bool = true
-    boundary_layer_drag::BoundaryLayerDrag{NF} = NoBoundaryLayerDrag(spectral_grid)
+    boundary_layer_drag::BoundaryLayerDrag{NF} = BulkRichardsonDrag(spectral_grid)
     temperature_relaxation::TemperatureRelaxation{NF} = NoTemperatureRelaxation(spectral_grid)
     static_energy_diffusion::VerticalDiffusion{NF} = NoVerticalDiffusion(spectral_grid)
     surface_thermodynamics::AbstractSurfaceThermodynamics{NF} = SurfaceThermodynamicsConstant(spectral_grid)
     surface_wind::AbstractSurfaceWind{NF} = SurfaceWind(spectral_grid)
     surface_heat_flux::AbstractSurfaceHeat{NF} = SurfaceSensibleHeat(spectral_grid)
-    
+    shortwave_radiation::AbstractShortwave{NF} = NoShortwave(spectral_grid)
+    longwave_radiation::AbstractLongwave{NF} = UniformCooling(spectral_grid)
+
     # NUMERICS
     time_stepping::TimeStepper{NF} = Leapfrog(spectral_grid)
     spectral_transform::SpectralTransform{NF} = SpectralTransform(spectral_grid)
@@ -215,6 +217,8 @@ function initialize!(model::PrimitiveDry;time::DateTime = DEFAULT_DATE)
     initialize!(model.boundary_layer_drag,model)
     initialize!(model.temperature_relaxation,model)
     initialize!(model.static_energy_diffusion,model)
+    initialize!(model.shortwave_radiation,model)
+    initialize!(model.longwave_radiation,model)
 
     # initial conditions
     prognostic_variables = PrognosticVariables(spectral_grid,model)
@@ -255,8 +259,8 @@ Base.@kwdef mutable struct PrimitiveWetModel{NF<:AbstractFloat, D<:AbstractDevic
 
     # PHYSICS/PARAMETERIZATIONS
     physics::Bool = true
-    clausis_clapeyron::AbstractClausiusClapeyron{NF} = ClausiusClapeyron(spectral_grid,atmosphere)
-    boundary_layer_drag::BoundaryLayerDrag{NF} = NoBoundaryLayerDrag(spectral_grid)
+    clausius_clapeyron::AbstractClausiusClapeyron{NF} = ClausiusClapeyron(spectral_grid,atmosphere)
+    boundary_layer_drag::BoundaryLayerDrag{NF} = BulkRichardsonDrag(spectral_grid)
     temperature_relaxation::TemperatureRelaxation{NF} = NoTemperatureRelaxation(spectral_grid)
     static_energy_diffusion::VerticalDiffusion{NF} = NoVerticalDiffusion(spectral_grid)
     humidity_diffusion::VerticalDiffusion{NF} = NoVerticalDiffusion(spectral_grid)
@@ -266,7 +270,9 @@ Base.@kwdef mutable struct PrimitiveWetModel{NF<:AbstractFloat, D<:AbstractDevic
     surface_heat_flux::AbstractSurfaceHeat{NF} = SurfaceSensibleHeat(spectral_grid)
     evaporation::AbstractEvaporation{NF} = SurfaceEvaporation(spectral_grid)
     convection::AbstractConvection{NF} = SimplifiedBettsMiller(spectral_grid)
-    
+    shortwave_radiation::AbstractShortwave{NF} = NoShortwave(spectral_grid)
+    longwave_radiation::AbstractLongwave{NF} = UniformCooling(spectral_grid)
+
     # NUMERICS
     time_stepping::TimeStepper{NF} = Leapfrog(spectral_grid)
     spectral_transform::SpectralTransform{NF} = SpectralTransform(spectral_grid)
@@ -318,6 +324,8 @@ function initialize!(model::PrimitiveWet;time::DateTime = DEFAULT_DATE)
     initialize!(model.humidity_diffusion,model)
     initialize!(model.large_scale_condensation,model)
     initialize!(model.convection,model)
+    initialize!(model.shortwave_radiation,model)
+    initialize!(model.longwave_radiation,model)
 
     # initial conditions
     prognostic_variables = PrognosticVariables(spectral_grid,model)
diff --git a/src/physics/column_variables.jl b/src/physics/column_variables.jl
index 8b25d49b1..ded89b696 100644
--- a/src/physics/column_variables.jl
+++ b/src/physics/column_variables.jl
@@ -146,15 +146,11 @@ function reset_column!(column::ColumnVariables{NF}) where NF
     column.flux_temp_upward .= 0
     column.flux_temp_downward .= 0
 
-    # Convection
+    # Convection and precipitation
     column.cloud_top = column.nlev+1            # also diagnostic from condensation
-    column.conditional_instability = false
-    column.activate_convection = false
-    column.excess_humid = 0
     column.precip_convection = 0
-
-    # Large-scale condensation
     column.precip_large_scale = 0
+
     return nothing
 end
 
diff --git a/src/physics/convection.jl b/src/physics/convection.jl
index 3f106f57b..5045c1f42 100644
--- a/src/physics/convection.jl
+++ b/src/physics/convection.jl
@@ -1,328 +1,5 @@
 abstract type AbstractConvection{NF} <: AbstractParameterization{NF} end
 
-Base.@kwdef struct SpeedyConvection{NF} <: AbstractConvection{NF}
-
-    "Number of vertical levels"
-    nlev::Int
-
-    "Minimum (normalised) surface pressure for the occurrence of convection [1]"
-    pres_threshold::NF = 0.8
-
-    "Relative humidity threshold for convection in PBL [1]"
-    humid_threshold_boundary::NF = 0.9
-
-    "Relative humidity threshold for convection in the troposphere [1]"
-    humid_threshold_troposphere::NF = 0.7
-
-    "Relaxation time for PBL humidity"
-    time_scale::Second = Hour(6)
-
-    "Maximum entrainment as a fraction of cloud-base mass flux"
-    max_entrainment::NF = 0.5
-
-    "Ratio between secondary and primary mass flux at cloud-base"
-    ratio_secondary_mass_flux::NF = 0.8
-
-    "Super saturation of relative humidity in the planteray boundar layer [1]"
-    super_saturation::NF = 1.01
-
-    "Mass flux limiter for humidity [kg/kg]"
-    mass_flux_limiter_humid::NF = 5e-3
-
-    # precomputed in initialize!
-    "Reference surface pressure [Pa]"
-    pres_ref::Base.RefValue{NF} = Base.Ref(zero(NF))
-
-    "latent heat of condensation [J/kg] for consistency with specific humidity [kg/kg], also called alhc"
-    latent_heat_condensation::Base.RefValue{NF} = Base.Ref(zero(NF))
-
-    "specific heat [J/kg/K] of air"
-    specific_heat::Base.RefValue{NF} = Base.Ref(zero(NF))
-
-    "Number of vertical levels for stratosphere"
-    n_stratosphere_levels::Base.RefValue{Int} = Base.Ref(0)
-
-    "Number of vertical levels for planetary boundary layer"
-    n_boundary_levels::Base.RefValue{Int} = Base.Ref(0)
-
-    "Mass flux entrainment profile in the vertical [1]"
-    entrainment_profile::Vector{NF} = zeros(NF,nlev)
-end
-
-SpeedyConvection(SG::SpectralGrid;kwargs...) = SpeedyConvection{SG.NF}(nlev=SG.nlev;kwargs...)
-
-function initialize!(convection::SpeedyConvection,model::PrimitiveWet)
-
-    (;σ_levels_full) = model.geometry
-    (;σ_tropopause, σ_boundary_layer) = model.atmosphere
-    (;entrainment_profile, nlev) = convection
-
-    # number of stratospheric levels, for nlev very small n can be nothing, then use 0
-    n = findlast(σ->σ<=σ_tropopause,σ_levels_full)
-    convection.n_stratosphere_levels[] = isnothing(n) ? 0 : n
-
-    # number of levels for the planetary boundary layer, same as above
-    n = findlast(σ->σ<=σ_boundary_layer,σ_levels_full)
-    convection.n_boundary_levels[] = isnothing(n) ? 0 : nlev - n
-    
-    # reference pressure
-    convection.pres_ref[] = model.atmosphere.pres_ref*100     # [hPa] -> [Pa]
-    convection.latent_heat_condensation[] = model.atmosphere.latent_heat_condensation
-    convection.specific_heat[] = model.atmosphere.cₚ
-
-    # Mass entrainment profile
-    entrainment_profile[1] = 0      # no entrainment in top layer
-    entrainment_profile[nlev] = 0   # no entrainment in bottom layer
-    for k = 2:nlev-1                # intermediate layers with minimum at σ=0
-        entrainment_profile[k] = max(0, (σ_levels_full[k] - 0.5))
-    end
-
-    # profile as fraction of cloud-base mass flux, normalise to max entrainment at nlev-1
-    entrainment_profile .*= convection.max_entrainment/entrainment_profile[nlev-1]
-end
-
-# dry model doesn't have convection TODO maybe it should?
-convection!(column::ColumnVariables,model::PrimitiveDry) = nothing
-
-# function barrier for all AbstractConvection
-function convection!(
-    column::ColumnVariables,
-    model::PrimitiveEquation,
-)
-    convection!(column,model.convection,model)
-end
-
-function convection!(
-    column::ColumnVariables,
-    scheme::SpeedyConvection,
-    model::PrimitiveEquation,
-)
-
-    moist_static_energy!(column, model.clausis_clapeyron)
-    vertical_interpolate!(column, model)                # Interpolate certain variables to half-levels
-
-    # always diagnose convection
-    diagnose_convection!(column, model.convection)
-
-    # but only execute if conditions are met
-    if column.conditional_instability && column.activate_convection
-        convection!(column,model.convection,model.constants,model.geometry,model.time_stepping)
-    end
-end
-
-"""
-$(TYPEDSIGNATURES)
-Check whether the convection scheme should be activated in the given atmospheric column.
-
-1. A conditional instability exists when the saturation moist energy (MSS) decreases with
-height, that is, there exists an atmospheric level k such that,
-
-    MSS(N) > MSS(k+h)
-
-where N is the planetary boundary layer (PBL) and k+h is the half-level at the lower
-boundary of the full level k.
-
-2. When a conditional instability exists, the convection scheme is activated when, either,
-
-    a. the actual moist static energy (MSE) at level N-1 (directly above the PBL) is greater
-       than the saturation moist static energy at some half-level k+h,
-
-            MSE(N-1) > MSS(k+h)
-
-    b. the humidity in both the PBL and one layer above exceeds a prescribed threshold,
-
-            Q(N)   > RH_cnv * Qˢᵃᵗ(N)
-            Q(N-1) > RH_cnv * Qˢᵃᵗ(N-1)
-
-The top-of-convection (TCN) layer, or cloud-top, is the largest value of k for which
-condition 1 is satisfied. The cloud-top layer may be subsequently adjusted upwards by the
-large-scale condensation parameterization, which is executed after this one."""
-function diagnose_convection!(column::ColumnVariables,convection::SpeedyConvection)
-
-    (; pres_ref, pres_threshold, humid_threshold_boundary) = convection
-    n_stratosphere_levels = convection.n_stratosphere_levels[]
-    n_boundary_levels = convection.n_boundary_levels[]
-    latent_heat = convection.latent_heat_condensation[]
-
-    (; nlev ) = column
-    (; humid, pres, sat_humid, moist_static_energy,
-    sat_moist_static_energy, sat_moist_static_energy_half) = column
-
-    # effectively disables convection over Himalaya/Tibet, Greenland and Antarctica
-    if pres[end] > pres_threshold*pres_ref[]
-
-        # First we pre-compute some values which we will need inside the loop
-        # 1. Saturation (or super-saturated) moist static energy in the PBL
-        sat_moist_static_energy_pbl =
-            max(moist_static_energy[nlev], sat_moist_static_energy[nlev])
-
-        # 2. Minimum of moist static energy in the lowest two levels
-        moist_static_energy_lower_trop =
-            min(moist_static_energy[nlev], moist_static_energy[nlev-1])
-
-        # 3. Humidity threshold for convection, defined in the PBL and one level above
-        humid_threshold_pbl = humid_threshold_boundary * sat_humid[nlev]
-        humid_threshold_above_pbl = humid_threshold_boundary * sat_humid[nlev-1]
-
-        # The range of this loop requires clarification, but in its current form it means
-        # that the top-of-convection level may be any tropospheric level, excluding the two
-        # layers directly above the PBL.
-        for k = (nlev-(n_boundary_levels+1)):-1:(n_stratosphere_levels+1)
-            # Condition 1: Conditional instability (MSS in PBL < MSS at this half-level)
-            if sat_moist_static_energy_pbl > sat_moist_static_energy_half[k]
-                column.conditional_instability = true
-                column.cloud_top = k
-            end
-
-            # Condition 2a: Gradient of actual moist static energy between lower and upper troposphere
-            if moist_static_energy_lower_trop > sat_moist_static_energy_half[k]
-                column.activate_convection = true
-                column.excess_humid = max(
-                    humid[nlev] - humid_threshold_pbl,
-                    (moist_static_energy[nlev] - sat_moist_static_energy_half[k]) / latent_heat,
-                )
-            end
-        end
-
-        if column.conditional_instability && column.activate_convection
-            return nothing  # Condition for convection already satisfied
-        end
-
-        # Condition 2b: Humidity exceeds threshold in both PBL and one layer above
-        if column.conditional_instability &&
-           (humid[nlev] > humid_threshold_pbl) &&
-           (humid[nlev-1] > humid_threshold_above_pbl)
-            column.activate_convection = true
-            column.excess_humid = humid[nlev] - humid_threshold_pbl
-        end
-    end
-    return nothing
-end
-
-"""
-$(TYPEDSIGNATURES)
-Compute fluxes and precipitation due to convection in the given atmospheric column.
-
-The scheme computes fluxes of mass, humidity and dry static energy. A part of the upward
-moisture flux at the lower boundary of the cloud-top (TCN) layer is converted into
-convective precipitation.
-
-For full details of the scheme see: http://users.ictp.it/~kucharsk/speedy_description/km_ver41_appendixA.pdf
-"""
-function convection!(
-    column::ColumnVariables{NF},
-    convection::SpeedyConvection,
-    C::DynamicsConstants,
-    G::Geometry,
-    T::TimeStepper,
-) where NF
-
-    (; gravity, pres_ref ) = C
-    (; nlev, pres, humid, humid_half, sat_humid, sat_humid_half, cloud_top) = column
-    (; dry_static_energy, dry_static_energy_half, cloud_top) = column
-    (; flux_temp_downward, flux_temp_upward, flux_humid_downward, flux_humid_upward) = column
-    (; temp_tend, humid_tend) = column
-    (; σ_levels_thick) = G
-
-    n_boundary_levels = convection.n_boundary_levels[]
-    time_scale = convection.time_scale.value
-    latent_heat_cₚ = convection.latent_heat_condensation[]/convection.specific_heat[]
-    (;pres_threshold, humid_threshold_troposphere, ratio_secondary_mass_flux) = convection
-    (;mass_flux_limiter_humid) = convection
-
-    # 1. Fluxes in the planetary boundary layer (PBL)
-    # Humidity at the upper boundary of the PBL
-    humid_top_of_pbl = min(humid_half[nlev-n_boundary_levels], humid[nlev-n_boundary_levels+1])
-    
-    # Maximum specific humidity in the PBL
-    max_humid_pbl = zero(NF)
-    for k in nlev-n_boundary_levels+1:nlev
-        max_humid_pbl = max(convection.super_saturation * humid[k], sat_humid[k])
-    end
-
-    # Cloud-base mass flux
-    pₛ = pres[end]                          # surface pressure
-    p_norm = pₛ/pres_ref                    # normalised surface pressure
-    Δp = pres_ref*σ_levels_thick[nlev]      # Pressure difference between bottom and top of PBL
-    
-    # excess humidity relative to humid difference across PBL
-    excess_humid = column.excess_humid / (max_humid_pbl - humid_top_of_pbl)
-
-    # Fortran SPEEDY documentation eq. (12) and (13)
-    mass_flux =
-        Δp / (gravity * time_scale) *
-        # Fortran SPEEDY extends this with some flux limiters:
-        # (original SPEEDY formulation doesn't make sense to me given it's supposed to be for
-        # "Minimum (normalised) surface pressure for the occurrence of convection")
-        # p_norm*min(1,  2 * (p_norm - pres_threshold) / (1 - pres_threshold)) *
-        max(0,(p_norm - pres_threshold) / (1 - pres_threshold)) *       # new version
-        min(mass_flux_limiter_humid, excess_humid)
-    
-    column.cloud_base_mass_flux = mass_flux
-
-    # Upward fluxes at upper boundary of PBL
-    flux_up_humid = mass_flux * max_humid_pbl                                       # eq. (10)
-    flux_up_static_energy = mass_flux * dry_static_energy[nlev-n_boundary_levels+1] # eq. (11)
-
-    # Downward fluxes at upper boundary of PBL
-    flux_down_humid = mass_flux * humid_top_of_pbl                                          # eq. (10) 
-    flux_down_static_energy = mass_flux * dry_static_energy_half[nlev-n_boundary_levels]    # eq. (11)
-
-    # Accumulate in fluxes for all parameterizations
-    flux_temp_upward[nlev-n_boundary_levels+1]   += flux_up_static_energy
-    flux_temp_downward[nlev-n_boundary_levels+1] += flux_down_static_energy
-    
-    flux_humid_upward[nlev-n_boundary_levels+1]   += flux_up_humid
-    flux_humid_downward[nlev-n_boundary_levels+1] += flux_down_humid
-
-    # 2. Fluxes for intermediate layers
-    for k = (nlev-n_boundary_levels):-1:(cloud_top+1)
-
-        # Mass entrainment
-        mass_entrainment = p_norm * convection.entrainment_profile[k] * column.cloud_base_mass_flux
-        mass_flux += mass_entrainment
-
-        # Upward fluxes at upper boundary
-        flux_up_static_energy += mass_entrainment * dry_static_energy[k]
-        flux_up_humid         += mass_entrainment * humid[k]
-
-        # Downward fluxes at upper boundary, _half vectors skip top (k=1/2)
-        flux_down_static_energy = mass_flux * dry_static_energy_half[k-1]
-        flux_down_humid         = mass_flux * humid_half[k-1]
-
-        # accumulate in fluxes that are translated to tendencies in tendencies.jl
-        flux_temp_upward[k]   += flux_up_static_energy
-        flux_temp_downward[k] += flux_down_static_energy
-        
-        flux_humid_upward[k]   += flux_up_humid
-        flux_humid_downward[k] += flux_down_humid
-
-        # Secondary moisture flux representing shallower, non-precipitating convective systems
-        # Occurs when RH in an intermediate layer falls below a threshold
-        Δhumid = humid_threshold_troposphere * sat_humid[k] - humid[k]
-        if Δhumid > 0
-            Δflux_humid = ratio_secondary_mass_flux * mass_flux * Δhumid
-
-            # a flux from bottom nlev to layer k, equivalent to net flux of Δflux_humid into k but out of nlev
-            for kk in k+1:nlev+1
-                flux_humid_upward[kk] += Δflux_humid
-            end
-        end
-    end
-
-    # 3. Convective precipitation in top-of-convection layer
-    precip_convection = max(flux_up_humid - mass_flux * sat_humid_half[cloud_top], 0)   # in [kg/m²/s]
-    column.precip_convection = precip_convection*T.Δt_sec/C.water_density               # convert to [m]
-
-    # Condensation of convectice precipiation creates a humidity and temperature (heating)
-    # tendency at the top of convection layer cloud_top
-    humid_tend_k = -precip_convection*gravity/(pₛ*σ_levels_thick[cloud_top])
-    humid_tend[cloud_top] += humid_tend_k
-    temp_tend[cloud_top] -= humid_tend_k*latent_heat_cₚ
-
-    return nothing
-end
-
 Base.@kwdef struct SimplifiedBettsMiller{NF} <: AbstractConvection{NF}
     "number of vertical layers/levels"
     nlev::Int
@@ -331,7 +8,7 @@ Base.@kwdef struct SimplifiedBettsMiller{NF} <: AbstractConvection{NF}
     time_scale::Second = Hour(4)
 
     "Relative humidity for reference profile"
-    relative_humidity::NF = 0.8
+    relative_humidity::NF = 0.7
 
     "temperature [K] reference profile to adjust to"
     temp_ref_profile::Vector{NF} = zeros(NF,nlev)
@@ -341,16 +18,26 @@ Base.@kwdef struct SimplifiedBettsMiller{NF} <: AbstractConvection{NF}
 end
 
 SimplifiedBettsMiller(SG::SpectralGrid;kwargs...) = SimplifiedBettsMiller{SG.NF}(nlev=SG.nlev;kwargs...)
-
 initialize!(::SimplifiedBettsMiller,::PrimitiveWet) = nothing
 
+# function barrier for all AbstractConvection
+function convection!(
+    column::ColumnVariables,
+    model::PrimitiveEquation,
+)
+    convection!(column,model.convection,model)
+end
+
+# TODO SimplifiedBettsMiller can be trimmed for PrimitiveDry
+convection!(column::ColumnVariables,model::PrimitiveDry) = nothing
+
 # function barrier to unpack model
 function convection!(
     column::ColumnVariables,
     scheme::SimplifiedBettsMiller,
     model::PrimitiveEquation,
 )
-    convection!(column, scheme, model.clausis_clapeyron, model.geometry, model.constants, model.time_stepping)
+    convection!(column, scheme, model.clausius_clapeyron, model.geometry, model.constants, model.time_stepping)
 end
 
 function convection!(
diff --git a/src/physics/define_column.jl b/src/physics/define_column.jl
index 01c97e68a..7a6c9ad66 100644
--- a/src/physics/define_column.jl
+++ b/src/physics/define_column.jl
@@ -64,84 +64,13 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV
     # THERMODYNAMICS
     surface_air_density::NF = 0
     const sat_humid::Vector{NF} = zeros(NF,nlev)                # Saturation specific humidity [kg/kg]
-    const rel_humid::Vector{NF} = zeros(NF,nlev)                # Relative humidity [1]
-    const dry_static_energy::Vector{NF} = zeros(NF,nlev)        # Dry static energy
-    const moist_static_energy::Vector{NF} = zeros(NF,nlev)      # Moist static energy
-    const sat_moist_static_energy::Vector{NF} = zeros(NF,nlev)  # Saturation moist static energy
+    const dry_static_energy::Vector{NF} = zeros(NF,nlev)        # Dry static energy [J/kg]
     const temp_virt::Vector{NF} = zeros(NF,nlev)                # virtual temperature [K]
     const bulk_richardson::Vector{NF} = zeros(NF,nlev)          # bulk richardson number [1]
     const geopot::Vector{NF} = zeros(NF,nlev)                   # gepotential height [m]
 
-    # and interpolated to half levels
-    const humid_half::Vector{NF} = zeros(NF,nlev)                    # Specific humidity interpolated to half-levels
-    const sat_humid_half::Vector{NF} = zeros(NF,nlev)                # Saturation specific humidity interpolated to half-levels
-    const dry_static_energy_half::Vector{NF} = zeros(NF,nlev)        # Dry static energy interpolated to half-levels
-    const sat_moist_static_energy_half::Vector{NF} = zeros(NF,nlev)  # Saturation moist static energy interpolated to half-levels
-
-    # 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_humid::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]
-
-    # Large-scale condensation
-    precip_large_scale::NF = 0                                  # precipitation [m]
-
-    # Longwave radiation
-    ## New vars in radlw_down!
-    wvi::Matrix{NF} = fill(NF(NaN), nlev, 2)  # Weights for vertical interpolation
-    tau2::Matrix{NF} = fill(NF(NaN), nlev, nband) # Transmissivity of atmospheric layers
-    dfabs::Vector{NF} = fill(NF(NaN), nlev)   # Flux of sw rad. absorbed by each atm. layer
-    fsfcd::NF = NaN                       # Downward-only flux of sw rad. at the surface
-    st4a::Matrix{NF} = fill(NF(NaN), nlev, 2) # Blackbody emission from full and half atmospheric levels
-    flux::Vector{NF} = fill(NF(NaN), nband)       # Radiative flux in different spectral bands
-
-    ## New vars in compute_bbe!
-    fsfcu::NF = NaN # surface blackbody emission (upward)
-    ts::NF = NaN    # surface temperature
-
-    ## New vars in radlw_up!
-    fsfc::NF = NaN # Net (downw.) flux of sw rad. at the surface
-    ftop::NF = NaN # Net (downw.) flux of sw rad. at the atm. top
-    stratc::Vector{NF} = fill(NF(NaN), n_stratosphere_levels) # Stratospheric correction term 
-    # Shortwave radiation: solar
-    tyear::NF = NF(NaN) # time as fraction of year (0-1, 0 = 1jan.h00)
-    csol::NF = NF(NaN)  # FIXME
-    topsr::NF = NF(NaN) # FIXME
-    # Shortwave radiation: solar_oz
-    fsol::NF = NF(NaN)   # Flux of incoming solar radiation
-    ozupp::NF = NF(NaN)  # Flux absorbed by ozone (upper stratos.)
-    ozone::NF = NF(NaN)  # Flux absorbed by ozone (lower stratos.)
-    zenit::NF = NF(NaN)  # Optical depth ratio (function of solar zenith angle)
-    stratz::NF = NF(NaN) # Stratospheric correction for polar night
-    # Shortwave radiation: radsw
-    albsfc::NF = NF(NaN) # Combined surface albedo (land + sea)
-    ssrd::NF = NF(NaN)   # Surface shortwave radiation (downward-only)
-    ssr::NF = NF(NaN)    # Surface shortwave radiation (net downward)
-    tsr::NF = NF(NaN)    # Top-of-atm. shortwave radiation (downward)
-    tend_t_rsw::Vector{NF} = fill(NF(NaN), nlev) # Tempterature tendency
-    norm_pres::NF = NF(NaN) # Normalized pressure (p/1000 hPa)
-    # Shortwave radiation: cloud
-    icltop::Int = typemax(Int) # Cloud top level (all clouds)
-    cloudc::NF = NF(NaN)       # Total cloud cover (fraction)
-    clstr::NF = NF(NaN)        # Stratiform cloud cover (fraction)
-    qcloud::NF = NF(NaN)       # Equivalent specific humidity of clouds
-    fmask::NF = NF(NaN)        # Fraction of land
-    # Shortwave radiation: shortwave_radiation
-    # rel_hum::Vector{NF} = fill(NF(NaN), nlev) # Relative humidity
-    grad_dry_static_energy::NF = NF(NaN)      # gradient of dry static energy
-end
-
-function plot(column::ColumnVariables,var::Symbol=:temp)
-    v = getfield(column,var)
-    z = 1:length(v)
-
-    x,y = column.lond, column.latd
-    title = "Column at $(y)˚N, $(x)˚E"
-    ylabel = "k"
-    yflip = true
-
-    UnicodePlots.lineplot(v,z;title,xlabel=string(var),ylabel,yflip)
+    # CONVECTION AND PRECIPITATION
+    cloud_top::Int = nlev + 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]
 end
\ No newline at end of file
diff --git a/src/physics/large_scale_condensation.jl b/src/physics/large_scale_condensation.jl
index f2713aee5..59d50ae15 100644
--- a/src/physics/large_scale_condensation.jl
+++ b/src/physics/large_scale_condensation.jl
@@ -81,7 +81,7 @@ function large_scale_condensation!(
     column::ColumnVariables,
     model::PrimitiveWet,
 )
-    saturation_humidity!(column, model.clausis_clapeyron)
+    saturation_humidity!(column, model.clausius_clapeyron)
     large_scale_condensation!(column,model.large_scale_condensation,model)
 end
 
@@ -178,7 +178,7 @@ function large_scale_condensation!(
     model::PrimitiveWet,
 )
     large_scale_condensation!(column,scheme,
-        model.clausis_clapeyron,model.geometry,model.constants,model.time_stepping)
+        model.clausius_clapeyron,model.geometry,model.constants,model.time_stepping)
 end
 
 """
diff --git a/src/physics/longwave_radiation.jl b/src/physics/longwave_radiation.jl
index 666e32274..d3bb9b981 100644
--- a/src/physics/longwave_radiation.jl
+++ b/src/physics/longwave_radiation.jl
@@ -1,241 +1,45 @@
-"""
-    function initialize_longwave_radiation!(
-        P::Parameters
-    )
+abstract type AbstractLongwave{NF} <: AbstractRadiation{NF} end
 
-Initialise variables and parameters used by the longwave radiation parametrization 
-"""
-function initialize_longwave_radiation!(K::ParameterizationConstants,
-                                        P::Parameters) 
-    radset!(K,P)
-end
-
-"""
-    function longwave_radiation!(
-        column::ColumnVariables{NF}, model::PrimitiveEquation
-    )
-
-Compute air temperature tendencies from longwave radiation for an atmospheric column.
-For more details see [http://users.ictp.it/~kucharsk/speedy_description/km_ver41_appendixA.pdf](http://users.ictp.it/~kucharsk/speedy_description/km_ver41_appendixA.pdf)
-"""
-function longwave_radiation!(   column::ColumnVariables,
-                                model::PrimitiveEquation)
+struct NoLongwave{NF} <: AbstractLongwave{NF} end
+NoLongwave(SG::SpectralGrid) = NoLongwave{SG.NF}()
+initialize!(::NoLongwave,::PrimitiveEquation) = nothing
 
-    radlw_down!(column, model)
-    compute_bbe!(column, model)
-    radlw_up!(column, model)
+# function barrier for all AbstractLongwave
+function longwave_radiation!(column::ColumnVariables,model::PrimitiveEquation)
+    longwave_radiation!(column,model.longwave_radiation,model)
 end
 
-"""
-    function radset!(model::PrimitiveEquation) where {NF<:AbstractFloat}
-        
-Compute energy fractions in four longwave bands as a function of temperature.
-
-"""
-function radset!(   K::ParameterizationConstants,
-                    P::Parameters)
-
-    (; NF, nband ) = P
-    (; epslw ) = P.radiation_coefs
-    (; fband ) = K
-
-    @assert nband == 4 "Only four bands are supported, given $nband"
+longwave_radiation!(::ColumnVariables,::NoLongwave,::PrimitiveEquation) = nothing
 
-    eps1 = 1 - epslw
-    for jtemp = 200:320
-        fband[jtemp, 2] = (NF(0.148) - NF(3.0e-6) * (jtemp - 247)^2) * eps1
-        fband[jtemp, 3] = (NF(0.356) - NF(5.2e-6) * (jtemp - 282)^2) * eps1
-        fband[jtemp, 4] = (NF(0.314) + NF(1.0e-5) * (jtemp - 315)^2) * eps1
-        fband[jtemp, 1] = eps1 - (fband[jtemp, 2] + fband[jtemp, 3] + fband[jtemp, 4])
-    end
-    for jb = 1:4
-        for jtemp = 100:199
-            fband[jtemp, jb] = fband[200, jb]
-        end
-        for jtemp = 321:400
-            fband[jtemp, jb] = fband[320, jb]
-        end
-    end
+Base.@kwdef struct UniformCooling{NF} <: AbstractLongwave{NF}
+    time_scale::Second = Hour(16)
+    temp_min::NF = 207.5
+    temp_stratosphere::NF = 200
+    time_scale_stratosphere::Second = Day(5)
 end
 
+UniformCooling(SG::SpectralGrid;kwargs...) = UniformCooling{SG.NF}(;kwargs...)
+initialize!(scheme::UniformCooling,model::PrimitiveEquation) = nothing
 
-"""
-    function radlw_down!(
-        column::ColumnVariables{NF}, model::PrimitiveEquation
-    ) where {NF<:AbstractFloat}
-
-Compute the downward flux of longwave radiation. Inputs variables are `temp``, `wvi`, `tau2`.
-Output column varables are `fsfcd`, `dfabs`, `flux`, `st4a`.
-"""
-function radlw_down!(
-    column::ColumnVariables{NF}, model::PrimitiveEquation
-) where {NF<:AbstractFloat}
-
-    (; nlev, temp, wvi, tau2 ) = column
-    (; nband, sbc ) = model.parameters
-    (; n_stratosphere_levels ) = model.geometry
-    (; fband ) = model.parameterization_constants
-    (; epslw, emisfc ) = model.parameters.radiation_coefs
-    
-    # 1. Blackbody emission from atmospheric levels.
-    #    The linearized gradient of the blakbody emission is computed
-    #    from temperatures at layer boundaries, which are interpolated 
-    #    assuming a linear dependence of T on log_sigma.
-    #    Above the first (top) level, the atmosphere is assumed isothermal.
-    # 
-    # Temperature at level boundaries
-    for k=1:nlev-1
-        column.st4a[k, 1] = temp[k] + wvi[k, 2] * (temp[k + 1] - temp[k])
-    end
-    
-    # Mean temperature in stratospheric layers
-    for k=1:n_stratosphere_levels-1
-        column.st4a[k, 2] = NF(0.75) * temp[k] + NF(0.25) * column.st4a[k, 1]
-        column.st4a[k+1, 2] = NF(0.50) * temp[k+1] + NF(0.25) * (column.st4a[k, 1] + column.st4a[k+1, 1])
-    end
-
-    # Temperature gradient in tropospheric layers
-    anis = NF(1.)
-    anish = NF(0.5) # 0.5 * anis
-
-    for k = n_stratosphere_levels+1:nlev-1
-        column.st4a[k, 2] = anish * max(column.st4a[k, 1] - column.st4a[k - 1, 1], NF(0.))
-    end
-    column.st4a[nlev, 2] = anis * max(temp[nlev] - column.st4a[nlev-1, 1], NF(0.))
-
-    # Blackbody emission in the stratosphere
-    for k = 1:n_stratosphere_levels
-        column.st4a[k, 1] = sbc * column.st4a[k, 2]^4
-        column.st4a[k, 2] = NF(0.)
-    end
-
-    # Blackbody emission in the troposphere
-    for k = n_stratosphere_levels+1:nlev
-            st3a = sbc * temp[k]^3
-            column.st4a[k, 1] = st3a * temp[k]
-            column.st4a[k, 2] = NF(4.) * st3a * column.st4a[k, 2]
-    end
-
-    # 2.0 Initialization of fluxes
-    column.fsfcd = 0.
-    column.dfabs .= 0.
-
-    # 3.0 Emission ad absorption of longwave downward flux.
-    #    For downward emission, a correction term depending on the
-    #    local temperature gradient and on the layer transmissivity is
-    #    added to the average (full - level) emission of each layer.
-
-    # 3.1  Stratosphere
-    k = n_stratosphere_levels-1
-    for jb = 1:nband-2
-        emis = 1. - tau2[k, jb]
-        brad = fband[convert(Int, temp[k]), jb] * (column.st4a[k, 1] + emis * column.st4a[k, 2])
-        column.flux[jb] = emis * brad
-        column.dfabs[k] = column.dfabs[k] - column.flux[jb]
-    end
-    column.flux[nband-1:nband] .= NF(0.)
-
-    # 3.2  Troposphere
-    for jb = 1:nband
-        for k = n_stratosphere_levels:nlev
-            emis = 1. - tau2[k, jb]
-            brad = fband[convert(Int, temp[k]), jb] * (column.st4a[k, 1] + emis * column.st4a[k, 2])
-            column.dfabs[k] = column.dfabs[k] + column.flux[jb]
-            column.flux[jb] = tau2[k, jb] * column.flux[jb] + emis * brad
-            column.dfabs[k] = column.dfabs[k] - column.flux[jb]
-        end
-    end
-
-    # 3.3 Surface downward flux
-    for jb = 1: nband
-        column.fsfcd = column.fsfcd + emisfc * column.flux[jb]
-    end
-    eps1 = epslw * emisfc
-    
-    # 3.4 Correction for "black" band (incl. surface reflection)
-    corlw = eps1 * column.st4a[nlev, 1]
-    column.dfabs[nlev] = column.dfabs[nlev] - corlw
-    column.fsfcd = column.fsfcd + corlw
-end
-
-"""
-    function compute_bbe!(
-        column::ColumnVariables{NF}, model::PrimitiveEquation
-    ) where {NF<:AbstractFloat}
-
-# Computes black-body (or grey-body) emissions.
-Input and output variables are `ts` and `fsfcu`, respectively.
-"""
-function compute_bbe!(
-    column::ColumnVariables{NF}, model::PrimitiveEquation
-) where {NF<:AbstractFloat}
-
-(; sbc ) = model.parameters
-(; emisfc ) = model.parameters.radiation_coefs
-(; ts ) = column
-
-    column.fsfcu = emisfc * sbc * ts^4
+function longwave_radiation!(
+    column::ColumnVariables,
+    scheme::UniformCooling,
+    model::PrimitiveEquation,
+)
+    longwave_radiation!(column,scheme)
 end
 
-"""
-    function radlw_up!(
-        column::ColumnVariables{NF}, model::PrimitiveEquation
-    ) where {NF<:AbstractFloat}
-
-# Computes the upward flux of longwave radiation.
-Input variables are  `nlev`, `temp`, `fsfcu`, `fsfcd`, `flux`, `ts`, `tau2`, `st4a`, `dfabs`, `stratc`, `σ_levels_thick`, `n_stratosphere_levels`. Output column variables are `fsfc` and `ftop`.
-"""
-function radlw_up!(
-    column::ColumnVariables{NF}, model::PrimitiveEquation
-) where {NF<:AbstractFloat}
-
-(; nband ) = model.parameters
-(; epslw, emisfc ) = model.parameters.radiation_coefs
-(; fband ) = model.parameterization_constants
-(; nlev, temp, fsfcu, fsfcd, flux, ts, tau2, st4a, dfabs, stratc ) = column
-(; σ_levels_thick, n_stratosphere_levels ) = model.geometry
-
-    column.fsfc = fsfcu - fsfcd
+function longwave_radiation!(
+    column::ColumnVariables{NF},
+    scheme::UniformCooling,
+) where NF
+    (;temp, temp_tend) = column
+    (;temp_min, temp_stratosphere) = scheme
     
-    refsfc = 1 - emisfc
-    for jb = 1: nband
-        column.flux[jb] = fband[convert(Int, ts), jb] * fsfcu + refsfc * flux[jb]
-    end
+    cooling = -inv(convert(NF,scheme.time_scale.value))
+    τ⁻¹ = inv(scheme.time_scale_stratosphere.value)
 
-    # 4.2  Troposphere
-    # Correction for "black" band
-    dfabs[nlev] = dfabs[nlev] + epslw * fsfcu
-    for jb = 1: nband
-        for k = nlev:-1:n_stratosphere_levels
-            emis = 1. - tau2[k, jb]
-            brad = fband[convert(Int, temp[k]), jb] * (st4a[k, 1] - emis * st4a[k, 2])
-            column.dfabs[k] = dfabs[k] + flux[jb]
-            column.flux[jb] = tau2[k, jb] * flux[jb] + emis * brad
-            column.dfabs[k] = dfabs[k] - flux[jb]
-        end
-    end
-    # 4.3  Stratosphere
-    k = n_stratosphere_levels - 1
-    for jb = 1: 2
-        emis = 1 - tau2[k, jb]
-        brad = fband[convert(Int, temp[k]), jb] * (st4a[k, 1] - emis * st4a[k, 2])
-        column.dfabs[k] = dfabs[k] + flux[jb]
-        column.flux[jb] = tau2[k, jb] * flux[jb] + emis * brad
-        column.dfabs[k] = dfabs[k] - flux[jb]
+    @inbounds for k in eachlayer(column)
+        temp_tend[k] += temp[k] > temp_min ? cooling : (temp_stratosphere - temp[k])*τ⁻¹
     end
-
-    # Correction for "black" band and polar night cooling
-    column.ftop = 0. # Init
-    for k=1:n_stratosphere_levels-1
-        corlw1 = σ_levels_thick[k] * stratc[k+1] * st4a[k, 1] + stratc[k]
-        corlw2 = σ_levels_thick[k+1] * stratc[k+1] * st4a[k+1, 1]
-        dfabs[k] -= corlw1
-        dfabs[k+1] -= corlw2
-        column.ftop += corlw1 + corlw2
-    end
-
-    # 4.4  Outgoing longwave radiation
-    for jb = 1:nband
-        column.ftop = column.ftop + flux[jb]
-    end
-end
+end
\ No newline at end of file
diff --git a/src/physics/shortwave_radiation.jl b/src/physics/shortwave_radiation.jl
index ed1c19f6e..fddd3dba7 100644
--- a/src/physics/shortwave_radiation.jl
+++ b/src/physics/shortwave_radiation.jl
@@ -1,389 +1,44 @@
-"""
-Parameters for radiation parameterizations.
-"""
-Base.@kwdef struct RadiationCoefs{NF<:Real} <: Coefficients
-    epslw::NF = 0.05    # Fraction of blackbody spectrum absorbed/emitted by PBL only
-    emisfc::NF = 0.98   # Longwave surface emissivity
+abstract type AbstractRadiation{NF} <: AbstractParameterization{NF} end
+abstract type AbstractShortwave{NF} <: AbstractRadiation{NF} end
 
-    p0::Real = 1e5      # Reference pressure (Pa)
+struct NoShortwave{NF} <: AbstractShortwave{NF} end
+NoShortwave(SG::SpectralGrid) = NoShortwave{SG.NF}()
+initialize!(::NoShortwave,::PrimitiveEquation) = nothing
 
-    # Shortwave radiation: sol_oz
-    solc::NF = 342.0    # Solar constant (area averaged) [W/m^2]
-    epssw::NF = 0.020   # Fraction of incoming solar radiation absorbed by ozone
-
-    # Shortwave radiation: cloud
-    rhcl1::NF = 0.30    # Relative humidity threshold corresponding to cloud cover = 0
-    rhcl2::NF = 1.00    # Relative humidity correponding to cloud cover = 1
-    rrcl::NF = 1 / (rhcl2 - rhcl1)
-    qcl::NF = 0.20      # Specific humidity threshold for cloud cover
-    pmaxcl::NF = 10.0   # Maximum value of precipitation (mm/day) contributing to cloud cover
-    wpcl::NF = 0.2      # Cloud cover weight for the square-root of precipitation (for p = 1 mm/day)
-    gse_s1::NF = 0.40   # Gradient of dry static energy corresponding to stratiform cloud cover = 1
-    gse_s0::NF = 0.25   # Gradient of dry static energy corresponding to stratiform cloud cover = 0
-    clsmax::NF = 0.60   # Maximum stratiform cloud cover
-    clsminl::NF = 0.15  # Minimum stratiform cloud cover over land (for RH = 1)
-
-    # Shortwave radiation: radsw
-    albcl::NF = 0.43    # Cloud albedo (for cloud cover = 1)
-    albcls::NF = 0.50   # Stratiform cloud albedo (for st. cloud cover = 1)
-    abscl1::NF = 0.015  # Absorptivity of clouds (visible band, maximum value)
-    abscl2::NF = 0.15   # Absorptivity of clouds (visible band, for dq_base = 1 g/kg)
-
-    absdry::NF = 0.033  # Absorptivity of dry air (visible band)
-    absaer::NF = 0.033  # Absorptivity of aerosols (visible band)
-    abswv1::NF = 0.022  # Absorptivity of water vapour (visible band, for dq = 1 g/kg)
-    abswv2::NF = 15.0   # Absorptivity of water vapour (near IR band, for dq = 1 g/kg)
-
-    ablwin::NF = 0.3    # Absorptivity of air in "window" band
-    ablco2::NF = 6.0    # Absorptivity of air in CO2 band
-    ablwv1::NF = 0.7    # Absorptivity of water vapour in H2O band 1 (weak), (for dq = 1 g/kg)
-    ablwv2::NF = 50.0   # Absorptivity of water vapour in H2O band 2 (strong), (for dq = 1 g/kg)
-    ablcl2::NF = 0.6    # Absorptivity of "thin" upper clouds in window and H2O bands
-    ablcl1::NF = 12.0   # Absorptivity of "thick" clouds in window band (below cloud top)
-end 
-
-"""
-    function shortwave_radiation!(
-        column::ColumnVariables{NF}, model::PrimitiveEquation
-    )
-
-Compute air temperature tendencies from shortwave radiation for an atmospheric column.
-For more details see http://users.ictp.it/~kucharsk/speedy_description/km_ver41_appendixA.pdf
-"""
-function shortwave_radiation!(
-    column::ColumnVariables{NF}, model::PrimitiveEquation
-) where {NF<:AbstractFloat}
-    (; humid, sat_vap_pres, dry_static_energy, geopot, norm_pres ) = column
-    (; cₚ ) = model.parameters
-    (; p0 ) = model.parameters.radiation_coefs
-    (; σ_levels_thick ) = model.geometry
-    (; gravity ) = model.constants
-
-    sol_oz!(column, model)
-
-    column.rel_hum .= humid ./ sat_vap_pres
-    column.grad_dry_static_energy = 
-        (dry_static_energy[end - 1] - dry_static_energy[end]) /
-        (geopot[end - 1] - geopot[end])
-    cloud!(column, model)
-
-    radsw!(column, model)
-
-    for k in eachlayer(column)
-        column.temp_tend[k] += 
-            column.tend_t_rsw[k] * inv(norm_pres) *
-            (gravity / (σ_levels_thick[k] * p0)) /cₚ
-    end
+# function barrier for all AbstractShortwave
+function shortwave_radiation!(column::ColumnVariables,model::PrimitiveEquation)
+    shortwave_radiation!(column,model.shortwave_radiation,model)
 end
 
-"""
-    function solar!(column::ColumnVariables{NF})
-
-Compute average daily flux of solar radiation for an atmospheric column,
-from Hartmann (1994).
-"""
-function solar!(column::ColumnVariables{NF}) where {NF<:AbstractFloat}
-    (; tyear, csol, latd ) = column
-
-    # Compute cosine and sine of latitude
-    clat = cos(latd)
-    slat = sin(latd)
-
-    # 1.0 Compute declination angle and Earth - Sun distance factor
-    pigr = 2 * asin(1)
-    α = 2 * pigr * tyear
-
-    ca1 = cos(α)
-    sa1 = sin(α)
-    ca2 = ca1 * ca1 - sa1 * sa1
-    sa2 = 2 * sa1 * ca1
-    ca3 = ca1 * ca2 - sa1 * sa2
-    sa3 = sa1 * ca2 + sa2 * ca1
-
-    decl =
-        NF(0.006918) - NF(0.399912) * ca1 + NF(0.070257) * sa1 - NF(0.006758) * ca2 + NF(0.000907) * sa2 -
-        NF(0.002697) * ca3 + NF(0.001480) * sa3
-
-    fdis = NF(1.000110) + NF(0.034221) * ca1 + NF(0.001280) * sa1 + NF(0.000719) * ca2 + NF(0.000077) * sa2
+shortwave_radiation!(::ColumnVariables,::NoShortwave,::PrimitiveEquation) = nothing
 
-    cdecl = cos(decl)
-    sdecl = sin(decl)
-    tdecl = sdecl / cdecl
-
-    # 2.0 Compute daily - average insolation at the atm. top
-    csolp = csol / pigr
-
-    ch0 = min(1, max(-1, -tdecl * slat / clat))
-    h0 = acos(ch0)
-    sh0 = sin(h0)
-
-    column.topsr = csolp * fdis * (h0 * slat * sdecl + sh0 * clat * cdecl)
+Base.@kwdef struct TransparentShortwave{NF} <: AbstractShortwave{NF}
+    albedo::NF = 0.3
+    S::Base.RefValue{NF} = Ref(zero(NF))
 end
 
-"""
-    function sol_oz!(
-        column::ColumnVariables{NF}, model::PrimitiveEquation
-    )
-
-Compute solar radiation parametres for an atmospheric column.
-"""
-function sol_oz!(
-    column::ColumnVariables{NF}, model::PrimitiveEquation
-) where {NF<:AbstractFloat}
-    (; tyear, latd ) = column
-    (; solc, epssw ) = model.parameters.radiation_coefs
-
-    # Compute cosine and sine of latitude
-    clat = cos(latd)
-    slat = sin(latd)
-
-    # α = year phase ( 0 - 2pi, 0 = winter solstice = 22dec.h00 )
-    α = 4 * asin(1) * (tyear + 10 / 365)
-    dα = NF(0.)
-
-    coz1 = 1 * max(0, cos(α - dα))
-    coz2 = NF(1.8)
-    azen = 1
-    nzen = 2
-
-    rzen = -cos(α) * NF(23.45) * asin(1.0) / 90
-    czen = cos(rzen)
-    szen = sin(rzen)
+TransparentShortwave(SG::SpectralGrid) = TransparentShortwave{SG.NF}()
 
-    fs0 = 6
-
-    # Solar radiation at the top
-    column.csol = 4 * solc
-    column.topsr = solar!(column)
-    column.fsol = column.topsr
-
-    flat2 = NF(1.5) * slat^2 - NF(0.5)
-
-    # Ozone depth in upper and lower stratosphere
-    column.ozupp = NF(0.5) * epssw
-    column.ozone = NF(0.4) * epssw * (1 + coz1 * slat + coz2 * flat2)
-
-    # Zenith angle correction to (downward) absorptivity
-    column.zenit = 1 + azen * (1 - (clat * czen + slat * szen))^nzen
-
-    # Ozone absorption in upper and lower stratosphere
-    column.ozupp = column.fsol * column.ozupp * column.zenit
-    column.ozone = column.fsol * column.ozone * column.zenit
-
-    # Polar night cooling in the stratosphere
-    column.stratz = max(fs0 - column.fsol, 0)
+function initialize!(scheme::TransparentShortwave,model::PrimitiveEquation)
+    (;solar_constant, gravity) = model.planet
+    (;cₚ, pres_ref) = model.atmosphere
+    Δσ = model.geometry.σ_levels_thick
+    scheme.S[] = (1 - scheme.albedo) * solar_constant * gravity / (Δσ[end] * (pres_ref*100) * cₚ)
 end
 
-"""
-    function cloud!(
-        column::ColumnVariables{NF}, model::PrimitiveEquation
-    )
-
-Compute shortwave radiation cloud contibutions for an atmospheric column.
-"""
-function cloud!(
-    column::ColumnVariables{NF}, model::PrimitiveEquation
-) where {NF<:AbstractFloat}
-    (; rhcl1, rhcl2, rrcl, qcl, pmaxcl ) = model.parameters.radiation_coefs
-    (; wpcl, gse_s1, gse_s0, clsmax, clsminl ) = model.parameters.radiation_coefs
-    (; humid, rel_hum, grad_dry_static_energy, precip_convection ) = column
-    (; precip_large_scale, cloud_top, nlev, fmask ) = column
-    (; n_stratosphere_levels ) = model.geometry
-
-    # 1.0 Cloud cover, defined as the sum of:
-    # - a term proportional to the square - root of precip. rate
-    # - a quadratic function of the max. relative humidity
-    #    in tropospheric layers above pbl where q > prog_qcl :
-    #     ( = 0 for rhmax < rhcl1, = 1 for rhmax > rhcl2 )
-    #    Cloud - top level: defined as the highest (i.e. least sigma)
-    #      between the top of convection / condensation and
-    #      the level of maximum relative humidity.
-    if (rel_hum[nlev - 1] > rhcl1)
-        column.cloudc = rel_hum[nlev - 1] - rhcl1
-        column.icltop = nlev - 1
-    else
-        column.cloudc = NF(0)
-        column.icltop = nlev + 1
-    end
-
-    for k in n_stratosphere_levels+1:(nlev - n_stratosphere_levels)
-        drh = rel_hum[k] - rhcl1
-        if (drh > column.cloudc) & (humid[k] > qcl)
-            column.cloudc = drh
-            column.icltop = k
-        end
-    end
-
-    pr1 = min(pmaxcl, NF(86.4) * (precip_convection + precip_large_scale))
-    column.cloudc = min(1, wpcl * sqrt(pr1) + min(1, column.cloudc * rrcl)^2)
-    column.icltop = min(cloud_top, column.icltop)
-
-    # 2.0  Equivalent specific humidity of clouds
-    column.qcloud = humid[nlev - 1]
-
-    # 3.0 Stratiform clouds at the top of pbl
-    clfact = NF(1.2)
-    rgse = 1 / (gse_s1 - gse_s0)
-
-    # Stratocumulus clouds over sea
-    fstab = max(0, min(1, rgse * (grad_dry_static_energy - gse_s0)))
-    column.clstr = fstab * max(clsmax - clfact * column.cloudc, 0)
-
-    # Stratocumulus clouds over land
-    clstrl = max(column.clstr, clsminl) * rel_hum[nlev]
-    column.clstr = column.clstr + fmask * (clstrl - column.clstr)
+function shortwave_radiation!(
+    column::ColumnVariables,
+    scheme::TransparentShortwave,
+    model::PrimitiveEquation,
+)
+    shortwave_radiation!(column,scheme,model.solar_zenith)
 end
 
-"""
-    function radsw!(
-        column::ColumnVariables{NF}, model::PrimitiveEquation
-    )
-
-Compute shortwave radiation fluxes for an atmospheric column.
-"""
-function radsw!(
-    column::ColumnVariables{NF}, model::PrimitiveEquation
-) where {NF<:AbstractFloat}
-    (; norm_pres, humid, icltop, cloudc, clstr, ozupp, ozone ) = column
-    (; zenit, stratz, fsol, qcloud, albsfc, nlev ) = column
-    (; σ_levels_full, σ_levels_thick, n_stratosphere_levels ) = model.geometry
-    (; albcl, albcls, abscl1, abscl2, absdry, absaer ) = model.parameters.radiation_coefs
-    (; abswv1, abswv2, ablwin, ablco2, ablwv1 ) = model.parameters.radiation_coefs
-    (; ablwv2, ablcl2, ablcl1, epslw ) = model.parameters.radiation_coefs
-
-    # Locals variables
-    sbands_flux = 2
-    flux::Vector{NF} = fill(NaN, sbands_flux)
-    stratc::Vector{NF} = fill(NaN, n_stratosphere_levels)
-
-    nl1 = nlev - 1
-
-    fband2 = NF(0.05)
-    fband1 = 1 - fband2
-
-    # 1.0 Initialization
-    column.tau2 = fill!(column.tau2, 0)
-
-    # Change to ensure only ICLTOP < = NLEV used
-    if (icltop <= nlev)
-        column.tau2[icltop, 3] = albcl * cloudc
-    end
-    column.tau2[nlev, 3] = albcls * clstr
-
-    # 2.0 Shortwave transmissivity:
-    # function of layer mass, ozone (in the statosphere),
-    # abs. humidity and cloud cover (in the troposphere)
-    psaz = norm_pres * zenit
-    acloud = cloudc * min(abscl1 * qcloud, abscl2)
-
-    deltap = psaz * σ_levels_thick[1]
-    column.tau2[1, 1] = exp(-deltap * absdry)
-
-    for k in n_stratosphere_levels:nl1
-        abs1 = absdry + absaer * σ_levels_full[k]^2
-        deltap = psaz * σ_levels_thick[k]
-        if (k >= icltop)
-            column.tau2[k, 1] = exp(-deltap * (abs1 + abswv1 * humid[k] + acloud))
-        else
-            column.tau2[k, 1] = exp(-deltap * (abs1 + abswv1 * humid[k]))
-        end
-    end
-
-    abs1 = absdry + absaer * σ_levels_full[nlev]^2
-    deltap = psaz * σ_levels_thick[nlev]
-    column.tau2[nlev, 1] = exp(-deltap * (abs1 + abswv1 * humid[nlev]))
-
-    for k in n_stratosphere_levels:nlev
-        deltap = psaz * σ_levels_thick[k]
-        column.tau2[k, 2] = exp(-deltap * abswv2 * humid[k])
-    end
-
-    # 3.0 Shortwave downward flux
-    # 3.1 Initialization of fluxes
-    column.tsr = fsol
-    flux[1] = fsol * fband1
-    flux[2] = fsol * fband2
-
-    # 3.2 Ozone and dry - air absorption in the stratosphere
-    for k in 1:n_stratosphere_levels
-        if k == 1 ozone_tmp = ozupp else ozone_tmp = ozone end
-        column.tend_t_rsw[k] = flux[1]
-        flux[1] = column.tau2[k, 1] * (flux[1] - ozone_tmp * norm_pres)
-        column.tend_t_rsw[k] = column.tend_t_rsw[k] - flux[1]
-    end
-    
-    # 3.3  Absorption and reflection in the troposphere
-    for k in n_stratosphere_levels+1:nlev
-        column.tau2[k, 3] = flux[1] * column.tau2[k, 3]
-        flux[1] = flux[1] - column.tau2[k, 3]
-        column.tend_t_rsw[k] = flux[1]
-        flux[1] = column.tau2[k, 1] * flux[1]
-        column.tend_t_rsw[k] = column.tend_t_rsw[k] - flux[1]
-    end
-
-    for k in n_stratosphere_levels:nlev
-        column.tend_t_rsw[k] = column.tend_t_rsw[k] + flux[2]
-        flux[2] = column.tau2[k, 2] * flux[2]
-        column.tend_t_rsw[k] = column.tend_t_rsw[k] - flux[2]
-    end
-
-    # 4.0 Shortwave upward flux
-    # 4.1  Absorption and reflection at the surface
-    column.ssrd = flux[1] + flux[2]
-    flux[1] = flux[1] * albsfc
-    column.ssr = column.ssrd - flux[1]
-
-    # 4.2  Absorption of upward flux
-    for k in nlev:-1:1
-        column.tend_t_rsw[k] = column.tend_t_rsw[k] + flux[1]
-        flux[1] = column.tau2[k, 1] * flux[1]
-        column.tend_t_rsw[k] = column.tend_t_rsw[k] - flux[1]
-        flux[1] = flux[1] + column.tau2[k, 3]
-    end
-
-    # 4.3  Net solar radiation = incoming - outgoing
-    column.tsr = column.tsr - flux[1]
-
-    # 5.0  Initialization of longwave radiation model
-    # 5.1  Longwave transmissivity:
-    #      function of layer mass, abs. humidity and cloud cover.
-
-    #    Cloud - free levels (stratosphere + PBL)
-    k = 1
-    deltap = norm_pres * σ_levels_thick[k]
-    column.tau2[k, 1] = exp(-deltap * ablwin)
-    column.tau2[k, 2] = exp(-deltap * ablco2)
-    column.tau2[k, 3] = 1.
-    column.tau2[k, 4] = 1.
-
-    for k in n_stratosphere_levels:(nlev - n_stratosphere_levels):nlev
-        deltap = norm_pres * σ_levels_thick[k]
-        column.tau2[k, 1] = exp(-deltap * ablwin)
-        column.tau2[k, 2] = exp(-deltap * ablco2)
-        column.tau2[k, 3] = exp(-deltap * ablwv1 * humid[k])
-        column.tau2[k, 4] = exp(-deltap * ablwv2 * humid[k])
-    end
-
-    # Cloudy layers (free troposphere)
-    acloud = cloudc * ablcl2
-
-    for k in n_stratosphere_levels+1:nl1
-        deltap = norm_pres * σ_levels_thick[k]
-        if (k < icltop)
-            acloud1 = acloud
-        else
-            acloud1 = ablcl1 * cloudc
-        end
-        column.tau2[k, 1] = exp(-deltap * (ablwin + acloud1))
-        column.tau2[k, 2] = exp(-deltap * ablco2)
-        column.tau2[k, 3] = exp(-deltap * max(ablwv1 * humid[k], acloud))
-        column.tau2[k, 4] = exp(-deltap * max(ablwv2 * humid[k], acloud))
-    end
-
-    # 5.2  Stratospheric correction terms
-    stratc[1] = stratz * norm_pres
-    for k in 2:n_stratosphere_levels
-        eps1 = epslw / (σ_levels_thick[k-1] + σ_levels_thick[k])
-        stratc[k] = eps1 * norm_pres
-    end
-end
+function shortwave_radiation!(
+    column::ColumnVariables,
+    scheme::TransparentShortwave,
+    solar_zenith::AbstractZenith,
+)
+    cos_zenith = solar_zenith.cos_zenith[column.ij]
+    column.temp_tend[end] += cos_zenith*scheme.S[]
+end
\ No newline at end of file
diff --git a/src/physics/surface_fluxes.jl b/src/physics/surface_fluxes.jl
index 12165fb02..09690be76 100644
--- a/src/physics/surface_fluxes.jl
+++ b/src/physics/surface_fluxes.jl
@@ -194,7 +194,7 @@ end
 # function barrier
 function evaporation!(  column::ColumnVariables,
                         model::PrimitiveWet)
-    evaporation!(column,model.evaporation,model.clausis_clapeyron)
+    evaporation!(column,model.evaporation,model.clausius_clapeyron)
 end
 
 function evaporation!(  column::ColumnVariables{NF},
diff --git a/src/physics/tendencies.jl b/src/physics/tendencies.jl
index abfddcdcc..e1eb51a27 100644
--- a/src/physics/tendencies.jl
+++ b/src/physics/tendencies.jl
@@ -53,12 +53,12 @@ function parameterization_tendencies!(
     static_energy_diffusion!(column,model)
     humidity_diffusion!(column,model)
 
-    # Calculate parametrizations (order of execution is important!)
+    # Calculate parametrizations
     convection!(column,model)
     large_scale_condensation!(column,model)
     # clouds!(column, model)
-    # shortwave_radiation!(column,model)
-    # longwave_radiation!(column,model)
+    shortwave_radiation!(column,model)
+    longwave_radiation!(column,model)
     surface_fluxes!(column,model)
 
     # sum fluxes on half levels up and down for every layer
diff --git a/src/physics/thermodynamics.jl b/src/physics/thermodynamics.jl
index 3d8564bee..33f6c6c95 100644
--- a/src/physics/thermodynamics.jl
+++ b/src/physics/thermodynamics.jl
@@ -94,15 +94,15 @@ end
 $(TYPEDSIGNATURES)
 Saturation humidity [kg/kg] from temperature [K], pressure [Pa] via
 
-    sat_vap_pres = clausis_clapeyron(temperature)
+    sat_vap_pres = clausius_clapeyron(temperature)
     saturation humidity = mol_ratio * sat_vap_pres / pressure"""
 function saturation_humidity(
     temp_kelvin::NF,
     pres::NF,
-    clausis_clapeyron::AbstractClausiusClapeyron{NF}
+    clausius_clapeyron::AbstractClausiusClapeyron{NF}
 ) where NF
-    sat_vap_pres = clausis_clapeyron(temp_kelvin)
-    return saturation_humidity(sat_vap_pres,pres;mol_ratio=clausis_clapeyron.mol_ratio)
+    sat_vap_pres = clausius_clapeyron(temp_kelvin)
+    return saturation_humidity(sat_vap_pres,pres;mol_ratio=clausius_clapeyron.mol_ratio)
 end
 
 """
@@ -183,11 +183,10 @@ function saturation_humidity!(
     column::ColumnVariables,
     clausius_clapeyron::AbstractClausiusClapeyron,
 )
-    (;sat_humid, rel_humid, pres, temp, humid) = column
+    (;sat_humid, pres, temp) = column
 
     for k in eachlayer(column)
         sat_humid[k] = saturation_humidity(temp[k],pres[k],clausius_clapeyron)
-        rel_humid[k] = humid[k]/sat_humid[k]
     end
 end
 
diff --git a/test/column_variables.jl b/test/column_variables.jl
index 030da6a88..a476422da 100644
--- a/test/column_variables.jl
+++ b/test/column_variables.jl
@@ -12,14 +12,10 @@
 
         # Convection
         @test column.cloud_top === column.nlev+1
-        @test column.conditional_instability === false
-        @test column.activate_convection === false
-        @test column.excess_humid === zero(NF)
-        @test column.precip_convection === zero(NF)
-        @test column.cloud_base_mass_flux === zero(NF)
 
         # Large-scale condensation
         @test column.precip_large_scale === zero(NF)
+        @test column.precip_convection === zero(NF)
     end
 end
 
diff --git a/test/convection.jl b/test/convection.jl
deleted file mode 100644
index f759dbd46..000000000
--- a/test/convection.jl
+++ /dev/null
@@ -1,45 +0,0 @@
-@testset "diagnose_convection!" begin
-    @testset for NF in (Float32, Float64)
-        _, diagn, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation)
-        (; nlev,σ_levels_full) = model.geometry
-        pres_thresh_cnv = model.constants.pres_thresh_cnv
-        
-        column = ColumnVariables{NF}(nlev = nlev)
-        column.humid .= 0 .+ 50 * rand(NF, nlev)         # Typical values between 0-50 g/kg
-        column.pres .= (10*rand(NF)+1e5)*vcat(σ_levels_full,1)   # Typical values around 1e5 Pa
-        column.sat_humid .= 0 .+ 50 * rand(NF, nlev)     # Typical values between 0-50 g/kg
-        column.dry_static_energy .= rand(NF, nlev)
-        column.moist_static_energy .= rand(NF, nlev)
-        column.sat_moist_static_energy .= rand(NF, nlev)
-        column.sat_moist_static_energy_half .= rand(NF, nlev)
-        
-        SpeedyWeather.diagnose_convection!(column, model)
-        
-        @test isfinite(column.excess_humidity)
-        @test column.excess_humidity >= 0
-    end
-end
-
-@testset "convection!" begin
-    @testset for NF in (Float32, Float64)
-        _, diagn, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation)
-        (; nlev,σ_levels_full) = model.geometry
-        pres_thresh_cnv = model.constants.pres_thresh_cnv
-        
-        column = ColumnVariables{NF}(nlev = nlev)
-        column.humid .= 0 .+ 50 * rand(NF, nlev)         # Typical values between 0-50 g/kg
-        column.pres .= (10*rand(NF)+1e5)*vcat(σ_levels_full,1)   # Typical values around 1e5 Pa
-        column.sat_humid .= 0 .+ 50 * rand(NF, nlev)     # Typical values between 0-50 g/kg
-        column.dry_static_energy .= rand(NF, nlev)
-        column.moist_static_energy .= rand(NF, nlev)
-        column.sat_moist_static_energy .= rand(NF, nlev)
-        column.sat_moist_static_energy_half .= rand(NF, nlev)
-        
-        SpeedyWeather.convection!(column, model)
-        
-        @test isfinite(column.cloud_base_mass_flux)
-        @test isfinite(column.precip_convection)
-        @test all(isfinite.(column.net_flux_humid))
-        @test all(isfinite.(column.net_flux_dry_static_energy))
-    end
-end
diff --git a/test/large_scale_condensation.jl b/test/large_scale_condensation.jl
deleted file mode 100644
index 4c1306e1c..000000000
--- a/test/large_scale_condensation.jl
+++ /dev/null
@@ -1,13 +0,0 @@
-@testset "Parametrization: large scale condensation" begin
-    @testset for NF in (Float32,Float64)
-        _, diagn, model = SpeedyWeather.initialize_speedy(NF,PrimitiveEquation)
-
-        column = ColumnVariables{NF}(nlev=diagn.nlev)
-
-        for ij in SpeedyWeather.eachgridpoint(diagn)
-            SpeedyWeather.reset_column!(column)
-            SpeedyWeather.get_column!(column,diagn,ij,model.geometry)
-            SpeedyWeather.large_scale_condensation!(column, model)
-        end
-    end
-end
diff --git a/test/longwave_radiation.jl b/test/longwave_radiation.jl
deleted file mode 100644
index 8ef85c8cb..000000000
--- a/test/longwave_radiation.jl
+++ /dev/null
@@ -1,95 +0,0 @@
-@testset "Longwave radiation: radset!" begin
-    @testset for NF in (Float32, Float64)
-        _, diagn, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation)
-        
-        # Just check the last band
-        @test isapprox(
-            model.parameterization_constants.fband[end, :],
-            [0.19498351, 0.12541235, 0.33106664, 0.2985375],
-            rtol=0.1
-            )
-    end
-end
-
-@testset "Longwave radiation: radlw_down!" begin
-    @testset for NF in (Float32, Float64)
-        _, diagn, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation, nlev=8)
-
-        nlev = model.parameters.nlev
-        nband = model.parameters.nband
-        n_stratosphere_levels = model.geometry.n_stratosphere_levels
-        column = ColumnVariables{NF}(nlev=nlev, nband=nband, n_stratosphere_levels=n_stratosphere_levels)
-        column.temp = fill(300., nlev)
-        column.wvi = fill(0.5, nlev, 2)
-        column.tau2 = fill(0.5, nlev, 4)
-        
-        SpeedyWeather.radlw_down!(column, model)
-
-        # Just check fsfcd and dfabs as used by radlw_up! 
-        @test column.fsfcd ≈ 447.29436216
-        @test column.dfabs ≈ [-71.86727228, -182.21961386, -91.10980693, -45.55490346, -22.77745173, -11.38872587, -5.69436293, -25.35141147]
-    end
-end
-
-@testset "Longwave radiation: compute_bbe!" begin
-    @testset for NF in (Float32, Float64)
-        _, diagn, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation, nlev=8)
-
-        nlev = model.parameters.nlev
-        nband = model.parameters.nband
-        n_stratosphere_levels = model.geometry.n_stratosphere_levels
-        column = ColumnVariables{NF}(nlev=nlev, nband=nband, n_stratosphere_levels=n_stratosphere_levels)
-        column.ts = 320.
-        
-        SpeedyWeather.compute_bbe!(column, model)
-
-        @test column.fsfcu ≈ 582.65174016
-    end
-end
-
-@testset "radlw_up!" begin
-    @testset for NF in (Float32, Float64)
-        _, diagn, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation, nlev=8)
-
-        nlev = model.parameters.nlev
-        nband = model.parameters.nband
-        n_stratosphere_levels = model.geometry.n_stratosphere_levels
-        column = ColumnVariables{NF}(nlev=nlev, nband=nband, n_stratosphere_levels=n_stratosphere_levels)
-        column.temp = fill(300., nlev)
-        column.wvi = fill(0.5, nlev, 2)
-        column.tau2 = fill(0.5, nlev, 4)
-        column.ts = 320.
-        column.stratc = fill(0.5, 2)
-        
-        SpeedyWeather.radlw_down!(column, model)
-        SpeedyWeather.compute_bbe!(column, model)
-        SpeedyWeather.radlw_up!(column, model)
-
-        # Just check what's needed
-        @test column.fsfc ≈ 135.357378
-        @test_skip column.ftop ≈ 474.76064406
-    end
-end
-
-@testset "longwave_radiation!" begin
-    @testset for NF in (Float32, Float64)
-        _, diagn, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation, nlev=8)
-
-        nlev = model.parameters.nlev
-        nband = model.parameters.nband
-        n_stratosphere_levels = model.geometry.n_stratosphere_levels
-        column = ColumnVariables{NF}(nlev=nlev, nband=nband, n_stratosphere_levels=n_stratosphere_levels)
-        column.temp = fill(300., nlev)
-        column.wvi = fill(0.5, nlev, 2)
-        column.tau2 = fill(0.5, nlev, 4)
-        column.ts = 320.
-        column.stratc = fill(0.5, 2)
-        
-        SpeedyWeather.longwave_radiation!(column, model)
-
-        # Just check what's needed
-        @test column.fsfc ≈ 135.357378
-        @test_skip column.ftop ≈ 474.76064406
-    end
-end
-
diff --git a/test/shortwave_radiation.jl b/test/shortwave_radiation.jl
deleted file mode 100644
index ff9d724a7..000000000
--- a/test/shortwave_radiation.jl
+++ /dev/null
@@ -1,156 +0,0 @@
-@testset "Shortwave radiation: solar!" begin
-    @testset for NF in (Float32, Float64)
-        topsr = [0., 283.06, 352.56537, 378.8209880964894]
-        latd = [-89., 0., 45., 89.]
-        @testset for i in 1:4
-            _, diagn, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation, nlev = 4)
-            nlev = diagn.nlev
-            nband = model.parameters.nband
-
-            column = ColumnVariables{NF}(nlev=nlev, nband=nband)
-            column.tyear = 0.5
-            column.csol = 1000.
-            column.latd = latd[i]
-
-            SpeedyWeather.solar!(column)
-
-            @test_skip column.topsr ≈ topsr[i]
-        end
-    end
-end
-
-@testset "Shortwave radiation: sol_oz!" begin
-    @testset for NF in (Float32, Float64)
-        _, diagn, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation, nlev = 4)
-        nlev = diagn.nlev
-        nband = model.parameters.nband
-
-        column = ColumnVariables{NF}(nlev=nlev, nband=nband)
-        column.tyear = 0.5
-        column.latd = 89.
-
-        SpeedyWeather.sol_oz!(column, model)
-
-        @test_skip column.topsr ≈ 518.2271117159975
-        @test_skip column.fsol ≈ 518.2271117159975
-        @test_skip column.ozupp ≈ 6.99611021887365
-        @test_skip column.ozone ≈ 15.666684
-        @test_skip column.zenit ≈ 1.3500085311452612
-        @test column.stratz ≈ 0.
-    end
-end
-
-@testset "Shortwave radiation: cloud!" begin
-    @testset for NF in (Float32, Float64)
-        _, diagn, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation, nlev = 4)
-        nlev = diagn.nlev
-        nband = model.parameters.nband
-
-        column = ColumnVariables{NF}(nlev=nlev, nband=nband)
-        column.humid .= [0., 0.03124937, 0.9748285, 6.7846994] # g/kg
-        column.rel_hum = [0.0, 1., 1., 0.8964701087948754]
-        column.grad_dry_static_energy = 0.4255393541723314
-        column.precip_convection = 1
-        column.precip_large_scale = 1
-        column.cloud_top = 1
-        column.fmask = 1.
-
-        SpeedyWeather.cloud!(column, model)
-
-        @test column.icltop ≈ 1.
-        @test column.cloudc ≈ 1.
-        @test column.clstr ≈ 0.13447051631923132
-        @test column.qcloud ≈ 0.9748285
-    end
-end
-
-@testset "Shortwave radiation: radsw!" begin
-    @testset for NF in (Float32, Float64)
-        _, diagn, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation, nlev = 4)
-        nlev = diagn.nlev
-        nband = model.parameters.nband
-
-        column = ColumnVariables{NF}(nlev=nlev, nband=nband)
-        column.norm_pres = 1.
-        column.humid .= [0., 0.03124937, 0.9748285, 6.7846994] # g/kg
-        column.albsfc = 0.5
-        # Same as returned from sol_oz and cloud
-        column.icltop = 1
-        column.cloudc = 0.90186596
-        column.clstr = 0.1485
-        column.ozupp = 6.99611021887365
-        column.ozone = 15.666684
-        column.zenit = 1.3500085311452612
-        column.stratz = 0.
-        column.fsol = 518.2271117159975
-        column.qcloud = 0.033334
-
-        SpeedyWeather.radsw!(column, model)
-
-        @test_skip column.ssrd ≈ 385.7997293028816
-        @test_skip column.ssr ≈ 192.8998646514408
-        @test_skip column.tsr ≈ 315.10016
-        @test_skip column.tau2 ≈ [0.9526258427031802 0.3788324531779398 1.0 1.0;
-                                0.9286716429916902 0.2276374390970355 0.994618802250959 0.6801722651144487;
-                                0.02148265489609933 0.12596241035550795 0.7900788064253661 4.906182803130207e-8;
-                                0.9287847242363858 0.22819245387064319 0.31050206404533354 5.234705430612321e-37]
-        @test_skip column.tend_t_rsw ≈ [11.947794979452055, 27.611640775148402, 42.565299927060124, 40.46336031165737]
-    end
-end
-
-@testset "shortwave_radiation!" begin
-    @testset for NF in (Float32, Float64)
-        _, diagn, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation, nlev = 4)
-        nlev = diagn.nlev
-        nband = model.parameters.nband
-
-        column = ColumnVariables{NF}(nlev=nlev, nband=nband)
-
-        # 1. Set variables for sol_oz
-        column.tyear = 0.5
-        column.latd = 89.
-
-        # 2. Compute sat_vap_pres and dry_static_energy
-        # and set remainig varables for cloud
-        column.pres .= 1. # nomalised
-        column.temp .= [208.40541, 219.8126 , 249.25502, 276.14264]
-        column.humid .= [0., 0.03124937, 0.9748285, 6.7846994] # g/kg
-        column.geopot .= [241932.19, 117422.14, 54618.79, 7626.5884]
-        SpeedyWeather.get_thermodynamics!(column, model)
-
-        # column.sat_vap_pres = [0.005265206274688095, 0.02494611200040165, 0.7012750147882735, 7.56823828640608]
-        # column.dry_static_energy = [451171.22164, 338113.9904, 304870.83008, 284873.79896]
-        # column.rel_hum = [0.0, 1.2526749659224194, 1.3900801817305823, 0.8964701087948754]
-        # column.grad_dry_static_energy = 0.4255393541723314
-
-        column.precip_convection = 1
-        column.precip_large_scale = 1
-        column.cloud_top = 1
-        column.fmask = 1.
-
-        # Set variables for radsw
-        # FIXME: pres is overloaded and will need to be
-        # fixed in other functions
-        column.norm_pres = column.pres[end]
-        column.albsfc = 0.5
-
-        SpeedyWeather.shortwave_radiation!(column, model)
-
-        # Results need to be close to those computed by radsw but
-        # are not the same as precision is maintained between functions
-        @test isapprox(column.ssrd, 385.7997293028816, rtol=0.1)
-        @test isapprox(column.ssr, 192.8998646514408, rtol=0.1)
-        @test_skip isapprox(column.tsr, 315.10016, rtol=0.1)
-        @test_skip isapprox(column.tau2, [
-            0.9526258427031802 0.3788324531779398 1.0 1.0;
-            0.9286716429916902 0.2276374390970355 0.994618802250959 0.6801722651144487;
-            0.02148265489609933 0.12596241035550795 0.7900788064253661 4.906182803130207e-8;
-            0.9287847242363858 0.22819245387064319 0.31050206404533354 5.234705430612321e-37
-            ], rtol=0.1)
-        @test_skip isapprox(column.tend_t_rsw, [
-            11.947794979452055, 27.611640775148402, 42.565299927060124, 40.46336031165737
-            ], rtol=0.1)
-        # Range is generally order 10^-5
-        @test_skip column.temp_tend ≈ [7.189150408440864e-6, 1.2141337616807728e-5, 1.3196339560843619e-5, 1.5994200814647754e-5]
-    end
-end
diff --git a/test/thermodynamics.jl b/test/thermodynamics.jl
deleted file mode 100644
index 435373d48..000000000
--- a/test/thermodynamics.jl
+++ /dev/null
@@ -1,117 +0,0 @@
-@testset "Thermodynamics: get_thermodynamics!" begin
-    @testset for NF in (Float32, Float64)
-        _, _, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation)
-        (; nlev,σ_levels_full) = model.geometry
-        
-        column = ColumnVariables{NF}(; nlev)
-        column.temp .= 200 .+ 150 * rand(NF, nlev)       # Typical values between 200-350K
-        column.humid .= 0 .+ 50 * rand(NF, nlev)         # Typical values between 0-50 g/kg
-        column.pres .= (10*rand(NF)+1e5)*vcat(σ_levels_full,1)   # Typical values around 1e5 Pa
-        column.geopot .= rand(NF, nlev)                  # What are typical values for geopotential?
-        
-        # For now, test that it runs with no errors
-        SpeedyWeather.get_thermodynamics!(column, model)
-    end
-end
-
-@testset "Thermodynamics: interpolate!" begin
-    @testset for NF in (Float32, Float64)
-        _, _, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation)
-        (; nlev) = model.geometry
-        
-        column = ColumnVariables{NF}(; nlev)
-        
-        A_full_level = rand(NF, nlev)
-        A_half_level = zeros(NF, nlev)
-        
-        SpeedyWeather.interpolate!(A_full_level, A_half_level, column, model)
-        
-        # Test that the half-level values lie between the enclosing full level values.
-        @test all(
-            (A_full_level[2:nlev] .> A_half_level[1:nlev-1] .> A_full_level[1:nlev-1])
-            .||
-            (A_full_level[2:nlev] .< A_half_level[1:nlev-1] .< A_full_level[1:nlev-1])
-            )
-    end
-end
-
-@testset "Saturation vapour pressure" begin
-    @testset for NF in (Float32, Float64)
-        _, diagn, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation)
-        (; nlev) = model.geometry
-        
-        column = ColumnVariables{NF}(; nlev)
-        column.temp .= 200 .+ 150 * rand(NF, nlev)  # Typical values between 200-350K
-        
-        SpeedyWeather.saturation_vapour_pressure!(column, model)
-        
-        @test all(column.sat_vap_pres .> 0.0)
-        @test all(column.sat_vap_pres .< 500.0)
-    end
-end
-
-@testset "Saturation specific humidity" begin
-    @testset for NF in (Float32, Float64)
-        _, diag, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation)
-        (; nlev,σ_levels_full) = model.geometry
-        
-        column = ColumnVariables{NF}(; nlev)
-        
-        column.temp = 200 .+ 150 * rand(NF, nlev)     # Typical values between 200-350 K
-        column.pres .= (10*rand(NF)+1e5)*vcat(σ_levels_full,1)   # Typical values around 1e5 Pa
-        
-        SpeedyWeather.saturation_vapour_pressure!(column, model)
-        SpeedyWeather.saturation_specific_humidity!(column, model)
-        
-        @test all(isfinite.(column.sat_humid))
-        @test !any(iszero.(column.sat_humid))
-    end
-end
-
-@testset "Dry static energy" begin
-    @testset for NF in (Float32, Float64)
-        _, _, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation)
-        (; nlev, σ_levels_full) = model.geometry
-        
-        column = ColumnVariables{NF}(; nlev)
-        column.temp .= 200 .+ 150 * rand(NF, nlev)  # Typical values between 200-350K
-        column.pres .= (10*rand(NF)+1e5)*vcat(σ_levels_full,1)   # Typical values around 1e5 Pa
-        
-        SpeedyWeather.dry_static_energy!(column, model)
-        
-        @test all(isfinite.(column.dry_static_energy))
-        @test !any(iszero.(column.dry_static_energy))
-    end
-end
-
-@testset "Moist static energy" begin
-    @testset for NF in (Float32, Float64)
-        _, _, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation)
-        (; nlev) = model.geometry
-        
-        column = ColumnVariables{NF}(; nlev)
-        column.dry_static_energy .= rand(NF, nlev)
-        column.humid .= 0 .+ 50 * rand(NF, nlev)         # Typical values between 0-50 g/kg
-        
-        SpeedyWeather.moist_static_energy!(column, model)
-        
-        @test all(isfinite.(column.moist_static_energy))
-        @test !any(iszero.(column.moist_static_energy))
-    end
-end
-
-@testset "Saturation moist static energy" begin
-    @testset for NF in (Float32, Float64)
-        _, _, model = SpeedyWeather.initialize_speedy(NF, PrimitiveEquation)
-        (; nlev) = model.geometry
-        
-        column = ColumnVariables{NF}(; nlev)
-        column.dry_static_energy .= rand(NF, nlev)
-        column.sat_humid .= 0 .+ 50 * rand(NF, nlev)         # Typical values between 0-50 g/kg
-        
-        SpeedyWeather.saturation_moist_static_energy!(column, model)
-        
-        @test all(isfinite.(column.sat_moist_static_energy))
-        @test !any(iszero.(column.sat_moist_static_energy))
-    end
-end
\ No newline at end of file