Merge pull request #433 from SpeedyWeather/mk/humidity_diffusion

Vertical diffusion of humidity

src/SpeedyWeather.jl

@@ -183,6 +183,7 @@ include("dynamics/vertical_advection.jl")
 include("dynamics/implicit.jl")
 include("dynamics/scaling.jl")
 include("dynamics/tendencies.jl")
+include("dynamics/hole_filling.jl")
 
 # PARAMETERIZATIONS
 include("physics/tendencies.jl")

src/abstract_types.jl

@@ -1,4 +1,5 @@
 # MODELS
+abstract type AbstractSimulation{Model} end
 abstract type ModelSetup end
 abstract type Barotropic <: ModelSetup end
 abstract type ShallowWater <: ModelSetup end

src/dynamics/hole_filling.jl

@@ -0,0 +1,22 @@
+abstract type AbstractHoleFilling{NF} end
+
+struct ClipNegatives{NF} <: AbstractHoleFilling{NF} end
+ClipNegatives(SG::SpectralGrid) = ClipNegatives{SG.NF}()
+
+# nothing to initialize
+initialize!(H::ClipNegatives,::PrimitiveWet) = nothing
+
+# function barrier
+function hole_filling!(
+    A::AbstractGrid,
+    H::ClipNegatives,
+    model::PrimitiveWet)
+
+    hole_filling!(A,H)
+end
+
+function hole_filling!(A::AbstractGrid,H::ClipNegatives)
+    @inbounds for ij in eachgridpoint(A)
+        A[ij] = max(A[ij],0)
+    end
+end

src/dynamics/initial_conditions.jl

