Skip to content

Commit d81e816

Browse files
gaelforgetgaelforget
authored
Merge pull request #171 from simone-silvestri/ss-gf-jmc/mitgcm-running-in-julia
Start interfacing MITgcm with julia via a dynamic library
2 parents 7581352 + 49c71da commit d81e816

11 files changed

Lines changed: 2616 additions & 0 deletions

‎Project.toml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
1111
Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
1212
FortranFiles = "c58ffaec-ab22-586d-bfc5-781a99fd0b10"
1313
Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
14+
Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
1415
MeshArrays = "cb8c808f-1acf-59a3-9d2b-6e38d009f683"
1516
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
1617
Scratch = "6c6a2e73-6563-6170-7368-637461726353"

‎examples/global_oce_latlon/run_global_oce_latlon.jl‎

Lines changed: 173 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,173 @@
1+
#=
2+
run_global_oce_latlon.jl
3+
4+
Drive the MITgcm global_oce_latlon experiment from Julia using MITgcm.jl.
5+
6+
This script:
7+
1. Loads the MITgcm shared library via MITgcmOceanSimulation
8+
2. Steps forward in time with prescribed atmospheric forcing
9+
3. Prints diagnostics at each step
10+
11+
Usage:
12+
Option A: Build from scratch (requires MITgcm source):
13+
julia -e 'using MITgcm; ocean = MITgcmOceanSimulation("/path/to/MITgcm")'
14+
15+
Option B: Use pre-built library:
16+
cd examples/global_oce_latlon
17+
julia run_global_oce_latlon.jl
18+
=#
19+
20+
using Printf
21+
using MITgcm
22+
23+
# ============================================================
24+
# Setup: use pre-built library or build from scratch
25+
# ============================================================
26+
27+
const SCRIPT_DIR = @__DIR__
28+
29+
# Check for pre-built library
30+
function find_prebuilt_library()
31+
for ext in (".dylib", ".so")
32+
path = joinpath(SCRIPT_DIR, "libmitgcm" * ext)
33+
isfile(path) && return path
34+
end
35+
return nothing
36+
end
37+
38+
prebuilt = find_prebuilt_library()
39+
run_dir = joinpath(SCRIPT_DIR, "run")
40+
41+
if prebuilt !== nothing && isdir(run_dir)
42+
println("Using pre-built library: $prebuilt")
43+
ocean = MITgcmOceanSimulation(prebuilt, run_dir)
44+
else
45+
mitgcm_dir = get(ENV, "MITGCM_DIR", joinpath(SCRIPT_DIR, "..", "..", "..", "MITgcm"))
46+
if !isdir(mitgcm_dir)
47+
error("MITgcm directory not found: $mitgcm_dir\n" *
48+
"Set MITGCM_DIR environment variable or run build_mitgcm_lib.sh first.")
49+
end
50+
println("Building MITgcm from: $mitgcm_dir")
51+
ocean = MITgcmOceanSimulation(mitgcm_dir; output_dir=SCRIPT_DIR)
52+
end
53+
54+
# ============================================================
55+
# Grid info
56+
# ============================================================
57+
58+
lib = ocean.library
59+
dims = lib.dims
60+
Nx, Ny, Nr = dims.Nx, dims.Ny, dims.Nr
61+
62+
println()
63+
println("=" ^ 60)
64+
println(" MITgcm global_oce_latlon - Julia Interface Demo")
65+
println("=" ^ 60)
66+
println()
67+
println("Grid dimensions: Nx=$Nx Ny=$Ny Nr=$Nr")
68+
@printf("Longitude range: %.1f to %.1f deg\n", minimum(ocean.xc), maximum(ocean.xc))
69+
@printf("Latitude range: %.1f to %.1f deg\n", minimum(ocean.yc), maximum(ocean.yc))
70+
println()
71+
72+
# Ocean mask from hFacC surface layer
73+
ocean_mask_2d = ocean.hfacc[:, :, 1]
74+
75+
# ============================================================
76+
# Prescribed atmospheric forcing (simple bulk formulas)
77+
# ============================================================
78+
79+
const ρ_air = 1.225 # air density (kg/m³)
80+
const c_p_air = 1004.0 # specific heat of air (J/kg/K)
81+
const C_d = 1.2e-3 # drag coefficient
82+
const C_h = 1.0e-3 # heat transfer coefficient
83+
84+
u_wind = zeros(Float64, Nx, Ny)
85+
T_air = zeros(Float64, Nx, Ny)
86+
87+
for j in 1:Ny, i in 1:Nx
88+
lat_rad = ocean.yc[i, j] * π / 180.0
89+
u_wind[i, j] = -10.0 * cos(3.0 * lat_rad)
90+
T_air[i, j] = 25.0 - 40.0 * sin(lat_rad)^2
91+
end
92+
93+
println("Prescribed atmosphere:")
94+
@printf(" u_wind range: %.2f to %.2f m/s\n", minimum(u_wind), maximum(u_wind))
95+
@printf(" T_air range: %.2f to %.2f °C\n", minimum(T_air), maximum(T_air))
96+
println()
97+
98+
function compute_bulk_fluxes!(fu, fv, qnet, sst, u_wind, T_air, mask)
99+
Nx, Ny = size(fu)
100+
for j in 1:Ny, i in 1:Nx
101+
if mask[i, j] > 0
102+
wind_speed = abs(u_wind[i, j])
103+
fu[i, j] = ρ_air * C_d * wind_speed * u_wind[i, j]
104+
fv[i, j] = 0.0
105+
qnet[i, j] = ρ_air * c_p_air * C_h * wind_speed * (sst[i, j] - T_air[i, j])
106+
else
107+
fu[i, j] = 0.0
108+
fv[i, j] = 0.0
109+
qnet[i, j] = 0.0
110+
end
111+
end
112+
end
113+
114+
# ============================================================
115+
# Helper: print field statistics
116+
# ============================================================
117+
118+
function field_stats(name::String, arr; mask=nothing)
119+
valid = mask !== nothing ? arr[mask .> 0] : arr[:]
120+
isempty(valid) && return @printf(" %-8s: no valid points\n", name)
121+
@printf(" %-8s: min=%12.5e max=%12.5e mean=%12.5e\n",
122+
name, minimum(valid), maximum(valid), sum(valid) / length(valid))
123+
end
124+
125+
# ============================================================
126+
# Time stepping loop
127+
# ============================================================
128+
129+
println("Initial state:")
130+
field_stats("theta", ocean.theta[:,:,1], mask=ocean_mask_2d)
131+
field_stats("salt", ocean.salt[:,:,1], mask=ocean_mask_2d)
132+
println()
133+
134+
nsteps = 20
135+
println("Running $nsteps time steps...")
136+
println("-" ^ 60)
137+
138+
for s in 1:nsteps
139+
# Compute bulk fluxes from current SST
140+
sst = @view ocean.theta[:, :, 1]
141+
compute_bulk_fluxes!(ocean.fu, ocean.fv, ocean.qnet, sst, u_wind, T_air, ocean_mask_2d)
142+
143+
# Set forcing in MITgcm
144+
set_fu!(lib, ocean.fu)
145+
set_fv!(lib, ocean.fv)
146+
set_qnet!(lib, ocean.qnet)
147+
set_empmr!(lib, ocean.empmr)
148+
set_saltflux!(lib, ocean.saltflux)
149+
150+
# Step forward
151+
step!(lib)
152+
refresh_state!(ocean)
153+
154+
niter = get_niter(lib)
155+
mtime = get_time(lib)
156+
157+
@printf("\nStep %3d | iter=%5d | time=%10.1f s (%.2f days)\n",
158+
s, niter, mtime, mtime / 86400.0)
159+
field_stats("theta", ocean.theta[:,:,1], mask=ocean_mask_2d)
160+
field_stats("salt", ocean.salt[:,:,1], mask=ocean_mask_2d)
161+
field_stats("uVel", ocean.uvel[:,:,1], mask=ocean_mask_2d)
162+
field_stats("etaN", ocean.etan, mask=ocean_mask_2d)
163+
end
164+
165+
println()
166+
println("-" ^ 60)
167+
println("Timestepping complete.")
168+
println()
169+
170+
# Finalize
171+
println("Finalizing MITgcm...")
172+
MITgcm.finalize!(lib)
173+
println("Done!")

0 commit comments

Comments
 (0)
⚡