All parameterizations at prev timestep, flux limiters removed, no sho…

…rtwave (#555)

* All parameterizations at prev timestep, flux limiters back to pre v0.10

* recalculate virtual temperature

* Take all flux limiters out

* Increase critical Richardson number to 10

* Remove all flux limiters

* TransparentShortwave back in

* surface heat flux max flux removed

* NoShortwave as default for primitive models

docs/src/surface_fluxes.md

@@ -191,11 +191,6 @@ with ``\rho_s = \frac{p_s}{R_d T_N}`` the surface air density calculated from
 surface pressure ``p_s`` and lowermost layer temperature ``T_N``.
 Better would be to extrapolate ``T_N`` to ``T_s`` a surface air temperature
 assuming adiabatic descent but that is currently not implemented.
-In practice we use a flux limiter for numerical stability which limits the magnitude
-(preserving the sign). Choosing ``F_{uv, max}^\uparrow = 0.5 Pa`` would mean
-that the drag for winds faster than about ``33~m/s`` (typical ``C = 5\cdot 10^{-4}``)
-does not further increase. This can help to prevent oscillations drag terms can produce
-for sudden strong wind gusts.
 
 ## Surface heat fluxes
 
@@ -212,9 +207,6 @@ land surface layer or the mixed layer in the ocean. We then compute
 F_T^\uparrow = \rho_s C V_0 c_p (T_{skin} - T_s)
 ```
 
-and apply a similar flux limiter as for the momentum flux to prevent a sudden strong heating
-or cooling. For ``F_{T, max}^\uparrow = 100~W/m^2`` 
-
 ## Surface evaporation
 
 The surface moisture flux, i.e. evaporation of soil moisture over land and evaporation of

src/dynamics/planet.jl

@@ -1,6 +1,6 @@
 abstract type AbstractPlanet <: AbstractModelComponent end
 
-const DEFAULT_ROTATION = 7.29e-5    # default angular frequency of Earth's rotation [1/s]
+const DEFAULT_ROTATION = 7.29e-5    # default angular frequency of Earth's rotation [1/s]
 const DEFAULT_GRAVITY = 9.81        # default gravitational acceleration on Earth [m/s²]
 
 export Earth
@@ -12,7 +12,7 @@ characteristics. Note that `radius` is not part of it as this should be chosen
 in `SpectralGrid`. Keyword arguments are
 $(TYPEDFIELDS)
 """
-Base.@kwdef mutable struct Earth{NF<:AbstractFloat} <: AbstractPlanet
+@kwdef mutable struct Earth{NF<:AbstractFloat} <: AbstractPlanet
 
     "angular frequency of Earth's rotation [rad/s]"
     rotation::NF = DEFAULT_ROTATION
@@ -39,7 +39,7 @@ Base.@kwdef mutable struct Earth{NF<:AbstractFloat} <: AbstractPlanet
     axial_tilt::NF = 23.4
 
     "Total solar irradiance at the distance of 1 AU [W/m²]"
-    solar_constant::NF = 1365/4     # for testing
+    solar_constant::NF = 1365
 end
 
 Earth(SG::SpectralGrid; kwargs...) = Earth{SG.NF}(; kwargs...)

src/dynamics/virtual_temperature.jl

@@ -50,6 +50,26 @@ function virtual_temperature!(  diagn::DiagnosticVariablesLayer,
     copyto!(temp_virt_grid, temp_grid)
 end
 
+function virtual_temperature!(
+    column::ColumnVariables,
+    model::PrimitiveEquation,
+)
+    (; temp, temp_virt, humid) = column
+    μ = model.atmosphere.μ_virt_temp
+
+    @. temp_virt = temp*(1 + μ*humid)
+    return nothing
+end
+
+function virtual_temperature!(
+    column::ColumnVariables,
+    model::PrimitiveDry,
+)
+    (; temp, temp_virt) = column
+    @. temp_virt = temp             # temp = temp_virt for PrimitiveDry
+    return nothing
+end
+
 """
 $(TYPEDSIGNATURES)
 Linear virtual temperature for `model::PrimitiveDry`: Just copy over

src/models/primitive_dry.jl

@@ -76,7 +76,7 @@ Base.@kwdef mutable struct PrimitiveDryModel{
     surface_wind::SUW = SurfaceWind(spectral_grid)
     surface_heat_flux::SH = SurfaceHeatFlux(spectral_grid)
     convection::CV = DryBettsMiller(spectral_grid)
-    shortwave_radiation::SW = TransparentShortwave(spectral_grid)
+    shortwave_radiation::SW = NoShortwave(spectral_grid)
     longwave_radiation::LW = JeevanjeeRadiation(spectral_grid)
 
     # NUMERICS

src/models/primitive_wet.jl

@@ -87,7 +87,7 @@ Base.@kwdef mutable struct PrimitiveWetModel{
     surface_evaporation::EV = SurfaceEvaporation(spectral_grid)
     large_scale_condensation::LSC = ImplicitCondensation(spectral_grid)
     convection::CV = SimplifiedBettsMiller(spectral_grid)
-    shortwave_radiation::SW = TransparentShortwave(spectral_grid)
+    shortwave_radiation::SW = NoShortwave(spectral_grid)
     longwave_radiation::LW = JeevanjeeRadiation(spectral_grid)
 
     # NUMERICS

src/physics/boundary_layer.jl

@@ -102,7 +102,7 @@ Base.@kwdef struct BulkRichardsonDrag{NF} <: AbstractBoundaryLayer
     z₀::NF = 3.21e-5
 
     "Critical Richardson number for stable mixing cutoff [1]"
-    Ri_c::NF = 1
+    Ri_c::NF = 10
 
     "Maximum drag coefficient, κ²/log(zₐ/z₀) for zₐ from reference temperature"
     drag_max::Base.RefValue{NF} = Ref(zero(NF))

src/physics/column_variables.jl

@@ -35,7 +35,7 @@ function get_column!(
 )
 
     (; σ_levels_full, ln_σ_levels_full) = geometry
-    (; temp_profile) = implicit     # reference temperature on this layer
+    (; temp_profile) = implicit     # reference temperature
 
     @boundscheck C.nlev == D.nlev || throw(BoundsError)
 
@@ -54,18 +54,14 @@ function get_column!(
     C.pres[1:end-1] .= σ_levels_full.*pₛ    # pressure on every level p = σ*pₛ
     C.pres[end] = pₛ                        # last element is surface pressure pₛ
 
-    @inbounds for (k, layer) = enumerate(D.layers)   # read out prognostic variables on grid
-        C.u[k] = layer.grid_variables.u_grid[ij]
-        C.v[k] = layer.grid_variables.v_grid[ij]
-        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] 
+    @inbounds for (k, layer) = enumerate(D.layers)
+        # read out prognostic variables on grid at previous time step
+        C.u[k] = layer.grid_variables.u_grid_prev[ij]
+        C.v[k] = layer.grid_variables.v_grid_prev[ij]
+
+        # add temp reference profile back in as temp_grid_prev is anomaly
+        C.temp[k] = layer.grid_variables.temp_grid_prev[ij] + temp_profile[k]
+        C.humid[k] = layer.grid_variables.humid_grid_prev[ij] 
     end
 
     # TODO skin = surface approximation for now

src/physics/define_column.jl

@@ -26,12 +26,6 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV
     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)]
     const pres::Vector{NF} = zeros(NF, nlev+1)       # pressure [Pa]

src/physics/large_scale_condensation.jl

@@ -10,12 +10,9 @@ export ImplicitCondensation
 """
 Large scale condensation as with implicit precipitation.
 $(TYPEDFIELDS)"""
-Base.@kwdef struct ImplicitCondensation{NF<:AbstractFloat} <: AbstractCondensation
+@kwdef struct ImplicitCondensation{NF<:AbstractFloat} <: AbstractCondensation
     "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 = 10           # currently not needed? maybe too high?
 
     "Time scale in multiples of time step Δt, the larger the less immediate"
     time_scale::NF = 3
@@ -70,9 +67,7 @@ function large_scale_condensation!(
     time_stepping::AbstractTimeStepper,
 )
 
-    (; 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
+    (; pres, temp, humid) = column          # prognostic vars (from previous time step for numerical stability)
     (; temp_tend, humid_tend) = column      # tendencies to write into
     (; sat_humid) = column                  # intermediate variable, calculated in thermodynamics!
 
@@ -84,7 +79,6 @@ 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, relative_humidity_threshold) = scheme
 
     @inbounds for k in eachindex(column)
@@ -97,9 +91,8 @@ function large_scale_condensation!(
             dqsat_dT = sat_humid[k] * Lᵥ_Rᵥ/temp[k]^2       # derivative of qsat wrt to temp
             δq /= ((1 + Lᵥ_cₚ*dqsat_dT) * time_scale*Δt_sec) 
 
-            # latent heat release with maximum heating limiter for stability
-            δT = min(max_heating, -Lᵥ_cₚ * δq)
-            δq = -δT/Lᵥ_cₚ                                  # also limit drying for enthalpy conservation
+            # latent heat release for enthalpy conservation
+            δT = -Lᵥ_cₚ * δq
 
             # If there is large-scale condensation at a level higher (i.e. smaller k) than
             # the cloud-top previously diagnosed due to convection, then increase the cloud-top

src/physics/shortwave_radiation.jl

@@ -34,5 +34,8 @@ function shortwave_radiation!(
 )
     (; cos_zenith, albedo) = column
     (; solar_constant) = planet
-    column.flux_temp_upward[end] += (1-albedo) * solar_constant * cos_zenith
+
+    # reduce strength by 1/4 as this currently just hits the air temperature in lowermost
+    # layer directly and not through a skin temperature
+    column.flux_temp_upward[end] += ((1-albedo) * solar_constant * cos_zenith)/4
 end

src/physics/surface_fluxes.jl

@@ -29,8 +29,8 @@ function surface_thermodynamics!(   column::ColumnVariables,
 
     # surface value is same as lowest model level, use previous time step
     # for numerical stability
-    column.surface_temp = column.temp_prev[end]      # todo use constant POTENTIAL temperature
-    column.surface_humid = column.humid_prev[end]    # humidity at surface is the same as 
+    column.surface_temp = column.temp[end]      # todo use constant POTENTIAL temperature
+    column.surface_humid = column.humid[end]    # humidity at surface is the same as 
 
     # surface air density via virtual temperature
     (; R_dry) = model.atmosphere
@@ -44,7 +44,7 @@ function surface_thermodynamics!(   column::ColumnVariables,
     (; R_dry) = model.atmosphere
     # surface value is same as lowest model level, but use previous
     # time step for numerical stability
-    column.surface_temp = column.temp_prev[end]   # todo use constant POTENTIAL temperature
+    column.surface_temp = column.temp[end]   # todo use constant POTENTIAL temperature
     column.surface_air_density = column.pres[end]/(R_dry*column.surface_temp)
 end
 
@@ -71,9 +71,6 @@ Base.@kwdef struct SurfaceWind{NF<:AbstractFloat} <: AbstractSurfaceWind
 
     "Otherwise, Drag coefficient over sea [1]"
     drag_sea::NF = 1.8e-3
-
-    "Flux limiter to cap the max of surface momentum fluxes [kg/m/s²]"
-    max_flux::NF = 10       # currently not needed? maybe too high?
 end
 
 SurfaceWind(SG::SpectralGrid; kwargs...) = SurfaceWind{SG.NF}(; kwargs...)
@@ -85,11 +82,10 @@ function surface_wind_stress!(  column::ColumnVariables,
 
     (; land_fraction) = column
     (; f_wind, V_gust, drag_land, drag_sea) = surface_wind
-    (; max_flux) = surface_wind
 
     # SPEEDY documentation eq. 49, but use previous time step for numerical stability
-    column.surface_u = f_wind*column.u_prev[end] 
-    column.surface_v = f_wind*column.v_prev[end]
+    column.surface_u = f_wind*column.u[end] 
+    column.surface_v = f_wind*column.v[end]
     (; surface_u, surface_v) = column
 
     # SPEEDY documentation eq. 50
@@ -105,15 +101,9 @@ function surface_wind_stress!(  column::ColumnVariables,
     V₀ = column.surface_wind_speed
     drag = land_fraction*drag_land + (1-land_fraction)*drag_sea
 
-    # add flux limiter to avoid heavy drag in (initial) shock
-    u_flux = ρ*drag*V₀*surface_u
-    v_flux = ρ*drag*V₀*surface_v
-    column.flux_u_upward[end] -= clamp(u_flux, -max_flux, max_flux)
-    column.flux_v_upward[end] -= clamp(v_flux, -max_flux, max_flux)
-
     # SPEEDY documentation eq. 52, 53, accumulate fluxes with +=
-    # column.flux_u_upward[end] -= ρ*drag*V₀*surface_u
-    # column.flux_v_upward[end] -= ρ*drag*V₀*surface_v
+    column.flux_u_upward[end] -= ρ*drag*V₀*surface_u
+    column.flux_v_upward[end] -= ρ*drag*V₀*surface_v
 
     return nothing
 end
@@ -136,9 +126,6 @@ Base.@kwdef struct SurfaceHeatFlux{NF<:AbstractFloat} <: AbstractSurfaceHeatFlux
 
     "Otherwise, use the following drag coefficient for heat fluxes over sea"
     heat_exchange_sea::NF = 0.9e-3
-
-    "Flux limiter for surface heat fluxes [W/m²]"
-    max_flux::NF = 100          # currently not needed? maybe too high?
 end
 
 SurfaceHeatFlux(SG::SpectralGrid; kwargs...) = SurfaceHeatFlux{SG.NF}(; kwargs...)
@@ -150,7 +137,7 @@ function surface_heat_flux!(
     model::PrimitiveEquation,
 )   
     cₚ = model.atmosphere.heat_capacity
-    (; heat_exchange_land, heat_exchange_sea, max_flux) = heat_flux
+    (; heat_exchange_land, heat_exchange_sea) = heat_flux
 
     ρ = column.surface_air_density
     V₀ = column.surface_wind_speed
@@ -168,10 +155,6 @@ function surface_heat_flux!(
     flux_land = ρ*drag_land*V₀*cₚ*(T_skin_land - T)
     flux_sea  = ρ*drag_sea*V₀*cₚ*(T_skin_sea  - T)
 
-    # flux limiter
-    flux_land = clamp(flux_land, -max_flux, max_flux)
-    flux_sea = clamp(flux_sea, -max_flux, max_flux)
-
     # mix fluxes for fractional land-sea mask
     land_available = isfinite(T_skin_land)
     sea_available = isfinite(T_skin_sea)

src/physics/thermodynamics.jl

@@ -121,6 +121,7 @@ Calculate geopotentiala and dry static energy for the primitive equation model."
 function get_thermodynamics!(column::ColumnVariables, model::PrimitiveEquation)
     geopotential!(column.geopot, column.temp, model.geopotential, column.surface_geopotential)
     dry_static_energy!(column, model.atmosphere)
+    virtual_temperature!(column, model)
 end
 
 """
@@ -148,13 +149,11 @@ function saturation_humidity!(
     column::ColumnVariables,
     clausius_clapeyron::AbstractClausiusClapeyron,
 )
-    (; 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
+    (; sat_humid, pres, temp) = column
 
-    for k in eachlayer(column)
+    @inbounds for k in eachlayer(column)
         sat_humid[k] = saturation_humidity(temp[k], pres[k], clausius_clapeyron)
     end
 end
@@ -175,7 +174,7 @@ function moist_static_energy!(
     (; sat_moist_static_energy, moist_static_energy, dry_static_energy) = column
     (; humid, sat_humid) = column
 
-    for k in eachlayer(column)
+    @inbounds for k in eachlayer(column)
         moist_static_energy[k] = dry_static_energy[k] + Lᵥ * humid[k]
         sat_moist_static_energy[k] = dry_static_energy[k] + Lᵥ * sat_humid[k]
     end

src/physics/vertical_diffusion.jl

@@ -16,7 +16,7 @@ initialize!(::NoVerticalDiffusion,::PrimitiveEquation) = nothing
 vertical_diffusion!(::ColumnVariables, ::NoVerticalDiffusion, ::PrimitiveEquation) = nothing
 
 export BulkRichardsonDiffusion
-Base.@kwdef struct BulkRichardsonDiffusion{NF} <: AbstractVerticalDiffusion
+@kwdef struct BulkRichardsonDiffusion{NF} <: AbstractVerticalDiffusion
     nlev::Int
 
     "[OPTION] von Kármán constant [1]"
@@ -26,7 +26,7 @@ Base.@kwdef struct BulkRichardsonDiffusion{NF} <: AbstractVerticalDiffusion
     z₀::NF = 3.21e-5
 
     "[OPTION] Critical Richardson number for stable mixing cutoff [1]"
-    Ri_c::NF = 1
+    Ri_c::NF = 10
 
     "[OPTION] Fraction of surface boundary layer"
     fb::NF = 0.1
@@ -216,7 +216,7 @@ function bulk_richardson!(
 )
     cₚ = atmosphere.heat_capacity
     (; u, v, geopot, temp_virt, nlev) = column
-    bulk_richardson = column.a      # reuse work array
+    bulk_richardson = column.d      # reuse work array
 
     # surface layer
     V² = u[nlev]^2 + v[nlev]^2