Skip to content

Commit 44992da

Browse files
bgroenks96bgroenks96
committed
Updates to coupling script
1 parent c29b790 commit 44992da

File tree

3 files changed

+114
-33
lines changed

3 files changed

+114
-33
lines changed

‎examples/Project.toml‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,10 @@
11
[deps]
2+
ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
3+
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
4+
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
25
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
36
CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
7+
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
48
GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
59
GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6"
610
Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
@@ -10,6 +14,7 @@ Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
1014
RingGrids = "d1845624-ad4f-453b-8ff4-a8db365bf3a7"
1115
SpeedyWeather = "9e226e20-d153-4fed-8a5b-493def4f21a9"
1216
SpeedyWeatherInternals = "34489162-d270-4603-8b96-37b04f830d73"
17+
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
1318
TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
1419
Terrarium = "80418d68-07fa-499d-ae2b-c12e531f5cd8"
1520

‎examples/speedy_dry_land.jl‎

Lines changed: 108 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -21,8 +21,8 @@ struct TerrariumDryLand{NF, TM<:ModelState{NF}} <: Speedy.AbstractDryLand
2121
model::TM
2222

2323
function TerrariumDryLand(model::ModelState{NF, Arch, Grid}; spectral_grid_kwargs...) where {NF, Arch, Grid<:ColumnRingGrid}
24-
spectral_grid = Speedy.SpectralGrid(model.grid.rings; NF, nlayers_soil=1, spectral_grid_kwargs...)
25-
return new{typeof(model)}(spectral_grid, model)
24+
spectral_grid = Speedy.SpectralGrid(model.grid.rings; NF, nlayers_soil=2, spectral_grid_kwargs...)
25+
return new{eltype(model), typeof(model)}(spectral_grid, model)
2626
end
2727
end
2828

@@ -45,51 +45,127 @@ function Speedy.timestep!(
4545
land::TerrariumDryLand,
4646
model::Speedy.PrimitiveEquation,
4747
) where {NF}
48-
Tair = diagn.grid.temp_grid[:,end] .- NF(273.15)
49-
Tair_land = get_input_field(land.model.inputs.fields, :air_temperature, XY())
50-
set!(Tair_land, Tair)
51-
# do timestep
48+
# get speedy state variables
49+
Tair = @view diagn.grid.temp_grid[:, end]
50+
pres = diagn.grid.pres_grid
51+
SwIn = diagn.physics.surface_shortwave_down
52+
LwIn = diagn.physics.surface_longwave_down
53+
SwOut = diagn.physics.surface_shortwave_up
54+
LwOut = diagn.physics.surface_longwave_up
55+
H = diagn.physics.sensible_heat_flux
56+
# update Terrarium input fields
57+
inputs = land.model.state.inputs
58+
# set!(inputs.air_temperature, Tair) # first set directly to avoid allocating new fields
59+
# set!(inputs.air_temperature, inputs.air_temperature - 273.15) # then convert to celsius
60+
# set!(inputs.air_pressure, pres) # similar with pressure
61+
# set!(inputs.air_pressure, exp(inputs.air_pressure)) # take exp to get pressure in Pa
62+
# set!(inputs.surface_shortwave_down, SwIn)
63+
# set!(inputs.surface_longwave_down, LwIn)
64+
# set!(inputs.surface_shortwave_up, SwOut)
65+
# set!(inputs.surface_longwave_up, LwOut)
66+
# set!(inputs.sensible_heat_flux, H)
67+
# do land timestep
5268
Terrarium.timestep!(land.model, Dates.second(progn.clock.Δt))
53-
# extract surface temperature and copy to speedy prognostic variable
69+
# SEB coupling:
70+
# Ts = land.model.state.skin_temperature
71+
# progn.land.soil_temperature .= Ts .+ NF(273.15)
72+
# progn.land.sensible_heat_flux .= land.model.state.sensible_heat_flux
73+
# Soil surface temperature coupling:
5474
Tsoil = land.model.state.temperature
55-
Tsurf = Tsoil[1:Tsoil.grid.Nx, 1, Tsoil.grid.Nz:Tsoil.grid.Nz]
56-
progn.land.soil_temperature .= Tsurf .+ NF(273.15)
75+
Nx, Nz = Tsoil.grid.Nx, Tsoil.grid.Nz
76+
Tsurf = @view Tsoil[1:Nx, 1, Nz-1:Nz]
77+
progn.land.soil_temperature .= reshape(Tsurf, :, 2) .+ NF(273.15)
78+
return nothing
5779
end
5880

