adding ClimaSeaIce

Author: juliasloan25

experiments/ClimaEarth/components/ocean/clima_seaice.jl

@@ -0,0 +1,70 @@
+import Oceananigans as OC
+import ClimaOcean as CO
+import ClimaCoupler: Checkpointer, FieldExchanger, FluxCalculator, Interfacer, Utilities
+import ClimaComms
+import ClimaCore as CC
+import Thermodynamics as TD
+import ClimaOcean.EN4: download_dataset
+using KernelAbstractions: @kernel, @index, @inbounds
+
+"""
+    ClimaSeaIceSimulation{SIM, A, OPROP, REMAP}
+
+The ClimaCoupler simulation object used to run with Oceananigans.
+This type is used by the coupler to indicate that this simulation
+is an surface/ocean simulation for dispatch.
+
+It contains the following objects:
+- `sea_ice::SIM`: The ClimaSeaIce simulation object.
+- `area_fraction::A`: A ClimaCore Field representing the surface area fraction of this component model on the exchange grid.
+- `melting_speed::MS`: An constant characteristic speed for melting/freezing.
+- `remapping::REMAP`: Objects needed to remap from the exchange (spectral) grid to Oceananigans spaces.
+"""
+struct ClimaSeaIceSimulation{SIM,A,MS,REMAP} <: Interfacer.SeaIceModelSimulation
+    sea_ice::SIM
+    area_fraction::A
+    melting_speed::MS
+    remapping::REMAP
+end
+
+"""
+    ClimaSeaIceSimulation()
+
+Creates an OceananigansSimulation object containing a model, an integrator, and
+a surface area fraction field.
+This type is used to indicate that this simulation is an ocean simulation for
+dispatch in coupling.
+
+Specific details about the default model configuration
+can be found in the documentation for `ClimaOcean.ocean_simulation`.
+"""
+function ClimaSeaIceSimulation(area_fraction, ocean; output_dir)
+    # Initialize the sea ice with the same grid as the ocean
+    # Initially no sea ice is present
+    grid = ocean.ocean.model.grid
+    advection = ocean.ocean.model.advection.T
+    sea_ice = CO.sea_ice_simulation(grid, ocean.ocean; advection)
+
+    melting_speed = 1e-4
+    remapping = ocean.remapping
+
+    # Before version 0.96.22, the NetCDFWriter was broken on GPU
+    if arch isa OC.CPU || pkgversion(OC) >= v"0.96.22"
+        # TODO maybe this is broken
+        # Save all tracers and velocities to a NetCDF file at daily frequency
+        outputs = prognostic_fields(sea_ice.model)
+        netcdf_writer = OC.NetCDFWriter(
+            sea_ice.model,
+            outputs;
+            schedule = OC.TimeInterval(86400), # Daily output
+            filename = joinpath(output_dir, "seaice_diagnostics.nc"),
+            indices = (:, :, grid.Nz),
+            overwrite_existing = true,
+            array_type = Array{Float32},
+        )
+        sea_ice.output_writers[:diagnostics] = netcdf_writer
+    end
+
+    sim = ClimaSeaIceSimulation(sea_ice, area_fraction, melting_speed, remapping)
+    return sim
+end

experiments/ClimaEarth/components/ocean/oceananigans.jl

@@ -88,6 +88,7 @@ function OceananigansSimulation(
         active_cells_map = true,
     )
 
+    # TODO use_restoring only if not using ClimaSeaIce
     use_restoring = start_date + Dates.Month(1) < stop_date
 
     if use_restoring
@@ -108,7 +109,21 @@ function OceananigansSimulation(
     end
 
     # Create ocean simulation
-    ocean = CO.ocean_simulation(grid; forcing)
+    free_surface = SplitExplicitFreeSurface(grid; substeps = 70)
+    momentum_advection = WENOVectorInvariant(order = 5)
+    tracer_advection = WENO(order = 5)
+    eddy_closure = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew = 2e3, κ_symmetric = 2e3)
+    vertical_mixing = ClimaOcean.OceanSimulations.default_ocean_closure()
+    horizontal_viscosity = HorizontalScalarDiffusivity(ν = 4000)
+
+    ocean = CO.ocean_simulation(
+        grid;
+        forcing,
+        momentum_advection,
+        tracer_advection,
+        free_surface,
+        closure = (eddy_closure, horizontal_viscosity, vertical_mixing),
+    )
 
     # Set initial condition to EN4 state estimate at start_date
     OC.set!(ocean.model, T = en4_temperature[1], S = en4_salinity[1])
@@ -147,8 +162,6 @@ function OceananigansSimulation(
 
     # Before version 0.96.22, the NetCDFWriter was broken on GPU
     if arch isa OC.CPU || pkgversion(OC) >= v"0.96.22"
-        # TODO: Add more diagnostics, make them dependent on simulation duration, take
-        # monthly averages
         # Save all tracers and velocities to a NetCDF file at daily frequency
         outputs = merge(ocean.model.tracers, ocean.model.velocities)
         netcdf_writer = OC.NetCDFWriter(
@@ -253,10 +266,8 @@ Interfacer.get_field(sim::OceananigansSimulation, ::Val{:roughness_momentum}) =
     Float32(5.8e-5)
 Interfacer.get_field(sim::OceananigansSimulation, ::Val{:beta}) = Float32(1)
 Interfacer.get_field(sim::OceananigansSimulation, ::Val{:emissivity}) = Float32(0.97)
-Interfacer.get_field(sim::OceananigansSimulation, ::Val{:surface_direct_albedo}) =
-    Float32(0.06)
-Interfacer.get_field(sim::OceananigansSimulation, ::Val{:surface_diffuse_albedo}) =
-    Float32(0.06)
+Interfacer.get_field(sim::OceananigansSimulation, ::Val{:surface_direct_albedo}) = Float32(0.011)
+Interfacer.get_field(sim::OceananigansSimulation, ::Val{:surface_diffuse_albedo}) = Float32(0.069)
 
 # NOTE: This is 3D, but it will be remapped to 2D
 Interfacer.get_field(sim::OceananigansSimulation, ::Val{:surface_temperature}) =