diff --git a/Project.toml b/Project.toml index 35bbe9e..872392f 100644 --- a/Project.toml +++ b/Project.toml @@ -13,8 +13,11 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" Git = "d7ba0133-e1db-5d97-8f8c-041e4b3a1eb2" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" @@ -26,6 +29,7 @@ Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d" IniFile = "83e8ac13-25f8-5344-8a64-a9f2b223428f" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9" +Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" Zarr = "0a941bbe-ad1d-11e8-39d9-ab76183a1d99" @@ -34,21 +38,24 @@ ClimateModelsCondaExt = ["Conda"] ClimateModelsIniFileExt = ["IniFile"] ClimateModelsMakieExt = ["Makie"] ClimateModelsNetCDFExt = ["NetCDF"] +ClimateModelsOceananigansExt = ["Oceananigans"] ClimateModelsPyCallExt = ["PyCall"] ClimateModelsZarrExt = ["Zarr"] [compat] +CFTime = "0.1" CSV = "0.6, 0.7, 0.8, 0.9, 0.10" Conda = "1" -CFTime = "0.1" DataDeps = "0.7" DataFrames = "1" Dataverse = "0.2" Git = "1" Glob = "1" IniFile = "0.5" +JLD2 = "0.4, 0.5" Makie = "0.21" NetCDF = "0.11, 0.12" +Oceananigans = "0.90" OrderedCollections = "1" PyCall = "1" Suppressor = "0.2" diff --git a/ext/ClimateModelsMakieExt.jl b/ext/ClimateModelsMakieExt.jl index a1a9155..ca81b94 100644 --- a/ext/ClimateModelsMakieExt.jl +++ b/ext/ClimateModelsMakieExt.jl @@ -22,6 +22,8 @@ function plot_examples(ID=Symbol,stuff...) Speedy_plot_input(stuff...) elseif ID==:Speedy_zm Speedy_plot_output_zm(stuff...) + elseif ID==:Speedy_xy + Speedy_plot_output_xy(stuff...) elseif ID==:Hector plot_Hector(stuff...) elseif ID==:Hector_scenarios diff --git a/ext/ClimateModelsNetCDFExt.jl b/ext/ClimateModelsNetCDFExt.jl index cb2c69d..7ef0680 100644 --- a/ext/ClimateModelsNetCDFExt.jl +++ b/ext/ClimateModelsNetCDFExt.jl @@ -1,9 +1,9 @@ module ClimateModelsNetCDFExt using NetCDF - import ClimateModels: write_CMIP6_mean, read_CMIP6_mean, ModelConfig + import ClimateModels: write_CMIP6_mean, read_CMIP6_mean, AbstractModelConfig, read_NetCDF - function write_CMIP6_mean(x::ModelConfig,mm,meta) + function write_CMIP6_mean(x::AbstractModelConfig,mm,meta) pth=joinpath(pathof(x),"output") filename = joinpath(pth,"MeanMaps.nc") varname = x.inputs["variable_id"] @@ -21,4 +21,9 @@ module ClimateModelsNetCDFExt return lon,lat,var end + function read_NetCDF(file) + ncfile = NetCDF.open(file) + ncfile.vars + end + end diff --git a/ext/ClimateModelsOceananigansExt.jl b/ext/ClimateModelsOceananigansExt.jl new file mode 100644 index 0000000..090b8e2 --- /dev/null +++ b/ext/ClimateModelsOceananigansExt.jl @@ -0,0 +1,109 @@ + +module ClimateModelsOceananigansExt + +using Oceananigans, ClimateModels +import ClimateModels: JLD2, @printf, Oceananigans_launch +import ClimateModels: Oceananigans_setup_grid, Oceananigans_setup_BC +import ClimateModels: Oceananigans_build_model, Oceananigans_build_simulation + +import Oceananigans.Units: minute, minutes, hour + +function Oceananigans_setup_grid(Nz=50,Lz=50) + #Nz = 50 # number of points in the vertical direction + #Lz = 50 # (m) domain depth + fz(k)=-Lz*(Nz+1-k)/Nz #fz.(1:Nz+1) gives the vertical grid for w points + return RectilinearGrid(size = (32, 32, Nz), x = (0, 64), y = (0, 64), z = fz) +end + +function Oceananigans_setup_BC(Qʰ,u₁₀,Ev) + dTdz = 0.1 # K m⁻¹, bottom boundary condition + ρₒ = 1026 # kg m⁻³, average density at the surface of the world ocean + cᴾ = 3991 # J K⁻¹ kg⁻¹, typical heat capacity for seawater + Qᵀ(x, y, t) = Qʰ(t) / (ρₒ * cᴾ) # K m s⁻¹, surface _temperature_ flux + T_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵀ), + bottom = GradientBoundaryCondition(dTdz)) + + cᴰ = 2.5e-3 # dimensionless drag coefficient + ρₐ = 1.225 # kg m⁻³, average density of air at sea-level + Qᵘ(x, y, t) = - ρₐ / ρₒ * cᴰ * u₁₀(t) * abs(u₁₀(t)) # m² s⁻² + u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ)) + + Qˢ(x, y, t, S) = - Ev(t) * S # [salinity unit] m s⁻¹ + evaporation_bc = FluxBoundaryCondition(Qˢ, field_dependencies=:S) + S_bcs = FieldBoundaryConditions(top=evaporation_bc) + + return (u=u_bcs,T=T_bcs,S=S_bcs) +end + +## + +function Oceananigans_build_model(grid,BC,IC) + + buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion=2e-4, haline_contraction=8e-4)) + + model = NonhydrostaticModel( + advection = UpwindBiasedFifthOrder(), + timestepper = :RungeKutta3, + grid = grid, + tracers = (:T, :S), + coriolis = FPlane(f=1e-4), + buoyancy = buoyancy, + closure = AnisotropicMinimumDissipation(), + boundary_conditions = (u=BC.u, T=BC.T, S=BC.S)) + + # initial conditions (as functions of x,y,z) + set!(model, u=IC.u, w=IC.u, T=IC.T, S=IC.S) + + return model +end + +function Oceananigans_build_simulation(model,Nh,rundir) + simulation = Simulation(model, Δt=10.0, stop_time=Nh*60minutes) + + wizard = TimeStepWizard(cfl=1.0, max_change=1.1, max_Δt=20.0) + simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) + + progress_message(sim) = + @printf("Iteration: %04d, time: %s, Δt: %s, max(|w|) = %.1e ms⁻¹, wall time: %s\n", + iteration(sim),prettytime(sim),prettytime(sim.Δt), + maximum(abs, sim.model.velocities.w),prettytime(sim.run_wall_time)) + simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(1minute)) + + eddy_viscosity = (; νₑ = model.diffusivity_fields.νₑ) + simulation.output_writers[:slices] = + JLD2OutputWriter(model, merge(model.velocities, model.tracers, eddy_viscosity), + dir = rundir, + filename = "daily_cycle.jld2", + indices = (:,Int(model.grid.Ny/2),:), + schedule = TimeInterval(1minute), + overwrite_existing = true) + + simulation.output_writers[:checkpointer] = + Checkpointer(model, schedule=TimeInterval(24hour), dir = rundir, prefix="model_checkpoint") + + ## + + fil=simulation.output_writers[:slices].filepath + + xw, yw, zw = nodes(model.velocities.w) + xT, yT, zT = nodes(model.tracers.T) + coords=(xw, yw, zw, xT, yT, zT) + + fil_coords=joinpath(rundir,"coords.jld2") + JLD2.jldopen(fil_coords, "w") do file + mygroup = JLD2.Group(file, "coords") + mygroup["xw"] = xw + mygroup["yw"] = yw + mygroup["zw"] = zw + mygroup["xT"] = xT + mygroup["yT"] = yT + mygroup["zT"] = zT + end + + return simulation +end + +Oceananigans_launch(x::OceananigansConfig) = run!(x.outputs["simulation"], pickup=true) + +end + diff --git a/src/ClimateModels.jl b/src/ClimateModels.jl index b3ac130..4899d63 100644 --- a/src/ClimateModels.jl +++ b/src/ClimateModels.jl @@ -1,13 +1,19 @@ module ClimateModels using UUIDs, Suppressor, OrderedCollections -using Pkg, Git, TOML +using Pkg, Git, TOML, Printf, JLD2 function plot_examples end function read_Zarr end function write_CMIP6_mean end function read_CMIP6_mean end function read_IniFile end +function read_NetCDF end +function Oceananigans_setup_grid end +function Oceananigans_setup_BC end +function Oceananigans_build_model end +function Oceananigans_build_simulation end +function Oceananigans_launch end include("interface.jl") include("notebooks.jl") @@ -15,20 +21,31 @@ include("toy_models.jl") include("files.jl") include("CMIP6.jl") include("Hector.jl") +include("FaIR.jl") +include("Speedy.jl") +include("Oceananigans.jl") import .notebooks: update import .downloads: add_datadep import .Hector: HectorConfig - -export AbstractModelConfig, ModelConfig, PlutoConfig, HectorConfig -export ModelRun, @ModelRun, PkgDevConfig, add_datadep, read_Zarr, read_IniFile +import .FaIR: FaIRConfig +import .Speedy: SpeedyConfig +import .Oceananigans: OceananigansConfig + +export AbstractModelConfig, ModelConfig, PlutoConfig +export HectorConfig, FaIRConfig, SpeedyConfig, OceananigansConfig +export ModelRun, @ModelRun, PkgDevConfig, add_datadep +export read_Zarr, read_IniFile, read_NetCDF export clean, build, compile, setup, launch, update, notebooks export put!, take!, pathof, readdir, log #export git_log_init, git_log_msg, git_log_fil #export git_log_prm, git_log_show #export monitor, help, pause #export train, compare, analyze -export RandomWalker, IPCC, CMIP6, Hector + +export RandomWalker, IPCC, CMIP6, Hector, FaIR, Speedy +#don't export Oceananigans module, which would conflict with Oceananigans.jl +#export Oceananigans #export OrderedDict, UUID, uuid4, @suppress diff --git a/src/FaIR.jl b/src/FaIR.jl new file mode 100644 index 0000000..409eb84 --- /dev/null +++ b/src/FaIR.jl @@ -0,0 +1,76 @@ + +module FaIR + +import ClimateModels +import ClimateModels: setup, AbstractModelConfig + +uuid4=ClimateModels.uuid4 +UUID=ClimateModels.UUID +OrderedDict=ClimateModels.OrderedDict + +function loop_over_scenarios() + scenarios=[:rcp26,:rcp45,:rcp60,:rcp85] + temperatures=[] + + (fair,forward,RCPs)=ClimateModels.pyimport(:fair) + + for i in scenarios + emissions=RCPs[i].Emissions.emissions + C,F,T = forward.fair_scm(emissions=emissions) + push!(temperatures,T) + end + + scenarios,temperatures +end + +""" + struct FaIRConfig <: AbstractModelConfig + +Concrete type of `AbstractModelConfig` for `FaIR` model. +""" +Base.@kwdef struct FaIRConfig <: AbstractModelConfig + model :: String = "FaIR" + configuration :: String = "rcp45" + inputs :: OrderedDict{Any,Any} = OrderedDict{Any,Any}() + outputs :: OrderedDict{Any,Any} = OrderedDict(:C=>Float64[],:F=>Float64[],:T=>Float64[]) + status :: OrderedDict{Any,Any} = OrderedDict{Any,Any}() + channel :: Channel{Any} = Channel{Any}(10) + folder :: String = tempdir() + ID :: UUID = uuid4() +end + +function setup(x :: FaIRConfig) + !isdir(x.folder) ? mkdir(x.folder) : nothing + !isdir(pathof(x)) ? mkdir(pathof(x)) : nothing + try + ClimateModels.pyimport(:fair) + catch + ClimateModels.conda(:fair) + end + !isdir(joinpath(pathof(x),"log")) ? log(x,"initial setup",init=true) : nothing + put!(x.channel,FaIR_launch) +end + +function FaIR_launch(x::FaIRConfig) + pth0=pwd() + cd(pathof(x)) + + (fair,forward,RCPs)=ClimateModels.pyimport(:fair) + + emissions=RCPs.rcp85.Emissions.emissions + C,F,T = forward.fair_scm(emissions=emissions) + + if isempty(x.outputs[:C]) + push!.(Ref(x.outputs[:C]),C[:]) + push!.(Ref(x.outputs[:F]),F[:]) + push!.(Ref(x.outputs[:T]),T[:]) + else + x.outputs[:C]=C[:] + x.outputs[:F]=F[:] + x.outputs[:T]=T[:] + end + + cd(pth0) +end + +end diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl new file mode 100644 index 0000000..85b3b21 --- /dev/null +++ b/src/Oceananigans.jl @@ -0,0 +1,203 @@ + +module Oceananigans + +using Printf, Random, JLD2, Statistics, Downloads + +using ClimateModels +import ClimateModels: build, setup, AbstractModelConfig, Oceananigans_launch +import ClimateModels: Oceananigans_setup_grid, Oceananigans_setup_BC +import ClimateModels: Oceananigans_build_model, Oceananigans_build_simulation + +OrderedDict=ClimateModels.OrderedDict +uuid4=ClimateModels.uuid4 +UUID=ClimateModels.UUID + +## + +function setup_initial_conditions(Tᵢ) + # Temperature initial condition: + T(x, y, z) = Tᵢ(z) + # Velocity initial condition: + u(x, y, z) = randn()*1e-5 + # Salinity initial condition: + S(x, y, z)=35.0 + + return (T=T,S=S,u=u) +end + +## + +function read_grid(MC) + fil_coords=joinpath(pathof(MC),"coords.jld2") + file = jldopen(fil_coords) + coords = load(fil_coords) + xw = file["coords"]["xw"] + yw = file["coords"]["yw"] + zw = file["coords"]["zw"] + xT = file["coords"]["xT"] + yT = file["coords"]["yT"] + zT = file["coords"]["zT"] + close(file) + return xw, yw, zw, xT, yT, zT +end + +# ╔═╡ e3453716-b8db-449a-a3bb-c918af91878e +function xz_read(fil,t) + # Open the file with our data + file = jldopen(fil) + + # Extract a vector of iterations + iterations = parse.(Int, keys(file["timeseries/t"])) + times = [file["timeseries/t/$iter"] for iter in iterations] + + iter=iterations[t] + + t = file["timeseries/t/$iter"] + w = file["timeseries/w/$iter"][:, 1, :] + T = file["timeseries/T/$iter"][:, 1, :] + S = file["timeseries/S/$iter"][:, 1, :] + νₑ = file["timeseries/νₑ/$iter"][:, 1, :] + + close(file) + + return t,w,T,S,νₑ +end + +function zt_read(fil,t) + # Open the file with our data + file = jldopen(fil) + + # Extract a vector of iterations + iterations = parse.(Int, keys(file["timeseries/t"])) + times = [file["timeseries/t/$iter"] for iter in iterations] + + iter=iterations[t] + + t = file["timeseries/t/$iter"] + w = sqrt.(mean(file["timeseries/w/$iter"][:, :, :].^2, dims=(1,2)))[:] + T = mean(file["timeseries/T/$iter"][:, :, :], dims=(1,2))[:] + S = mean(file["timeseries/S/$iter"][:, :, :], dims=(1,2))[:] + νₑ = sqrt.(mean(file["timeseries/νₑ/$iter"][:, :, :].^2, dims=(1,2)))[:] + + close(file) + + return t,w,T,S,νₑ +end + +function xz_plot_prep(MC,i) + fil=joinpath(pathof(MC),"daily_cycle.jld2") + t,w,T,S,νₑ=xz_read(fil,i) + xw, yw, zw, xT, yT, zT=read_grid(MC) + tt="$(round(t/86400)) days" + #tt=prettytime(t) + (tt,w,T,S,νₑ,xw, yw, zw, xT, yT, zT) +end + +function tz_slice(MC;nt=1,wli=missing,Tli=missing,Sli=missing,νli=missing) + xw, yw, zw, xT, yT, zT=read_grid(MC) + + fil=joinpath(pathof(MC),"daily_cycle.jld2") + Tall=Matrix{Float64}(undef,length(zT),nt) + Sall=Matrix{Float64}(undef,length(zT),nt) + wall=Matrix{Float64}(undef,length(zw),nt) + νₑall=Matrix{Float64}(undef,length(zT),nt) + for tt in 1:nt + t,w,T,S,νₑ=zt_read(fil,tt) + Tall[:,tt]=T + Sall[:,tt]=S + wall[:,tt]=w + νₑall[:,tt]=νₑ + end + + permutedims(Tall),permutedims(Sall),permutedims(wall),permutedims(νₑall) +end + +function nt_from_jld2(MC) + fil=joinpath(pathof(MC),"daily_cycle.jld2") + file = jldopen(fil) + iterations = parse.(Int, keys(file["timeseries/t"])) + times = [file["timeseries/t/$iter"] for iter in iterations] + close(file) + nt=(length(times)) + + return nt +end + +## + +""" + OceananigansConfig() + +Concrete type of `AbstractModelConfig` for `Oceananigans.jl` +""" +Base.@kwdef struct OceananigansConfig <: AbstractModelConfig + model :: String = "Oceananigans" + configuration :: String = "daily_cycle" + inputs :: OrderedDict{Any,Any} = OrderedDict{Any,Any}() + outputs :: OrderedDict{Any,Any} = OrderedDict{Any,Any}() + status :: OrderedDict{Any,Any} = OrderedDict{Any,Any}() + channel :: Channel{Any} = Channel{Any}(10) + folder :: String = tempdir() + ID :: UUID = uuid4() +end + +function build(x::OceananigansConfig) + rundir=pathof(x) + model=Oceananigans_build_model(x.outputs["grid"],x.outputs["BC"],x.outputs["IC"]) + simulation=Oceananigans_build_simulation(model,x.inputs["Nh"],rundir) + + x.outputs["model"]=model + x.outputs["simulation"]=simulation + + return true +end + +function rerun(x::OceananigansConfig) + simulation=demo.build_simulation(x.outputs["model"],x.inputs["Nh"],pathof(x)) + x.outputs["simulation"]=simulation + Oceananigans_launch(x) +end + +function setup(x::OceananigansConfig) + + if x.configuration=="daily_cycle" + Qʰ(t) = 200.0 * (1.0-2.0*(mod(t,86400.0)>43200.0)) # W m⁻², surface heat flux (>0 means ocean cooling) + u₁₀(t) = 4.0 * (1.0-0.9*(mod(t,86400.0)>43200.0)) # m s⁻¹, wind speed 10 meters above ocean surface + Ev(t) = 1e-7 * (1.0-2.0*(mod(t,86400.0)>43200.0)) # m s⁻¹, evaporation rate + Tᵢ(z) = 20 + 0.1 * z #initial temperature condition (function of z=-depth) + + grid=Oceananigans_setup_grid() + IC=setup_initial_conditions(Tᵢ) + BC=Oceananigans_setup_BC(Qʰ,u₁₀,Ev) + else + grid=missing + IC=missing + BC=missing + error("unnknown model configuration") + end + + !isdir(joinpath(pathof(x),"log")) ? log(x,"initial setup",init=true) : nothing + put!(x.channel,Oceananigans_launch) + + x.outputs["grid"]=grid + x.outputs["IC"]=IC + x.outputs["BC"]=BC + + if haskey(x.inputs,"checkpoint") + checkpoint_file=joinpath(x,basename(x.inputs["checkpoint"])) + if occursin("http",x.inputs["checkpoint"]) + Downloads.download(x.inputs["checkpoint"],checkpoint_file) + else + cp(x.inputs["checkpoint"],checkpoint_file) + end + end + + println("Oceananigans run directory is \n "*pathof(x)) + + return true +end + +## + +end + diff --git a/src/Speedy.jl b/src/Speedy.jl new file mode 100644 index 0000000..0ff544d --- /dev/null +++ b/src/Speedy.jl @@ -0,0 +1,152 @@ + +module Speedy + +import ClimateModels +import ClimateModels: build, setup, AbstractModelConfig, read_NetCDF + +using Suppressor +git=ClimateModels.git +DataFrame=ClimateModels.DataFrame +uuid4=ClimateModels.uuid4 +UUID=ClimateModels.UUID +OrderedDict=ClimateModels.OrderedDict + +""" + struct SpeedyConfig <: AbstractModelConfig + +Concrete type of `AbstractModelConfig` for `SPEEDY` model. +""" +Base.@kwdef struct SpeedyConfig <: AbstractModelConfig + model :: String = "speedy" + configuration :: String = "default" + inputs :: OrderedDict{Any,Any} = OrderedDict{Any,Any}() + outputs :: OrderedDict{Any,Any} = OrderedDict{Any,Any}() + status :: OrderedDict{Any,Any} = OrderedDict{Any,Any}() + channel :: Channel{Any} = Channel{Any}(10) + folder :: String = tempdir() + ID :: UUID = uuid4() +end + +function build(x :: SpeedyConfig) + pth0=pwd() + pth=pathof(x) + + cd(pth) + if Sys.isapple() +# ENV["NETCDF"] = "/usr/local/" + ENV["NETCDF"]="/opt/homebrew/" + else + ENV["NETCDF"] = "/usr/" + end + @suppress run(`bash build.sh`) + cd(pth0) +end + +function SPEEDY_launch(x::SpeedyConfig) + pth0=pwd() + pth=pathof(x) + cd(pth) + @suppress run(`bash run.sh`) + cd(pth0) +end + +function setup(x :: SpeedyConfig) + !isdir(joinpath(x.folder)) ? mkdir(joinpath(x.folder)) : nothing + pth=pathof(x) + !isdir(pth) ? mkdir(pth) : nothing + + url="https://github.com/gaelforget/speedy.f90" + @suppress run(`$(git()) clone -b more_diags $url $pth`) + + !isdir(joinpath(pth,"log")) ? log(x,"initial setup",init=true) : nothing + + put!(x.channel,SPEEDY_launch) +end + +function write_ini_file(MC,nmonths) + (ny,nm)=divrem(nmonths+1,12) + nm==0 ? (ny,nm)=(ny-1,12) : nothing + (ny,nm) + + nml_ini="! nsteps_out = model variables are output every nsteps_out timesteps +! nstdia = model diagnostics are printed every nstdia timesteps +¶ms +nsteps_out = 180 +nstdia = 180 +/ +! start_datetime = the start datetime of the model integration +! end_datetime = the end datetime of the model integration +&date +start_datetime%year = 1982 +start_datetime%month = 1 +start_datetime%day = 1 +start_datetime%hour = 0 +start_datetime%minute = 0 +end_datetime%year = $(1982+ny) +end_datetime%month = $(nm) +end_datetime%day = 1 +end_datetime%hour = 0 +end_datetime%minute = 0 +/ +" + + write(joinpath(pathof(MC),"namelist.nml"),nml_ini) + + "Done with parameter file" +end + +function list_files_output(x::SpeedyConfig) + pth=pathof(x) + tmp=readdir(joinpath(pth,"rundir")) + files=tmp[findall(occursin.(".nc",tmp))[2:end]] + nt=length(files) + [joinpath(pth,"rundir",t) for t in files] +end + +function read_output(files,varname="t",time=1) + nt=length(files) + isnothing(time) ? t=1 : t=mod(time,Base.OneTo(nt)) + + lon = read_NetCDF(files[t])["lon"][:] + lat = read_NetCDF(files[t])["lat"][:] + lev = read_NetCDF(files[t])["lev"][:] + tmp = read_NetCDF(files[t])[varname] + + return (lon=lon,lat=lat,lev,values=tmp,fil=files[t]) +end + +function get_msk(rundir) + file=joinpath(rundir,"surface.nc") + lsm=read_NetCDF(file)["lsm"][:,:] + msk=Float64.(reverse(lsm,dims=2)) + msk[findall(msk[:,:].==1.0)].=NaN + msk[findall(msk[:,:].<1.0)].=1.0 + msk +end + +function read_namelist(x:: SpeedyConfig) + + fil=joinpath(pathof(x),"rundir/namelist.nml") + nml=readlines(fil) + i=findall([nml[i][1]!=='!' for i in 1:length(nml)]) + nml=nml[i] + + i0=findall([nml[i][1]=='&' for i in 1:length(nml)]) + i1=findall([nml[i][1]=='/' for i in 1:length(nml)]) + + tmp0=OrderedDict() + for i in 1:length(i0) + tmp1=Symbol(nml[i0[i]][2:end]) + tmp2=OrderedDict() + for j in i0[i]+1:i1[i]-1 + tmp3=split(nml[j],'=') + tmp2[tmp3[1]]=parse(Int,tmp3[2]) + end + tmp0[tmp1]=tmp2 + end + + return tmp0 +end + +end + diff --git a/test/Project.toml b/test/Project.toml index e66a5c7..9f5dacf 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -6,9 +6,14 @@ IniFile = "83e8ac13-25f8-5344-8a64-a9f2b223428f" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9" +Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Zarr = "0a941bbe-ad1d-11e8-39d9-ab76183a1d99" + +[compat] +Oceananigans = "0.90" + diff --git a/test/runtests.jl b/test/runtests.jl index 912736b..a110529 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,9 +1,55 @@ using ClimateModels, Documenter, Test, PyCall, Conda, CairoMakie -import Zarr, NetCDF, IniFile +import Zarr, NetCDF, IniFile, Oceananigans + +@testset "Oceananigans" begin + Nhours=1 + checkpoint_url="https://zenodo.org/record/8322234/files/model_checkpoint_iteration42423.jld2" + inputs=Dict("Nh" => 144+Nhours, "checkpoint" => checkpoint_url) + MC=OceananigansConfig(configuration="daily_cycle",inputs=inputs) + run(MC) + + nt=ClimateModels.Oceananigans.nt_from_jld2(MC) + XZ=ClimateModels.Oceananigans.xz_plot_prep(MC,1) + xz_fig=ClimateModels.plot_examples(:Oceananigans_xz,XZ...) + + T,S,w,νₑ=ClimateModels.Oceananigans.tz_slice(MC,nt=nt) + xw, yw, zw, xT, yT, zT=ClimateModels.Oceananigans.read_grid(MC) + tz_fig=ClimateModels.plot_examples(:Oceananigans_tz,xw, yw, zw, xT, yT, zT,T,S,w,νₑ) + @test isa(tz_fig,Figure) +end ClimateModels.conda(:fair) ClimateModels.pyimport(:fair) +@testset "FaIR" begin + MC=FaIRConfig() + run(MC) + scenarios,temperatures=FaIR.loop_over_scenarios() + f=ClimateModels.plot_examples(:FaIR,scenarios,temperatures) + @test isa(f,Figure) +end + +@testset "Speedy" begin + MC=SpeedyConfig() + run(MC) + + nml=Speedy.read_namelist(MC) + files=Speedy.list_files_output(MC) + + myvar="t"; ti=1 + tmp=Speedy.read_output(files,myvar,ti) + f_xy=ClimateModels.plot_examples(:Speedy_xy,tmp,myvar,ti,8) + f_zm=ClimateModels.plot_examples(:Speedy_zm,tmp,myvar,ti) + + rundir=joinpath(MC,"rundir") + msk=Speedy.get_msk(rundir) + myvar="sst"; to=1 + ncfile = NetCDF.open(joinpath(rundir,"sea_surface_temperature.nc")) + f=ClimateModels.plot_examples(:Speedy_input,ncfile,myvar,to,msk) + + @test isa(f,Figure) +end + @testset "JuliaClimate/Notebooks" begin nbs=notebooks.list() path=joinpath(tempdir(),"nbs")