From 3d76db97fd8d872a24e614b12b01538bfefdf475 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Wed, 25 Dec 2024 06:18:59 +1100 Subject: [PATCH 1/4] Add JRA55 in docs (#319) * disable ecco2 example * fix docstring rendering * add JRA55 module in docs library * fix rendering --- docs/make.jl | 4 +- docs/src/library/internals.md | 7 ++++ docs/src/library/public.md | 7 ++++ src/DataWrangling/JRA55.jl | 79 +++++++++++++++++------------------ 4 files changed, 55 insertions(+), 42 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 4f9164d5..317e81e6 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -13,7 +13,7 @@ const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples") const OUTPUT_DIR = joinpath(@__DIR__, "src/literated") to_be_literated = [ - "ecco_inspect_temperature_salinity.jl", + # "ecco_inspect_temperature_salinity.jl", "generate_bathymetry.jl", "generate_surface_fluxes.jl", "single_column_os_papa_simulation.jl", @@ -41,7 +41,7 @@ pages = [ "Home" => "index.md", "Examples" => [ - "Inspect ECCO2 data" => "literated/ecco_inspect_temperature_salinity.md", + # "Inspect ECCO2 data" => "literated/ecco_inspect_temperature_salinity.md", "Generate bathymetry" => "literated/generate_bathymetry.md", "Surface fluxes" => "literated/generate_surface_fluxes.md", "Single-column simulation" => "literated/single_column_os_papa_simulation.md", diff --git a/docs/src/library/internals.md b/docs/src/library/internals.md index 629b0216..503acc88 100644 --- a/docs/src/library/internals.md +++ b/docs/src/library/internals.md @@ -38,6 +38,13 @@ Modules = [ClimaOcean.ECCO] Public = false ``` +## JRA55 + +```@autodocs +Modules = [ClimaOcean.JRA55] +Public = false +``` + ## Bathymetry ```@autodocs diff --git a/docs/src/library/public.md b/docs/src/library/public.md index 13f94c53..31a77c2b 100644 --- a/docs/src/library/public.md +++ b/docs/src/library/public.md @@ -39,6 +39,13 @@ Modules = [ClimaOcean.ECCO] Private = false ``` +## JRA55 + +```@autodocs +Modules = [ClimaOcean.JRA55] +Private = false +``` + ## Bathymetry ```@autodocs diff --git a/src/DataWrangling/JRA55.jl b/src/DataWrangling/JRA55.jl index 96ee8f9e..e1774ad2 100644 --- a/src/DataWrangling/JRA55.jl +++ b/src/DataWrangling/JRA55.jl @@ -23,7 +23,7 @@ using ClimaOcean.OceanSeaIceModels: using CUDA: @allowscalar using NCDatasets -using JLD2 +using JLD2 using Dates using Scratch @@ -286,40 +286,40 @@ new_backend(::JRA55NetCDFBackend, start, length) = JRA55NetCDFBackend(start, len """ JRA55_field_time_series(variable_name; architecture = CPU(), - time_indices = nothing, - latitude = nothing, - longitude = nothing, + grid = nothing, location = nothing, url = nothing, + dir = download_jra55_cache, filename = nothing, shortname = nothing, + latitude = nothing, + longitude = nothing, backend = InMemory(), + time_indexing = Cyclical(), preprocess_chunk_size = 10, - preprocess_architecture = CPU()) + preprocess_architecture = CPU(), + time_indices = nothing) Return a `FieldTimeSeries` containing atmospheric reanalysis data for `variable_name`, which describes one of the variables in the "repeat year forcing" dataset derived -from the Japanese 55-year atmospheric reanalysis for driving ocean-sea-ice models (JRA55-do). -For more information about the derivation of the repeat year forcing dataset, see - -"Stewart et al., JRA55-do-based repeat year forcing datasets for driving ocean–sea-ice models", -Ocean Modelling, 2020, https://doi.org/10.1016/j.ocemod.2019.101557. - -The `variable_name`s (and their `shortname`s used in NetCDF files) -available from the JRA55-do are: - - - `:river_freshwater_flux` ("friver") - - `:rain_freshwater_flux` ("prra") - - `:snow_freshwater_flux` ("prsn") - - `:iceberg_freshwater_flux` ("licalvf") - - `:specific_humidity` ("huss") - - `:sea_level_pressure` ("psl") - - `:relative_humidity` ("rhuss") - - `:downwelling_longwave_radiation` ("rlds") - - `:downwelling_shortwave_radiation` ("rsds") - - `:temperature` ("ras") - - `:eastward_velocity` ("uas") - - `:northward_velocity` ("vas") +from the Japanese 55-year atmospheric reanalysis for driving ocean-sea ice models (JRA55-do). +For more information about the derivation of the repeat-year forcing dataset, see + +> Stewart et al. (2020). JRA55-do-based repeat year forcing datasets for driving ocean–sea-ice models, _Ocean Modelling_, **147**, 101557, https://doi.org/10.1016/j.ocemod.2019.101557. + +The `variable_name`s (and their `shortname`s used in NetCDF files) available from the JRA55-do are: +- `:river_freshwater_flux` ("friver") +- `:rain_freshwater_flux` ("prra") +- `:snow_freshwater_flux` ("prsn") +- `:iceberg_freshwater_flux` ("licalvf") +- `:specific_humidity` ("huss") +- `:sea_level_pressure` ("psl") +- `:relative_humidity` ("rhuss") +- `:downwelling_longwave_radiation` ("rlds") +- `:downwelling_shortwave_radiation` ("rsds") +- `:temperature` ("ras") +- `:eastward_velocity` ("uas") +- `:northward_velocity` ("vas") Keyword arguments ================= @@ -352,8 +352,8 @@ Keyword arguments - `interpolated_file`: file holding an Oceananigans compatible version of the JRA55 data. If it does not exist it will be generated. -- `time_chunks_in_memory`: number of fields held in memory. If `nothing` the whole timeseries is - loaded (not recommended). +- `time_chunks_in_memory`: number of fields held in memory. If `nothing` then the whole timeseries + is loaded (not recommended). """ function JRA55_field_time_series(variable_name; architecture = CPU(), @@ -639,20 +639,20 @@ JRA55PrescribedAtmosphere(arch::Distributed, time_indices=Colon(); kw...) = # TODO: allow the user to pass dates """ JRA55PrescribedAtmosphere(architecture::AA, time_indices=Colon(); - backend = nothing, - time_indexing = Cyclical(), - reference_height = 10, # meters - include_rivers_and_icebergs = false, - other_kw...) + backend = nothing, + time_indexing = Cyclical(), + reference_height = 10, # meters + include_rivers_and_icebergs = false, + other_kw...) Return a `PrescribedAtmosphere` representing JRA55 reanalysis data. """ function JRA55PrescribedAtmosphere(architecture::AA, time_indices=Colon(); - backend = nothing, - time_indexing = Cyclical(), - reference_height = 10, # meters - include_rivers_and_icebergs = false, - other_kw...) + backend = nothing, + time_indexing = Cyclical(), + reference_height = 10, # meters + include_rivers_and_icebergs = false, + other_kw...) if isnothing(backend) # apply a default Ni = try @@ -667,7 +667,7 @@ function JRA55PrescribedAtmosphere(architecture::AA, time_indices=Colon(); end kw = (; time_indices, time_indexing, backend, architecture) - kw = merge(kw, other_kw) + kw = merge(kw, other_kw) ua = JRA55_field_time_series(:eastward_velocity; kw...) va = JRA55_field_time_series(:northward_velocity; kw...) @@ -719,4 +719,3 @@ function JRA55PrescribedAtmosphere(architecture::AA, time_indices=Colon(); end end # module - From 0baa94f0d5fac23de11c08f99bb65f8b871a2210 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 26 Dec 2024 05:45:56 +1100 Subject: [PATCH 2/4] Fix axis labels in near-global example + Docstring rendering + Few typos (#315) * fix docstring rendering + few typos * add coords in plots + phrasing tweaks * drop inspect ecco * interpolate s -> (Center, Center) * markdown formatting requires 2 spaces for list within a list * there is no default in remove_minor_basins! * fix math rendering * Update near_global_ocean_simulation.jl * fix bug + add title with time * Update make.jl * Update make.jl * Update src/Bathymetry.jl Co-authored-by: Simone Silvestri * homogenize font label size * homogenize font label size --------- Co-authored-by: Simone Silvestri --- docs/make.jl | 4 +- examples/near_global_ocean_simulation.jl | 66 +++++++++------ src/Bathymetry.jl | 83 ++++++++++--------- .../CrossRealmFluxes/surface_temperature.jl | 35 ++++---- 4 files changed, 104 insertions(+), 84 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 317e81e6..0964c0f6 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -13,7 +13,6 @@ const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples") const OUTPUT_DIR = joinpath(@__DIR__, "src/literated") to_be_literated = [ - # "ecco_inspect_temperature_salinity.jl", "generate_bathymetry.jl", "generate_surface_fluxes.jl", "single_column_os_papa_simulation.jl", @@ -41,7 +40,6 @@ pages = [ "Home" => "index.md", "Examples" => [ - # "Inspect ECCO2 data" => "literated/ecco_inspect_temperature_salinity.md", "Generate bathymetry" => "literated/generate_bathymetry.md", "Surface fluxes" => "literated/generate_surface_fluxes.md", "Single-column simulation" => "literated/single_column_os_papa_simulation.md", @@ -59,7 +57,7 @@ pages = [ makedocs(sitename = "ClimaOcean.jl"; format, - pages, + pages, modules = [ClimaOcean], doctest = true, clean = true, diff --git a/examples/near_global_ocean_simulation.jl b/examples/near_global_ocean_simulation.jl index 6db964bc..ae7ce604 100644 --- a/examples/near_global_ocean_simulation.jl +++ b/examples/near_global_ocean_simulation.jl @@ -25,11 +25,11 @@ using Printf # # We define a global grid with a horizontal resolution of 1/4 degree and 40 vertical levels. # The grid is a `LatitudeLongitudeGrid` spanning latitudes from 75°S to 75°N. -# We use an exponential vertical spacing to better resolve the upper ocean layers. +# We use an exponential vertical spacing to better resolve the upper-ocean layers. # The total depth of the domain is set to 6000 meters. # Finally, we specify the architecture for the simulation, which in this case is a GPU. -arch = GPU() +arch = GPU() Nx = 1440 Ny = 600 @@ -49,10 +49,10 @@ grid = LatitudeLongitudeGrid(arch; # # We use `regrid_bathymetry` to derive the bottom height from ETOPO1 data. # To smooth the interpolated data we use 5 interpolation passes. We also fill in -# all the minor enclosed basins except the 3 largest `major_basins`, as well as regions -# that are shallower than `minimum_depth`. +# (i) all the minor enclosed basins except the 3 largest `major_basins`, as well as +# (ii) regions that are shallower than `minimum_depth`. -bottom_height = regrid_bathymetry(grid; +bottom_height = regrid_bathymetry(grid; minimum_depth = 10meters, interpolation_passes = 5, major_basins = 3) @@ -63,7 +63,9 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_ h = grid.immersed_boundary.bottom_height -fig, ax, hm = heatmap(h, colormap=:deep, colorrange=(-depth, 0)) +fig = Figure(size = (800, 450)) +ax = Axis(fig[1, 1]) +hm = heatmap!(ax, h, colormap=:deep, colorrange=(-depth, 0)) cb = Colorbar(fig[0, 1], hm, label="Bottom height (m)", vertical=false) hidedecorations!(ax) save("bathymetry.png", fig) @@ -81,7 +83,7 @@ ocean = ocean_simulation(grid) ocean.model -# We initialize the ocean model to ECCO2 temperature and salinity for January 1, 1993. +# We initialize the ocean model with ECCO2 temperature and salinity for January 1, 1993. date = DateTimeProlepticGregorian(1993, 1, 1) set!(ocean.model, T=ECCOMetadata(:temperature; dates=date), @@ -90,8 +92,7 @@ set!(ocean.model, T=ECCOMetadata(:temperature; dates=date), # ### Prescribed atmosphere and radiation # # Next we build a prescribed atmosphere state and radiation model, -# which will drive the development of the ocean simulation. -# We use the default `Radiation` model, +# which will drive the ocean simulation. We use the default `Radiation` model, ## The radiation model specifies an ocean albedo emissivity to compute the net radiative ## fluxes. The default ocean albedo is based on Payne (1982) and depends on cloud cover @@ -101,9 +102,8 @@ set!(ocean.model, T=ECCOMetadata(:temperature; dates=date), radiation = Radiation(arch) # The atmospheric data is prescribed using the JRA55 dataset. -# The JRA55 dataset provides atmospheric -# data such as temperature, humidity, and wind fields to calculate turbulent fluxes -# using bulk formulae, see [`CrossRealmFluxes`](@ref). +# The JRA55 dataset provides atmospheric data such as temperature, humidity, and winds +# to calculate turbulent fluxes using bulk formulae, see [`CrossRealmFluxes`](@ref). # The number of snapshots that are loaded into memory is determined by # the `backend`. Here, we load 41 snapshots at a time into memory. @@ -116,8 +116,8 @@ atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(41)) coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) -# We then create a coupled simulation, starting with a time step of 90 seconds -# and running the simulation for 10 days. +# We then create a coupled simulation. We start with a small-ish time step of 90 seconds. +# We run the simulation for 10 days with this small-ish time step. simulation = Simulation(coupled_model; Δt=90, stop_time=10days) @@ -125,7 +125,7 @@ simulation = Simulation(coupled_model; Δt=90, stop_time=10days) wall_time = Ref(time_ns()) -function progress(sim) +function progress(sim) ocean = sim.model.ocean u, v, w = ocean.model.velocities T = ocean.model.tracers.T @@ -154,7 +154,8 @@ simulation.callbacks[:progress] = Callback(progress, TimeInterval(5days)) # # We define output writers to save the simulation data at regular intervals. # In this case, we save the surface fluxes and surface fields at a relatively high frequency (every day). -# The `indices` keyword argument allows us to save down a slice at the surface, which is located at `k = grid.Nz` +# The `indices` keyword argument allows us to save only a slice of the three dimensional variable. +# Below, we use `indices` to save only the values of the variables at the surface, which corresponds to `k = grid.Nz` outputs = merge(ocean.model.tracers, ocean.model.velocities) ocean.output_writers[:surface] = JLD2OutputWriter(ocean.model, outputs; @@ -167,21 +168,24 @@ ocean.output_writers[:surface] = JLD2OutputWriter(ocean.model, outputs; # ### Spinning up the simulation # -# We spin up the simulation with a very short time-step to resolve the "initialization shock" -# associated with starting from ECCO initial conditions that are both interpolated and also +# We spin up the simulation with a small-ish time-step to resolve the "initialization shock" +# associated with starting from ECCO2 initial conditions that are both interpolated and also # satisfy a different dynamical balance than our simulation. run!(simulation) # ### Running the simulation for real +# After the initial spin up of 10 days, we can increase the time-step and run for longer. + simulation.stop_time = 60days simulation.Δt = 10minutes run!(simulation) # ## A pretty movie # -# It's time to make a pretty movie of the simulation. First we plot a snapshot: +# It's time to make a pretty movie of the simulation. First we load the output we've been saving on +# disk and plot the final snapshot: u = FieldTimeSeries("near_global_surface_fields.jld2", "u"; backend = OnDisk()) v = FieldTimeSeries("near_global_surface_fields.jld2", "v"; backend = OnDisk()) @@ -209,7 +213,9 @@ end un = Field{Face, Center, Nothing}(u.grid) vn = Field{Center, Face, Nothing}(v.grid) -s = Field(sqrt(un^2 + vn^2)) + +s = @at (Center, Center, Nothing) sqrt(un^2 + vn^2) # compute √(u²+v²) and interpolate back to Center, Center +s = Field(s) sn = @lift begin parent(un) .= parent(u[$n]) @@ -220,26 +226,34 @@ sn = @lift begin view(sn, :, :, 1) end -fig = Figure(size = (800, 1200)) +title = @lift string("Near-global 1/4 degree ocean simulation after ", + prettytime(times[$n] - times[1])) + +λ, φ, _ = nodes(T) # T, e, and s all live on the same grid locations + +fig = Figure(size = (1000, 1500)) axs = Axis(fig[1, 1], xlabel="Longitude (deg)", ylabel="Latitude (deg)") axT = Axis(fig[2, 1], xlabel="Longitude (deg)", ylabel="Latitude (deg)") axe = Axis(fig[3, 1], xlabel="Longitude (deg)", ylabel="Latitude (deg)") -hm = heatmap!(axs, sn, colorrange = (0, 0.5), colormap = :deep, nan_color=:lightgray) -Colorbar(fig[1, 2], hm, label = "Surface speed (m s⁻¹)") +hm = heatmap!(axs, λ, φ, sn, colorrange = (0, 0.5), colormap = :deep, nan_color=:lightgray) +Colorbar(fig[1, 2], hm, label = "Surface Speed (m s⁻¹)") -hm = heatmap!(axT, Tn, colorrange = (-1, 30), colormap = :magma, nan_color=:lightgray) +hm = heatmap!(axT, λ, φ, Tn, colorrange = (-1, 30), colormap = :magma, nan_color=:lightgray) Colorbar(fig[2, 2], hm, label = "Surface Temperature (ᵒC)") -hm = heatmap!(axe, en, colorrange = (0, 1e-3), colormap = :solar, nan_color=:lightgray) +hm = heatmap!(axe, λ, φ, en, colorrange = (0, 1e-3), colormap = :solar, nan_color=:lightgray) Colorbar(fig[3, 2], hm, label = "Turbulent Kinetic Energy (m² s⁻²)") + +Label(fig[0, :], title) + save("snapshot.png", fig) nothing #hide # ![](snapshot.png) -# And now a movie: +# And now we make a movie: record(fig, "near_global_ocean_surface.mp4", 1:Nt, framerate = 8) do nn n[] = nn diff --git a/src/Bathymetry.jl b/src/Bathymetry.jl index bf9ee264..b9148d06 100644 --- a/src/Bathymetry.jl +++ b/src/Bathymetry.jl @@ -32,65 +32,70 @@ end """ regrid_bathymetry(target_grid; - url = "https://www.ngdc.noaa.gov/thredds/fileServer/global/ETOPO2022/60s/60s_surface_elev_netcdf", - height_above_water = , + height_above_water = nothing, minimum_depth = 0, - dir = download_cache, - filename = "ETOPO_2022_v1_60s_N90W180_surface.nc") + dir = download_bathymetry_cache, + url = "https://www.ngdc.noaa.gov/thredds/fileServer/global/ETOPO2022/60s/60s_surface_elev_netcdf", + filename = "ETOPO_2022_v1_60s_N90W180_surface.nc", + interpolation_passes = 1, + major_basins = 1) -Regrid bathymetry associated with the NetCDF file at `path = joinpath(dir, filename)` to `target_grid`. +Return bathymetry associated with the NetCDF file at `path = joinpath(dir, filename)` regridded onto `target_grid`. If `path` does not exist, then a download is attempted from `joinpath(url, filename)`. Arguments ========= -- `target_grid`: grid to interpolate onto +- `target_grid`: grid to interpolate the bathymetry onto. Keyword Arguments ================= -- `height_above_water`: limits the maximum height of above-water topography (where h > 0). If - `nothing` the original topography is retained +- `height_above_water`: limits the maximum height of above-water topography (where ``h > 0``) before inetrpolating. + Default: `nothing`, which implies that the original topography is retained. - `minimum_depth`: minimum depth for the shallow regions, defined as a positive value. - `h > - minimum_depth` will be considered land + `h > - minimum_depth` is considered land. Default: 0. -- `dir`: directory of the bathymetry-containing file +- `dir`: directory of the bathymetry-containing file. Default: `download_bathymetry_cache`. -- `filename`: file containing bathymetric data. Must be netcdf with fields: - (1) `lat` vector of latitude nodes - (2) `lon` vector of longitude nodes - (3) `z` matrix of depth values +- `filename`: file containing bathymetric data. Must be netCDF with fields: + 1. `lat` vector of latitude nodes + 2. `lon` vector of longitude nodes + 3. `z` matrix of depth values - `interpolation_passes`: regridding/interpolation passes. The bathymetry is interpolated in - `interpolation_passes - 1` intermediate steps. With more steps the - final bathymetry will be smoother. - - Example: interpolating from a 400x200 grid to a 100x100 grid in 4 passes involves: + `interpolation_passes - 1` intermediate steps. The more the interpolation + steps the smoother the final bathymetry becomes. - * 400x200 → 325x175 - * 325x175 → 250x150 - * 250x150 → 175x125 - * 175x125 → 100x100 + Example + ======= + + Interpolating from a 400 x 200 grid to a 100 x 100 grid in 4 passes involves: + + * 400 x 200 → 325 x 175 + * 325 x 175 → 250 x 150 + * 250 x 150 → 175 x 125 + * 175 x 125 → 100 x 100 If _coarsening_ the original grid, linear interpolation in passes is equivalent to applying a smoothing filter, with more passes increasing the strength of the filter. - If _refining_ the original grid, additional passes will not help and no intermediate - steps will be performed. + If _refining_ the original grid, additional passes do not help and no intermediate + steps are performed. - `major_basins`: Number of "independent major basins", or fluid regions fully encompassed by land, that are retained by [`remove_minor_basins!`](@ref). Basins are removed by order of size: - the smallest basins are removed first. `major_basins=1` will retain only the largest basin. - Default: `Inf`, which does not remove any basins. + the smallest basins are removed first. `major_basins = 1` retains only the largest basin. + If `Inf` then no basins are removed. Default: 1. """ function regrid_bathymetry(target_grid; height_above_water = nothing, minimum_depth = 0, dir = download_bathymetry_cache, - url = "https://www.ngdc.noaa.gov/thredds/fileServer/global/ETOPO2022/60s/60s_surface_elev_netcdf", + url = "https://www.ngdc.noaa.gov/thredds/fileServer/global/ETOPO2022/60s/60s_surface_elev_netcdf", filename = "ETOPO_2022_v1_60s_N90W180_surface.nc", interpolation_passes = 1, - major_basins = 1) # Allow an `Inf` number of ``lakes'' + major_basins = 1) # Allow an `Inf` number of "lakes" filepath = joinpath(dir, filename) fileurl = url * "/" * filename # joinpath on windows creates the wrong url @@ -191,7 +196,7 @@ function regrid_bathymetry(target_grid; end # Here we can either use `regrid!` (three dimensional version) or `interpolate` -function interpolate_bathymetry_in_passes(native_z, target_grid; +function interpolate_bathymetry_in_passes(native_z, target_grid; passes = 10) Nλt, Nφt = Nt = size(target_grid) Nλn, Nφn = Nn = size(native_z) @@ -201,9 +206,9 @@ function interpolate_bathymetry_in_passes(native_z, target_grid; if Nλt > Nλn || Nφt > Nφn target_z = Field{Center, Center, Nothing}(target_grid) interpolate!(target_z, native_z) - @info string("Skipping passes for interplating bathymetry of size $Nn ", '\n', + @info string("Skipping passes for interpolating bathymetry of size $Nn ", '\n', "to target grid of size $Nt. Interpolation passes may only ", '\n', - "be used to refine bathymetryand requires that the bathymetry ", '\n', + "be used to refine bathymetry and require that the bathymetry ", '\n', "is larger than the target grid in both horizontal directions.") return target_z end @@ -230,9 +235,9 @@ function interpolate_bathymetry_in_passes(native_z, target_grid; @debug "Bathymetry interpolation pass $pass with size $new_size" new_grid = LatitudeLongitudeGrid(architecture(target_grid), - size = new_size, - latitude = (latitude[1], latitude[2]), - longitude = (longitude[1], longitude[2]), + size = new_size, + latitude = (latitude[1], latitude[2]), + longitude = (longitude[1], longitude[2]), z = (0, 1), topology = (TX, TY, Bounded)) @@ -276,7 +281,7 @@ function regrid_bathymetry(target_grid::DistributedGrid; kw...) local_bottom_height = Field{Center, Center, Nothing}(target_grid) set!(local_bottom_height, bottom_height) fill_halo_regions!(local_bottom_height) - + return local_bottom_height end @@ -291,14 +296,14 @@ Arguments ========= - `z_data`: A 2D array representing the bathymetry data. -- `keep_major_basins`: The maximum number of connected regions to keep. - Default is `Inf`, which means all connected regions are kept. +- `keep_major_basins`: The maximum number of connected regions to keep. + If `Inf` is provided then all connected regions are kept. """ function remove_minor_basins!(zb::Field, keep_major_basins) zb_cpu = on_architecture(CPU(), zb) TX = topology(zb_cpu.grid, 1) - + Nx, Ny, _ = size(zb_cpu.grid) zb_data = maybe_extend_longitude(zb_cpu, TX()) # Outputs a 2D AbstractArray @@ -326,7 +331,7 @@ function maybe_extend_longitude(zb_cpu, ::Periodic) # Update offsets yoffsets = zb_cpu.data.offsets[2] xoffsets = - nx - + return OffsetArray(zb_parent, xoffsets, yoffsets) end diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index 3da5b500..46c6e5b7 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -12,11 +12,12 @@ import Thermodynamics as AtmosphericThermodynamics #### Bulk surface temperature (the easiest case) #### -""" - struct BulkTemperature end +""" + struct BulkTemperature -Represents the surface temperature used in fixed-point iteration for surface fluxes following similarity theory. -The surface temperature is not calculated but provided by either the ocean or the sea ice model. +A type to represent the surface temperature used in fixed-point iteration for surface +fluxes following similarity theory. The surface temperature is not calculated but instead +provided by either the ocean or the sea ice model. """ struct BulkTemperature end @@ -29,17 +30,19 @@ struct BulkTemperature end """ struct SkinTemperature - internal_flux :: I - end A type to represent the surface temperature used in the flux calculation. The surface temperature is calculated from the flux balance at the surface. In particular, the surface temperature ``θₛ`` is the root of: -F(θₛ) - Jᵀ = 0 (all fluxes positive upwards) +```math +F(θₛ) - Jᵀ = 0 +``` + +where ``Jᵀ`` are the fluxes at the top of the surface (turbulent + radiative), and +``F`` is the internal diffusive flux dependent on the surface temperature itself. -where Jᵀ are the fluxes at the top of the surface (turbulent + radiative), and F is the internal diffusive flux -dependent on the surface temperature itself. +Note that all fluxes positive upwards. """ struct SkinTemperature{I} internal_flux :: I @@ -77,19 +80,19 @@ DiffusiveFlux(FT; κ = 1e-2, δ = 1.0) = DiffusiveFlux(convert(FT, δ), convert( # that can be explored in the future. @inline flux_balance_temperature(F::DiffusiveFlux, θₒ, Jᵀ) = θₒ - Jᵀ / F.κ * F.δ -# he flaw here is that the ocean emissivity and albedo are fixed, but they might be a function of the +# the flaw here is that the ocean emissivity and albedo are fixed, but they might be a function of the # surface temperature, so we might need to pass the radiation and the albedo and emissivity as arguments. -@inline function compute_surface_temperature(st::SkinTemperature, θₛ, ℂ, 𝒬₀, - ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, g, - prescribed_heat_fluxes, +@inline function compute_surface_temperature(st::SkinTemperature, θₛ, ℂ, 𝒬₀, + ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, g, + prescribed_heat_fluxes, radiation_properties) Rd = prescribed_heat_fluxes # net downwelling radiation (positive out of the ocean) - + # upwelling radiation is calculated explicitly Ru = upwelling_radiation(θₛ, radiation_properties) Rn = Rd + Ru # Net radiation (positive out of the ocean) - + u★ = Σ★.momentum θ★ = Σ★.temperature q★ = Σ★.water_vapor @@ -104,4 +107,4 @@ DiffusiveFlux(FT; κ = 1e-2, δ = 1.0) = DiffusiveFlux(convert(FT, δ), convert( θₛ = flux_balance_temperature(st.internal_flux, θₒ, Jᵀ) return θₛ -end \ No newline at end of file +end From c6bffe74e38f0f4afb2dcc57eb08c975e68b6ffa Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 27 Dec 2024 14:53:11 +1100 Subject: [PATCH 3/4] hidedecorations was making the colormap labels go out of fig (#320) --- examples/near_global_ocean_simulation.jl | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/examples/near_global_ocean_simulation.jl b/examples/near_global_ocean_simulation.jl index ae7ce604..f869526a 100644 --- a/examples/near_global_ocean_simulation.jl +++ b/examples/near_global_ocean_simulation.jl @@ -63,11 +63,8 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_ h = grid.immersed_boundary.bottom_height -fig = Figure(size = (800, 450)) -ax = Axis(fig[1, 1]) -hm = heatmap!(ax, h, colormap=:deep, colorrange=(-depth, 0)) -cb = Colorbar(fig[0, 1], hm, label="Bottom height (m)", vertical=false) -hidedecorations!(ax) +fig, ax, hm = heatmap(h, colormap=:deep, colorrange=(-depth, 0)) +Colorbar(fig[0, 1], hm, label="Bottom height (m)", vertical=false) save("bathymetry.png", fig) nothing #hide From 048c893c5cf60c0df64b049cd4efb7bef6b4f2f5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Sun, 5 Jan 2025 17:26:41 +0100 Subject: [PATCH 4/4] correct (#321) --- .../atmosphere_ocean_fluxes.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl index 1d8a0f65..afdb7b96 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl @@ -353,11 +353,11 @@ end @inbounds begin # +0: cooling, -0: heating - Qv[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.latent_heat) - Qc[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.sensible_heat) - Fv[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.water_vapor) - ρτx[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.x_momentum) - ρτy[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.y_momentum) + Qv[i, j, 1] = ifelse(inactive, zero(grid), turbulent_fluxes.latent_heat) + Qc[i, j, 1] = ifelse(inactive, zero(grid), turbulent_fluxes.sensible_heat) + Fv[i, j, 1] = ifelse(inactive, zero(grid), turbulent_fluxes.water_vapor) + ρτx[i, j, 1] = ifelse(inactive, zero(grid), turbulent_fluxes.x_momentum) + ρτy[i, j, 1] = ifelse(inactive, zero(grid), turbulent_fluxes.y_momentum) Ts[i, j, 1] = surface_temperature end end @@ -433,10 +433,10 @@ end inactive = inactive_node(i, j, kᴺ, grid, c, c, c) @inbounds begin - τx[i, j, 1] = ifelse(inactive, 0, atmos_ocean_τx) - τy[i, j, 1] = ifelse(inactive, 0, atmos_ocean_τy) - Jᵀ[i, j, 1] = ifelse(inactive, 0, atmos_ocean_Jᵀ) - Jˢ[i, j, 1] = ifelse(inactive, 0, atmos_ocean_Jˢ) + τx[i, j, 1] = ifelse(inactive, zero(grid), atmos_ocean_τx) + τy[i, j, 1] = ifelse(inactive, zero(grid), atmos_ocean_τy) + Jᵀ[i, j, 1] = ifelse(inactive, zero(grid), atmos_ocean_Jᵀ) + Jˢ[i, j, 1] = ifelse(inactive, zero(grid), atmos_ocean_Jˢ) end end