diff --git a/experiments/omip_prototype/download_data.jl b/experiments/omip_prototype/download_data.jl new file mode 100644 index 00000000..2881bcd7 --- /dev/null +++ b/experiments/omip_prototype/download_data.jl @@ -0,0 +1,8 @@ +using ClimaOcean +using ClimaOcean.JRA55 +using ClimaOcean.DataWrangling: download_dataset + +atmosphere = JRA55PrescribedAtmosphere(; dir="forcing_data/", + dataset=JRA55MultipleYears(), + backend=JRA55NetCDFBackend(10), + include_rivers_and_icebergs=true) \ No newline at end of file diff --git a/experiments/omip_prototype/sixth_degree_omip.jl b/experiments/omip_prototype/sixth_degree_omip.jl new file mode 100644 index 00000000..6a482bb0 --- /dev/null +++ b/experiments/omip_prototype/sixth_degree_omip.jl @@ -0,0 +1,160 @@ +using ClimaOcean +using ClimaSeaIce +using Oceananigans +using Oceananigans.Grids +using Oceananigans.Units +using Oceananigans.OrthogonalSphericalShellGrids +using ClimaOcean.OceanSimulations +using ClimaOcean.ECCO +using ClimaOcean.JRA55 +using ClimaOcean.DataWrangling +using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium +using Printf +using Dates +using CUDA + +import Oceananigans.OutputWriters: checkpointer_address + +function synch!(clock1::Clock, clock2) + # Synchronize the clocks + clock1.time = clock2.time + clock1.iteration = clock2.iteration + clock1.last_Δt = clock2.last_Δt +end + +synch!(model1, model2) = synch!(model1.clock, model2.clock) + +arch = GPU() +r_faces = ClimaOcean.exponential_z_faces(; Nz=60, depth=6200) +z_faces = MutableVerticalDiscretization(r_faces) + +Nx = 2160 # longitudinal direction +Ny = 1080 # meridional direction +Nz = length(r_faces) - 1 + +grid = TripolarGrid(arch; + size = (Nx, Ny, Nz), + z = z_faces, + halo = (7, 7, 7)) + +bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1) +grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true) + +##### +##### A Propgnostic Ocean model +##### + +using Oceananigans.TurbulenceClosures: ExplicitTimeDiscretization +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEVerticalDiffusivity, CATKEMixingLength, CATKEEquation + +momentum_advection = WENOVectorInvariant() +tracer_advection = WENO(order=7) + +free_surface = SplitExplicitFreeSurface(grid; substeps=70) + +mixing_length = CATKEMixingLength(Cᵇ=0.01) +turbulent_kinetic_energy_equation = CATKEEquation(Cᵂϵ=1.0) + +catke_closure = CATKEVerticalDiffusivity(ExplicitTimeDiscretization(); mixing_length, turbulent_kinetic_energy_equation) +closure = (catke_closure, VerticalScalarDiffusivity(κ=1e-5, ν=1e-5)) + +ocean = ocean_simulation(grid; Δt=1minutes, + momentum_advection, + tracer_advection, + free_surface, + closure) + +dataset = ECCO4Monthly() + +set!(ocean.model, T=Metadatum(:temperature; dataset), + S=Metadatum(:salinity; dataset)) + +##### +##### A Prognostic Sea-ice model +##### + +# Remember to pass the SSS as a bottom bc to the sea ice! +SSS = view(ocean.model.tracers.S.data, :, :, grid.Nz) +bottom_heat_boundary_condition = IceWaterThermalEquilibrium(SSS) + +# Default sea-ice dynamics +sea_ice_dynamics = ClimaOcean.SeaIceSimulations.default_sea_ice_dynamics(grid; ocean) + +sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition, + dynamics = sea_ice_dynamics, + advection=WENO(order=7)) + +set!(sea_ice.model, h=Metadatum(:sea_ice_thickness; dataset), + ℵ=Metadatum(:sea_ice_concentration; dataset)) + +##### +##### A Prescribed Atmosphere model +##### + +dir = "./forcing_data" +dataset = MultiYearJRA55() +backend = JRA55NetCDFBackend(40) + +atmosphere = JRA55PrescribedAtmosphere(arch; dir, dataset, backend, include_rivers_and_icebergs=true) +radiation = Radiation() + +##### +##### An ocean-sea ice coupled model +##### + +omip = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +omip = Simulation(omip, Δt=20, stop_time=60days) + +# Figure out the outputs.... + +checkpointer_address(::SeaIceModel) = "SeaIceModel" + +ocean.output_writers[:checkpointer] = Checkpointer(ocean.model, + schedule = IterationInterval(10000), + prefix = "ocean_checkpoint", + overwrite_existing = true) + +sea_ice.output_writers[:checkpointer] = Checkpointer(sea_ice.model, + schedule = IterationInterval(10000), + prefix = "sea_ice_checkpoint", + overwrite_existing = true) + +wall_time = Ref(time_ns()) + +using Statistics + +function progress(sim) + sea_ice = sim.model.sea_ice + ocean = sim.model.ocean + hmax = maximum(sea_ice.model.ice_thickness) + ℵmax = maximum(sea_ice.model.ice_concentration) + Tmax = maximum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + Tmin = minimum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature) + umax = maximum(ocean.model.velocities.u) + vmax = maximum(ocean.model.velocities.v) + wmax = maximum(ocean.model.velocities.w) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt)) + msg2 = @sprintf("max(h): %.2e m, max(ℵ): %.2e ", hmax, ℵmax) + msg4 = @sprintf("extrema(T): (%.2f, %.2f) ᵒC, ", Tmax, Tmin) + msg5 = @sprintf("maximum(u): (%.2f, %.2f, %.2f) m/s, ", umax, vmax, wmax) + msg6 = @sprintf("wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg4 * msg5 * msg6 + + wall_time[] = time_ns() + + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(omip, progress, IterationInterval(50)) + +run!(omip) + +omip.Δt = 10minutes +omip.stop_time = 58 * 365days + +run!(omip) diff --git a/src/DataWrangling/JRA55/JRA55_field_time_series.jl b/src/DataWrangling/JRA55/JRA55_field_time_series.jl index e7bbc78b..3a346382 100644 --- a/src/DataWrangling/JRA55/JRA55_field_time_series.jl +++ b/src/DataWrangling/JRA55/JRA55_field_time_series.jl @@ -189,7 +189,6 @@ function set!(fts::JRA55NetCDFFTSMultipleYears, backend=fts.backend) LX, LY, LZ = location(fts) i₁, i₂, j₁, j₂, TX = compute_bounding_indices(nothing, nothing, fts.grid, LX, LY, λc, φc) - if issorted(nn) data = ds[name][i₁:i₂, j₁:j₂, nn] else @@ -458,7 +457,7 @@ function JRA55FieldTimeSeries(metadata::JRA55Metadata, architecture=CPU(), FT=Fl ds = Dataset(filepath) data = ds[shortname][i₁:i₂, j₁:j₂, time_indices_in_memory] close(ds) - + copyto!(interior(fts, :, :, 1, :), data) fill_halo_regions!(fts) diff --git a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl index 1532a3be..dff50271 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl @@ -261,8 +261,8 @@ end @inbounds begin Ts = surface_temperature[i, j, kᴺ] Ts = convert_to_kelvin(sea_ice_properties.temperature_units, Ts) - ℵi = ice_concentration[i, j, 1] - + ℵi = ice_concentration[i, j, kᴺ] + Qs = downwelling_radiation.Qs[i, j, 1] Qℓ = downwelling_radiation.Qℓ[i, j, 1] Qc = atmosphere_sea_ice_fluxes.sensible_heat[i, j, 1] # sensible or "conductive" heat flux diff --git a/src/SeaIceSimulations.jl b/src/SeaIceSimulations.jl index 08f3de9e..d3650501 100644 --- a/src/SeaIceSimulations.jl +++ b/src/SeaIceSimulations.jl @@ -17,6 +17,7 @@ using Oceananigans.Operators using ClimaSeaIce using ClimaSeaIce: SeaIceModel, SlabSeaIceThermodynamics, PhaseTransitions, ConductiveFlux using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium +using ClimaSeaIce.Rheologies using ClimaSeaIce.SeaIceDynamics: SplitExplicitSolver, SemiImplicitStress, SeaIceMomentumEquation, StressBalanceFreeDrift using ClimaSeaIce.Rheologies: IceStrength, ElastoViscoPlasticRheology diff --git a/test/test_jra55.jl b/test/test_jra55.jl index 5d45ba47..481fdd05 100644 --- a/test/test_jra55.jl +++ b/test/test_jra55.jl @@ -1,5 +1,6 @@ include("runtests_setup.jl") +using ClimaOcean.JRA55 using ClimaOcean.JRA55: download_JRA55_cache using ClimaOcean.OceanSeaIceModels: PrescribedAtmosphere