Clean up ECCO functionality (#169)

* Clean up ECCO stuff

* Update for ECCO

* Import propagate_inbounds

* ecco -> ECCO in make.jl

* Fix ecco -> ECCO

* Update src/ClimaOcean.jl

Co-authored-by: Simone Silvestri <silvestri.simone0@gmail.com>

* Update src/DataWrangling/ECCO/ECCO_metadata.jl

Co-authored-by: Simone Silvestri <silvestri.simone0@gmail.com>

* Add a README to ECCO

* Update README.md

* Update ECCO_metadata.jl

* Fix capitalization in tests

* Case sensitivity

* Capitalization again

---------

Co-authored-by: Simone Silvestri <silvestri.simone0@gmail.com>

Author: glwagner

docs/make.jl

@@ -15,11 +15,11 @@ const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
 const OUTPUT_DIR   = joinpath(@__DIR__, "src/literated")
 
 to_be_literated = [
-    # "inspect_ecco_data.jl",
+    # "inspect_ECCO_data.jl",
     "generate_bathymetry.jl",
     "generate_surface_fluxes.jl",
     # "single_column_simulation.jl",
-    # "mediterranean_simulation_with_ecco_restoring.jl",
+    # "mediterranean_simulation_with_ECCO_restoring.jl",
     "near_global_ocean_simulation.jl"
 ]
 
@@ -52,11 +52,11 @@ pages = [
         ],
 
     "Examples" => [
-        # "Inspect ECCO2 data" => "literated/inspect_ecco_data.md",
+        # "Inspect ECCO2 data" => "literated/inspect_ECCO_data.md",
         "Generate bathymetry" => "literated/generate_bathymetry.md",
         "Surface fluxes" => "literated/generate_surface_fluxes.md",
         # "Single column simulation" => "literated/single_column_simulation.md",
-        # "Mediterranean simulation with ECCO restoring" => "literated/mediterranean_simulation_with_ecco_restoring.md",
+        # "Mediterranean simulation with ECCO restoring" => "literated/mediterranean_simulation_with_ECCO_restoring.md",
         "Near-global Ocean simulation" => "literated/near_global_ocean_simulation.md",
         ]
 ]

examples/generate_surface_fluxes.jl

@@ -16,25 +16,21 @@ using ClimaOcean.OceanSimulations
 using Oceananigans
 using CairoMakie
 
-# We start by defining a grid. The ECCO2 grid is a good starting point.
-# The ECCO2 grid is not "immersed" by default, but we can use the ECCO mask
-# to define the "bathymetry" of the ECCO fields.
-# `ecco2_center_mask` produces a field with `true` values where data is missing (i.e., in immersed cells).
-# We can use this mask as an immersed boundary for our grid.
-# Let's create the grid and visualize the mask.
-
-mask = ecco_mask()
+# # Computing fluxes on the ECCO2 grid
+#
+# We start by building the ECCO2 grid, using `ECCO_mask`
+# to identify missing cells.
+
+mask = ECCO_mask()
 grid = mask.grid
 grid = ImmersedBoundaryGrid(grid, GridFittedBoundary(mask))
 
 fig = Figure()
 ax  = Axis(fig[1, 1])
 heatmap!(ax, interior(grid.immersed_boundary.mask, :, :, grid.Nz))
+save("ECCO_continents.png", fig) # hide
 
-save("ecco_continents.png", fig)
-nothing #hide
-
-# ![](ecco_continents.png)
+# ![](ECCO_continents.png)
 
 # Next, we construct our atmosphere and ocean.
 #
@@ -50,103 +46,50 @@ nothing #hide
 #
 # We invoke the constructor with only the first two time indices, corresponding to 
 # January 1st (at 00:00 AM and 03:00 AM).
