From 0c44b113eb301460075e721ea9daf74833df14bc Mon Sep 17 00:00:00 2001 From: Milan Date: Thu, 19 Oct 2023 15:10:00 -0400 Subject: [PATCH] soil moisture and vegetation --- src/dynamics/diagnostic_variables.jl | 1 + src/dynamics/models.jl | 8 +- src/dynamics/ocean.jl | 34 +++- src/dynamics/prognostic_variables.jl | 8 +- src/dynamics/time_integration.jl | 17 +- src/physics/column_variables.jl | 1 + src/physics/define_column.jl | 1 + src/physics/land.jl | 290 +++++++++++++++++++++++++++ src/physics/surface_fluxes.jl | 63 +----- 9 files changed, 343 insertions(+), 80 deletions(-) create mode 100644 src/physics/land.jl diff --git a/src/dynamics/diagnostic_variables.jl b/src/dynamics/diagnostic_variables.jl index e24bb1593..d48657756 100644 --- a/src/dynamics/diagnostic_variables.jl +++ b/src/dynamics/diagnostic_variables.jl @@ -138,6 +138,7 @@ Base.@kwdef struct SurfaceVariables{NF<:AbstractFloat,Grid<:AbstractGrid{NF}} precip_large_scale::Grid = zeros(Grid,nlat_half) # large scale precipitation (for output) precip_convection::Grid = zeros(Grid,nlat_half) # convective precipitation (for output) + soil_moisture_availability::Grid = zeros(Grid,nlat_half) end # generator function based on a SpectralGrid diff --git a/src/dynamics/models.jl b/src/dynamics/models.jl index d9c864759..875d49be5 100644 --- a/src/dynamics/models.jl +++ b/src/dynamics/models.jl @@ -153,7 +153,7 @@ Base.@kwdef struct PrimitiveDryModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr # BOUNDARY CONDITIONS land_sea_mask::AbstractLandSeaMask{NF} = LandSeaMask(spectral_grid) ocean::AbstractOcean{NF} = SeasonalOceanClimatology(spectral_grid) - land::AbstractLand{NF} = SeasonalLandClimatology(spectral_grid) + land::AbstractLand{NF} = SeasonalLandTemperature(spectral_grid) # PHYSICS/PARAMETERIZATIONS physics::Bool = true @@ -237,7 +237,9 @@ Base.@kwdef struct PrimitiveWetModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr # BOUNDARY CONDITIONS land_sea_mask::AbstractLandSeaMask{NF} = LandSeaMask(spectral_grid) ocean::AbstractOcean{NF} = SeasonalOceanClimatology(spectral_grid) - land::AbstractLand{NF} = SeasonalLandClimatology(spectral_grid) + land::AbstractLand{NF} = SeasonalLandTemperature(spectral_grid) + soil::AbstractSoil{NF} = SeasonalSoilMoisture(spectral_grid) + vegetation::AbstractVegetation{NF} = VegetationClimatology(spectral_grid) # PHYSICS/PARAMETERIZATIONS physics::Bool = true @@ -287,6 +289,8 @@ function initialize!(model::PrimitiveWet) initialize!(model.land_sea_mask) initialize!(model.ocean) initialize!(model.land) + initialize!(model.soil) + initialize!(model.vegetation) # parameterizations initialize!(model.boundary_layer_drag,model) diff --git a/src/dynamics/ocean.jl b/src/dynamics/ocean.jl index f759e5383..5fa5d418e 100644 --- a/src/dynamics/ocean.jl +++ b/src/dynamics/ocean.jl @@ -21,6 +21,9 @@ Base.@kwdef struct SeasonalOceanClimatology{NF,Grid<:AbstractGrid{NF}} <: Abstra "filename of sea surface temperatures" file::String = "sea_surface_temperature.nc" + "Variable name in netcdf file" + varname::String = "sst" + "Grid the sea surface temperature file comes on" file_Grid::Type{<:AbstractGrid} = FullGaussianGrid @@ -39,33 +42,42 @@ function SeasonalOceanClimatology(SG::SpectralGrid;kwargs...) end function initialize!(ocean::SeasonalOceanClimatology{NF,Grid}) where {NF,Grid} + load_monthly_climatology!(ocean.monthly_temperature,ocean) +end + +function load_monthly_climatology!( + monthly::Vector{Grid}, + scheme; + varname::String = scheme.varname +) where {Grid<:AbstractGrid} + # LOAD NETCDF FILE - if ocean.path == "SpeedyWeather.jl/input_data" - path = joinpath(@__DIR__,"../../input_data",ocean.file) + if scheme.path == "SpeedyWeather.jl/input_data" + path = joinpath(@__DIR__,"../../input_data",scheme.file) else - path = joinpath(ocean.path,ocean.file) + path = joinpath(scheme.path,scheme.file) end ncfile = NCDataset(path) # create interpolator from grid in file to grid used in model nx, ny = ncfile.dim["lon"], ncfile.dim["lat"] npoints = nx*ny - NF_file = typeof(ncfile["sst"].attrib["_FillValue"]) - sst = ocean.file_Grid(zeros(NF_file,npoints)) - interp = RingGrids.interpolator(NF,ocean.monthly_temperature[1],sst) + NF_file = typeof(ncfile[varname].attrib["_FillValue"]) + grid = scheme.file_Grid(zeros(NF_file,npoints)) + interp = RingGrids.interpolator(Float32,monthly[1],grid) - # interpolate and store in ocean + # interpolate and store in monthly for month in 1:12 - sst_this_month = ncfile["sst"][:,:,month] + this_month = ncfile[varname][:,:,month] ij = 0 for j in 1:ny for i in 1:nx ij += 1 - x = sst_this_month[i,j] - sst[ij] = ismissing(x) ? ocean.missing_value : x + x = this_month[i,j] + grid[ij] = ismissing(x) ? scheme.missing_value : x end end - interpolate!(ocean.monthly_temperature[month],sst,interp) + interpolate!(monthly[month],grid,interp) end return nothing end diff --git a/src/dynamics/prognostic_variables.jl b/src/dynamics/prognostic_variables.jl index 694f1c46d..d7d6c66ce 100644 --- a/src/dynamics/prognostic_variables.jl +++ b/src/dynamics/prognostic_variables.jl @@ -136,11 +136,11 @@ Base.@kwdef mutable struct PrognosticVariablesLand{NF<:AbstractFloat,Grid<:Abstr "Snow depth [m]" const snow_depth::Grid = zeros(Grid,nlat_half) - "Soil wetness layer 1, volume fraction [1]" - const soil_wetness_layer1::Grid = zeros(Grid,nlat_half) + "Soil moisture layer 1, volume fraction [1]" + const soil_moisture_layer1::Grid = zeros(Grid,nlat_half) - "Soil wetness layer 2, volume fraction [1]" - const soil_wetness_layer2::Grid = zeros(Grid,nlat_half) + "Soil moisture layer 2, volume fraction [1]" + const soil_moisture_layer2::Grid = zeros(Grid,nlat_half) end # generator function based on a SpectralGrid diff --git a/src/dynamics/time_integration.jl b/src/dynamics/time_integration.jl index eb39e517a..c668d17fc 100644 --- a/src/dynamics/time_integration.jl +++ b/src/dynamics/time_integration.jl @@ -221,12 +221,17 @@ function timestep!( progn::PrognosticVariables{NF}, # all prognostic variables model.feedback.nars_detected && return nothing # exit immediately if NaRs already present (;time) = progn.clock # current time - ocean_timestep!(progn.ocean,time,model) - land_timestep!(progn.land,time,model) - - # switch on/off all physics - model.physics && parameterization_tendencies!(diagn,progn,time,model) - model.physics || zero_tendencies!(diagn) # set tendencies to zero otherwise + if model.physics # switch on/off all physics parameterizations + # time step ocean (temperature and TODO sea ice) and land (temperature and soil moisture) + ocean_timestep!(progn.ocean,time,model) + land_timestep!(progn.land,time,model) + soil_moisture_availability!(diagn.surface,progn.land,model) + + # calculate all parameterizations + parameterization_tendencies!(diagn,progn,time,model) + else # set tendencies to zero otherwise for accumulators + zero_tendencies!(diagn) + end dynamics_tendencies!(diagn,progn,model,lf2) # dynamical core implicit_correction!(diagn,model.implicit,progn) # semi-implicit time stepping corrections diff --git a/src/physics/column_variables.jl b/src/physics/column_variables.jl index d4a485e66..9a1989ce7 100644 --- a/src/physics/column_variables.jl +++ b/src/physics/column_variables.jl @@ -40,6 +40,7 @@ function get_column!( C::ColumnVariables, # TODO skin = surface approximation for now C.skin_temperature_sea = P.ocean.sea_surface_temperature[ij] C.skin_temperature_land = P.land.land_surface_temperature[ij] + C.soil_moisture_availability = D.surface.soil_moisture_availability[ij] end """Recalculate ring index if not provided.""" diff --git a/src/physics/define_column.jl b/src/physics/define_column.jl index 77968dcef..efa5a10bc 100644 --- a/src/physics/define_column.jl +++ b/src/physics/define_column.jl @@ -53,6 +53,7 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV surface_wind_speed::NF = 0 skin_temperature_sea::NF = 0 skin_temperature_land::NF = 0 + soil_moisture_availability::NF = 0 # THERMODYNAMICS surface_air_density::NF = 0 diff --git a/src/physics/land.jl b/src/physics/land.jl new file mode 100644 index 000000000..e5053c719 --- /dev/null +++ b/src/physics/land.jl @@ -0,0 +1,290 @@ +## LAND TEMPERATURE +abstract type AbstractLand{NF,Grid} end + +function Base.show(io::IO,L::AbstractLand) + println(io,"$(typeof(L)) <: AbstractLand") + keys = propertynames(L) + print_fields(io,L,keys) +end + +Base.@kwdef struct SeasonalLandTemperature{NF,Grid<:AbstractGrid{NF}} <: AbstractLand{NF,Grid} + + "number of latitudes on one hemisphere, Equator included" + nlat_half::Int + + # OPTIONS + "Time step used to update land surface temperatures" + Δt::Dates.Day = Dates.Day(3) + + "path to the folder containing the land temperature file, pkg path default" + path::String = "SpeedyWeather.jl/input_data" + + "filename of land surface temperatures" + file::String = "land_surface_temperature.nc" + + "variable name in netcdf file" + varname::String = "lst" + + "Grid the land surface temperature file comes on" + file_Grid::Type{<:AbstractGrid} = FullGaussianGrid + + "The missing value in the data respresenting ocean" + missing_value::NF = NF(NaN) + + # to be filled from file + "Monthly land surface temperatures [K], interpolated onto Grid" + monthly_temperature::Vector{Grid} = [zeros(Grid,nlat_half) for _ in 1:12] +end + +# generator function +function SeasonalLandTemperature(SG::SpectralGrid;kwargs...) + (;NF,Grid,nlat_half) = SG + return SeasonalLandTemperature{NF,Grid{NF}}(;nlat_half,kwargs...) +end + +function initialize!(land::SeasonalLandTemperature{NF,Grid}) where {NF,Grid} + load_monthly_climatology!(land.monthly_temperature,land) +end + +function initialize!( land::PrognosticVariablesLand, + time::DateTime, + model::PrimitiveEquation) + land_timestep!(land,time,model.land,initialize=true) + model isa PrimitiveWet && soil_timestep!(land,time,model.soil) +end + +# function barrier +function land_timestep!(land::PrognosticVariablesLand, + time::DateTime, + model::PrimitiveEquation; + initialize::Bool = false) + # the time step is dictated by the land "model" + executed = land_timestep!(land,time,model.land;initialize) + executed && model isa PrimitiveWet && soil_timestep!(land,time,model.soil) +end + +function land_timestep!(land::PrognosticVariablesLand{NF}, + time::DateTime, + land_model::SeasonalLandTemperature; + initialize::Bool = false) where NF + + # escape immediately if Δt of land model hasn't passed yet + # unless the land hasn't been initialized yet + initialize || (time - land.time) < land_model.Δt && return false # = executed + + # otherwise update land prognostic variables: + land.time = time + this_month = Dates.month(time) + next_month = (this_month % 12) + 1 # mod for dec 12 -> jan 1 + + # linear interpolation weight between the two months + # TODO check whether this shifts the climatology by 1/2 a month + weight = convert(NF,Dates.days(time-Dates.firstdayofmonth(time))/Dates.daysinmonth(time)) + + (;monthly_temperature) = land_model + @. land.land_surface_temperature = (1-weight) * monthly_temperature[this_month] + + weight * monthly_temperature[next_month] + + return true # = executed +end + +## SOIL MOISTURE +abstract type AbstractSoil{NF,Grid} end + +function Base.show(io::IO,S::AbstractSoil) + println(io,"$(typeof(S)) <: AbstractSoil") + keys = propertynames(S) + print_fields(io,S,keys) +end + +Base.@kwdef struct SeasonalSoilMoisture{NF,Grid<:AbstractGrid{NF}} <: AbstractSoil{NF,Grid} + + "number of latitudes on one hemisphere, Equator included" + nlat_half::Int + + # OPTIONS + "Depth of top soil layer [m]" + D_top::Float64 = 0.07 + + "Depth of root layer [m]" + D_root::Float64 = 0.21 + + "Soil wetness at field capacity [volume fraction]" + W_cap::Float64 = 0.3 + + "Soil wetness at wilting point [volume fraction]" + W_wilt::Float64 = 0.17 + + # READ CLIMATOLOGY FROM FILE + "path to the folder containing the soil moisture file, pkg path default" + path::String = "SpeedyWeather.jl/input_data" + + "filename of soil moisture" + file::String = "soil_moisture.nc" + + "variable name in netcdf file" + varname_layer1::String = "swl1" + varname_layer2::String = "swl2" + + "Grid the soil moisture file comes on" + file_Grid::Type{<:AbstractGrid} = FullGaussianGrid + + "The missing value in the data respresenting ocean" + missing_value::NF = NF(NaN) + + # to be filled from file + "Monthly soil moisture volume fraction [1], top layer, interpolated onto Grid" + monthly_soil_moisture_layer1::Vector{Grid} = [zeros(Grid,nlat_half) for _ in 1:12] + + "Monthly soil moisture volume fraction [1], 2nd layer, interpolated onto Grid" + monthly_soil_moisture_layer2::Vector{Grid} = [zeros(Grid,nlat_half) for _ in 1:12] +end + +# generator function +function SeasonalSoilMoisture(SG::SpectralGrid;kwargs...) + (;NF,Grid,nlat_half) = SG + return SeasonalSoilMoisture{NF,Grid{NF}}(;nlat_half,kwargs...) +end + +function initialize!(soil::SeasonalSoilMoisture{NF,Grid}) where {NF,Grid} + load_monthly_climatology!(soil.monthly_soil_moisture_layer1,soil,varname=soil.varname_layer1) + load_monthly_climatology!(soil.monthly_soil_moisture_layer2,soil,varname=soil.varname_layer2) +end + +function soil_timestep!(land::PrognosticVariablesLand{NF}, + time::DateTime, + soil_model::SeasonalSoilMoisture) where NF + + this_month = Dates.month(time) + next_month = (this_month % 12) + 1 # mod for dec 12 -> jan 1 + + # linear interpolation weight between the two months + # TODO check whether this shifts the climatology by 1/2 a month + weight = convert(NF,Dates.days(time-Dates.firstdayofmonth(time))/Dates.daysinmonth(time)) + + (;monthly_soil_moisture_layer1, monthly_soil_moisture_layer2) = soil_model + @. land.soil_moisture_layer1 = (1-weight) * monthly_soil_moisture_layer1[this_month] + + weight * monthly_soil_moisture_layer1[next_month] + @. land.soil_moisture_layer2 = (1-weight) * monthly_soil_moisture_layer2[this_month] + + weight * monthly_soil_moisture_layer2[next_month] + + return nothing +end + +## SOIL MOISTURE +abstract type AbstractVegetation{NF,Grid} end + +function Base.show(io::IO,A::AbstractVegetation) + println(io,"$(typeof(A)) <: AbstractVegetation") + keys = propertynames(A) + print_fields(io,A,keys) +end + +Base.@kwdef struct VegetationClimatology{NF,Grid<:AbstractGrid{NF}} <: AbstractVegetation{NF,Grid} + + "number of latitudes on one hemisphere, Equator included" + nlat_half::Int + + # OPTIONS + "Combine high and low vegetation factor, a in high + a*low [1]" + low_veg_factor::Float64 = 0.8 + + "path to the folder containing the soil moisture file, pkg path default" + path::String = "SpeedyWeather.jl/input_data" + + "filename of soil moisture" + file::String = "vegetation.nc" + + "variable name in netcdf file" + varname_vegh::String = "vegh" + varname_vegl::String = "vegl" + + "Grid the soil moisture file comes on" + file_Grid::Type{<:AbstractGrid} = FullGaussianGrid + + "The missing value in the data respresenting ocean" + missing_value::NF = NF(NaN) + + # to be filled from file + "High vegetation cover [1], interpolated onto Grid" + high_cover::Grid = zeros(Grid,nlat_half) + + "Low vegetation cover [1], interpolated onto Grid" + low_cover::Grid = zeros(Grid,nlat_half) +end + +# generator function +function VegetationClimatology(SG::SpectralGrid;kwargs...) + (;NF,Grid,nlat_half) = SG + return VegetationClimatology{NF,Grid{NF}}(;nlat_half,kwargs...) +end + +function initialize!(vegetation::VegetationClimatology) + + # LOAD NETCDF FILE + if vegetation.path == "SpeedyWeather.jl/input_data" + path = joinpath(@__DIR__,"../../input_data",vegetation.file) + else + path = joinpath(vegetation.path,vegetation.file) + end + ncfile = NCDataset(path) + + # high and low vegetation cover + vegh = vegetation.file_Grid(ncfile[vegetation.varname_vegh][:,:]) + vegl = vegetation.file_Grid(ncfile[vegetation.varname_vegl][:,:]) + + # interpolate onto grid + high_vegetation_cover = vegetation.high_cover + low_vegetation_cover = vegetation.low_cover + interpolator = RingGrids.interpolator(Float32,high_vegetation_cover,vegh) + interpolate!(high_vegetation_cover,vegh,interpolator) + interpolate!(low_vegetation_cover,vegl,interpolator) +end + +# function barrier +function soil_moisture_availability!( + diagn::SurfaceVariables, + land::PrognosticVariablesLand, + model::PrimitiveWet) + soil_moisture_availability!(diagn,land,model.soil,model.vegetation) +end + +function soil_moisture_availability!( + ::SurfaceVariables, + ::PrognosticVariablesLand, + ::PrimitiveDry) + return nothing +end + +function soil_moisture_availability!( + diagn::SurfaceVariables{NF}, + land::PrognosticVariablesLand, + soil::AbstractSoil, + vegetation::AbstractVegetation) where NF + + (;soil_moisture_availability) = diagn + (;soil_moisture_layer1, soil_moisture_layer2) = land + (;high_cover, low_cover) = vegetation + + D_top = convert(NF,soil.D_top) + D_root = convert(NF,soil.D_root) + W_cap = convert(NF,soil.W_cap) + W_wilt = convert(NF,soil.W_wilt) + low_veg_factor = convert(NF,vegetation.low_veg_factor) + + # precalculate + r = 1/(D_top*W_cap + D_root*(W_cap - W_wilt)) + + @inbounds for ij in eachgridpoint(soil_moisture_availability, + soil_moisture_layer1, + soil_moisture_layer2, + high_cover,low_cover) + + # Fortran SPEEDY source/land_model.f90 line 111 origin unclear + veg = max(0,high_cover[ij] + low_veg_factor*low_cover[ij]) + + # Fortran SPEEDY documentation eq. 51 + soil_moisture_availability[ij] = r*(D_top*soil_moisture_layer1[ij] + + veg*D_root*max(soil_moisture_layer2[ij]- W_wilt,0)) + end +end diff --git a/src/physics/surface_fluxes.jl b/src/physics/surface_fluxes.jl index 3f13a74f7..af67f9033 100644 --- a/src/physics/surface_fluxes.jl +++ b/src/physics/surface_fluxes.jl @@ -7,7 +7,6 @@ function surface_fluxes!(column::ColumnVariables,model::PrimitiveEquation) surface_wind_stress!(column,model.surface_wind) # now call other heat and humidity fluxes - # soil_moisture!(column,model.vegetation) sensible_heat_flux!(column,model.surface_heat_flux,model.constants) evaporation!(column,model) end @@ -116,8 +115,8 @@ function sensible_heat_flux!( column::ColumnVariables{NF}, flux_sea = ρ*heat_exchange_sea*V₀*cₚ*(T_skin_sea - T) # mix fluxes for fractional land-sea mask - land_available = ~isnan(T_skin_land) - sea_available = ~isnan(T_skin_sea) + land_available = isfinite(T_skin_land) + sea_available = isfinite(T_skin_sea) if land_available && sea_available column.flux_temp_upward[end] += land_fraction*flux_land + (1-land_fraction)*flux_sea @@ -137,59 +136,11 @@ function sensible_heat_flux!( column::ColumnVariables{NF}, end Base.@kwdef struct SurfaceEvaporation{NF<:AbstractFloat} <: AbstractEvaporation{NF} - - nlat_half::Int - moisture_exchange_land::Float64 = 1.2e-3 # for neutral stability moisture_exchange_sea::Float64 = 0.9e-3 - - "Depth of top soil layer [m]" - D_top::Float64 = 0.07 - - "Depth of root layer [m]" - D_root::Float64 = 0.21 - - "Soil wetness at field capacity [volume fraction]" - W_cap::Float64 = 0.3 - - "Soil wetness at wilting point [volume fraction]" - W_wilt::Float64 = 0.17 - - # FILE OPTIONS - "path to the folder containing the soil moisture and vegetation file, pkg path default" - path::String = "SpeedyWeather.jl/input_data" - - "filename of that netcdf file" - file::String = "soil_moisture.nc" - - # "Grid the soil and vegetation file comes on" - # file_Grid::Type{<:AbstractGrid} = FullGaussianGrid - - # # arrays to be initialised - # "Soil moisture availability index" - # soil_moisture_availability::Grid = zeros(Grid{NF},nlat_half) end -SurfaceEvaporation(SG::SpectralGrid;kwargs...) = SurfaceEvaporation{SG.NF}(nlat_half=SG.nlat_half;kwargs...) - -function initialize!(evaporation::SurfaceEvaporation) - - # LOAD NETCDF FILE - if evaporation.path == "SpeedyWeather.jl/input_data" - path = joinpath(@__DIR__,"../../input_data",evaporation.file) - else - path = joinpath(evaporation.path,evaporation.file) - end - ncfile = NCDataset(path) - - # high resolution land-sea mask - lsm_highres = evaporation.file_Grid(ncfile["lsm"][:,:]) - - # average onto grid cells of the model - RingGrids.grid_cell_average!(land_sea_mask.land_sea_mask,lsm_highres) -end - -SurfaceEvaporation(SG::SpectralGrid;kwargs...) = SurfaceEvaporation{SG.NF}(nlat_half=SG.nlat_half;kwargs...) +SurfaceEvaporation(SG::SpectralGrid;kwargs...) = SurfaceEvaporation{SG.NF}(;kwargs...) # don't do anything for dry core function evaporation!( column::ColumnVariables, @@ -210,9 +161,7 @@ function evaporation!( column::ColumnVariables{NF}, (;skin_temperature_sea, skin_temperature_land, pres, humid) = column moisture_exchange_land = convert(NF,evaporation.moisture_exchange_land) moisture_exchange_sea = convert(NF,evaporation.moisture_exchange_sea) - - # α = soil_moisture.availability[column.ij] - α = convert(NF,0.2) + α = column.soil_moisture_availability # SATURATION HUMIDITY OVER LAND AND OCEAN surface_pressure = pres[end] @@ -226,8 +175,8 @@ function evaporation!( column::ColumnVariables{NF}, flux_land = ρ*moisture_exchange_land*V₀*α*max(sat_humid_land - humid[end],0) # mix fluxes for fractional land-sea mask - land_available = ~isnan(skin_temperature_land) - sea_available = ~isnan(skin_temperature_sea) + land_available = isfinite(skin_temperature_land) && isfinite(α) + sea_available = isfinite(skin_temperature_sea) if land_available && sea_available column.flux_humid_upward[end] += land_fraction*flux_land + (1-land_fraction)*flux_sea