5981
ring_grid = RingGrids.FullGaussianGrid(24)
60-
grid = ColumnRingGrid(CPU(), Float32, ExponentialSpacing(N=30), ring_grid)
61-
model = SoilModel(grid)
82+
Nz = 30
83+
Δz_min = 0.1 # currently the coupling is only stable with a large surface layer
84+
grid = ColumnRingGrid(CPU(), Float32, ExponentialSpacing(;N=Nz, Δz_min), ring_grid)
85+
# grid = ColumnGrid(CPU(), Float32, ExponentialSpacing(N=30))
6286
# Initial conditions
63-
initializer = FieldInitializers(
87+
soil_initializer = FieldInitializers(
6488
# steady-ish state initial condition for soil temperature
65-
temperature = (x,z) -> 10 + 0.02*z,
89+
temperature = (x,z) -> 10 - 0.02f0*z,
6690
# fully saturated soil
67-
saturation_water_ice = 1.0,
91+
saturation_water_ice = 1.0f0,
92+
)
93+
94+
# Coupled model with prescribed surface energy balance
95+
# currently not working :(
96+
# surface_energy_balance = SurfaceEnergyBalance(
97+
# eltype(grid);
98+
# radiative_fluxes = PrescribedRadiativeFluxes(),
99+
# turbulent_fluxes = PrescribedTurbulentFluxes(),
100+
# )
101+
# soil_model = SoilModel(grid, initializer=soil_initializer)
102+
# model = CoupledSoilAtmosphereModel(soil_model; surface_energy_balance)
103+
# # Quick check to make sure everything works OK
104+
# land_state = initialize(model)
105+
# set!(land_state.state.surface_shortwave_down, 200.0)
106+
# set!(land_state.state.surface_longwave_down, 100.0)
107+
# set!(land_state.state.air_temperature, 11.0)
108+
# set!(land_state.state.air_pressure, 101_325)
109+
# set!(land_state.state.specific_humidity, 0.001)
110+
# timestep!(land_state)
111+
# run!(land_state, period=Day(1))
112+
# @assert all(isfinite.(land_state.state.skin_temperature))
113+
# @assert all(isfinite.(land_state.state.latent_heat_flux))
114+
# @assert all(isfinite.(land_state.state.temperature))
115+
116+
# Soil model with prescribed surface temperautre BC
117+
soil_bcs = SoilBoundaryConditions(
118+
eltype(grid);
119+
top=PrescribedValue(:temperature, Input(:air_temperature, units=u"°C")),
120+
bottom=DefaultBoundaryConditions()
68121
)
69-
# Periodic surface temperature with annual cycle
70-
T_ub = PrescribedValue(:temperature, Input(:air_temperature, units=u"°C"))
71-
boundary_conditions = SoilBoundaryConditions(eltype(grid), top=T_ub)
72-
model = SoilModel(grid; initializer, boundary_conditions)
73-
state = initialize(model)
74-
75-
# Set up Speedy land and atmosphere model
76-
land = TerrariumDryLand(state)
77-
spectral_grid = land.spectral_grid #Speedy.SpectralGrid(grid.rings)
78-
land_sea_mask = Speedy.RockyPlanetMask(spectral_grid)
79-
output = Speedy.NetCDFOutput(spectral_grid, Speedy.PrimitiveDryModel, path="outputs/")
80-
primitive_dry = Speedy.PrimitiveDryModel(spectral_grid; land, land_sea_mask, output)
81-
Speedy.add!(primitive_dry.output, Speedy.SoilTemperatureOutput())
82-
sim = Speedy.initialize!(primitive_dry)
83-
Speedy.run!(sim, period=Day(30), output=true)
122+
model = SoilModel(grid, initializer=soil_initializer, boundary_conditions=soil_bcs)
123+
124+
# Initialize Terrarium model
125+
land_state = initialize(model)
126+
land = TerrariumDryLand(land_state)
127+
128+
# First set up and run a decoupled Speedy model for spinup
129+
land_sea_mask = Speedy.RockyPlanetMask(land.spectral_grid)
130+
primitive_dry_spinup = Speedy.PrimitiveDryModel(land.spectral_grid; land_sea_mask)
131+
sim_spinup = Speedy.initialize!(primitive_dry_spinup)
132+
Speedy.run!(sim_spinup, period=Day(5))
133+
134+
# Then run coupled model
135+
# surface_heat_flux = Speedy.SurfaceHeatFlux(land=Speedy.PrescribedLandHeatFlux())
136+
output = Speedy.NetCDFOutput(land.spectral_grid, Speedy.PrimitiveDryModel, path="outputs/")
137+
primitive_dry_coupled = Speedy.PrimitiveDryModel(land.spectral_grid; land, land_sea_mask, output)
138+
# add soil temperature as output variable for Speedy simulation
139+
Speedy.add!(primitive_dry_coupled.output, Speedy.SoilTemperatureOutput())
140+
# initialize coupled simulation
141+
sim_coupled = Speedy.Simulation(sim_spinup.prognostic_variables, sim_spinup.diagnostic_variables, primitive_dry_coupled)
142+
# Speedy.timestep!(sim)
143+
# spin up
144+
Speedy.run!(sim_coupled, period=Day(30), output=true)
145+
heatmap(sim.diagnostic_variables.grid.temp_grid[:,8])
146+
heatmap(sim.diagnostic_variables)
147+
148+
spectral_grid2 = Speedy.SpectralGrid(trunc=31)
149+
primitive_dry = Speedy.PrimitiveDryModel(spectral_grid2)
150+
sim2 = Speedy.initialize!(primitive_dry)
151+
Speedy.run!(sim2)
152+
84153
# animate surface layer air temperature
85154
Speedy.animate(sim, variable="temp", coastlines=false, level=spectral_grid.nlayers, output_file="plots/speedy_terrarium_dry_air_temperature.mp4")
86155
Speedy.animate(sim, variable="st", coastlines=false, level=1, output_file="plots/speedy_terrarium_dry_soil_temperature.mp4")
87156

88-
@show state.state.air_temperature
89-
@show state.state.temperature
157+
@show land_state.state.air_temperature
158+
@show land_state.state.skin_temperature
90159
@show sim.prognostic_variables.land.soil_temperature
91160

92161
heatmap(sim.diagnostic_variables.grid.temp_grid[:,end])
93162
heatmap(sim.prognostic_variables.land.soil_temperature[:,end])
94163

95-
@profview @btime Speedy.timestep!(sim.prognostic_variables, sim.diagnostic_variables, land, sim.model)
164+
using BenchmarkTools
165+
166+
@profview @btime Speedy.timestep!($sim)
167+
168+
default_spectral_grid = Speedy.SpectralGrid(grid.rings)
169+
default_primitive_dry = Speedy.PrimitiveDryModel(default_spectral_grid)
170+
default_sim = Speedy.initialize!(default_primitive_dry)
171+
Speedy.run!(default_sim, period=Day(1))

‎src/processes/prescribed_atmosphere.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -66,7 +66,7 @@ end
6666

6767
function PrescribedAtmosphere(
6868
::Type{NF};
69-
altitude::NF = 10.0, # Default to 10 m
69+
altitude::NF = NF(10), # Default to 10 m
7070
humidity::AbstractHumidity = SpecificHumidity(),
7171
precip::AbstractPrecipitation = TwoPhasePrecipitation(),
7272
radiation::AbstractIncomingRadiation = LongShortWaveRadiation(),

0 commit comments

Comments
 (0)