Slab ocean and net surface and top-of-atmosphere fluxes (#613)

* Sum of surface fluxes as diagnostic variable

* Slab ocean

* update changelog

* ocean mixed-layer heat capacity like in Frierson 2006

* transparent shortwave with fluxes

* export SlabOcean

* Jeevanjee extended to surface

* surface fluxes algorithm simplified

* land/ocean interfaces generalised

* heat capacity via mixed layer depth, specific heat capacity and density

* reset tendencies depending on model

* fix initialize! for ocean variables

* outgoing longwave radiation

* test ocean models

* add outgoing shortwave radiation netCDF variable

* typo osr not olr

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Unreleased
 
+- Slab ocean and net surface fluxes [#613](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/613)
 - added compat entries for FiniteDifferences.jl and Enzyme.jl [#620](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/620) [#622](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/622)
 - update to GPUArrays v11, JLArrays v0.2 and remove CUDA from tests [#590](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/590)
 - Power spectrum for n-dim LowerTriangularArrays [#618](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/618)

src/dynamics/diagnostic_variables.jl

@@ -262,6 +262,18 @@ $(TYPEDFIELDS)"""
     "Availability of soil moisture to evaporation [1]"
     soil_moisture_availability::GridVariable2D = zeros(GridVariable2D, nlat_half)
 
+    "Surface flux of heat [W/m^2]"
+    surface_flux_heat::GridVariable2D = zeros(GridVariable2D, nlat_half)
+
+    "Surface flux of humidity [?]"
+    surface_flux_humid::GridVariable2D = zeros(GridVariable2D, nlat_half)
+
+    "Outgoing shortwave radiation [W/m^2]"
+    outgoing_shortwave_radiation::GridVariable2D = zeros(GridVariable2D, nlat_half)
+
+    "Outgoing longwave radiation [W/m^2]"
+    outgoing_longwave_radiation::GridVariable2D = zeros(GridVariable2D, nlat_half)
+
     "Cosine of solar zenith angle [1]"
     cos_zenith::GridVariable2D = zeros(GridVariable2D, nlat_half)           
 end

src/dynamics/time_integration.jl

@@ -300,16 +300,17 @@ function timestep!(
     (; time) = progn.clock                           # current time
 
     # set the tendencies back to zero for accumulation
-    fill!(diagn.tendencies, 0, PrimitiveWet)
+    fill!(diagn.tendencies, 0, typeof(model))
 
     if model.physics                                # switch on/off all physics parameterizations
+        # calculate all parameterizations
+        parameterization_tendencies!(diagn, progn, time, model)
+
         # time step ocean (temperature and TODO sea ice) and land (temperature and soil moisture)
+        # with fluxes from parameterizations
         ocean_timestep!(progn, diagn, model)
         land_timestep!(progn, diagn, model)
         soil_moisture_availability!(diagn, progn, model)
-
-        # calculate all parameterizations
-        parameterization_tendencies!(diagn, progn, time, model)
     end
 
     if model.dynamics                                           # switch on/off all dynamics

src/models/primitive_dry.jl

@@ -150,8 +150,9 @@ function initialize!(model::PrimitiveDry; time::DateTime = DEFAULT_DATE)
     initialize!(prognostic_variables.particles, model)
 
     # initialize ocean and land
-    initialize!(prognostic_variables.ocean, time, model)
+    initialize!(prognostic_variables.ocean, prognostic_variables, diagnostic_variables, model)
     initialize!(prognostic_variables.land,  prognostic_variables, diagnostic_variables, model)
 
+    # pack prognostic, diagnostic variables and model into a simulation
     return Simulation(prognostic_variables, diagnostic_variables, model)
 end

src/models/primitive_wet.jl

@@ -154,8 +154,8 @@ function initialize!(model::PrimitiveWet; time::DateTime = DEFAULT_DATE)
 
     # initial conditions
     prognostic_variables = PrognosticVariables(spectral_grid, model)
-    initialize!(prognostic_variables,model.initial_conditions, model)
-    (;clock) = prognostic_variables
+    initialize!(prognostic_variables, model.initial_conditions, model)
+    (; clock) = prognostic_variables
     clock.time = time       # set the current time
     clock.start = time      # and store the start time
 
@@ -166,8 +166,9 @@ function initialize!(model::PrimitiveWet; time::DateTime = DEFAULT_DATE)
     initialize!(prognostic_variables.particles, model)
 
     # initialize ocean and land
-    initialize!(prognostic_variables.ocean, time, model)
+    initialize!(prognostic_variables.ocean, prognostic_variables, diagnostic_variables, model)
     initialize!(prognostic_variables.land,  prognostic_variables, diagnostic_variables, model)
 
+    # pack prognostic, diagnostic variables and model into a simulation
     return Simulation(prognostic_variables, diagnostic_variables, model)
 end

src/output/netcdf_output.jl

@@ -895,6 +895,136 @@ function output!(
     return nothing
 end
 
+## SURFACE FLUXES
+
+"""Defines netCDF output for a specific variables, see `VorticityOutput` for details.
+Fields are $(TYPEDFIELDS)"""
+@kwdef mutable struct SurfaceFluxHeatOutput <: AbstractOutputVariable
+    name::String = "surface_flux_heat"
+    unit::String = "W/m^2"
+    long_name::String = "Surface heat fluxes (positive down)"
+    dims_xyzt::NTuple{4, Bool} = (true, true, false, true)
+    missing_value::Float64 = NaN
+    compression_level::Int = 3
+    shuffle::Bool = true
+    keepbits::Int = 7
+end
+
+"""$(TYPEDSIGNATURES)
+`output!` method for `variable`, see `output!(::NetCDFOutput, ::VorticityOutput, ...)` for details."""
+function output!(
+    output::NetCDFOutput,
+    variable::SurfaceFluxHeatOutput,
+    progn::PrognosticVariables,
+    diagn::DiagnosticVariables,
+    model::AbstractModel,
+)
+    flux = output.grid2D
+    (; surface_flux_heat) = diagn.physics
+    RingGrids.interpolate!(flux, surface_flux_heat, output.interpolator)
+
+    round!(flux, variable.keepbits)
+    i = output.output_counter   # output time step to write
+    output.netcdf_file[variable.name][:, :, i] = flux
+    return nothing
+end
+
+"""Defines netCDF output for a specific variables, see `VorticityOutput` for details.
+Fields are $(TYPEDFIELDS)"""
+@kwdef mutable struct OutgoingLongwaveRadiationOutput <: AbstractOutputVariable
+    name::String = "olr"
+    unit::String = "W/m^2"
+    long_name::String = "Outgoing longwave radiation"
+    dims_xyzt::NTuple{4, Bool} = (true, true, false, true)
+    missing_value::Float64 = NaN
+    compression_level::Int = 3
+    shuffle::Bool = true
+    keepbits::Int = 7
+end
+
+"""$(TYPEDSIGNATURES)
+`output!` method for `variable`, see `output!(::NetCDFOutput, ::VorticityOutput, ...)` for details."""
+function output!(
+    output::NetCDFOutput,
+    variable::OutgoingLongwaveRadiationOutput,
+    progn::PrognosticVariables,
+    diagn::DiagnosticVariables,
+    model::AbstractModel,
+)
+    olr = output.grid2D
+    (; outgoing_longwave_radiation) = diagn.physics
+    RingGrids.interpolate!(olr, outgoing_longwave_radiation, output.interpolator)
+
+    round!(olr, variable.keepbits)
+    i = output.output_counter   # output time step to write
+    output.netcdf_file[variable.name][:, :, i] = olr
+    return nothing
+end
+
+"""Defines netCDF output for a specific variables, see `VorticityOutput` for details.
+Fields are $(TYPEDFIELDS)"""
+@kwdef mutable struct OutgoingShortwaveRadiationOutput <: AbstractOutputVariable
+    name::String = "osr"
+    unit::String = "W/m^2"
+    long_name::String = "Outgoing shortwave radiation"
+    dims_xyzt::NTuple{4, Bool} = (true, true, false, true)
+    missing_value::Float64 = NaN
+    compression_level::Int = 3
+    shuffle::Bool = true
+    keepbits::Int = 7
+end
+
+"""$(TYPEDSIGNATURES)
+`output!` method for `variable`, see `output!(::NetCDFOutput, ::VorticityOutput, ...)` for details."""
+function output!(
+    output::NetCDFOutput,
+    variable::OutgoingShortwaveRadiationOutput,
+    progn::PrognosticVariables,
+    diagn::DiagnosticVariables,
+    model::AbstractModel,
+)
+    osr = output.grid2D
+    (; outgoing_shortwave_radiation) = diagn.physics
+    RingGrids.interpolate!(osr, outgoing_shortwave_radiation, output.interpolator)
+
+    round!(osr, variable.keepbits)
+    i = output.output_counter   # output time step to write
+    output.netcdf_file[variable.name][:, :, i] = osr
+    return nothing
+end
+
+"""Defines netCDF output for a specific variables, see `VorticityOutput` for details.
+Fields are $(TYPEDFIELDS)"""
+@kwdef mutable struct SurfaceFluxHumidOutput <: AbstractOutputVariable
+    name::String = "surface_flux_humid"
+    unit::String = "kg/s/m^2"
+    long_name::String = "Surface humidity fluxes (positive down)"
+    dims_xyzt::NTuple{4, Bool} = (true, true, false, true)
+    missing_value::Float64 = NaN
+    compression_level::Int = 3
+    shuffle::Bool = true
+    keepbits::Int = 7
+end
+
+"""$(TYPEDSIGNATURES)
+`output!` method for `variable`, see `output!(::NetCDFOutput, ::VorticityOutput, ...)` for details."""
+function output!(
+    output::NetCDFOutput,
+    variable::SurfaceFluxHumidOutput,
+    progn::PrognosticVariables,
+    diagn::DiagnosticVariables,
+    model::AbstractModel,
+)
+    flux = output.grid2D
+    (; surface_flux_humid) = diagn.physics
+    RingGrids.interpolate!(flux, surface_flux_humid, output.interpolator)
+
+    round!(flux, variable.keepbits)
+    i = output.output_counter   # output time step to write
+    output.netcdf_file[variable.name][:, :, i] = flux
+    return nothing
+end
+
 ## RANDOM PATTERN
 
 """Defines netCDF output for a specific variables, see `VorticityOutput` for details.
@@ -929,6 +1059,42 @@ function output!(
     return nothing
 end
 
+## RANDOM PATTERN
+
+"""Defines netCDF output for a specific variables, see `VorticityOutput` for details.
+Fields are $(TYPEDFIELDS)"""
+@kwdef mutable struct SeaSurfaceTemperatureOutput <: AbstractOutputVariable
+    name::String = "sst"
+    unit::String = "degC"
+    long_name::String = "sea surface temperature"
+    dims_xyzt::NTuple{4, Bool} = (true, true, false, true)
+    missing_value::Float64 = NaN
+    compression_level::Int = 3
+    shuffle::Bool = true
+    keepbits::Int = 10
+end
+
+"""$(TYPEDSIGNATURES)
+`output!` method for `variable`, see `output!(::NetCDFOutput, ::VorticityOutput, ...)` for details."""
+function output!(
+    output::NetCDFOutput,
+    variable::SeaSurfaceTemperatureOutput,
+    progn::PrognosticVariables,
+    diagn::DiagnosticVariables,
+    model::AbstractModel,
+)
+    sst = output.grid2D
+    sst_grid = progn.ocean.sea_surface_temperature
+    RingGrids.interpolate!(sst, sst_grid, output.interpolator)
+
+    sst .-= 273.15  # convert to ˚C
+
+    round!(sst, variable.keepbits)
+    i = output.output_counter   # output time step to write
+    output.netcdf_file[variable.name][:, :, i] = sst
+    return nothing
+end
+
 """
 $(TYPEDSIGNATURES)
 Checks existing `run_????` folders in `path` to determine a 4-digit `id` number

src/physics/column_variables.jl

@@ -132,17 +132,27 @@ function write_column_tendencies!(
         diagn.tendencies.humid_tend_grid[ij, k] = column.humid_tend[k]
     end
 
-    # accumulate (set back to zero when netcdf output)
+    # accumulate precipitation [m]
     diagn.physics.precip_large_scale[ij] += column.precip_large_scale
     diagn.physics.precip_convection[ij] += column.precip_convection
 
-    # Output cloud top in height [m] from geopotential height divided by gravity,
-    # but NaN for no clouds
+    # Cloud top in height [m] from geopotential height divided by gravity, 0 for no clouds
     diagn.physics.cloud_top[ij] = column.cloud_top == nlayers+1 ? 0 : column.geopot[column.cloud_top]
     diagn.physics.cloud_top[ij] /= planet.gravity
 
     # just use layer index 1 (top) to nlayers (surface) for analysis, but 0 for no clouds
     # diagn.physics.cloud_top[ij] = column.cloud_top == nlayers+1 ? 0 : column.cloud_top
+
+    # net surface fluxes of humidity and temperature, defined as positive downward (i.e. into ocean/land)
+    diagn.physics.surface_flux_heat[ij] = column.flux_temp_downward[nlayers+1] -
+                                            column.flux_temp_upward[nlayers+1]
+    diagn.physics.surface_flux_humid[ij] = column.flux_humid_downward[nlayers+1] -
+                                            column.flux_humid_upward[nlayers+1]
+
+    # outgoing radiation
+    diagn.physics.outgoing_longwave_radiation[ij] = column.outgoing_longwave_radiation
+    diagn.physics.outgoing_shortwave_radiation[ij] = column.outgoing_shortwave_radiation
+
     return nothing
 end
 

src/physics/define_column.jl

@@ -84,9 +84,12 @@ $(TYPEDFIELDS)"""
     precip_large_scale::NF = 0              # precipitation due to large-scale condensation [m]
 
     # RADIATION
-    cos_zenith::NF = 0                                      # cosine of solar zenith angle
-    albedo::NF = 0                                          # surface albedo
+    cos_zenith::NF = 0                      # cosine of solar zenith angle
+    albedo::NF = 0                          # surface albedo
 
+    outgoing_longwave_radiation::NF = 0     # OLR [W/m^2]
+    outgoing_shortwave_radiation::NF = 0    # same for shortwave reflection [W/m^2]
+
     # optical depth of the atmosphere, on half levels, for shortwave and longwave radiation
     const optical_depth_shortwave::MatrixType = zeros(NF, nlayers+1, nbands_shortwave)
     const optical_depth_longwave::MatrixType = zeros(NF, nlayers+1, nbands_longwave)

src/physics/land.jl

@@ -20,7 +20,6 @@ function land_timestep!(
     diagn::DiagnosticVariables,
     model::PrimitiveEquation,
 )
-    # the time step is dictated by the land "model" 
     land_timestep!(progn, diagn, model.land, model)
     model isa PrimitiveWet && soil_timestep!(progn, diagn, model.soil, model)
 end

src/physics/longwave_radiation.jl

@@ -117,7 +117,29 @@ function longwave_radiation!(
     (; α, time_scale) = scheme
     Tₜ = scheme.temp_tropopause
 
-    Fₖ::NF = 0                      # flux into lowermost layer from below
+    (; skin_temperature_sea, skin_temperature_land, land_fraction) = column
+
+    # mix radiative fluxes for fractional land-sea mask if both available 
+    land_available = isfinite(skin_temperature_land)
+    sea_available = isfinite(skin_temperature_sea)
+
+    # land fraction-weighted average of surface temperatures from land/sea
+    local surface_temperature::NF = 0
+
+    if land_available && sea_available
+        surface_temperature = land_fraction*skin_temperature_land + (1-land_fraction)*skin_temperature_sea
+    elseif land_available
+        surface_temperature = skin_temperature_land
+    elseif sea_available
+        surface_temperature = skin_temperature_sea
+    else    # fallback: use surface air temperature to have zero flux
+        surface_temperature = T[end]
+    end
+
+    # extension to Jeevanjee: Include temperature flux between surface and lowermost air temperature
+    local Fₖ::NF    # flux into lowermost layer from surface land/sea below
+    Fₖ = (T[end] - surface_temperature) * α * (Tₜ - surface_temperature)
+    F[end] += Fₖ
 
     # integrate from surface up, skipping surface (k=nlayers+1) and top-of-atmosphere flux (k=1)
     @inbounds for k in nlayers:-1:2
@@ -128,6 +150,10 @@ function longwave_radiation!(
 
     # Relax the uppermost level towards prescribed "tropopause temperature"
     temp_tend[1] += (Tₜ - T[1])/time_scale.value
+
+    # for diagnostic, use Fₖ as the outgoing longwave radiation although it's technically into the
+    # uppermost layer from below (not out of it)
+    column.outgoing_longwave_radiation = Fₖ
 end
 
 # dummy one band radiation for now