Skip to content

Commit

Permalink
start import udpates
Browse files Browse the repository at this point in the history
  • Loading branch information
juliasloan25 committed Apr 8, 2024
1 parent 3b97838 commit a86ae8b
Show file tree
Hide file tree
Showing 12 changed files with 102 additions and 94 deletions.
54 changes: 32 additions & 22 deletions experiments/AMIP/components/atmosphere/climaatmos.jl
Original file line number Diff line number Diff line change
@@ -1,13 +1,11 @@
# atmos_init: for ClimaAtmos pre-AMIP interface
using StaticArrays
using Statistics: mean
using LinearAlgebra: norm
# atmos_init: for ClimaAtmos interface
import StaticArrays as SA
import Statistics as Stats
import LinearAlgebra as LinAlg

import ClimaAtmos as CA
import ClimaAtmos: CT1, CT2, CT12, CT3, C3, C12, unit_basis_vector_data,
import ClimaCore
import SurfaceFluxes as SF
using ClimaCore
using ClimaCore.Utilities: half

import ClimaCoupler.Interfacer: AtmosModelSimulation
import ClimaCoupler.FluxCalculator:
Expand Down Expand Up @@ -103,10 +101,22 @@ function get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:radiative_energy_flux
(; face_lw_flux_dn, face_lw_flux_up, face_sw_flux_dn, face_sw_flux_up) =
atmos_sim.integrator.p.radiation.radiation_model

LWd_TOA = ClimaCore.Fields.level(CA.RRTMGPI.array2field(FT.(face_lw_flux_dn), face_space), nz_faces - half)
LWu_TOA = ClimaCore.Fields.level(CA.RRTMGPI.array2field(FT.(face_lw_flux_up), face_space), nz_faces - half)
SWd_TOA = ClimaCore.Fields.level(CA.RRTMGPI.array2field(FT.(face_sw_flux_dn), face_space), nz_faces - half)
SWu_TOA = ClimaCore.Fields.level(CA.RRTMGPI.array2field(FT.(face_sw_flux_up), face_space), nz_faces - half)
LWd_TOA = ClimaCore.Fields.level(
CA.RRTMGPI.array2field(FT.(face_lw_flux_dn), face_space),
nz_faces - ClimaCore.Utilities.half,
)
LWu_TOA = ClimaCore.Fields.level(
CA.RRTMGPI.array2field(FT.(face_lw_flux_up), face_space),
nz_faces - ClimaCore.Utilities.half,
)
SWd_TOA = ClimaCore.Fields.level(
CA.RRTMGPI.array2field(FT.(face_sw_flux_dn), face_space),
nz_faces - ClimaCore.Utilities.half,
)
SWu_TOA = ClimaCore.Fields.level(
CA.RRTMGPI.array2field(FT.(face_sw_flux_up), face_space),
nz_faces - ClimaCore.Utilities.half,
)