@@ -463,7 +463,6 @@ function initialize_humidity!(  progn::PrognosticVariables,
     spectral_truncation!(humid_surf)
 
     # Specific humidity at tropospheric levels (stratospheric humidity remains zero)
-    a = model.spectral_transform.norm_sphere
     for k in n_stratosphere_levels+1:nlev
         for lm in eachharmonic(humid_surf)
             progn.layers[k].timesteps[1].humid[lm] = humid_surf[lm]*σ_levels_full[k]^scale_height_ratio

src/dynamics/models.jl

@@ -3,7 +3,7 @@ $(TYPEDSIGNATURES)
 Simulation is a container struct to be used with `run!(::Simulation)`.
 It contains
 $(TYPEDFIELDS)"""
-struct Simulation{Model<:ModelSetup}
+struct Simulation{Model<:ModelSetup} <: AbstractSimulation{Model}
     "define the current state of the model"
     prognostic_variables::PrognosticVariables
 
@@ -19,7 +19,7 @@ $(SIGNATURES)
 The BarotropicModel struct holds all other structs that contain precalculated constants,
 whether scalars or arrays that do not change throughout model integration.
 $(TYPEDFIELDS)"""
-Base.@kwdef struct BarotropicModel{NF<:AbstractFloat, D<:AbstractDevice} <: Barotropic
+Base.@kwdef mutable struct BarotropicModel{NF<:AbstractFloat, D<:AbstractDevice} <: Barotropic
     spectral_grid::SpectralGrid = SpectralGrid(nlev=1)
 
     # DYNAMICS
@@ -81,7 +81,7 @@ $(SIGNATURES)
 The ShallowWaterModel struct holds all other structs that contain precalculated constants,
 whether scalars or arrays that do not change throughout model integration.
 $(TYPEDFIELDS)"""
-Base.@kwdef struct ShallowWaterModel{NF<:AbstractFloat, D<:AbstractDevice} <: ShallowWater
+Base.@kwdef mutable struct ShallowWaterModel{NF<:AbstractFloat, D<:AbstractDevice} <: ShallowWater
     spectral_grid::SpectralGrid = SpectralGrid(nlev=1)
 
     # DYNAMICS
@@ -145,7 +145,7 @@ $(SIGNATURES)
 The PrimitiveDryModel struct holds all other structs that contain precalculated constants,
 whether scalars or arrays that do not change throughout model integration.
 $(TYPEDFIELDS)"""
-Base.@kwdef struct PrimitiveDryModel{NF<:AbstractFloat, D<:AbstractDevice} <: PrimitiveDry
+Base.@kwdef mutable struct PrimitiveDryModel{NF<:AbstractFloat, D<:AbstractDevice} <: PrimitiveDry
     spectral_grid::SpectralGrid = SpectralGrid()
 
     # DYNAMICS
@@ -162,8 +162,8 @@ Base.@kwdef struct PrimitiveDryModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr
 
     # PHYSICS/PARAMETERIZATIONS
     physics::Bool = true
-    boundary_layer_drag::BoundaryLayerDrag{NF} = NoBoundaryLayerDrag(spectral_grid)
-    temperature_relaxation::TemperatureRelaxation{NF} = NoTemperatureRelaxation(spectral_grid)
+    boundary_layer_drag::BoundaryLayerDrag{NF} = LinearDrag(spectral_grid)
+    temperature_relaxation::TemperatureRelaxation{NF} = HeldSuarez(spectral_grid)
     static_energy_diffusion::VerticalDiffusion{NF} = StaticEnergyDiffusion(spectral_grid)
     surface_thermodynamics::AbstractSurfaceThermodynamics{NF} = SurfaceThermodynamicsConstant(spectral_grid)
     surface_wind::AbstractSurfaceWind{NF} = SurfaceWind(spectral_grid)
@@ -233,7 +233,7 @@ $(SIGNATURES)
 The PrimitiveDryModel struct holds all other structs that contain precalculated constants,
 whether scalars or arrays that do not change throughout model integration.
 $(TYPEDFIELDS)"""
-Base.@kwdef struct PrimitiveWetModel{NF<:AbstractFloat, D<:AbstractDevice} <: PrimitiveWet
+Base.@kwdef mutable struct PrimitiveWetModel{NF<:AbstractFloat, D<:AbstractDevice} <: PrimitiveWet
     spectral_grid::SpectralGrid = SpectralGrid()
 
     # DYNAMICS
@@ -256,6 +256,7 @@ Base.@kwdef struct PrimitiveWetModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr
     boundary_layer_drag::BoundaryLayerDrag{NF} = NoBoundaryLayerDrag(spectral_grid)
     temperature_relaxation::TemperatureRelaxation{NF} = NoTemperatureRelaxation(spectral_grid)
     static_energy_diffusion::VerticalDiffusion{NF} = StaticEnergyDiffusion(spectral_grid)
+    humidity_diffusion::VerticalDiffusion{NF} = HumidityDiffusion(spectral_grid)
     large_scale_condensation::AbstractCondensation{NF} = SpeedyCondensation(spectral_grid)
     surface_thermodynamics::AbstractSurfaceThermodynamics{NF} = SurfaceThermodynamicsConstant(spectral_grid)
     surface_wind::AbstractSurfaceWind{NF} = SurfaceWind(spectral_grid)
@@ -269,7 +270,8 @@ Base.@kwdef struct PrimitiveWetModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr
     horizontal_diffusion::HorizontalDiffusion{NF} = HyperDiffusion(spectral_grid)
     implicit::AbstractImplicit{NF} = ImplicitPrimitiveEq(spectral_grid)
     vertical_advection::VerticalAdvection{NF} = CenteredVerticalAdvection(spectral_grid)
-
+    hole_filling::AbstractHoleFilling{NF} = ClipNegatives(spectral_grid)
+
     # INTERNALS
     geometry::Geometry{NF} = Geometry(spectral_grid)
     constants::DynamicsConstants{NF} = DynamicsConstants(spectral_grid,planet,atmosphere,geometry)
@@ -309,6 +311,7 @@ function initialize!(model::PrimitiveWet;time::DateTime = DEFAULT_DATE)
     initialize!(model.boundary_layer_drag,model)
     initialize!(model.temperature_relaxation,model)
     initialize!(model.static_energy_diffusion,model)
+    initialize!(model.humidity_diffusion,model)
     initialize!(model.large_scale_condensation,model)
     initialize!(model.convection,model)
 

src/dynamics/tendencies.jl

@@ -377,17 +377,23 @@ function vordiv_tendencies!(
     return nothing
 end
 
+# function barrier
 function tendencies_physics_only!(
     diagn::DiagnosticVariablesLayer,
     model::PrimitiveEquation
 )
-    tendencies_physics_only!(diagn,model.geometry,model.spectral_transform)
+    wet_core = model isa PrimitiveWet
+    tendencies_physics_only!(diagn,model.geometry,model.spectral_transform,wet_core)
 end
 
+"""For dynamics=false, after calling parameterization_tendencies! call this function
+to transform the physics tendencies from grid-point to spectral space including the
+necessary coslat⁻¹ scaling."""
 function tendencies_physics_only!(
     diagn::DiagnosticVariablesLayer,
     G::Geometry,
     S::SpectralTransform,
+    wet_core::Bool = true
 )
     (;coslat⁻¹) = G
     (;u_tend_grid, v_tend_grid, temp_tend_grid, humid_tend_grid) = diagn.tendencies  # already contains parameterizations
@@ -411,16 +417,13 @@ function tendencies_physics_only!(
     spectral!(u_tend,u_tend_grid,S)
     spectral!(v_tend,v_tend_grid,S)
     spectral!(temp_tend,temp_tend_grid,S)
-    spectral!(humid_tend,humid_tend_grid,S)
+    wet_core && spectral!(humid_tend,humid_tend_grid,S)
 
     curl!(vor_tend,u_tend,v_tend,S)         # ∂ζ/∂t = ∇×(u_tend,v_tend)
     divergence!(div_tend,u_tend,v_tend,S)   # ∂D/∂t = ∇⋅(u_tend,v_tend)
     return nothing
 end
 
-
-
-
 """
 $(TYPEDSIGNATURES)
 Function barrier to unpack `model`."""
@@ -847,6 +850,7 @@ function SpeedyTransforms.gridded!( diagn::DiagnosticVariablesLayer,
     U = diagn.dynamics_variables.a      # reuse work arrays for velocities spectral
     V = diagn.dynamics_variables.b      # U = u*coslat, V=v*coslat
 
+    # retain previous time step for vertical advection
     @. temp_grid_prev = temp_grid
     @. u_grid_prev = u_grid
     @. v_grid_prev = v_grid
@@ -862,10 +866,14 @@ function SpeedyTransforms.gridded!( diagn::DiagnosticVariablesLayer,
     gridded!(vor_grid,vor,S)                # get vorticity on grid from spectral vor
     gridded!(div_grid,div,S)                # get divergence on grid from spectral div
     gridded!(temp_grid,temp,S)              # (absolute) temperature
-    wet_core && gridded!(humid_grid,humid,S)# specific humidity (wet core only)
+
+    if wet_core                             # specific humidity (wet core only)
+        gridded!(humid_grid,humid,S)        
+        hole_filling!(humid_grid,model.hole_filling,model)  # remove negative humidity
+    end
 
     # include humidity effect into temp for everything stability-related
-    temperature_average!(diagn,temp,S)      # TODO: do at frequency of reinitialize implicit?
+    temperature_average!(diagn,temp,S)
     virtual_temperature!(diagn,temp,model)  # temp = virt temp for dry core
 
     # transform from U,V in spectral to u,v on grid (U,V = u,v*coslat)

src/dynamics/time_integration.jl

@@ -339,10 +339,10 @@ function timestep!( progn::PrognosticVariables{NF}, # all prognostic variables
             progn_layer = progn.layers[k]
             progn_layer_lf = progn_layer.timesteps[lf2]
 
-            horizontal_diffusion!(progn_layer,diagn_layer,model)    # implicit diffusion of vor, div, temp
-            leapfrog!(progn_layer,diagn_layer,dt,lf1,model)         # time step forward for vor, div, temp
+            horizontal_diffusion!(progn_layer,diagn_layer,model)    # for vor, div, temp, humid
+            leapfrog!(progn_layer,diagn_layer,dt,lf1,model)         # time step forward for vor, div, temp, humid
             gridded!(diagn_layer,progn_layer_lf,model)              # propagate spectral state to grid
-        else                                                        # surface level
+        else                                                        # surface level, time step for pressure
             leapfrog!(progn.surface,diagn.surface,dt,lf1,model)
             (;pres_grid) = diagn.surface
             pres_lf = progn.surface.timesteps[lf2].pres

src/output/feedback.jl

@@ -161,7 +161,7 @@ function nar_detection!(feedback::Feedback,progn::PrognosticVariables)
     if ~nars_detected_here
         nars_vor = ~isfinite(vor[1])    # just check first mode
                                         # spectral transform propagates NaRs globally anyway
-        nars_vor && @warn "NaR detected at time step $i"
+        nars_vor && @warn "NaN or Inf detected at time step $i"
         nars_detected_here |= nars_vor
     end
 

src/physics/boundary_layer.jl

@@ -24,12 +24,15 @@ end
 """Linear boundary layer drag Following Held and Suarez, 1996 BAMS
 $(TYPEDFIELDS)"""
 Base.@kwdef struct LinearDrag{NF<:AbstractFloat} <: BoundaryLayerDrag{NF}
+
+    # DIMENSIONS
+    nlev::Int
+
     # PARAMETERS
-    σb::Float64 = 0.7           # sigma coordinate below which linear drag is applied
-    time_scale::Float64 = 24    # [hours] time scale for linear drag coefficient at σ=1 (=1/kf in HS96)
+    σb::NF = 0.7                    # sigma coordinate below which linear drag is applied
+    time_scale::Second = Hour(24)   # time scale for linear drag coefficient at σ=1 (=1/kf in HS96)
 
     # PRECOMPUTED CONSTANTS
-    nlev::Int = 0
     drag_coefs::Vector{NF} = zeros(NF,nlev)
 end
 
@@ -45,9 +48,9 @@ function initialize!(   scheme::LinearDrag,
                         model::PrimitiveEquation)
 
     (;σ_levels_full) = model.geometry
-    (;σb,time_scale,drag_coefs) = scheme
+    (;σb, time_scale, drag_coefs) = scheme
 
-    kf = 1/(time_scale*3600)                        # scale with radius as ∂ₜu is; hrs -> sec
+    kf = 1/time_scale.value
 
     for (k,σ) in enumerate(σ_levels_full)
         drag_coefs[k] = -kf*max(0,(σ-σb)/(1-σb))    # drag only below σb, lin increasing to kf at σ=1
@@ -60,7 +63,7 @@ Compute tendency for boundary layer drag of a `column` and add to its tendencies
 function boundary_layer_drag!(  column::ColumnVariables,
                                 scheme::LinearDrag)
 
-    (;u,v,u_tend,v_tend) = column
+    (;u, v, u_tend, v_tend) = column
     (;drag_coefs) = scheme
 
     @inbounds for k in eachlayer(column)

src/physics/column_variables.jl

@@ -31,7 +31,7 @@ function get_column!(   C::ColumnVariables,
         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.humid[k] = layer.grid_variables.humid_grid[ij]
+        C.humid[k] = layer.grid_variables.humid_grid[ij] 
     end
 
     # TODO skin = surface approximation for now
@@ -53,6 +53,28 @@ function get_column!(   C::ColumnVariables,
     get_column!(C,D,P,ij,jring,G,L)
 end
 
+function get_column(    S::AbstractSimulation,
+                        ij::Integer)
+    (;prognostic_variables, diagnostic_variables) = S
+    (;geometry, land_sea_mask) = S.model
+
+    column = deepcopy(S.diagnostic_variables.columns[1])
+    reset_column!(column)
+
+    get_column!(column,
+                diagnostic_variables,
+                prognostic_variables,
+                ij,
+                geometry,
+                land_sea_mask)
+
+    # execute all parameterizations for this column to return a consistent state
+    parameterization_tendencies!(column,S.model)
+
+    @info "Receiving column at $(column.latd)˚N, $(column.lond)˚E."
+    return column
+end
+
 """
 $(TYPEDSIGNATURES)
 Write the parametrization tendencies from `C::ColumnVariables` into the horizontal fields

src/physics/define_column.jl

@@ -18,11 +18,11 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV
     latd::NF = 0                            # latitude, needed for shortwave radiation
     land_fraction::NF = 0                   # fraction of the column that is over land
 
-    # PROGNOSTIC VARIABLES, last element is surface
-    const u::Vector{NF} = zeros(NF,nlev+1)          # zonal velocity [m/s]
-    const v::Vector{NF} = zeros(NF,nlev+1)          # meridional velocity [m/s]
-    const temp::Vector{NF} = zeros(NF,nlev+1)       # temperature [K]
-    const humid::Vector{NF} = zeros(NF,nlev+1)      # specific humidity [kg/kg]
+    # PROGNOSTIC VARIABLES
+    const u::Vector{NF} = zeros(NF,nlev)        # zonal velocity [m/s]
+    const v::Vector{NF} = zeros(NF,nlev)        # meridional velocity [m/s]
+    const temp::Vector{NF} = zeros(NF,nlev)     # temperature [K]
+    const humid::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)]
@@ -50,6 +50,10 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV
     const flux_humid_upward::Vector{NF} = zeros(NF,nlev+1)
     const flux_humid_downward::Vector{NF} = zeros(NF,nlev+1)
 
