diff --git a/Project.toml b/Project.toml index b9816ea7b..fa00c3a66 100644 --- a/Project.toml +++ b/Project.toml @@ -25,6 +25,7 @@ Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" Scratch = "6c6a2e73-6563-6170-7368-637461726353" SeawaterPolynomials = "d496a93d-167e-4197-9f49-d3af4ff8fe40" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" @@ -33,6 +34,7 @@ Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" [weakdeps] +MITgcm = "dce5fa8e-68ce-4431-a242-9469c69627a0" ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3" Breeze = "660aa2fb-d4c8-4359-a52c-9c057bc511da" CDSAPI = "8a7b9de3-9c00-473e-88b4-7eccd7ef2fea" @@ -49,6 +51,7 @@ NumericalEarthArchGDALExt = "ArchGDAL" NumericalEarthBreezeExt = "Breeze" NumericalEarthCDSAPIExt = "CDSAPI" NumericalEarthCopernicusMarineExt = "CopernicusMarine" +NumericalEarthMITgcmExt = "MITgcm" NumericalEarthReactantExt = "Reactant" NumericalEarthSpeedyWeatherExt = ["SpeedyWeather", "ConservativeRegridding"] NumericalEarthVerosExt = ["PythonCall", "CondaPkg"] @@ -75,6 +78,7 @@ Glob = "1" ImageMorphology = "0.4" JLD2 = "0.4, 0.5, 0.6" KernelAbstractions = "0.9" +MITgcm = "0.5" MeshArrays = "0.3, 0.4, 0.5" NCDatasets = "0.12, 0.13, 0.14" Oceananigans = "0.107, 0.108" diff --git a/examples/mitgcm/.gitignore b/examples/mitgcm/.gitignore new file mode 100644 index 000000000..914de2476 --- /dev/null +++ b/examples/mitgcm/.gitignore @@ -0,0 +1,2 @@ +# MITgcm build artifacts +build/ diff --git a/examples/mitgcm/mitgcm_config/code/SIZE.h b/examples/mitgcm/mitgcm_config/code/SIZE.h new file mode 100644 index 000000000..695fce260 --- /dev/null +++ b/examples/mitgcm/mitgcm_config/code/SIZE.h @@ -0,0 +1,63 @@ +CBOP +C !ROUTINE: SIZE.h +C !INTERFACE: +C include SIZE.h +C !DESCRIPTION: \bv +C *==========================================================* +C | SIZE.h Declare size of underlying computational grid. +C *==========================================================* +C | The design here supports a three-dimensional model grid +C | with indices I,J and K. The three-dimensional domain +C | is comprised of nPx*nSx blocks (or tiles) of size sNx +C | along the first (left-most index) axis, nPy*nSy blocks +C | of size sNy along the second axis and one block of size +C | Nr along the vertical (third) axis. +C | Blocks/tiles have overlap regions of size OLx and OLy +C | along the dimensions that are subdivided. +C *==========================================================* +C \ev +C +C Voodoo numbers controlling data layout: +C sNx :: Number of X points in tile. +C sNy :: Number of Y points in tile. +C OLx :: Tile overlap extent in X. +C OLy :: Tile overlap extent in Y. +C nSx :: Number of tiles per process in X. +C nSy :: Number of tiles per process in Y. +C nPx :: Number of processes to use in X. +C nPy :: Number of processes to use in Y. +C Nx :: Number of points in X for the full domain. +C Ny :: Number of points in Y for the full domain. +C Nr :: Number of points in vertical direction. +CEOP + INTEGER sNx + INTEGER sNy + INTEGER OLx + INTEGER OLy + INTEGER nSx + INTEGER nSy + INTEGER nPx + INTEGER nPy + INTEGER Nx + INTEGER Ny + INTEGER Nr + PARAMETER ( + & sNx = 30, + & sNy = 20, + & OLx = 2, + & OLy = 2, + & nSx = 3, + & nSy = 2, + & nPx = 1, + & nPy = 1, + & Nx = sNx*nSx*nPx, + & Ny = sNy*nSy*nPy, + & Nr = 15) + +C MAX_OLX :: Set to the maximum overlap region size of any array +C MAX_OLY that will be exchanged. Controls the sizing of exch +C routine buffers. + INTEGER MAX_OLX + INTEGER MAX_OLY + PARAMETER ( MAX_OLX = OLx, + & MAX_OLY = OLy ) diff --git a/examples/mitgcm/mitgcm_config/code/packages.conf b/examples/mitgcm/mitgcm_config/code/packages.conf new file mode 100644 index 000000000..fd5c8bdcf --- /dev/null +++ b/examples/mitgcm/mitgcm_config/code/packages.conf @@ -0,0 +1,7 @@ +#-- packages for coupled library mode with KPP vertical mixing +# EXF, CAL, PROFILES are disabled: surface forcing comes from the host (Julia). +gfd +cd_code +gmredi +kpp +diagnostics diff --git a/examples/mitgcm/mitgcm_config/input/POLY3.COEFFS b/examples/mitgcm/mitgcm_config/input/POLY3.COEFFS new file mode 100644 index 000000000..44c919c8a --- /dev/null +++ b/examples/mitgcm/mitgcm_config/input/POLY3.COEFFS @@ -0,0 +1,31 @@ +15 + 13.4965485754395 32.6000000000000 24.5481924347446 + 13.4882386805755 32.6000000000000 24.8185506216191 + 13.4764025555903 32.6000000000000 25.2009783814515 + 8.4690806223310 35.1500000000000 28.6417699272743 + 8.4510089195050 35.1500000000000 29.3928269010235 + 5.9397585609655 34.9000000000000 30.5632984773008 + 4.4248880992503 34.9000000000000 31.9847214990702 + 2.9114111861653 34.7500000000000 33.5047690245101 + 1.8932862598040 34.8000000000000 35.3412270497367 + 1.3669344372630 34.8000000000000 37.2987636266591 + 1.3272007124547 34.8000000000000 39.4052882836211 + 0.7905919628270 34.8000000000000 41.7881156235655 + 0.7358877560078 34.7500000000000 44.2594853649941 + 0.1858449458101 34.7500000000000 47.0435345001690 + 0.5950242216911 34.7500000000000 49.8352041101475 +-2.0248441577781e-01 7.7327632363592e-01 -4.6235907431162e-03 -1.7855486745386e-03 6.7161214818854e-05 1.9780689793512e-05 1.8596621330081e-05 2.7193175334807e-06 7.7992579107390e-06 +-2.0356900098784e-01 7.7287763354454e-01 -4.6078685547211e-03 -1.7797434187354e-03 6.9849851237856e-05 1.9700367626079e-05 1.8541107779796e-05 2.6814780158476e-06 7.6950638206688e-06 +-2.0510288558783e-01 7.7231413929174e-01 -4.5856232199577e-03 -1.7715237666829e-03 7.3641940080245e-05 1.9586693371392e-05 1.8462890416846e-05 2.6279856728835e-06 7.5480688984286e-06 +-1.6315278510835e-01 7.8362068555215e-01 -4.9742804881192e-03 -2.2358797318897e-03 -4.1623282245313e-04 2.3781382246111e-05 2.8390410822447e-05 1.8232580281910e-05 8.7007573824530e-05 +-1.6668703762104e-01 7.8244679887675e-01 -4.9199187906716e-03 -2.2150360379047e-03 -4.0253524411189e-04 2.3370269790246e-05 2.8117133690396e-05 1.7907936881655e-05 8.5353460010715e-05 +-1.4434702349368e-01 7.8636558552016e-01 -5.2437823005010e-03 -2.4430632367892e-03 -5.0604107382534e-04 3.0094061345352e-05 3.4995169663676e-05 2.6906045499581e-05 1.5782380566825e-04 +-1.3422530714233e-01 7.8774544403399e-01 -5.4057628884322e-03 -2.5425581389254e-03 -2.6671549575359e-04 3.4876625315646e-05 3.8742849620680e-05 2.2351400258486e-05 1.1390437277313e-04 +-1.2494903729136e-01 7.8886002954110e-01 -5.5923621375103e-03 -2.6509369714478e-03 -1.8969464504401e-04 4.2909291159896e-05 4.4780770266509e-05 2.8883428644997e-05 2.1006556257136e-04 +-1.2342388397754e-01 7.8859195674364e-01 -5.6185001971557e-03 -2.6849010357721e-03 4.1405949512546e-06 4.6245953907873e-05 4.6728604344334e-05 1.9161965321620e-05 1.1115190475704e-04 +-1.2890117744188e-01 7.8668953911886e-01 -5.5282092351178e-03 -2.6686338397672e-03 4.5996647260460e-06 4.6776611042599e-05 4.7296302315891e-05 2.2513552111924e-05 1.7041983515548e-04 +-1.4114412009280e-01 7.8313783149598e-01 -5.3256878521832e-03 -2.5926220481035e-03 2.2428825626142e-05 4.4369484204344e-05 4.5782525710578e-05 2.1123578875789e-05 1.6000623056750e-04 +-1.4945904203962e-01 7.8048541390059e-01 -5.1851051549971e-03 -2.5533751190741e-03 1.1992432484732e-04 4.4094804560539e-05 4.5370934503017e-05 1.2282632879782e-05 6.5241387954551e-05 +-1.6392270705111e-01 7.7627008453271e-01 -4.9448680074238e-03 -2.4654573880316e-03 9.3717129874410e-05 4.1174124268602e-05 4.3790593946149e-05 1.6928793605944e-05 1.3505403387493e-04 +-1.7490789271395e-01 7.7294291979081e-01 -4.7570017089526e-03 -2.4106158149213e-03 1.6880513678366e-04 4.0088827917669e-05 4.2910040951280e-05 8.3333879620506e-06 4.1916874947312e-05 +-1.9573834386012e-01 7.6706502751275e-01 -4.4224744309248e-03 -2.2656921102643e-03 2.0327244553310e-04 3.6548670987625e-05 4.0510685127073e-05 3.9178994487299e-06 9.2592973542915e-06 diff --git a/examples/mitgcm/mitgcm_config/input/data b/examples/mitgcm/mitgcm_config/input/data new file mode 100644 index 000000000..f4afa8848 --- /dev/null +++ b/examples/mitgcm/mitgcm_config/input/data @@ -0,0 +1,79 @@ +# ==================== +# | Model parameters | +# ==================== +# +# Based on global_oce_latlon, modified for coupled library mode: +# - EXF disabled (forcing from Julia) +# - KPP vertical mixing enabled +# - Timestep set by Julia at runtime via set_timestep! +# +# Continuous equation parameters + &PARM01 + tRef = 15*20., + sRef = 15*35., + viscAr=1.E-3, + viscAh=5.E5, + diffKhT=1.E3, + diffKrT=3.E-5, + diffKhS=1.E3, + diffKrS=3.E-5, + rhoConst=1035., + rotationPeriod=86400., + gravity=9.81, + eosType = 'POLY3', + ivdc_kappa=0.0, + implicitDiffusion=.TRUE., + implicitViscosity=.TRUE., + allowFreezing=.TRUE., + useRealFreshWaterFlux=.FALSE., + useCDscheme=.FALSE., + staggerTimeStep=.FALSE., +# useJamartWetPoints=.TRUE., + useEnergyConservingCoriolis=.TRUE., +# vectorInvariantMomentum=.TRUE., +# tempAdvScheme=33, +# saltAdvScheme=33, +# turn on looped cells + hFacMin=.05, +# hFacMindr=50., +# set precision of data files + readBinaryPrec=32, + usesinglecpuio=.TRUE., + & + +# Elliptic solver parameters + &PARM02 + cg2dMaxIters=500, + cg2dTargetResidual=1.E-13, + & + +# Time stepping parameters + &PARM03 + nIter0 = 0, + nTimeSteps = 2000000000, + deltaTmom = 1200.0, + deltaTtracer= 1200.0, + deltaTClock = 1200.0, + abEps = 0.1, + pChkptFreq = 311040000., + dumpFreq = 311040000., + monitorFreq = 864000., + & + +# Gridding parameters + &PARM04 + usingSphericalPolarGrid=.TRUE., + delR= 50., 70., 100., 140., 190., + 240., 290., 340., 390., 440., + 490., 540., 590., 640., 690., + ygOrigin=-80., + dySpacing=4., + dxSpacing=4., + & + +# Input datasets + &PARM05 + bathyFile= 'bathymetry.bin', + hydrogThetaFile='lev_t.bin', + hydrogSaltFile= 'lev_s.bin', + & diff --git a/examples/mitgcm/mitgcm_config/input/data.diagnostics b/examples/mitgcm/mitgcm_config/input/data.diagnostics new file mode 100644 index 000000000..7e11e925c --- /dev/null +++ b/examples/mitgcm/mitgcm_config/input/data.diagnostics @@ -0,0 +1,5 @@ +# Diagnostics disabled for shared-library mode (timestep may change at runtime) + &DIAGNOSTICS_LIST + & + &DIAG_STATIS_PARMS + & diff --git a/examples/mitgcm/mitgcm_config/input/data.gmredi b/examples/mitgcm/mitgcm_config/input/data.gmredi new file mode 100644 index 000000000..802719f25 --- /dev/null +++ b/examples/mitgcm/mitgcm_config/input/data.gmredi @@ -0,0 +1,9 @@ +# GM+Redi package parameters: + &GM_PARM01 + GM_background_K = 1.e+3, + GM_taper_scheme = 'gkw91', + GM_maxSlope = 1.e-2, + GM_Kmin_horiz = 50., + GM_Scrit = 4.e-3, + GM_Sd = 1.e-3, + & diff --git a/examples/mitgcm/mitgcm_config/input/data.kpp b/examples/mitgcm/mitgcm_config/input/data.kpp new file mode 100644 index 000000000..63249f681 --- /dev/null +++ b/examples/mitgcm/mitgcm_config/input/data.kpp @@ -0,0 +1,3 @@ +# KPP vertical mixing parameters + &KPP_PARM01 + & diff --git a/examples/mitgcm/mitgcm_config/input/data.pkg b/examples/mitgcm/mitgcm_config/input/data.pkg new file mode 100644 index 000000000..d5f77f64e --- /dev/null +++ b/examples/mitgcm/mitgcm_config/input/data.pkg @@ -0,0 +1,8 @@ +# Package configuration for coupled library mode. +# EXF is disabled: surface forcing comes from the host (Julia). +# KPP is enabled for vertical mixing. + &PACKAGES + useGMRedi = .FALSE., + useKPP = .TRUE., + useDiagnostics = .TRUE., + & diff --git a/examples/mitgcm/mitgcm_ocean_forced_simulation.jl b/examples/mitgcm/mitgcm_ocean_forced_simulation.jl new file mode 100644 index 000000000..6b22ee9ec --- /dev/null +++ b/examples/mitgcm/mitgcm_ocean_forced_simulation.jl @@ -0,0 +1,93 @@ +# # MITgcm Ocean Forced Simulation +# +# This script runs MITgcm (as a shared library) on the `global_oce_latlon` +# grid: 90×40×15, 4°×4° lat-lon, 80S–80N, with JRA55 atmospheric forcing +# applied through NumericalEarth's coupling framework. +# +# Initial conditions and bathymetry come from the MITgcm tutorial binary files. +# Physics: KPP vertical mixing, GM/Redi, CD scheme. +# +# Run for 2 years to compare against `numerical_earth_example.jl`. + +using NumericalEarth +using MITgcm +using Oceananigans +using Oceananigans.Units +using Printf +using Statistics + +# ## Build and load MITgcm + +example_dir = @__DIR__ +config_dir = joinpath(example_dir, "mitgcm_config") +code_dir = joinpath(config_dir, "code") +input_dir = joinpath(config_dir, "input") +build_dir = joinpath(example_dir, "build_jra55") + +mitgcm_dir = get(ENV, "MITGCM_DIR", "") +if !isdir(mitgcm_dir) + @info "Downloading MITgcm source..." + mitgcm_dir = MITgcm.download_mitgcm_source() +end + +@info "Building MITgcm..." mitgcm_dir code_dir input_dir +ocean = MITgcmOceanSimulation(mitgcm_dir; + output_dir = build_dir, + code_dir, + input_dir, + verbose = false) + +lib = ocean.library +Nx, Ny, Nr = lib.dims.Nx, lib.dims.Ny, lib.dims.Nr +@info "MITgcm loaded" Nx Ny Nr + +# ## Build coupled model with JRA55 forcing + +atmos = JRA55PrescribedAtmosphere() +radiation = Radiation(ocean_emissivity=0.0, sea_ice_emissivity=0.0) + +coupled_model = OceanSeaIceModel(ocean, nothing; atmosphere=atmos, radiation) + +Δt = 1200 +stop_time = 360days +simulation = Simulation(coupled_model; Δt, stop_time) + +# ## Progress callback + +wall_time = Ref(time_ns()) + +function progress(sim) + ocean = sim.model.ocean + lib = ocean.library + mtime = MITgcm.get_time(lib) + + refresh_state!(ocean) + mask = view(ocean.hfacc, :, :, 1) + ocean_sst = ocean.theta[:, :, 1][mask .> 0] + ocean_sss = ocean.salt[:, :, 1][mask .> 0] + eta = ocean.etan[mask .> 0] + + elapsed = 1e-9 * (time_ns() - wall_time[]) + + @printf("iter %6d | day %7.1f | SST: [%6.2f, %5.2f] | SSS: [%5.2f, %5.2f] | η: [%.4f, %.4f] mean=%.4e | wall: %s\n", + iteration(sim), mtime / 86400, + minimum(ocean_sst), maximum(ocean_sst), + minimum(ocean_sss), maximum(ocean_sss), + minimum(eta), maximum(eta), mean(eta), + prettytime(elapsed)) + + wall_time[] = time_ns() + return nothing +end + +add_callback!(simulation, progress, TimeInterval(10days)) + +# ## Run + +@info "Running MITgcm + JRA55 simulation..." Δt stop_time +wall_time[] = time_ns() +run!(simulation) + +@info "Simulation complete. Finalizing..." +MITgcm.finalize!(lib) +@info "Done!" diff --git a/ext/NumericalEarthMITgcmExt/NumericalEarthMITgcmExt.jl b/ext/NumericalEarthMITgcmExt/NumericalEarthMITgcmExt.jl new file mode 100644 index 000000000..d843d367b --- /dev/null +++ b/ext/NumericalEarthMITgcmExt/NumericalEarthMITgcmExt.jl @@ -0,0 +1,10 @@ +module NumericalEarthMITgcmExt + +using NumericalEarth +using MITgcm +using Oceananigans + +include("mitgcm_ocean_coupling.jl") +include("mitgcm_state_exchanger.jl") + +end # module NumericalEarthMITgcmExt diff --git a/ext/NumericalEarthMITgcmExt/mitgcm_ocean_coupling.jl b/ext/NumericalEarthMITgcmExt/mitgcm_ocean_coupling.jl new file mode 100644 index 000000000..49f78b880 --- /dev/null +++ b/ext/NumericalEarthMITgcmExt/mitgcm_ocean_coupling.jl @@ -0,0 +1,96 @@ +using Oceananigans: LatitudeLongitudeGrid, Periodic, Bounded, Flat, CPU, Field, + Face, Center + +import Oceananigans.TimeSteppers: time_step!, initialize! +import Oceananigans.Architectures: architecture +import Oceananigans.Fields: set! + +using NumericalEarth.EarthSystemModels: EarthSystemModel +import NumericalEarth.EarthSystemModels: default_nan_checker, + reference_density, + heat_capacity, + ocean_temperature, + ocean_salinity, + ocean_surface_salinity, + ocean_surface_temperature, + ocean_surface_velocities + +import Base: eltype + +using MITgcm: MITgcmOceanSimulation, MITgcmLibrary, + set_timestep!, get_timestep, step!, refresh_state!, + set_fu!, set_fv!, set_qnet!, set_empmr!, set_qsw!, set_saltflux! + +default_nan_checker(model::EarthSystemModel{<:Any, <:Any, <:MITgcmOceanSimulation}) = nothing +initialize!(::MITgcmOceanSimulation) = nothing + +architecture(model::EarthSystemModel{<:Any, <:Any, <:MITgcmOceanSimulation}) = CPU() +eltype(model::EarthSystemModel{<:Any, <:Any, <:MITgcmOceanSimulation}) = Float64 + +##### +##### Time stepping +##### + +""" + time_step!(ocean::MITgcmOceanSimulation, Δt) + +Advance the MITgcm ocean by `Δt` seconds. Sets all MITgcm timestep parameters +(deltaT, deltaTMom, deltaTClock, deltaTFreeSurf, dTtracerLev) to `Δt`, then +takes a single forward step. +""" +function time_step!(ocean::MITgcmOceanSimulation, Δt) + set_timestep!(ocean.library, Float64(Δt)) + step!(ocean.library) + refresh_state!(ocean) + return nothing +end + +reference_density(ocean::MITgcmOceanSimulation) = ocean.reference_density +heat_capacity(ocean::MITgcmOceanSimulation) = ocean.heat_capacity + +ocean_temperature(ocean::MITgcmOceanSimulation) = ocean.theta +ocean_salinity(ocean::MITgcmOceanSimulation) = ocean.salt + +ocean_surface_temperature(ocean::MITgcmOceanSimulation) = view(ocean.theta, :, :, 1) +ocean_surface_salinity(ocean::MITgcmOceanSimulation) = view(ocean.salt, :, :, 1) + +function ocean_surface_velocities(ocean::MITgcmOceanSimulation) + u_surface = view(ocean.uvel, :, :, 1) + v_surface = view(ocean.vvel, :, :, 1) + return u_surface, v_surface +end + +""" + surface_grid(ocean::MITgcmOceanSimulation) + +Construct an Oceananigans `LatitudeLongitudeGrid` from MITgcm's grid arrays. +Assumes the MITgcm grid is a regular lat-lon grid. +""" +function surface_grid(ocean::MITgcmOceanSimulation) + xc = ocean.xc + yc = ocean.yc + Nx = size(xc, 1) + Ny = size(xc, 2) + + # Extract 1D coordinate vectors (MITgcm stores 2D, but for lat-lon they are separable) + λc = xc[:, 1] # longitude centers + φc = yc[1, :] # latitude centers + + # Compute face (edge) coordinates from centers + # Assume uniform spacing (can be extended for stretched grids) + Δλ = length(λc) > 1 ? λc[2] - λc[1] : 1.0 + Δφ = length(φc) > 1 ? φc[2] - φc[1] : 1.0 + + λf = vcat(λc[1] - Δλ/2, [λc[i] + Δλ/2 for i in 1:Nx]) + φf = vcat(φc[1] - Δφ/2, [φc[j] + Δφ/2 for j in 1:Ny]) + + # Determine topology: periodic in longitude if grid spans ~360 degrees + lon_span = λf[end] - λf[1] + TX = lon_span >= 359.0 ? Periodic : Bounded + + return LatitudeLongitudeGrid(size = (Nx, Ny), + longitude = λf, + latitude = φf, + topology = (TX, Bounded, Flat), + halo = (2, 2)) +end diff --git a/ext/NumericalEarthMITgcmExt/mitgcm_state_exchanger.jl b/ext/NumericalEarthMITgcmExt/mitgcm_state_exchanger.jl new file mode 100644 index 000000000..28e45e311 --- /dev/null +++ b/ext/NumericalEarthMITgcmExt/mitgcm_state_exchanger.jl @@ -0,0 +1,154 @@ +using NumericalEarth.Oceans + +import NumericalEarth.EarthSystemModels.InterfaceComputations: + net_fluxes, + initialize!, + ComponentExchanger, + exchange_grid, + compute_sea_ice_ocean_fluxes! + +import NumericalEarth.EarthSystemModels: + interpolate_state!, + update_net_fluxes! + +import NumericalEarth.Oceans: get_radiative_forcing + +using Oceananigans.Fields: interior + +using MITgcm: MITgcmOceanSimulation, + set_fu!, set_fv!, set_qnet!, set_empmr!, set_qsw!, set_saltflux! + +##### +##### ComponentExchanger +##### + +function ComponentExchanger(ocean::MITgcmOceanSimulation, grid) + # MITgcm uses Arakawa C-grid: u at (Face, Center), v at (Center, Face) + state = (; u = Field{Face, Center, Nothing}(grid), + v = Field{Center, Face, Nothing}(grid), + T = Field{Center, Center, Nothing}(grid), + S = Field{Center, Center, Nothing}(grid)) + + return ComponentExchanger(state, nothing) +end + +exchange_grid(atmosphere, ocean::MITgcmOceanSimulation, sea_ice) = surface_grid(ocean) + +##### +##### Net fluxes container +##### + +@inline function net_fluxes(ocean::MITgcmOceanSimulation) + grid = surface_grid(ocean) + u = Field{Face, Center, Nothing}(grid) + v = Field{Center, Face, Nothing}(grid) + T = Field{Center, Center, Nothing}(grid) + S = Field{Center, Center, Nothing}(grid) + + return (; u, v, T, S) +end + +# ============================================================ +# State interpolation +# ============================================================ + +function interpolate_state!(exchanger, exchange_grid, ocean::MITgcmOceanSimulation, coupled_model) + u = exchanger.state.u + v = exchanger.state.v + T = exchanger.state.T + S = exchanger.state.S + + u_surf, v_surf = ocean_surface_velocities(ocean) + T_surf = ocean_surface_temperature(ocean) + S_surf = ocean_surface_salinity(ocean) + + Nx = size(ocean.xc, 1) + Ny = size(ocean.xc, 2) + + # u at (Face, Center): with periodic longitude, interior is (Nx, Ny) — direct match + u_interior = interior(u, :, :, 1) + u_interior .= u_surf + + # v at (Center, Face): with bounded latitude, interior is (Nx, Ny+1). + # MITgcm vVel(i,j) is at the southern face of cell (i,j), so maps to face indices 1:Ny. + # The Ny+1-th face (northern boundary) is set to zero. + v_interior = interior(v, :, :, 1) + v_interior[:, 1:Ny] .= v_surf + v_interior[:, Ny+1] .= 0 + + # T, S at (Center, Center): interior is (Nx, Ny) — direct match + interior(T, :, :, 1) .= T_surf + interior(S, :, :, 1) .= S_surf + + return nothing +end + +initialize!(exchanger::ComponentExchanger, grid, ::MITgcmOceanSimulation) = nothing + +get_radiative_forcing(ocean::MITgcmOceanSimulation) = nothing + +# MITgcm handles its own freezing/temperature limiting internally, +# so no-op the FreezingLimitedOceanTemperature sea-ice-ocean flux computation +using NumericalEarth.EarthSystemModels: EarthSystemModel, NoSeaIceInterface +using NumericalEarth.SeaIces: FreezingLimitedOceanTemperature + +const MITgcmFreezingLimited = EarthSystemModel{<:FreezingLimitedOceanTemperature, <:Any, <:MITgcmOceanSimulation, <:NoSeaIceInterface} +compute_sea_ice_ocean_fluxes!(::MITgcmFreezingLimited) = nothing + +# ============================================================ +# Flux update: coupled model → MITgcm surface forcing +# ============================================================ + +function update_net_fluxes!(coupled_model, ocean::MITgcmOceanSimulation) + + # Use the generic ocean flux assembler to compute net fluxes + # on the exchange grid (momentum, heat, freshwater) + exchange_grid = coupled_model.interfaces.exchanger.grid + Oceans.update_net_ocean_fluxes!(coupled_model, ocean, exchange_grid) + + net_ocean_fluxes = coupled_model.interfaces.net_fluxes.ocean + + Nx = size(ocean.xc, 1) + Ny = size(ocean.xc, 2) + + ρ₀ = ocean.reference_density + cₚ = ocean.heat_capacity + + # Extract assembled fluxes. + # NumericalEarth computes: + # τx, τy = kinematic stress (m²/s² = N/m² / ρ₀) + # Jᵀ = temperature tendency (K·m/s = W/m² / (ρ₀·cₚ)) + # Jˢ = salinity tendency (PSU·m/s) + # MITgcm expects: + # fu, fv = wind stress (N/m²) → multiply by ρ₀ + # Qnet = net heat flux (W/m², positive = ocean cooling) → multiply by ρ₀·cₚ + # EmPmR = freshwater flux (kg/m²/s) → derive from Jˢ (or set to 0 if not available) + # Qsw = shortwave (W/m²) → set to 0 (already included in Qnet) + # saltFlux = salt flux (g/m²/s) → derive from Jˢ + + τx = interior(net_ocean_fluxes.u, :, :, 1) + τy = interior(net_ocean_fluxes.v, :, :, 1) + JT = interior(net_ocean_fluxes.T, :, :, 1) + JS = interior(net_ocean_fluxes.S, :, :, 1) + + for j in 1:Ny, i in 1:Nx + ocean.fu[i, j] = - τx[i, j] * ρ₀ + ocean.fv[i, j] = - τy[i, j] * ρ₀ + ocean.qnet[i, j] = JT[i, j] * ρ₀ * cₚ + ocean.saltflux[i, j] = JS[i, j] * ρ₀ + end + + # These are not needed + fill!(ocean.empmr, 0) + fill!(ocean.qsw, 0) + + # Push to MITgcm + set_fu!(ocean.library, ocean.fu) + set_fv!(ocean.library, ocean.fv) + set_qnet!(ocean.library, ocean.qnet) + set_empmr!(ocean.library, ocean.empmr) + set_qsw!(ocean.library, ocean.qsw) + set_saltflux!(ocean.library, ocean.saltflux) + + return nothing +end diff --git a/src/Oceans/ocean_simulation.jl b/src/Oceans/ocean_simulation.jl index b2876b97c..9a40abbb6 100644 --- a/src/Oceans/ocean_simulation.jl +++ b/src/Oceans/ocean_simulation.jl @@ -123,6 +123,7 @@ end equation_of_state = TEOS10EquationOfState(; reference_density), boundary_conditions::NamedTuple = NamedTuple(), radiative_forcing = default_radiative_forcing(grid), + materialize_buoyancy_gradients = true, warn = true, verbose = false) @@ -206,6 +207,7 @@ function ocean_simulation(grid; tracer_advection = WENO(order=7), equation_of_state = TEOS10EquationOfState(; reference_density), boundary_conditions::NamedTuple = NamedTuple(), + materialize_buoyancy_gradients = true, radiative_forcing = default_radiative_forcing(grid), warn = true, verbose = false) @@ -302,6 +304,7 @@ function ocean_simulation(grid; # conditions even when a user-bc is supplied). boundary_conditions = merge(default_boundary_conditions, boundary_conditions) buoyancy = SeawaterBuoyancy(; gravitational_acceleration, equation_of_state) + bouyancy = Oceananigans.BuoyancyFormulations.BuoyancyForce(grid, buoyancy; materialize_gradients = materialize_buoyancy_gradients) if tracer_advection isa NamedTuple tracer_advection = with_tracers(tracers, tracer_advection, default_tracer_advection()) diff --git a/test/test_mitgcm.jl b/test/test_mitgcm.jl new file mode 100644 index 000000000..9aab27d1a --- /dev/null +++ b/test/test_mitgcm.jl @@ -0,0 +1,97 @@ +include("runtests_setup.jl") + +using Test +using NumericalEarth +using Oceananigans +using MITgcm + +# Path to pre-built library from MITgcm.jl examples +const MITGCM_JL_DIR = pkgdir(MITgcm) +const EXAMPLE_DIR = joinpath(MITGCM_JL_DIR, "examples", "global_oce_latlon") +const LIB_PATH = joinpath(EXAMPLE_DIR, Sys.isapple() ? "libmitgcm.dylib" : "libmitgcm.so") +const RUN_DIR = joinpath(EXAMPLE_DIR, "run") + +if !isfile(LIB_PATH) || !isdir(RUN_DIR) + @warn "Pre-built MITgcm library not found at $EXAMPLE_DIR. Skipping MITgcm tests." +else + +@testset "MITgcm ocean model interface" begin + # Test that the extension is available + MITgcmExt = Base.get_extension(NumericalEarth, :NumericalEarthMITgcmExt) + @test !isnothing(MITgcmExt) + + # Test creating a MITgcmOceanSimulation from pre-built library + ocean = MITgcmOceanSimulation(LIB_PATH, RUN_DIR) + @test ocean isa MITgcmOceanSimulation + + # Test dimensions + @test size(ocean.theta) == (90, 40, 15) + @test size(ocean.salt) == (90, 40, 15) + @test size(ocean.uvel) == (90, 40, 15) + @test size(ocean.vvel) == (90, 40, 15) + @test size(ocean.etan) == (90, 40) + @test size(ocean.xc) == (90, 40) + @test size(ocean.yc) == (90, 40) + @test length(ocean.rc) == 15 + @test length(ocean.drf) == 15 + + # Test basic interface functions + ρₒ = NumericalEarth.EarthSystemModels.reference_density(ocean) + cₚ = NumericalEarth.EarthSystemModels.heat_capacity(ocean) + @test ρₒ ≈ 1029.0 + @test cₚ ≈ 3994.0 + + # Test ocean state accessors + T = NumericalEarth.EarthSystemModels.ocean_temperature(ocean) + S = NumericalEarth.EarthSystemModels.ocean_salinity(ocean) + @test T === ocean.theta + @test S === ocean.salt + @test ndims(T) == 3 + @test ndims(S) == 3 + + T_surf = NumericalEarth.EarthSystemModels.ocean_surface_temperature(ocean) + S_surf = NumericalEarth.EarthSystemModels.ocean_surface_salinity(ocean) + @test size(T_surf) == (90, 40) + @test size(S_surf) == (90, 40) + + u_surf, v_surf = NumericalEarth.EarthSystemModels.ocean_surface_velocities(ocean) + @test size(u_surf) == (90, 40) + @test size(v_surf) == (90, 40) + + # Test surface grid construction + grid = MITgcmExt.surface_grid(ocean) + @test grid isa Oceananigans.Grids.LatitudeLongitudeGrid + + # Test time stepping + theta_before = copy(ocean.theta) + MITgcmExt.time_step!(ocean, 43200.0) # 12 hours + @test ocean.theta != theta_before # State should change + + # Test a few more steps + for _ in 1:3 + MITgcmExt.time_step!(ocean, 43200.0) + end + @test MITgcm.get_niter(ocean.library) == 4 + + # Test with PrescribedAtmosphere and EarthSystemModel + @testset "Coupled model with PrescribedAtmosphere" begin + grid = MITgcmExt.surface_grid(ocean) + Nx, Ny = size(ocean.xc) + + atmosphere = NumericalEarth.PrescribedAtmosphere(grid) + + coupled_model = NumericalEarth.EarthSystemModels.OceanOnlyModel(ocean; atmosphere) + @test coupled_model isa NumericalEarth.EarthSystemModels.EarthSystemModel + + # Initialize and step the coupled model + NumericalEarth.EarthSystemModels.initialize!(coupled_model) + NumericalEarth.EarthSystemModels.time_step!(coupled_model, 43200.0) + + @test MITgcm.get_niter(ocean.library) == 5 + end + + # Finalize + MITgcm.finalize!(ocean.library) +end + +end # if pre-built library exists