return @. -(LWd_TOA + SWd_TOA - LWu_TOA - SWu_TOA)
else
Expand Down Expand Up @@ -141,7 +151,7 @@ get_field(sim::ClimaAtmosSimulation, ::Val{:air_temperature}) =
TD.air_temperature.(thermo_params, sim.integrator.p.precomputed.ᶜts)
get_field(sim::ClimaAtmosSimulation, ::Val{:liquid_precipitation}) = sim.integrator.p.precipitation.col_integrated_rain
get_field(sim::ClimaAtmosSimulation, ::Val{:radiative_energy_flux_sfc}) =
ClimaCore.Fields.level(sim.integrator.p.radiation.ᶠradiation_flux, half)
ClimaCore.Fields.level(sim.integrator.p.radiation.ᶠradiation_flux, ClimaCore.Utilities.half)
get_field(sim::ClimaAtmosSimulation, ::Val{:snow_precipitation}) = sim.integrator.p.precipitation.col_integrated_snow
get_field(sim::ClimaAtmosSimulation, ::Val{:turbulent_energy_flux}) =
ClimaCore.Geometry.WVector.(sim.integrator.p.precomputed.sfc_conditions.ρ_flux_h_tot)
Expand All @@ -155,18 +165,18 @@ get_field(atmos_sim::ClimaAtmosSimulation, ::Val{:water}) = atmos_sim.integrator
get_field(sim::ClimaAtmosSimulation, ::Val{:height_int}) =
ClimaCore.Spaces.level(ClimaCore.Fields.coordinate_field(sim.integrator.u.c).z, 1)
get_field(sim::ClimaAtmosSimulation, ::Val{:height_sfc}) =
ClimaCore.Spaces.level(ClimaCore.Fields.coordinate_field(sim.integrator.u.f).z, half)
ClimaCore.Spaces.level(ClimaCore.Fields.coordinate_field(sim.integrator.u.f).z, ClimaCore.Utilities.half)
function get_field(sim::ClimaAtmosSimulation, ::Val{:uv_int})
uₕ_int = ClimaCore.Geometry.UVVector.(ClimaCore.Spaces.level(sim.integrator.u.c.uₕ, 1))
return @. StaticArrays.SVector(uₕ_int.components.data.:1, uₕ_int.components.data.:2)
return @. SA.SVector(uₕ_int.components.data.:1, uₕ_int.components.data.:2)
end