+    surface_u::NF = 0
+    surface_v::NF = 0
+    surface_temp::NF = 0
+    surface_humid::NF = 0
     surface_wind_speed::NF = 0
     skin_temperature_sea::NF = 0
     skin_temperature_land::NF = 0
@@ -58,6 +62,7 @@ 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 sat_vap_pres::Vector{NF} = zeros(NF,nlev)             # Saturation vapour pressure [Pa]
     const dry_static_energy::Vector{NF} = zeros(NF,nlev)        # Dry static energy
     const moist_static_energy::Vector{NF} = zeros(NF,nlev)      # Moist static energy
@@ -121,6 +126,18 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV
     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
+    # 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)
 end

src/physics/large_scale_condensation.jl

@@ -19,7 +19,10 @@ Base.@kwdef struct SpeedyCondensation{NF<:AbstractFloat} <: AbstractCondensation
     time_scale::Second = Hour(4)
 
     "Flux limiter for humidity tendency"
-    max_flux::NF = 5
+    max_flux::NF = 0.1
+
+    "Flux limiter for condensation heating"
+    max_heating::NF = 0.001
 
     # precomputed arrays
     n_stratosphere_levels::Base.RefValue{Int} = Ref(0)
@@ -106,6 +109,7 @@ function large_scale_condensation!(
     (;nlev) = column
     pₛ = pres[end]                          # surface pressure
     pₛ_norm² = (pₛ/constants.pres_ref)^2
+    max_heating = -scheme.max_heating*time_scale⁻¹
 
     # precompute scaling constant, undo radius scaling here as directly used for output
     (;gravity, water_density) = constants
@@ -123,12 +127,12 @@ function large_scale_condensation!(
 
             # accumulate in tendencies (nothing is added if humidity not above threshold)
             humid_tend_k = -(humid[k] - humid_threshold) * time_scale⁻¹                   # eq. 22
-            
+
             # with flux limiter, use max and - as humid_tend_k is always negative
             humid_tend_k = max(humid_tend_k, -pₛ_norm²*humid_tend_max[k])
 
             # condensation heating, eq. 23 
-            temp_tend[k] -= scheme.latent_heat_cₚ[] * humid_tend_k
+            temp_tend[k] -= max(max_heating, scheme.latent_heat_cₚ[] * humid_tend_k)
 
             # 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