Merge pull request #502 from SpeedyWeather/mk/jeevanjee

Jeevanjee longwave radiation

src/dynamics/diagnostic_variables.jl

@@ -129,18 +129,19 @@ Base.@kwdef struct SurfaceVariables{NF<:AbstractFloat, Grid<:AbstractGrid{NF}}
     pres_tend::LTM{Complex{NF}} = zeros(LTM{Complex{NF}}, trunc+2, trunc+1)
     pres_tend_grid::Grid = zeros(Grid, nlat_half)
 
-    ∇lnp_x::Grid = zeros(Grid, nlat_half)        # zonal gradient of log surf pressure
-    ∇lnp_y::Grid = zeros(Grid, nlat_half)        # meridional gradient of log surf pres
+    ∇lnp_x::Grid = zeros(Grid, nlat_half)           # zonal gradient of log surf pressure
+    ∇lnp_y::Grid = zeros(Grid, nlat_half)           # meridional gradient of log surf pres
 
-    u_mean_grid::Grid = zeros(Grid, nlat_half)   # vertical average of: zonal velocity *coslat
-    v_mean_grid::Grid = zeros(Grid, nlat_half)   # meridional velocity *coslat
-    div_mean_grid::Grid = zeros(Grid, nlat_half) # divergence
+    u_mean_grid::Grid = zeros(Grid, nlat_half)      # vertical average of: zonal velocity *coslat
+    v_mean_grid::Grid = zeros(Grid, nlat_half)      # meridional velocity *coslat
+    div_mean_grid::Grid = zeros(Grid, nlat_half)    # divergence
     div_mean::LTM{Complex{NF}} = zeros(LTM{Complex{NF}}, trunc+2, trunc+1)    # divergence (in spectral though)
 
-    precip_large_scale::Grid = zeros(Grid, nlat_half)    # large scale precipitation (for output)
-    precip_convection::Grid = zeros(Grid, nlat_half)     # convective precipitation (for output)
-    cloud_top::Grid = zeros(Grid, nlat_half)             # cloud top [hPa]
+    precip_large_scale::Grid = zeros(Grid, nlat_half)   # large scale precipitation (for output)
+    precip_convection::Grid = zeros(Grid, nlat_half)    # convective precipitation (for output)
+    cloud_top::Grid = zeros(Grid, nlat_half)            # cloud top [hPa]
     soil_moisture_availability::Grid = zeros(Grid, nlat_half)
+    cos_zenith::Grid = zeros(Grid, nlat_half)           # cosine of solar zenith angle
 end
 
 # generator function based on a SpectralGrid

src/dynamics/orography.jl

