Skip to content

Commit 6538f49

Browse files
erikvansebilleerikvansebille
authored
Merge pull request #2481 from Parcels-code/convert_croco
Create convert.croco_to_sgrid function
2 parents 8d2c14e + f48aa8f commit 6538f49

16 files changed

Lines changed: 640 additions & 506 deletions

‎docs/user_guide/examples/tutorial_croco_3D.ipynb‎

Lines changed: 324 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,324 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# 🖥️ CROCO tutorial"
8+
]
9+
},
10+
{
11+
"cell_type": "markdown",
12+
"metadata": {},
13+
"source": [
14+
"This tutorial will show how to run a 3D simulation with output from the CROCO model."
15+
]
16+
},
17+
{
18+
"cell_type": "markdown",
19+
"metadata": {},
20+
"source": [
21+
"## Example setup"
22+
]
23+
},
24+
{
25+
"cell_type": "markdown",
26+
"metadata": {},
27+
"source": [
28+
"We start with loading the relevant modules and the data."
29+
]
30+
},
31+
{
32+
"cell_type": "code",
33+
"execution_count": null,
34+
"metadata": {},
35+
"outputs": [],
36+
"source": [
37+
"import matplotlib.pyplot as plt\n",
38+
"import numpy as np\n",
39+
"import xarray as xr\n",
40+
"\n",
41+
"import parcels\n",
42+
"\n",
43+
"data_folder = parcels.download_example_dataset(\"CROCOidealized_data\")\n",
44+
"ds_fields = xr.open_dataset(data_folder / \"CROCO_idealized.nc\")\n",
45+
"\n",
46+
"ds_fields.load(); # Preload data to speed up access"
47+
]
48+
},
49+
{
50+
"cell_type": "markdown",
51+
"metadata": {},
52+
"source": [
53+
"The simulation will be a very simple, idealised flow: a purely zonal flow over a sloping bottom. This flow (which is somewhat unrealistic of course) nicely showcases that particles stay on their initial depth levels, even though the sigma-layers slope down. \n",
54+
"\n",
55+
"This flow has been created by first running the example from the [Shelf front example on the CROCO website](https://croco-ocean.gitlabpages.inria.fr/croco_doc/model/model.test_cases.shelfront.html). Then, we took the restart file are manually replaced all the `u`-velocities with `1` m/s and all the `v`-velocities with `0` m/s. This way we get a purely zonal flow. We then started a new simulation from the restart file, and CROCO then automatically calculated the `w` velocities to match the new zonal field. We saved the `time=0` snapshot from this new run and use it below."
56+
]
57+
},
58+
{
59+
"cell_type": "markdown",
60+
"metadata": {},
61+
"source": [
62+
"Now we create a FieldSet object using the `FieldSet.from_croco()` method. Note that CROCO is a C-grid (with similar indexing at MITgcm)."
63+
]
64+
},
65+
{
66+
"cell_type": "markdown",
67+
"metadata": {},
68+
"source": [
69+
"```{note}\n",
70+
"In the code below, we use the `w` velocity field for vertical velocity. However, it is unclear whether this is always the right choice. CROCO (and ROMS) also output an `omega` field, which may be more appropriate to use. The idealised simulation below only works when using `w`, though. In other simulations, it is recommended to test whether `omega` provides more realistic results. See https://github.com/OceanParcels/Parcels/discussions/1728 for more information.\n",
71+
"```"
72+
]
73+
},
74+
{
75+
"cell_type": "code",
76+
"execution_count": null,
77+
"metadata": {},
78+
"outputs": [],
79+
"source": [
80+
"fields = {\n",
81+
" \"U\": ds_fields[\"u\"],\n",
82+
" \"V\": ds_fields[\"v\"],\n",
83+
" \"W\": ds_fields[\"w\"],\n",
84+
" \"h\": ds_fields[\"h\"],\n",
85+
" \"zeta\": ds_fields[\"zeta\"],\n",
86+
" \"Cs_w\": ds_fields[\"Cs_w\"],\n",
87+
"}\n",
88+
"\n",
89+
"coords = ds_fields[[\"x_rho\", \"y_rho\", \"s_w\", \"time\"]]\n",
90+
"ds_fset = parcels.convert.croco_to_sgrid(fields=fields, coords=coords)\n",
91+
"\n",
92+
"fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)\n",
93+
"\n",
94+
"# Add the critical depth (`hc`) as a constant to the fieldset\n",
95+
"fieldset.add_constant(\"hc\", ds_fields.hc.item())"
96+
]
97+
},
98+
{
99+
"cell_type": "markdown",
100+
"metadata": {},
101+
"source": [
102+
"Now we can use this Fieldset to advect particles as we would normally do. Note that the particle depths should be provided in (negative) meters, not in sigma-coordinates."
103+
]
104+
},
105+
{
106+
"cell_type": "code",
107+
"execution_count": null,
108+
"metadata": {
109+
"tags": [
110+
"hide-output"
111+
]
112+
},
113+
"outputs": [],
114+
"source": [
115+
"X, Z = np.meshgrid(\n",
116+
" [40e3, 80e3, 120e3],\n",
117+
" [100, -10, -130, -250, -400, -850, -1400, -1550],\n",
118+
")\n",
119+
"Y = np.ones(X.size) * 98000\n",
120+
"\n",
121+
"\n",
122+
"def DeleteParticle(particles, fieldset):\n",
123+
" any_error = particles.state >= 50\n",
124+
" particles[any_error].state = parcels.StatusCode.Delete\n",
125+
"\n",
126+
"\n",
127+
"pset = parcels.ParticleSet(\n",
128+
" fieldset=fieldset, pclass=parcels.Particle, lon=X, lat=Y, z=Z\n",
129+
")\n",
130+
"\n",
131+
"outputfile = parcels.ParticleFile(\n",
132+
" store=\"croco_particles3D.zarr\",\n",
133+
" outputdt=np.timedelta64(5000, \"s\"),\n",
134+
")\n",
135+
"\n",
136+
"pset.execute(\n",
137+
" [parcels.kernels.AdvectionRK4_3D_CROCO, DeleteParticle],\n",
138+
" runtime=np.timedelta64(50_000, \"s\"),\n",
139+
" dt=np.timedelta64(100, \"s\"),\n",
140+
" output_file=outputfile,\n",
141+
")"
142+
]
143+
},
144+
{
145+
"cell_type": "markdown",
146+
"metadata": {},
147+
"source": [
148+
"Now we plot the particle trajectories below. Note that the particles stay on their initial depth levels, even though the sigma-layers slope down. Also note that particles released above the surface (where depth >0) or below the bathymetry are not advected (due to the `DeleteParticle` kernel)."
149+
]
150+
},
151+
{
152+
"cell_type": "code",
153+
"execution_count": null,
154+
"metadata": {
155+
"tags": [
156+
"nbsphinx-thumbnail"
157+
]
158+
},
159+
"outputs": [],
160+
"source": [
161+
"fig, ax = plt.subplots(1, 1, figsize=(6, 4))\n",
162+
"ds = xr.open_zarr(\"croco_particles3D.zarr\")\n",
163+
"\n",
164+
"ax.plot(X / 1e3, Z, \"k.\", label=\"Initial positions\")\n",
165+
"ax.plot(ds.lon.T / 1e3, ds.z.T, \".-\")\n",
166+
"\n",
167+
"for z in ds_fields.s_w.values:\n",
168+
" ax.plot(\n",
169+
" fieldset.h.grid.lon[0, :] / 1e3,\n",
170+
" fieldset.h.data[0, 0, 25, :] * z,\n",
171+
" \"k\",\n",
172+
" linewidth=0.5,\n",
173+
" )\n",
174+
"ax.set_xlabel(\"X [km]\")\n",
175+
"ax.set_xlim(30, 170)\n",
176+
"ax.set_ylabel(\"Depth [m]\")\n",
177+
"ax.set_title(\"Particles in idealized CROCO velocity field using AdvectionRK4_3D_CROCO\")\n",
178+
"plt.tight_layout()\n",
179+
"plt.show()"
180+
]
181+
},
182+
{
183+
"cell_type": "markdown",
184+
"metadata": {},
185+
"source": [
186+
"### A CROCO simulation with no vertical velocities"
187+
]
188+
},
189+
{
190+
"cell_type": "markdown",
191+
"metadata": {},
192+
"source": [
193+
"It may be insightful to compare this 3D run with the `AdvectionRK4_3D_CROCO` kernel with a run where the vertical velocity (`W`) is set to zero. In that case, the particles will not stay on their initial depth levels but instead follow sigma-layers."
194+
]
195+
},
196+
{
197+
"cell_type": "code",
198+
"execution_count": null,
199+
"metadata": {
200+
"tags": [
201+
"hide-output"
202+
]
203+
},
204+
"outputs": [],
205+
"source": [
206+
"fieldset_noW = parcels.FieldSet.from_sgrid_conventions(ds_fset)\n",
207+
"fieldset_noW.W.data[:] = 0.0\n",
208+
"fieldset_noW.add_constant(\"hc\", ds_fields.hc.item())\n",
209+
"\n",
210+
"X, Z = np.meshgrid(\n",
211+
" [40e3, 80e3, 120e3],\n",
212+
" [100, -10, -130, -250, -400, -850, -1400, -1550],\n",
213+
")\n",
214+
"Y = np.ones(X.size) * 100\n",
215+
"\n",
216+
"pset_noW = parcels.ParticleSet(\n",
217+
" fieldset=fieldset_noW, pclass=parcels.Particle, lon=X, lat=Y, z=Z\n",
218+
")\n",
219+
"\n",
220+
"outputfile = parcels.ParticleFile(\n",
221+
" store=\"croco_particles_noW.zarr\", outputdt=np.timedelta64(5000, \"s\")\n",
222+
")\n",
223+
"\n",
224+
"pset_noW.execute(\n",
225+
" [parcels.kernels.AdvectionRK4_3D_CROCO, DeleteParticle],\n",
226+
" runtime=np.timedelta64(50_000, \"s\"),\n",
227+
" dt=np.timedelta64(100, \"s\"),\n",
228+
" output_file=outputfile,\n",
229+
")\n",
230+
"\n",
231+
"fig, ax = plt.subplots(1, 1, figsize=(6, 4))\n",
232+
"ds = xr.open_zarr(\"croco_particles_noW.zarr\")\n",
233+
"\n",
234+
"ax.plot(X / 1e3, Z, \"k.\", label=\"Initial positions\")\n",
235+
"ax.plot(ds.lon.T / 1e3, ds.z.T, \".-\")\n",
236+
"\n",
237+
"for z in ds_fields.s_w.values:\n",
238+
" ax.plot(\n",
239+
" fieldset.h.grid.lon[0, :] / 1e3,\n",
240+
" fieldset.h.data[0, 0, 25, :] * z,\n",
241+
" \"k\",\n",
242+
" linewidth=0.5,\n",
243+
" )\n",
244+
"ax.set_xlabel(\"X [km]\")\n",
245+
"ax.set_xlim(30, 170)\n",
246+
"ax.set_ylabel(\"Depth [m]\")\n",
247+
"ax.set_title(\"Particles in idealized CROCO velocity field with W=0\")\n",
248+
"plt.tight_layout()\n",
249+
"plt.show()"
250+
]
251+
},
252+
{
253+
"cell_type": "markdown",
254+
"metadata": {},
255+
"source": [
256+
"## The algorithms used"
257+
]
258+
},
259+
{
260+
"cell_type": "markdown",
261+
"metadata": {},
262+
"source": [
263+
"When using Croco data, Parcels needs to convert the particles.depth to sigma-coordinates, before doing any interpolation. This is done with the `convert_z_to_sigma_croco()` function. Interpolating onto a Field is then done like:\n",
264+
"\n",
265+
"```python\n",
266+
"def SampleTempCroco(particles, fieldset):\n",
267+
" \"\"\"Sample temperature field on a CROCO sigma grid by first converting z to sigma levels.\"\"\"\n",
268+
" sigma = convert_z_to_sigma_croco(fieldset, particles.time, particles.z, particles.lat, particles.lon, particles)\n",
269+
" particles.temp = fieldset.T[particles.time, sigma, particles.lat, particles.lon, particles]\n",
270+
"```"
271+
]
272+
},
273+
{
274+
"cell_type": "markdown",
275+
"metadata": {},
276+
"source": [
277+
"For Advection, you will need to use the `AdvectionRK4_3D_CROCO` kernel, which works slightly different from the normal 3D advection kernel because it converts the vertical velocity in sigma-units. The conversion from depth to sigma is done at every time step, using the `convert_z_to_sigma_croco()` function.\n",
278+
"\n",
279+
"In particular, the following algorithm is used (note that the RK4 version is slightly more complex than this Euler-Forward version, but the idea is identical). Also note that the vertical velocity is linearly interpolated here, which gives much better results than the default C-grid interpolation.\n",
280+
"\n",
281+
"```python\n",
282+
"sigma = particles.z / fieldset.h[particles.time, 0, particles.lat, particles.lon]\n",
283+
"\n",
284+
"sigma_levels = convert_z_to_sigma_croco(fieldset, particles.time, particles.z, particles.lat, particles.lon, particles)\n",
285+
"(u, v) = fieldset.UV[time, sigma_levels, particle.lat, particle.lon, particle]\n",
286+
"w = fieldset.W[time, sigma_levels, particle.lat, particle.lon, particle]\n",
287+
"\n",
288+
"# scaling the w with the sigma level of the particle\n",
289+
"w_sigma = w * sigma / fieldset.h[particles.time, 0, particles.lat, particles.lon]\n",
290+
"\n",
291+
"lon_new = particles.lon + u*particles.dt\n",
292+
"lat_new = particles.lat + v*particles.dt\n",
293+
"\n",
294+
"# calculating new sigma level\n",
295+
"sigma_new = sigma + w_sigma*particles.dt \n",
296+
"\n",
297+
"# converting back from sigma to depth, at _new_ location\n",
298+
"depth_new = sigma_new * fieldset.h[particles.time, 0, lat_new, lon_new]\n",
299+
"```"
300+
]
301+
}
302+
],
303+
"metadata": {
304+
"kernelspec": {
305+
"display_name": "docs",
306+
"language": "python",
307+
"name": "python3"
308+
},
309+
"language_info": {
310+
"codemirror_mode": {
311+
"name": "ipython",
312+
"version": 3
313+
},
314+
"file_extension": ".py",
315+
"mimetype": "text/x-python",
316+
"name": "python",
317+
"nbconvert_exporter": "python",
318+
"pygments_lexer": "ipython3",
319+
"version": "3.14.2"
320+
}
321+
},
322+
"nbformat": 4,
323+
"nbformat_minor": 2
324+
}

0 commit comments

Comments
 (0)