function update_field!(atmos_sim::ClimaAtmosSimulation, ::Val{:co2}, field)
if atmos_sim.integrator.p.atmos.radiation_mode isa CA.RRTMGPI.GrayRadiation
@warn("Gray radiation model initialized, skipping CO2 update", maxlog = 1)
return
else
atmos_sim.integrator.p.radiation.radiation_model.volume_mixing_ratio_co2 .= mean(parent(field))
atmos_sim.integrator.p.radiation.radiation_model.volume_mixing_ratio_co2 .= Stats.mean(parent(field))
end
end
# extensions required by the Interfacer
Expand All @@ -189,15 +199,15 @@ function update_field!(sim::ClimaAtmosSimulation, ::Val{:turbulent_fluxes}, fiel

Y = sim.integrator.u
surface_local_geometry = ClimaCore.Fields.level(ClimaCore.Fields.local_geometry_field(Y.f), ClimaCore.Fields.half)
surface_normal = @. C3(unit_basis_vector_data(C3, surface_local_geometry))
surface_normal = @. CA.C3(CA.unit_basis_vector_data(CA.C3, surface_local_geometry))

# get template objects for the contravariant components of the momentum fluxes (required by Atmos boundary conditions)
vec_ct12_ct1 = @. CT12(CT2(unit_basis_vector_data(CT1, surface_local_geometry)), surface_local_geometry)
vec_ct12_ct2 = @. CT12(CT2(unit_basis_vector_data(CT2, surface_local_geometry)), surface_local_geometry)
vec_ct12_ct1 = @. CA.CT12(CA.CT2(CA.unit_basis_vector_data(CA.CT1, surface_local_geometry)), surface_local_geometry)
vec_ct12_ct2 = @. CA.CT12(CA.CT2(CA.unit_basis_vector_data(CA.CT2, surface_local_geometry)), surface_local_geometry)

sim.integrator.p.precomputed.sfc_conditions.ρ_flux_uₕ .= (
surface_normal .⊗
C12.(
surface_normal.CA .⊗
CA.C12.(
swap_space!(ones(axes(vec_ct12_ct1)), F_turb_ρτxz) .* vec_ct12_ct1 .+
swap_space!(ones(axes(vec_ct12_ct2)), F_turb_ρτyz) .* vec_ct12_ct2,
surface_local_geometry,
Expand Down Expand Up @@ -421,7 +431,7 @@ function water_albedo_from_atmosphere!(


# set the direct and diffuse surface albedos
direct_albedo .= CA.surface_albedo_direct(α_model).(λ, μ, norm.(ClimaCore.Fields.level(Y.c.uₕ, 1)))
diffuse_albedo .= CA.surface_albedo_diffuse(α_model).(λ, μ, norm.(ClimaCore.Fields.level(Y.c.uₕ, 1)))
direct_albedo .= CA.surface_albedo_direct(α_model).(λ, μ, LinAlg.norm.(ClimaCore.Fields.level(Y.c.uₕ, 1)))
diffuse_albedo .= CA.surface_albedo_diffuse(α_model).(λ, μ, LinAlg.norm.(ClimaCore.Fields.level(Y.c.uₕ, 1)))

end
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@
# these extensions add extra diagnostics to the atmos model output
import ClimaAtmos.Diagnostics as CAD

import ClimaAtmos.Diagnostics: add_diagnostic_variable!

"""
Expand Down
69 changes: 33 additions & 36 deletions experiments/AMIP/components/land/climaland_bucket.jl
Original file line number Diff line number Diff line change
@@ -1,22 +1,19 @@
# slab_rhs!
using ClimaCore
import ClimaTimeSteppers as CTS
import Dates
import Thermodynamics as TD
using Dates: DateTime
using ClimaComms: AbstractCommsContext

import ClimaCore
import ClimaTimeSteppers as CTS
import ClimaParams

import ClimaLand
using ClimaLand.Bucket: BucketModel, BucketModelParameters, AbstractAtmosphericDrivers, AbstractRadiativeDrivers
import ClimaLand.Bucket: PrescribedSurfaceAlbedo, PrescribedBaregroundAlbedo
using ClimaLand:
make_exp_tendency,
initialize,
make_set_initial_cache,
surface_evaporative_scaling,
CoupledRadiativeFluxes,
CoupledAtmosphere
import ClimaLand.Parameters as LP
import ClimaLand as CL
import CL.Bucket:
BucketModel,
BucketModelParameters,
AbstractAtmosphericDrivers,
AbstractRadiativeDrivers,
PrescribedSurfaceAlbedo,
PrescribedBaregroundAlbedo
import CL.Parameters as LP

import ClimaCoupler.Interfacer: LandModelSimulation, get_field, update_field!, name, step!, reinit!
import ClimaCoupler.FluxCalculator: update_turbulent_fluxes_point!, surface_thermo_state
Expand Down Expand Up @@ -56,7 +53,7 @@ function bucket_init(
saveat::Float64,
area_fraction,
stepper = CTS.RK4(),
date_ref::DateTime,
date_ref::Dates.DateTime,
t_start::Float64,
) where {FT}
if config != "sphere"
Expand Down Expand Up @@ -100,11 +97,11 @@ function bucket_init(
# Note that this does not take into account topography of the surface, which is OK for this land model.
# But it must be taken into account when computing surface fluxes, for Δz.
domain = make_land_domain(space, (-d_soil, FT(0.0)), n_vertical_elements)
args = (params, CoupledAtmosphere{FT}(), CoupledRadiativeFluxes{FT}(), domain)
args = (params, CL.CoupledAtmosphere{FT}(), CL.CoupledRadiativeFluxes{FT}(), domain)
model = BucketModel{FT, typeof.(args)...}(args...)

# Initial conditions with no moisture
Y, p, coords = initialize(model)
Y, p, coords = CL.initialize(model)

# Get temperature anomaly function
T_functions = Dict("aquaplanet" => temp_anomaly_aquaplanet, "amip" => temp_anomaly_amip)
Expand All @@ -121,12 +118,12 @@ function bucket_init(
Y.bucket.σS .= 0.0

# Set initial aux variable values
set_initial_cache! = make_set_initial_cache(model)
set_initial_cache! = CL.make_set_initial_cache(model)
set_initial_cache!(p, Y, tspan[1])

exp_tendency! = make_exp_tendency(model)
exp_tendency! = CL.make_exp_tendency(model)
ode_algo = CTS.ExplicitAlgorithm(stepper)
bucket_ode_function = CTS.ClimaODEFunction(T_exp! = exp_tendency!, dss! = ClimaLand.dss!)
bucket_ode_function = CTS.ClimaODEFunction(T_exp! = exp_tendency!, dss! = CL.dss!)
prob = ODEProblem(bucket_ode_function, Y, tspan, p)
integrator = init(prob, ode_algo; dt = dt, saveat = saveat, adaptive = false)

Expand All @@ -141,17 +138,17 @@ end
get_field(sim::BucketSimulation, ::Val{:air_density}) = sim.integrator.p.bucket.ρ_sfc
get_field(sim::BucketSimulation, ::Val{:area_fraction}) = sim.area_fraction
get_field(sim::BucketSimulation, ::Val{:beta}) =
ClimaLand.surface_evaporative_scaling(sim.model, sim.integrator.u, sim.integrator.p)
CL.surface_evaporative_scaling(sim.model, sim.integrator.u, sim.integrator.p)
get_field(sim::BucketSimulation, ::Val{:roughness_buoyancy}) = sim.model.parameters.z_0b
get_field(sim::BucketSimulation, ::Val{:roughness_momentum}) = sim.model.parameters.z_0m
get_field(sim::BucketSimulation, ::Val{:surface_direct_albedo}) =
ClimaLand.surface_albedo(sim.model, sim.integrator.u, sim.integrator.p)
CL.surface_albedo(sim.model, sim.integrator.u, sim.integrator.p)
get_field(sim::BucketSimulation, ::Val{:surface_diffuse_albedo}) =
ClimaLand.surface_albedo(sim.model, sim.integrator.u, sim.integrator.p)
CL.surface_albedo(sim.model, sim.integrator.u, sim.integrator.p)
get_field(sim::BucketSimulation, ::Val{:surface_humidity}) =
ClimaLand.surface_specific_humidity(sim.model, sim.integrator.u, sim.integrator.p, sim.integrator.t)
CL.surface_specific_humidity(sim.model, sim.integrator.u, sim.integrator.p, sim.integrator.t)
get_field(sim::BucketSimulation, ::Val{:surface_temperature}) =
ClimaLand.surface_temperature(sim.model, sim.integrator.u, sim.integrator.p, sim.integrator.t)
CL.surface_temperature(sim.model, sim.integrator.u, sim.integrator.p, sim.integrator.t)

"""
get_field(bucket_sim::BucketSimulation, ::Val{:energy})
Expand Down Expand Up @@ -179,7 +176,7 @@ end
Extension of Interfacer.get_field that provides the total water contained in the bucket, including the liquid water in snow.
"""
function get_field(bucket_sim::BucketSimulation, ::Val{:water})
ρ_cloud_liq = ClimaLand.LP.ρ_cloud_liq(bucket_sim.model.parameters.earth_param_set)
ρ_cloud_liq = CL.LP.ρ_cloud_liq(bucket_sim.model.parameters.earth_param_set)
return
@. (bucket_sim.integrator.u.bucket.σS + bucket_sim.integrator.u.bucket.W + bucket_sim.integrator.u.bucket.Ws) *
ρ_cloud_liq # kg water / m2
Expand Down Expand Up @@ -248,27 +245,27 @@ function get_model_prog_state(sim::BucketSimulation)
end

###
### ClimaLand.jl bucket model-specific functions (not explicitly required by ClimaCoupler.jl)
### CL.jl bucket model-specific functions (not explicitly required by ClimaCoupler.jl)
###

# TODO remove this function after ClimaLand v0.8.1 update
function ClimaLand.turbulent_fluxes(atmos::CoupledAtmosphere, model::BucketModel, Y, p, t)
function CL.turbulent_fluxes(atmos::CL.CoupledAtmosphere, model::BucketModel, Y, p, t)
# coupler has done its thing behind the scenes already
model_name = ClimaLand.name(model)
model_name = CL.name(model)
model_cache = getproperty(p, model_name)
return model_cache.turbulent_fluxes
end


function ClimaLand.initialize_drivers(a::CoupledAtmosphere{FT}, coords) where {FT}
function CL.initialize_drivers(a::CL.CoupledAtmosphere{FT}, coords) where {FT}
keys = (:P_liq, :P_snow)
types = ([FT for k in keys]...,)
domain_names = ([:surface for k in keys]...,)
model_name = :drivers
# intialize_vars packages the variables as a named tuple,
# as part of a named tuple with `model_name` as the key.
# Here we just want the variable named tuple itself
vars = ClimaLand.initialize_vars(keys, types, domain_names, coords, model_name)
vars = CL.initialize_vars(keys, types, domain_names, coords, model_name)
return vars.drivers
end

Expand Down Expand Up @@ -327,7 +324,7 @@ function make_land_domain(
subsurface_space = ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace(atmos_boundary_space, vert_center_space)
space = (; surface = atmos_boundary_space, subsurface = subsurface_space)

return ClimaLand.Domains.SphericalShell{FT}(radius, depth, nothing, nelements, npolynomial, space)
return CL.Domains.SphericalShell{FT}(radius, depth, nothing, nelements, npolynomial, space)
end

"""
Expand All @@ -336,7 +333,7 @@ Returns the surface temperature of the earth, computed from the state u.
"""
function get_land_temp_from_state(land_sim, u)
# required by viz_explorer.jl
return ClimaLand.surface_temperature(land_sim.model, u, land_sim.integrator.p, land_sim.integrator.t)
return CL.surface_temperature(land_sim.model, u, land_sim.integrator.p, land_sim.integrator.t)
end

"""
Expand All @@ -349,5 +346,5 @@ or 3D dss buffer stored in the cache depending on space of each variable in
`sim.integrator.u`.
"""
function dss_state!(sim::BucketSimulation)
ClimaLand.dss!(sim.integrator.u, sim.integrator.p, sim.integrator.t)
CL.dss!(sim.integrator.u, sim.integrator.p, sim.integrator.t)
end
9 changes: 4 additions & 5 deletions experiments/AMIP/components/ocean/eisenman_seaice.jl
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
import SciMLBase: ODEProblem, init
import SciMLBase

using ClimaCore
using ClimaCore.Fields: getindex
import ClimaCore
import ClimaTimeSteppers as CTS

import ClimaCoupler: FluxCalculator
Expand Down Expand Up @@ -82,8 +81,8 @@ function eisenman_seaice_init(
thermo_params = thermo_params,
dss_buffer = ClimaCore.Spaces.create_dss_buffer(ClimaCore.Fields.zeros(space)),
)
problem = ODEProblem(ode_function, Y, Float64.(tspan), cache)
integrator = init(problem, ode_algo, dt = Float64(dt), saveat = Float64(saveat), adaptive = false)
problem = SciMLBase.ODEProblem(ode_function, Y, Float64.(tspan), cache)
integrator = SciMLBase.init(problem, ode_algo, dt = Float64(dt), saveat = Float64(saveat), adaptive = false)

sim = EisenmanIceSimulation(params, Y, space, integrator)
@warn name(sim) *
Expand Down
10 changes: 5 additions & 5 deletions experiments/AMIP/components/ocean/prescr_seaice.jl
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
import SciMLBase: ODEProblem, init
import SciMLBase

using ClimaCore
import ClimaCore
import ClimaTimeSteppers as CTS
import Thermodynamics as TD

import ClimaCoupler.Interfacer: SeaIceModelSimulation, get_field, update_field!, name, step!, reinit!
import ClimaCoupler.FluxCalculator: update_turbulent_fluxes_point!
using ClimaCoupler: Regridder
import ClimaCoupler: Regridder
import ClimaCoupler.Utilities: swap_space!
import ClimaCoupler.BCReader: float_type_bcf

Expand Down Expand Up @@ -88,8 +88,8 @@ function ice_init(::Type{FT}; tspan, saveat, dt, space, area_fraction, thermo_pa
ode_algo = CTS.ExplicitAlgorithm(stepper)
ode_function = CTS.ClimaODEFunction(T_exp! = ice_rhs!, dss! = weighted_dss_slab!)

problem = ODEProblem(ode_function, Y, Float64.(tspan), (; additional_cache..., params = params))
integrator = init(problem, ode_algo, dt = Float64(dt), saveat = Float64(saveat), adaptive = false)
problem = SciMLBase.ODEProblem(ode_function, Y, Float64.(tspan), (; additional_cache..., params = params))
integrator = SciMLBase.init(problem, ode_algo, dt = Float64(dt), saveat = Float64(saveat), adaptive = false)

sim = PrescribedIceSimulation(params, Y, space, integrator)

Expand Down
8 changes: 4 additions & 4 deletions experiments/AMIP/components/ocean/slab_ocean.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
import SciMLBase: ODEProblem, init
import SciMLBase

using ClimaCore
import ClimaCore
import ClimaTimeSteppers as CTS
import ClimaCoupler.Interfacer: OceanModelSimulation, get_field, update_field!, name, step!, reinit!
import ClimaCoupler.FluxCalculator: update_turbulent_fluxes_point!
Expand Down Expand Up @@ -98,8 +98,8 @@ function ocean_init(
ode_algo = CTS.ExplicitAlgorithm(stepper)
ode_function = CTS.ClimaODEFunction(T_exp! = slab_ocean_rhs!, dss! = weighted_dss_slab!)

problem = ODEProblem(ode_function, Y, Float64.(tspan), cache)
integrator = init(problem, ode_algo, dt = Float64(dt), saveat = Float64(saveat), adaptive = false)
problem = SciMLBase.ODEProblem(ode_function, Y, Float64.(tspan), cache)
integrator = SciMLBase.init(problem, ode_algo, dt = Float64(dt), saveat = Float64(saveat), adaptive = false)

sim = SlabOceanSimulation(params, Y, space, integrator)

Expand Down
3 changes: 2 additions & 1 deletion experiments/AMIP/user_io/amip_visualizer.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
import ClimaComms
using ClimaCore
import ClimaCore

import ClimaCoupler.PostProcessor: postprocess

include("plot_helper.jl")
Expand Down
7 changes: 3 additions & 4 deletions experiments/AMIP/user_io/debug_plots.jl
Original file line number Diff line number Diff line change
@@ -1,8 +1,7 @@
using Plots
using ClimaCorePlots
using Printf
import ClimaCore
import Plots

using ClimaCoupler.Interfacer: ComponentModelSimulation, SurfaceModelSimulation
using ClimaCore

# plotting functions for the coupled simulation
"""
Expand Down
7 changes: 4 additions & 3 deletions experiments/AMIP/user_io/ncep_visualizer.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
using Downloads
using NCDatasets
import Downloads
import NCDatasets as NCD

import ClimaCoupler.Diagnostics: get_var
import ClimaCoupler.PostProcessor: postprocess

Expand Down Expand Up @@ -88,7 +89,7 @@ Downloads and reads nc datafile of a specified NCEP variable
function download_read_nc(data_source::NCEPMonthlyDataSource, https::String, ncep_vname::String)
local_file = joinpath(data_source.tmp_dir, ncep_vname * ".nc")
Downloads.download(https, local_file)
NCDataset(local_file) do ds
NCD.NCDataset(local_file) do ds
t_i = findall(x -> Dates.yearmonth(x) == Dates.yearmonth(data_source.month_date[1]), Array(ds["time"])) # time index of month in file
d_i = length(size(Array(ds[ncep_vname]))) # index of time in the dimension list
lev = "level" in keys(ds) ? Array(ds["level"]) : [Float64(-999)]
Expand Down
Loading

0 comments on commit a86ae8b

Please sign in to comment.