Force the MITgcm using JRA55

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"

examples/mitgcm/.gitignore

@@ -0,0 +1,2 @@
+# MITgcm build artifacts
+build/

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 )

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

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

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',
+ &

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
+ &

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,
+ &

examples/mitgcm/mitgcm_config/input/data.kpp

@@ -0,0 +1,3 @@
+# KPP vertical mixing parameters
+ &KPP_PARM01
+ &

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.,
+ &

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!"

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