-# By passing the ECCO grid, we automatically interpolate the atmospheric data onto the grid.
-# Note that this is recommended only for small simulations (in terms of grid size).
-# By omitting the grid, the interpolation will be done on the fly.
-#
-# We construct the ocean simulation without considering advection, closures, or Coriolis effects since
-# we will not time-step the ocean but only use it to construct the fluxes.
-
-atmosphere  = JRA55_prescribed_atmosphere(1:2; backend = InMemory(), grid = grid.underlying_grid)
-
-ocean = ocean_simulation(grid; momentum_advection = nothing,
-                                 tracer_advection = nothing,
-                                          closure = nothing,
-                                         coriolis = nothing)
-
-# Now that we have an atmosphere and a container for the ocean, we need to populate
-# our ocean with initial conditions. To do this, we can use the ECCO2 dataset by
-# `set!`ting the model with the `ECCOMetadata`. If no date is specified,
-# the fields corresponding to January 1st, 1992 (the first available date in
-# ECCO2 dataset) are used.
-# This command will download the fields to the local machine.
-
-set!(ocean.model;
-      T = ECCOMetadata(:temperature), 
-      S = ECCOMetadata(:salinity))
-
-# The final step is to construct a coupled model.
-# The coupled model requires an ocean, which we have just constructed and initialized,
-# an atmosphere, which we have downloaded from the JRA55 dataset, a sea ice model 
-# (in this case we do not account for sea ice by defining `sea_ice = nothing`), 
-# and a radiation model. The default radiation model assumes two spectral bands: 
-# a shortwave band modeling visible and UV light, and a longwave band that accounts for
-# near, mid and far infrared (mostly far infrared given the low emission temperature). 
-# By constructing the coupled model, the `update_state!` function, which calculates the fluxes,
-# will be triggered.
-
-radiation = Radiation()
-sea_ice   = nothing
-coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
+
+atmosphere  = JRA55_prescribed_atmosphere(1:2; backend = InMemory())
+ocean = ocean_simulation(grid)
+
+# Now that we have an atmosphere and ocean, we `set!` the ocean temperature and salinity
+# to the ECCO2 data by first creating T, S metadata objects,
+
+T_metadata = ECCOMetadata(:temperature)
+S_metadata = ECCOMetadata(:salinity)
+
+# Note that if a date is not provided to `ECCOMetadata`, then the default Jan 1st, 1992 is used.
+# To copy the ECCO state into `ocean.model`, we use `set!`,
+
+set!(ocean.model; T=T_metadata, S=S_metadata)
+
+# Finally, we construct a coupled model, which will compute fluxes during construction.
+# We omit `sea_ice` so the model is ocean-only, and use the default `Radiation()` that
+# uses the two-band shortwave (visible and UV) + longwave (mid and far infrared)
+# decomposition of the radiation spectrum.
+
+coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation=Radiation())
 
 # Now that the surface fluxes are computed, we can extract and visualize them.
 # The turbulent fluxes are stored in `coupled_model.fluxes.turbulent`.
-# 
-# Qs = coupled_model.fluxes.turbulent.fields.sensible_heat : the sensible heat flux (in Wm⁻²)
-# Ql = coupled_model.fluxes.turbulent.fields.latent_heat   : the latent heat flux  (in Wm⁻²)
-# τx = coupled_model.fluxes.turbulent.fields.x_momentum    : the zonal wind stress (in Nm)
-# τy = coupled_model.fluxes.turbulent.fields.y_momentum    : the meridional wind stress (in Nm)
-# Mv = coupled_model.fluxes.turbulent.fields.water_vapor   : evaporation (in kg m⁻²s⁻¹)
-#
-# They are 2D fields (3D data structures with one point in the vertical). To extract the data, we use the 
-# `interior` functionality from Oceananigans.
-
-turbulent_fluxes = coupled_model.fluxes.turbulent.fields
 
-Qs = interior(turbulent_fluxes.sensible_heat, :, :, 1)
-Ql = interior(turbulent_fluxes.latent_heat,   :, :, 1)
-τx = interior(turbulent_fluxes.x_momentum,    :, :, 1)
-τy = interior(turbulent_fluxes.y_momentum,    :, :, 1)
-Mv = interior(turbulent_fluxes.water_vapor,   :, :, 1)
-nothing #hide
+fluxes = coupled_model.fluxes.turbulent.fields
 
 fig = Figure(size = (800, 400))
-ax = Axis(fig[1, 1], title = "Sensible heat flux (Wm⁻²)")
-hm = heatmap!(ax, Qs; colormap = :bwr)
-hidedecorations!(ax)
-save("sensible_heat_flux.png", fig)
-nothing #hide
-# ![](sensible_heat_flux.png)
 
-fig = Figure(size = (800, 400))
-ax = Axis(fig[1, 2], title = "Latent heat flux (Wm⁻²)")
-heatmap!(ax, Ql; colormap = :bwr)
-hidedecorations!(ax)
-save("latent_heat_flux.png", fig)
-nothing #hide
-# ![](latent_heat_flux.png)
+ax = Axis(fig[1, 1], title = "Sensible heat flux (W m⁻²)")
+heatmap!(ax, fluxes.sensible_heat; colormap = :bwr)
 
-fig = Figure(size = (800, 400))
-ax = Axis(fig[2, 1], title = "Zonal wind stress (Nm)")
-heatmap!(ax, τx; colormap = :bwr)
-hidedecorations!(ax)
-save("zonal_wind_stress.png", fig)
-nothing #hide
-# ![](zonal_wind_stress.png)
+ax = Axis(fig[1, 2], title = "Latent heat flux (W m⁻²)")
+heatmap!(ax, fluxes.latent_heat; colormap = :bwr)
 
