|
| 1 | +import math |
| 2 | + |
| 3 | +import numpy as np |
| 4 | +import pandas as pd |
| 5 | +import uxarray as ux |
| 6 | + |
| 7 | +__all__ = ["Nx", "datasets"] |
| 8 | + |
| 9 | +Nx = 20 |
| 10 | +vmax = 1.0 |
| 11 | +delta = 0.1 |
| 12 | + |
| 13 | + |
| 14 | +def _stommel_gyre_delaunay(): |
| 15 | + lon, lat = np.meshgrid(np.linspace(0, 60.0, Nx, dtype=np.float32), np.linspace(0, 60.0, Nx, dtype=np.float32)) |
| 16 | + lon_flat = lon.ravel() |
| 17 | + lat_flat = lat.ravel() |
| 18 | + |
| 19 | + # mask any point on one of the boundaries |
| 20 | + mask = ( |
| 21 | + np.isclose(lon_flat, 0.0) | np.isclose(lon_flat, 60.0) | np.isclose(lat_flat, 0.0) | np.isclose(lat_flat, 60.0) |
| 22 | + ) |
| 23 | + |
| 24 | + boundary_points = np.flatnonzero(mask) |
| 25 | + |
| 26 | + uxgrid = ux.Grid.from_points( |
| 27 | + (lon_flat, lat_flat), |
| 28 | + method="regional_delaunay", |
| 29 | + boundary_points=boundary_points, |
| 30 | + ) |
| 31 | + |
| 32 | + # Define arrays U (zonal), V (meridional) and P (sea surface height) |
| 33 | + U = np.zeros((1, 1, lat.size), dtype=np.float64) |
| 34 | + V = np.zeros((1, 1, lat.size), dtype=np.float64) |
| 35 | + P = np.zeros((1, 1, lat.size), dtype=np.float64) |
| 36 | + |
| 37 | + for i, (x, y) in enumerate(zip(lon_flat, lat_flat, strict=False)): |
| 38 | + xi = x / 60.0 |
| 39 | + yi = y / 60.0 |
| 40 | + |
| 41 | + P[0, 0, i] = -vmax * delta * (1 - xi) * (math.exp(-xi / delta) - 1) * np.sin(math.pi * yi) |
| 42 | + U[0, 0, i] = -vmax * (1 - math.exp(-xi / delta) - xi) * np.cos(math.pi * yi) |
| 43 | + V[0, 0, i] = vmax * ((2.0 - xi) * math.exp(-xi / delta) - 1) * np.sin(math.pi * yi) |
| 44 | + |
| 45 | + u = ux.UxDataArray( |
| 46 | + data=U, |
| 47 | + name="U", |
| 48 | + uxgrid=uxgrid, |
| 49 | + dims=["time", "nz1", "n_node"], |
| 50 | + coords=dict( |
| 51 | + time=(["time"], pd.to_datetime(["2000-01-01"])), |
| 52 | + nz1=(["nz1"], [0]), |
| 53 | + ), |
| 54 | + attrs=dict( |
| 55 | + description="zonal velocity", units="m/s", location="node", mesh="delaunay", Conventions="UGRID-1.0" |
| 56 | + ), |
| 57 | + ) |
| 58 | + v = ux.UxDataArray( |
| 59 | + data=V, |
| 60 | + name="V", |
| 61 | + uxgrid=uxgrid, |
| 62 | + dims=["time", "nz1", "n_node"], |
| 63 | + coords=dict( |
| 64 | + time=(["time"], pd.to_datetime(["2000-01-01"])), |
| 65 | + nz1=(["nz1"], [0]), |
| 66 | + ), |
| 67 | + attrs=dict( |
| 68 | + description="meridional velocity", units="m/s", location="node", mesh="delaunay", Conventions="UGRID-1.0" |
| 69 | + ), |
| 70 | + ) |
| 71 | + p = ux.UxDataArray( |
| 72 | + data=P, |
| 73 | + name="p", |
| 74 | + uxgrid=uxgrid, |
| 75 | + dims=["time", "nz1", "n_node"], |
| 76 | + coords=dict( |
| 77 | + time=(["time"], pd.to_datetime(["2000-01-01"])), |
| 78 | + nz1=(["nz1"], [0]), |
| 79 | + ), |
| 80 | + attrs=dict(description="pressure", units="N/m^2", location="node", mesh="delaunay", Conventions="UGRID-1.0"), |
| 81 | + ) |
| 82 | + |
| 83 | + return ux.UxDataset({"U": u, "V": v, "p": p}, uxgrid=uxgrid) |
| 84 | + |
| 85 | + |
| 86 | +datasets = { |
| 87 | + "stommel_gyre_delaunay": _stommel_gyre_delaunay(), |
| 88 | +} |
0 commit comments