@@ -78,7 +78,7 @@ function initialize!(   orog::ZonalRidge,
 
     ηᵥ = (1-η₀)*π/2                     # ηᵥ-coordinate of the surface [1]
     A = u₀*cos(ηᵥ)^(3/2)                # amplitude [m/s]
-    RΩ = radius*rotation                # [m/s]
+    RΩ = radius*rotation                # [m/s]
     g⁻¹ = inv(gravity)                  # inverse gravity [s²/m]
 
     for ij in eachindex(φ, orography)
@@ -114,7 +114,7 @@ Base.@kwdef mutable struct EarthOrography{NF<:AbstractFloat, Grid<:AbstractGrid{
     scale::Float64 = 1
 
     "smooth the orography field?"
-    smoothing::Bool = true
+    smoothing::Bool = false
 
     "power of Laplacian for smoothing"
     smoothing_power::Float64 = 1.0

src/dynamics/planet.jl

@@ -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
+    solar_constant::NF = 1365/4     # for testing
 end
 
 Earth(SG::SpectralGrid; kwargs...) = Earth{SG.NF}(; kwargs...)

src/dynamics/time_integration.jl

@@ -395,7 +395,7 @@ function time_stepping!(
 
     # FIRST TIMESTEPS: EULER FORWARD THEN 1x LEAPFROG
     # considered part of the model initialisation
-    first_timesteps!(progn,diagn, model, output)
+    first_timesteps!(progn, diagn, model, output)
 
     # only now initialise feedback for benchmark accuracy
     initialize!(feedback, clock, model)
@@ -405,7 +405,7 @@ function time_stepping!(
         timestep!(progn, diagn, 2Δt, model) # calculate tendencies and leapfrog forward
         timestep!(clock, Δt_millisec)       # time of lf=2 and diagn after timestep!
 
-        progress!(feedback, progn)           # updates the progress meter bar
+        progress!(feedback, progn)          # updates the progress meter bar
         write_output!(output, clock.time, diagn)
         callback!(model.callbacks, progn, diagn, model)
     end

src/models/abstract_models.jl

@@ -30,7 +30,7 @@ model_class(::Type{<:PrimitiveDry}) = PrimitiveDry
 model_class(::Type{<:PrimitiveWet}) = PrimitiveWet
 model_class(model::ModelSetup) = model_class(typeof(model))
 
-# model type is the parameter-free type of a model
+# model type is the parameter-free type of a model
 # TODO what happens if we have several concrete types under each abstract type?
 model_type(::Type{<:Barotropic}) = BarotropicModel
 model_type(::Type{<:ShallowWater}) = ShallowWaterModel
@@ -44,7 +44,7 @@ function Base.show(io::IO, M::ModelSetup)
     n = length(properties)
     for (i, key) in enumerate(properties)
         val = getfield(M, key)
-        s = i == n ? "└" : "├"  # choose ending └ for last property
+        s = i == n ? "└" : "├"  # choose ending └ for last property
         p = i == n ? print : println
         p(io, "$s $key: $(typeof(val))")
     end

src/models/barotropic.jl

@@ -79,8 +79,8 @@ function initialize!(model::Barotropic; time::DateTime = DEFAULT_DATE)
     # initial conditions
     prognostic_variables = PrognosticVariables(spectral_grid, model)
     initialize!(prognostic_variables, model.initial_conditions, model)
-    prognostic_variables.clock.time = time       # set the current time
-    prognostic_variables.clock.start = time      # and store the start time
+    prognostic_variables.clock.time = time       # set the current time
+    prognostic_variables.clock.start = time      # and store the start time
 
     # particle advection
     initialize!(model.particle_advection, model)

src/models/primitive_dry.jl

@@ -1,9 +1,14 @@
 export PrimitiveDryModel
 
 """
-$(SIGNATURES)
-The PrimitiveDryModel struct holds all other structs that contain precalculated constants,
-whether scalars or arrays that do not change throughout model integration.
+The PrimitiveDryModel contains all model components (themselves structs) needed for the
+simulation of the primitive equations without humidity. To be constructed like
+
+    model = PrimitiveDryModel(; spectral_grid, kwargs...)
+
+with `spectral_grid::SpectralGrid` used to initalize all non-default components
+passed on as keyword arguments, e.g. `planet=Earth(spectral_grid)`. Fields, representing
+model components, are
 $(TYPEDFIELDS)"""
 Base.@kwdef mutable struct PrimitiveDryModel{
     NF<:AbstractFloat,
@@ -20,6 +25,7 @@ Base.@kwdef mutable struct PrimitiveDryModel{
     OC<:AbstractOcean,
     LA<:AbstractLand,
     ZE<:AbstractZenith,
+    AL<:AbstractAlbedo,
     BL<:AbstractBoundaryLayer,
     TR<:AbstractTemperatureRelaxation,
     VD<:AbstractVerticalDiffusion,
@@ -58,6 +64,7 @@ Base.@kwdef mutable struct PrimitiveDryModel{
     ocean::OC = SeasonalOceanClimatology(spectral_grid)
     land::LA = SeasonalLandTemperature(spectral_grid)
     solar_zenith::ZE = WhichZenith(spectral_grid, planet)
+    albedo::AL = AlbedoClimatology(spectral_grid)
 
     # PHYSICS/PARAMETERIZATIONS
     physics::Bool = true
@@ -68,8 +75,8 @@ Base.@kwdef mutable struct PrimitiveDryModel{
     surface_wind::SUW = SurfaceWind(spectral_grid)
     surface_heat_flux::SH = SurfaceSensibleHeat(spectral_grid)
     convection::CV = DryBettsMiller(spectral_grid)
-    shortwave_radiation::SW = NoShortwave(spectral_grid)
-    longwave_radiation::LW = NoLongwave(spectral_grid)
+    shortwave_radiation::SW = TransparentShortwave(spectral_grid)
+    longwave_radiation::LW = JeevanjeeRadiation(spectral_grid)
 
     # NUMERICS
     device_setup::DS = DeviceSetup(CPUDevice())
@@ -111,6 +118,7 @@ function initialize!(model::PrimitiveDry; time::DateTime = DEFAULT_DATE)
     initialize!(model.ocean, model)
     initialize!(model.land, model)
     initialize!(model.solar_zenith, time, model)
+    initialize!(model.albedo, model)
 
     # parameterizations
     initialize!(model.boundary_layer_drag, model)

src/models/primitive_wet.jl

@@ -1,9 +1,14 @@
 export PrimitiveWetModel
 
 """
-$(SIGNATURES)
-The PrimitiveDryModel struct holds all other structs that contain precalculated constants,
-whether scalars or arrays that do not change throughout model integration.
+The PrimitiveWetModel contains all model components (themselves structs) needed for the
+simulation of the primitive equations with humidity. To be constructed like
+
+    model = PrimitiveWetModel(; spectral_grid, kwargs...)
+
+with `spectral_grid::SpectralGrid` used to initalize all non-default components
+passed on as keyword arguments, e.g. `planet=Earth(spectral_grid)`. Fields, representing
+model components, are
 $(TYPEDFIELDS)"""
 Base.@kwdef mutable struct PrimitiveWetModel{
     NF<:AbstractFloat,
@@ -20,6 +25,7 @@ Base.@kwdef mutable struct PrimitiveWetModel{
     OC<:AbstractOcean,
     LA<:AbstractLand,
     ZE<:AbstractZenith,
+    AL<:AbstractAlbedo,
     SO<:AbstractSoil,
     VG<:AbstractVegetation,
     CC<:AbstractClausiusClapeyron,
@@ -64,6 +70,7 @@ Base.@kwdef mutable struct PrimitiveWetModel{
     ocean::OC = SeasonalOceanClimatology(spectral_grid)
     land::LA = SeasonalLandTemperature(spectral_grid)
     solar_zenith::ZE = WhichZenith(spectral_grid, planet)
+    albedo::AL = AlbedoClimatology(spectral_grid)
     soil::SO = SeasonalSoilMoisture(spectral_grid)
     vegetation::VG = VegetationClimatology(spectral_grid)
 
@@ -79,8 +86,8 @@ Base.@kwdef mutable struct PrimitiveWetModel{
     evaporation::EV = SurfaceEvaporation(spectral_grid)
     large_scale_condensation::LSC = ImplicitCondensation(spectral_grid)
     convection::CV = SimplifiedBettsMiller(spectral_grid)
-    shortwave_radiation::SW = NoShortwave(spectral_grid)
-    longwave_radiation::LW = UniformCooling(spectral_grid)
+    shortwave_radiation::SW = TransparentShortwave(spectral_grid)
+    longwave_radiation::LW = JeevanjeeRadiation(spectral_grid)
 
     # NUMERICS
     device_setup::DS = DeviceSetup(CPUDevice())
@@ -125,6 +132,7 @@ function initialize!(model::PrimitiveWet; time::DateTime = DEFAULT_DATE)
     initialize!(model.soil, model)
     initialize!(model.vegetation, model)
     initialize!(model.solar_zenith, time, model)
+    initialize!(model.albedo, model)
 
     # parameterizations
     initialize!(model.boundary_layer_drag, model)

src/models/shallow_water.jl

@@ -82,8 +82,8 @@ function initialize!(model::ShallowWater; time::DateTime = DEFAULT_DATE)
     # initial conditions
     prognostic_variables = PrognosticVariables(spectral_grid, model)
     initialize!(prognostic_variables, model.initial_conditions, model)
-    prognostic_variables.clock.time = time       # set the current time
-    prognostic_variables.clock.start = time      # and store the start time
+    prognostic_variables.clock.time = time       # set the current time
+    prognostic_variables.clock.start = time      # and store the start time
 
     # particle advection
     initialize!(model.particle_advection, model)

src/models/tree.jl

@@ -4,7 +4,8 @@ export tree
 $(TYPEDSIGNATURES)
 Create a tree of fields inside a Simulation instance and fields within these fields
 as long as they are defined within the modules argument (default SpeedyWeather).
-Other keyword arguments are `max_level::Integer=10`, `with_types::Bool=false`."""
+Other keyword arguments are `max_level::Integer=10`, `with_types::Bool=false`
+or `with_size::Bool=false."""
 function tree(
     S::Simulation{M};
     modules=SpeedyWeather,
@@ -21,7 +22,8 @@ end
 $(TYPEDSIGNATURES)
 Create a tree of fields inside a model and fields within these fields
 as long as they are defined within the modules argument (default SpeedyWeather).
-Other keyword arguments are `max_level::Integer=10`, `with_types::Bool=false`."""
+Other keyword arguments are `max_level::Integer=10`, `with_types::Bool=false`
+or `with_size::Bool=false."""
 function tree(
     M::ModelSetup;
     modules=SpeedyWeather,
@@ -38,7 +40,8 @@ end
 $(TYPEDSIGNATURES)
 Create a tree of fields inside S and fields within these fields
 as long as they are defined within the modules argument (default SpeedyWeather).
-Other keyword arguments are `max_level::Integer=10`, `with_types::Bool=false`."""
+Other keyword arguments are `max_level::Integer=10`, `with_types::Bool=false`
+or `with_size::Bool=false."""
 function tree(S; modules=SpeedyWeather, with_size::Bool = false, kwargs...)
     s = "$(typeof(S))"
     s = ~with_size ? s : s*" ($(prettymemory(Base.summarysize(S))))"

src/output/feedback.jl

@@ -72,7 +72,7 @@ $(TYPEDSIGNATURES)
 Initializes the a `Feedback` struct."""
 function initialize!(feedback::Feedback, clock::Clock, model::ModelSetup)
 
-    # set to false to recheck for NaRs
+    # set to false to recheck for NaRs
     feedback.nars_detected = false
 
     # hack: redefine element in global constant dt_in_sec
@@ -158,17 +158,12 @@ function nar_detection!(feedback::Feedback, progn::PrognosticVariables)
 
     feedback.nars_detected && return nothing    # escape immediately if nans already detected
     i = feedback.progress_meter.counter         # time step
-    nars_detected_here = false
-    (; vor ) = progn.layers[end].timesteps[2] # only check for surface vorticity
-
-    if ~nars_detected_here
-        nars_vor = ~isfinite(vor[1])    # just check first mode
-                                        # spectral transform propagates NaRs globally anyway
-        nars_vor && @warn "NaN or Inf detected at time step $i"
-        nars_detected_here |= nars_vor
-    end
+    (; vor ) = progn.layers[end].timesteps[2]   # only check for surface vorticity
 
-    feedback.nars_detected |= nars_detected_here
+    # just check first harmonic, spectral transform propagates NaRs globally anyway
+    nars_detected_here = ~isfinite(vor[1])
+    nars_detected_here && @warn "NaN or Inf detected at time step $i"
+    feedback.nars_detected = nars_detected_here
 end
 
 """

src/output/plot.jl

@@ -1,5 +1,3 @@
 # to avoid ambiguity with exporting plot
 plot(L::LowerTriangularMatrix{T}; kwargs...) where T = LowerTriangularMatrices.plot(L; kwargs...)
-plot(G::AbstractGrid{T}; kwargs...) where T = RingGrids.plot(G; kwargs...)
-
-;
+plot(G::AbstractGrid{T}; kwargs...) where T = RingGrids.plot(G; kwargs...)

src/output/schedule.jl

@@ -51,7 +51,7 @@ happen only once if they coincide on a given time step."""
 function initialize!(scheduler::Schedule, clock::Clock)
     schedule = falses(clock.n_timesteps)    # initialize schedule as BitVector
 
-    # PERIODIC SCHEDULE, always AFTER scheduler.every time period has passed
+    # PERIODIC SCHEDULE, always AFTER scheduler.every time period has passed
     if scheduler.every.value < typemax(Int)
         every_n_timesteps = max(1,round(Int, scheduler.every/clock.Δt))
         schedule[every_n_timesteps:every_n_timesteps:end] .= true

src/physics/albedo.jl

@@ -1 +1,61 @@
-abstract type AbstractAlbedo end
+abstract type AbstractAlbedo <: AbstractModelComponent end
+
+## GLOBAL CONSTANT ALBEDO
+
+export GlobalConstantAlbedo
+Base.@kwdef struct GlobalConstantAlbedo{NF,Grid} <: AbstractAlbedo
+    nlat_half::Int
+    α::NF = 0.8
+    albedo::Grid = zeros(Grid, nlat_half)
+end
+
+function GlobalConstantAlbedo(SG::SpectralGrid; kwargs...)
+    (;NF, Grid, nlat_half) = SG
+    GlobalConstantAlbedo{NF, Grid{NF}}(; nlat_half, kwargs...)
+end
+
+function initialize!(albedo::GlobalConstantAlbedo,::PrimitiveEquation)
+    albedo.albedo .= albedo.α
+end
+
+## ALEBDO CLIMATOLOGY
+
+export AlbedoClimatology
+Base.@kwdef struct AlbedoClimatology{NF,Grid} <: AbstractAlbedo
+    "number of latitudes on one hemisphere, Equator included"
+    nlat_half::Int
+
+    "[OPTION] path to the folder containing the albedo file, pkg path default"
+    path::String = "SpeedyWeather.jl/input_data"
+
+    "[OPTION] filename of albedo"
+    file::String = "albedo.nc"
+
+    "[OPTION] variable name in netcdf file"
+    varname::String = "alb"
+
+    "[OPTION] Grid the albedo file comes on"
+    file_Grid::Type{<:AbstractGrid} = FullGaussianGrid
+
+    "Albedo climatology"
+    albedo::Grid = zeros(Grid, nlat_half)
+end
+
+function AlbedoClimatology(SG::SpectralGrid; kwargs...)
+    (; NF, Grid, nlat_half) = SG
+    AlbedoClimatology{NF,Grid{NF}}(;nlat_half, kwargs...)
+end
+
+function initialize!(albedo::AlbedoClimatology, model::PrimitiveEquation)
+
+    # LOAD NETCDF FILE
+    if albedo.path == "SpeedyWeather.jl/input_data"
+        path = joinpath(@__DIR__, "../../input_data", albedo.file)
+    else
+        path = joinpath(albedo.path, albedo.file)
+    end
+    ncfile = NCDataset(path)
+
+    a = albedo.file_Grid(ncfile[albedo.varname][:, :])
+    interpolate!(albedo.albedo, a)
+end