-fig = Figure(size = (800, 400))
-ax = Axis(fig[2, 2], title = "Meridional wind stress (Nm)")
-heatmap!(ax, τy; colormap = :bwr)
-hidedecorations!(ax)
-save("meridional_wind_stress.png", fig)
-nothing #hide
-# ![](meridional_wind_stress.png)
+ax = Axis(fig[2, 1], title = "Zonal wind stress (N m)")
+heatmap!(ax, fluxes.x_momentum; colormap = :bwr)
 
-fig = Figure(size = (800, 400))
-ax = Axis(fig[3, 1], title = "Water vapor flux (kg m⁻²s⁻¹)")
+ax = Axis(fig[2, 2], title = "Meridional wind stress (N m)")
+heatmap!(ax, fluxes.y_momentum; colormap = :bwr)
+
+ax = Axis(fig[3, 1], title = "Water vapor flux (kg m⁻² s⁻¹)")
 heatmap!(ax, Mv; colormap = :bwr)
-hidedecorations!(ax)
-save("water_vapor_flux.png", fig)
-nothing #hide
-# ![](water_vapor_flux.png)
+
+save("fluxes.png", fig)
+# ![](fluxes.png)
 

examples/inspect_ecco_data.jl

@@ -14,13 +14,13 @@ using Oceananigans
 using CairoMakie
 using Printf
 
-using ClimaOcean: ECCO
+using ClimaOcean.DataWrangling.ECCO: ECCO_field
 
-# The function `ecco_field` provided by `ClimaOcean.DataWrangling.ECCO` will automatically
+# The function `ECCO_field` provided by `ClimaOcean.DataWrangling.ECCO` will automatically
 # download ECCO data if it doesn't already exist at the default location.
 
-T = ECCO.ecco_field(:temperature)
-S = ECCO.ecco_field(:salinity)
+T = ECCO_field(:temperature)
+S = ECCO_field(:salinity)
 
 # Next, we massage the ECCO data by inserting NaNs in "land cells", which
 # are diagnosed by having an unphysically low temperature.
@@ -50,14 +50,14 @@ grid = T.grid
 Nz = size(grid, 3)
 k = Observable(Nz)
 
-Tk = @lift interior(T, :, :, $k)
-Sk = @lift interior(S, :, :, $k)
+Tk = @lift view(T, :, :, $k)
+Sk = @lift view(S, :, :, $k)
 
 # Finally, we make a nice plot with a label that displays depth, colorbars,
 # and light gray within land cells.
 
-hmT = heatmap!(axT, λ, φ, Tk, nan_color=:lightgray, colorrange=(-2, 30), colormap=:thermal)
-hmS = heatmap!(axS, λ, φ, Sk, nan_color=:lightgray, colorrange=(31, 37), colormap=:haline)
+hmT = heatmap!(axT, Tk, nan_color=:lightgray, colorrange=(-2, 30), colormap=:thermal)
+hmS = heatmap!(axS, Sk, nan_color=:lightgray, colorrange=(31, 37), colormap=:haline)
 
 Colorbar(fig[1, 2], hmT, label="Temperature (ᵒC)")
 Colorbar(fig[2, 2], hmS, label="Salinity (psu)")

src/ClimaOcean.jl

@@ -9,15 +9,14 @@ export
     JRA55_prescribed_atmosphere,
     JRA55NetCDFBackend,
     ECCOMetadata,
-    ecco2_field,
     regrid_bathymetry,
     retrieve_bathymetry,
     stretched_vertical_faces,
     exponential_z_faces,
     PowerLawStretching, LinearStretching,
     exponential_z_faces,
     JRA55_field_time_series,
-    ecco_field, ECCOMetadata,
+    ECCO_field, ECCOMetadata,
     ocean_simulation,
     initialize!
 
@@ -72,7 +71,7 @@ using .OceanSeaIceModels
 using .OceanSimulations
 using .DataWrangling: JRA55, ECCO
 using ClimaOcean.DataWrangling.JRA55: JRA55_prescribed_atmosphere, JRA55NetCDFBackend
-using ClimaOcean.DataWrangling.ECCO: ecco_field, ECCOMetadata
+using ClimaOcean.DataWrangling.ECCO
 
 end # module
 

src/DataWrangling/DataWrangling.jl

@@ -73,7 +73,7 @@ end
 
 include("inpaint_mask.jl")
 include("JRA55.jl")
-include("ECCO.jl")
+include("ECCO/ECCO.jl")
 
 using .JRA55
 using .ECCO