diff --git a/Makefile b/Makefile index 816f198..c87fcdc 100644 --- a/Makefile +++ b/Makefile @@ -183,6 +183,7 @@ COMMON_SRCS = $(SRCDIR)/kdtree.c \ $(SRCDIR)/mesh.c \ $(SRCDIR)/regrid.c \ $(SRCDIR)/file_netcdf.c \ + $(SRCDIR)/file_mitgcm.c \ $(SRCDIR)/colormaps.c \ $(SRCDIR)/view.c @@ -274,6 +275,7 @@ $(OBJDIR)/mesh.o: $(SRCDIR)/mesh.c $(SRCDIR)/mesh.h $(SRCDIR)/ushow.defines.h $(OBJDIR)/regrid.o: $(SRCDIR)/regrid.c $(SRCDIR)/regrid.h $(SRCDIR)/mesh.h \ $(SRCDIR)/kdtree.h $(SRCDIR)/ushow.defines.h $(OBJDIR)/file_netcdf.o: $(SRCDIR)/file_netcdf.c $(SRCDIR)/file_netcdf.h $(SRCDIR)/ushow.defines.h +$(OBJDIR)/file_mitgcm.o: $(SRCDIR)/file_mitgcm.c $(SRCDIR)/file_mitgcm.h $(SRCDIR)/mesh.h $(SRCDIR)/ushow.defines.h $(OBJDIR)/colormaps.o: $(SRCDIR)/colormaps.c $(SRCDIR)/colormaps.h $(SRCDIR)/ushow.defines.h $(OBJDIR)/view.o: $(SRCDIR)/view.c $(SRCDIR)/view.h $(SRCDIR)/file_netcdf.h \ $(SRCDIR)/regrid.h $(SRCDIR)/colormaps.h $(SRCDIR)/ushow.defines.h diff --git a/README.md b/README.md index 5dfae9a..9d68e8d 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ A fast, ncview‑inspired visualization tool for structured and unstructured geo ## Features -- **Multiple formats**: Supports netCDF, zarr and GRIB +- **Multiple formats**: Supports netCDF, zarr, GRIB, and MITgcm binary (MDS) - **Unified data handling**: Treats all data as collections of points with lon/lat coordinates - **Fast visualization**: KDTree-based nearest-neighbor interpolation to regular grid - **X11/Xaw interface**: Works over SSH with X forwarding @@ -189,7 +189,7 @@ No libraries should show as "not found". ## Usage ```bash -./ushow [options] [file2 ...] +./ushow [options] [file2 ...] Options: -m, --mesh Mesh file with coordinates (for unstructured data) @@ -258,6 +258,12 @@ GRIB file (requires `make WITH_GRIB=1`): ./uterm data.grib --color ``` +MITgcm binary data (MDS format, always compiled — no extra libraries): +```bash +./ushow /path/to/mitgcm/run # Open directory with .data/.meta files +./ushow /path/to/diags3D.0000008760.data # Open specific .data file (uses parent dir) +``` + MPAS unstructured data: ```bash ./ushow data.nc grid.nc # MPAS UGRID (face_lon/face_lat) @@ -310,6 +316,7 @@ The test suite includes: - **test_range_popup**: Range popup logic (symmetric computation, value parsing) - **test_timeseries**: Time series reading, multi-file concatenation, and CF time unit conversion - **test_file_netcdf**: NetCDF file I/O +- **test_file_mitgcm**: MITgcm MDS binary file I/O - **test_file_zarr**: Zarr file I/O (when built with `WITH_ZARR=1`) - **test_yac_3d**: Per-depth masked YAC interpolation (when built with `WITH_YAC=1`) - **test_integration**: End-to-end workflow tests @@ -389,6 +396,13 @@ The first visit to each new depth level builds a masked interpolation (timing de ## Supported Data Formats - **NetCDF**: Full support (NetCDF-3 and NetCDF-4) +- **MITgcm MDS**: Binary `.data`/`.meta` file pairs (always compiled, no extra libraries) + - Reads grid coordinates from XC/YC binary files + - Supports 2D and 3D diagnostic fields with multiple variables per file + - Automatic land masking via hFacC grid file + - Depth levels from RC grid file + - Multi-timestep iteration discovery + - Automatic velocity rotation on LLC/cubed-sphere grids: if AngleCS/AngleSN (or CS/SN) grid files are present, UVEL/VVEL and other paired velocity fields are rotated from local grid coordinates to geographic (eastward/northward) components - **Zarr**: v2 format (requires `make WITH_ZARR=1`) - Compression: LZ4, Blosc (with various inner codecs), or uncompressed - Data types: Float32, Float64, Int64 diff --git a/desighn/mitgcm_velocity_rotation.md b/desighn/mitgcm_velocity_rotation.md new file mode 100644 index 0000000..ccdfb6f --- /dev/null +++ b/desighn/mitgcm_velocity_rotation.md @@ -0,0 +1,107 @@ +# MITgcm Velocity Rotation on LLC/Cubed-Sphere Grids + +## The Problem + +On MITgcm LLC (Lat-Lon-Cap) and cubed-sphere grids, UVEL and VVEL are stored in **local grid coordinates**, not geographic (eastward/northward). Each face of the cubed sphere has its own local coordinate system, so the i-direction and j-direction point in different geographic directions on different faces. + +When visualized on a lon-lat map without rotation, this causes visible discontinuities at face boundaries — velocities appear to jump abruptly where one face meets another, because the local axes change orientation. + +Tracer fields (THETA, SALT, ETAN, etc.) are scalar values and are unaffected by grid orientation — they display correctly as-is. + +### Grid staggering (secondary issue) + +UVEL lives on the XG/YC grid (western cell face), VVEL lives on the XC/YG grid (southern cell face), while tracers live on XC/YC (cell centers). This creates a small spatial offset but is NOT the primary cause of the visible discontinuities. The dominant effect is the face rotation. + +## The Solution: CS/SN Rotation + +MITgcm can output two grid files: +- `AngleCS.data` / `CS.data` — cosine of the angle between grid-i and geographic-east +- `AngleSN.data` / `SN.data` — sine of the angle between grid-i and geographic-east + +With these, geographic velocities are computed as: + +``` +u_east = CS * UVEL - SN * VVEL +v_north = SN * UVEL + CS * VVEL +``` + +This transforms the local (i,j) velocity components to proper (east, north) components. + +## Implementation Plan + +### 1. Load CS/SN during `mitgcm_open()` + +In `MitgcmStore`, add: +```c +float *angle_cs; /* AngleCS or CS [ny*nx], NULL if unavailable */ +float *angle_sn; /* AngleSN or SN [ny*nx], NULL if unavailable */ +``` + +During `mitgcm_open()`, try loading in order: +1. `AngleCS.data` / `AngleSN.data` +2. `CS.data` / `SN.data` + +Both must be present; if either is missing, set both to NULL (no rotation). + +### 2. Mark velocity variables + +In `MitgcmVarData`, add: +```c +int velocity_component; /* 0=none, 1=U-component, 2=V-component */ +``` + +During `mitgcm_scan_variables()`, detect velocity fields by name: +- U-component: names starting with "U" or containing "UVEL" +- V-component: names starting with "V" or containing "VVEL" + +This heuristic covers UVEL, VVEL, and diagnostic velocity fields. + +### 3. Pair U/V variables + +When scanning variables, if both a U-component and V-component exist in the same diagnostic group (same prefix, same dimensions), link them: +```c +USVar *paired_velocity; /* Pointer to the other component */ +``` + +### 4. Apply rotation during `mitgcm_read_slice()` + +After reading the raw slab, if `angle_cs` is available and the variable is a velocity component: + +```c +if (store->angle_cs && vd->velocity_component != 0 && vd->paired_velocity) { + /* Need to read the paired component too */ + float *other = malloc(slab_size * sizeof(float)); + /* read paired component at same time/depth */ + + if (vd->velocity_component == 1) { + /* This is U, other is V */ + for (i = 0; i < slab_size; i++) { + float u = data[i], v = other[i]; + data[i] = store->angle_cs[i] * u - store->angle_sn[i] * v; + } + } else { + /* This is V, other is U */ + for (i = 0; i < slab_size; i++) { + float u = other[i], v = data[i]; + data[i] = store->angle_sn[i] * u + store->angle_cs[i] * v; + } + } + free(other); +} +``` + +### 5. Handle missing CS/SN gracefully + +If CS/SN files are not present, velocities display in raw local coordinates (current behavior). No error, no warning — the data is still valid, just not geographically oriented. + +Optionally print an info message: +``` +MITgcm: AngleCS/AngleSN not found, velocity fields shown in local grid coordinates +``` + +### Considerations + +- **Performance**: rotation requires reading two fields per slice instead of one. For interactive use this should be fine (binary reads are fast). +- **Timeseries**: the `mitgcm_read_timeseries()` function would also need to read both components at each timestep for rotated output. This doubles the I/O for velocity time series. +- **Grid staggering**: for full correctness, UVEL should be interpolated from XG to XC before rotation. This is a C-grid averaging: `u_at_center[i] = 0.5 * (uvel[i] + uvel[i-1])` in the i-direction (with wrapping within each face). This is a secondary refinement. +- **Variable naming**: rotated variables could be renamed to indicate geographic orientation (e.g., "UVEL_east", "VVEL_north"), or the original names kept with a visual indicator. diff --git a/src/file_mitgcm.c b/src/file_mitgcm.c new file mode 100644 index 0000000..e4f8827 --- /dev/null +++ b/src/file_mitgcm.c @@ -0,0 +1,1163 @@ +/* + * file_mitgcm.c - MITgcm MDS binary data reading + * + * Reads MITgcm .data/.meta file pairs. The .meta files are ASCII with + * MATLAB-style syntax; .data files are raw big-endian IEEE float arrays. + */ + +#include "file_mitgcm.h" +#include "mesh.h" +#include +#include +#include +#include +#include +#include +#include + +/* Internal structures */ + +typedef struct { + int n_dims; + int dims[MAX_DIMS]; /* C-order dimensions (reversed from Fortran) */ + char dataprec[32]; /* "float32" or "float64" */ + int n_records; + int n_fields; + char fields[MAX_VARS][16]; /* field names (stripped) */ + float missing_value; + int has_missing; + int timestep_number; + int has_timestep; +} MetaInfo; + +typedef struct { + char directory[MAX_NAME_LEN]; + int nx, ny, nz; + float *depth_values; /* RC[nz] */ + float *hfac; /* hFacC[nz*ny*nx], NULL if unavailable */ + float *angle_cs; /* AngleCS[ny*nx], NULL if unavailable */ + float *angle_sn; /* AngleSN[ny*nx], NULL if unavailable */ + float missing_value; +} MitgcmStore; + +typedef struct { + char prefix[MAX_NAME_LEN]; /* e.g. "diags3D" */ + int field_index; /* index within multi-field file */ + int *iterations; /* sorted iteration numbers */ + int n_iterations; + int is_3d; + int nx, ny, nz; + float missing_value; + int has_missing; + int velocity_component; /* 0=none, 1=U-component, 2=V-component */ + USVar *paired_velocity; /* paired U/V variable, NULL if unpaired */ +} MitgcmVarData; + +/* ---- Byte swapping ---- */ + +static int is_little_endian(void) { + uint32_t x = 1; + return *(unsigned char *)&x; +} + +static void swap_endian_float(float *data, size_t n) { + uint32_t *p = (uint32_t *)data; + for (size_t i = 0; i < n; i++) { + uint32_t v = p[i]; + p[i] = ((v >> 24) & 0xFF) | ((v >> 8) & 0xFF00) | + ((v << 8) & 0xFF0000) | ((v << 24) & 0xFF000000); + } +} + +/* ---- Meta file parsing ---- */ + +/* Skip whitespace */ +static const char *skip_ws(const char *p) { + while (*p == ' ' || *p == '\t' || *p == '\n' || *p == '\r') p++; + return p; +} + +/* Find value after "key = " in buffer. Returns pointer after '=' or NULL. */ +static const char *find_key(const char *buf, const char *key) { + const char *p = buf; + size_t klen = strlen(key); + while ((p = strstr(p, key)) != NULL) { + const char *after = p + klen; + after = skip_ws(after); + if (*after == '=') { + return skip_ws(after + 1); + } + p = after; + } + return NULL; +} + +/* Parse integer from "[ NNN ]" or "[ NNN ;" */ +static int parse_bracket_int(const char *p, int *val) { + p = skip_ws(p); + if (*p != '[') return 0; + p = skip_ws(p + 1); + char *end; + long v = strtol(p, &end, 10); + if (end == p) return 0; + *val = (int)v; + return 1; +} + +/* Parse float from "[ NNN ]" */ +static int parse_bracket_float(const char *p, float *val) { + p = skip_ws(p); + if (*p != '[') return 0; + p = skip_ws(p + 1); + char *end; + double v = strtod(p, &end); + if (end == p) return 0; + *val = (float)v; + return 1; +} + +/* Parse a string from "[ 'value' ]" */ +static int parse_bracket_string(const char *p, char *out, size_t outlen) { + p = skip_ws(p); + if (*p != '[') return 0; + p = skip_ws(p + 1); + if (*p != '\'') return 0; + p++; + size_t i = 0; + while (*p && *p != '\'' && i < outlen - 1) { + out[i++] = *p++; + } + out[i] = '\0'; + return 1; +} + +static int parse_meta(const char *meta_path, MetaInfo *info) { + memset(info, 0, sizeof(*info)); + info->missing_value = -999.0f; + + FILE *fp = fopen(meta_path, "r"); + if (!fp) return -1; + + /* Read entire file */ + fseek(fp, 0, SEEK_END); + long fsize = ftell(fp); + fseek(fp, 0, SEEK_SET); + if (fsize <= 0 || fsize > 100000) { fclose(fp); return -1; } + + char *buf = malloc(fsize + 1); + if (!buf) { fclose(fp); return -1; } + size_t nread = fread(buf, 1, fsize, fp); + fclose(fp); + buf[nread] = '\0'; + + /* nDims */ + const char *p = find_key(buf, "nDims"); + if (p) parse_bracket_int(p, &info->n_dims); + + /* dimList - parse triplets, take first value of each triplet */ + p = find_key(buf, "dimList"); + if (p) { + p = skip_ws(p); + if (*p == '[') p++; + int fortran_dims[MAX_DIMS]; + int nd = 0; + while (nd < info->n_dims && nd < MAX_DIMS) { + p = skip_ws(p); + char *end; + long v1 = strtol(p, &end, 10); + if (end == p) break; + fortran_dims[nd] = (int)v1; + /* Skip the other two values in the triplet */ + p = end; + for (int skip = 0; skip < 2; skip++) { + while (*p && (*p == ',' || *p == ' ' || *p == '\t' || *p == '\n' || *p == '\r')) p++; + strtol(p, &end, 10); + p = end; + } + while (*p && (*p == ',' || *p == ' ' || *p == '\t' || *p == '\n' || *p == '\r')) p++; + nd++; + } + /* Reverse Fortran order to C order */ + for (int i = 0; i < nd; i++) { + info->dims[i] = fortran_dims[nd - 1 - i]; + } + } + + /* dataprec */ + p = find_key(buf, "dataprec"); + if (p) parse_bracket_string(p, info->dataprec, sizeof(info->dataprec)); + + /* nrecords */ + p = find_key(buf, "nrecords"); + if (p) parse_bracket_int(p, &info->n_records); + + /* timeStepNumber */ + p = find_key(buf, "timeStepNumber"); + if (p && parse_bracket_int(p, &info->timestep_number)) { + info->has_timestep = 1; + } + + /* missingValue */ + p = find_key(buf, "missingValue"); + if (p && parse_bracket_float(p, &info->missing_value)) { + info->has_missing = 1; + } + + /* nFlds */ + p = find_key(buf, "nFlds"); + if (p) parse_bracket_int(p, &info->n_fields); + + /* fldList */ + p = find_key(buf, "fldList"); + if (p) { + /* Skip to opening brace */ + while (*p && *p != '{') p++; + if (*p == '{') p++; + + int fi = 0; + while (fi < info->n_fields && fi < MAX_VARS) { + p = skip_ws(p); + if (*p == '}' || *p == '\0') break; + if (*p == '\'') { + p++; + size_t j = 0; + while (*p && *p != '\'' && j < 15) { + info->fields[fi][j++] = *p++; + } + info->fields[fi][j] = '\0'; + if (*p == '\'') p++; + /* Trim trailing spaces */ + while (j > 0 && info->fields[fi][j - 1] == ' ') { + info->fields[fi][--j] = '\0'; + } + fi++; + } else { + p++; + } + } + } + + free(buf); + return 0; +} + +/* ---- Binary reading helpers ---- */ + +static float *read_binary_floats(const char *path, size_t offset, size_t count) { + FILE *fp = fopen(path, "rb"); + if (!fp) return NULL; + + float *data = malloc(count * sizeof(float)); + if (!data) { fclose(fp); return NULL; } + + if (offset > 0) fseek(fp, (long)offset, SEEK_SET); + + size_t nread = fread(data, sizeof(float), count, fp); + fclose(fp); + + if (nread != count) { + free(data); + return NULL; + } + + if (is_little_endian()) { + swap_endian_float(data, count); + } + + return data; +} + +/* ---- Path helpers ---- */ + +static int is_directory(const char *path) { + struct stat st; + return (stat(path, &st) == 0 && S_ISDIR(st.st_mode)); +} + +static int file_exists(const char *path) { + struct stat st; + return (stat(path, &st) == 0); +} + +/* Build companion .meta path from .data path */ +static void data_to_meta_path(const char *data_path, char *meta_path, size_t len) { + strncpy(meta_path, data_path, len - 1); + meta_path[len - 1] = '\0'; + size_t slen = strlen(meta_path); + if (slen > 5 && strcmp(meta_path + slen - 5, ".data") == 0) { + strcpy(meta_path + slen - 5, ".meta"); + } +} + +/* ---- Public API ---- */ + +int mitgcm_is_mitgcm(const char *path) { + if (!path) return 0; + + /* If it's a .data file, check for companion .meta */ + size_t plen = strlen(path); + if (plen > 5 && strcmp(path + plen - 5, ".data") == 0) { + char meta[MAX_NAME_LEN]; + data_to_meta_path(path, meta, sizeof(meta)); + return file_exists(meta); + } + + /* If it's a directory, check for grid meta files or any .meta */ + if (is_directory(path)) { + char test[MAX_NAME_LEN]; + snprintf(test, sizeof(test), "%s/XC.meta", path); + if (file_exists(test)) return 1; + + /* Check for any .meta file with iteration number pattern */ + DIR *dir = opendir(path); + if (!dir) return 0; + struct dirent *entry; + while ((entry = readdir(dir)) != NULL) { + size_t nlen = strlen(entry->d_name); + if (nlen > 5 && strcmp(entry->d_name + nlen - 5, ".meta") == 0) { + closedir(dir); + return 1; + } + } + closedir(dir); + } + + return 0; +} + +USFile *mitgcm_open(const char *path) { + if (!path) return NULL; + + char directory[MAX_NAME_LEN]; + + /* Resolve directory */ + if (is_directory(path)) { + strncpy(directory, path, MAX_NAME_LEN - 1); + directory[MAX_NAME_LEN - 1] = '\0'; + } else { + /* Extract parent directory from .data file path */ + strncpy(directory, path, MAX_NAME_LEN - 1); + directory[MAX_NAME_LEN - 1] = '\0'; + char *slash = strrchr(directory, '/'); + if (slash) { + *slash = '\0'; + } else { + strcpy(directory, "."); + } + } + + /* Remove trailing slash */ + size_t dlen = strlen(directory); + if (dlen > 1 && directory[dlen - 1] == '/') directory[dlen - 1] = '\0'; + + /* Determine grid dimensions from meta files */ + int nx = 0, ny = 0, nz = 0; + + /* Try XC.meta first (2D: nx, ny) */ + char meta_path[MAX_NAME_LEN]; + snprintf(meta_path, sizeof(meta_path), "%s/XC.meta", directory); + MetaInfo minfo; + if (parse_meta(meta_path, &minfo) == 0 && minfo.n_dims >= 2) { + /* C-order: dims are [ny, nx] for 2D */ + ny = minfo.dims[0]; + nx = minfo.dims[1]; + } + + /* Try hFacC.meta for nz (3D: nz, ny, nx) */ + snprintf(meta_path, sizeof(meta_path), "%s/hFacC.meta", directory); + if (parse_meta(meta_path, &minfo) == 0 && minfo.n_dims >= 3) { + nz = minfo.dims[0]; + if (nx == 0) { ny = minfo.dims[1]; nx = minfo.dims[2]; } + } + + /* If still no grid dims, scan for any diagnostic meta with dimList */ + if (nx == 0 || ny == 0) { + DIR *dir = opendir(directory); + if (dir) { + struct dirent *entry; + while ((entry = readdir(dir)) != NULL) { + size_t nlen = strlen(entry->d_name); + if (nlen > 5 && strcmp(entry->d_name + nlen - 5, ".meta") == 0 && + strchr(entry->d_name, '.') != entry->d_name + nlen - 5) { + /* Has a dot before .meta → likely a diagnostic file */ + snprintf(meta_path, sizeof(meta_path), "%s/%s", directory, entry->d_name); + if (parse_meta(meta_path, &minfo) == 0 && minfo.n_dims >= 2) { + if (minfo.n_dims == 2) { + ny = minfo.dims[0]; nx = minfo.dims[1]; + } else if (minfo.n_dims == 3) { + nz = minfo.dims[0]; ny = minfo.dims[1]; nx = minfo.dims[2]; + } + break; + } + } + } + closedir(dir); + } + } + + if (nx == 0 || ny == 0) { + fprintf(stderr, "MITgcm: Cannot determine grid dimensions from meta files in %s\n", directory); + return NULL; + } + + printf("MITgcm grid: %d x %d", nx, ny); + if (nz > 0) printf(" x %d", nz); + printf("\n"); + + /* Allocate store */ + MitgcmStore *store = calloc(1, sizeof(MitgcmStore)); + if (!store) return NULL; + strncpy(store->directory, directory, MAX_NAME_LEN - 1); + store->nx = nx; + store->ny = ny; + store->nz = nz; + store->missing_value = -999.0f; + + /* Read RC.data for depth values */ + if (nz > 0) { + char rc_path[MAX_NAME_LEN]; + snprintf(rc_path, sizeof(rc_path), "%s/RC.data", directory); + if (file_exists(rc_path)) { + store->depth_values = read_binary_floats(rc_path, 0, nz); + if (store->depth_values) { + printf("MITgcm: loaded %d depth levels (%.1f to %.1f m)\n", + nz, store->depth_values[0], store->depth_values[nz - 1]); + } + } + } + + /* Read hFacC.data for land masking */ + if (nz > 0) { + char hfac_path[MAX_NAME_LEN]; + snprintf(hfac_path, sizeof(hfac_path), "%s/hFacC.data", directory); + if (file_exists(hfac_path)) { + size_t total = (size_t)nz * ny * nx; + store->hfac = read_binary_floats(hfac_path, 0, total); + if (store->hfac) { + printf("MITgcm: loaded hFacC mask (%zu cells)\n", total); + } + } + } + + /* Load AngleCS/AngleSN (or CS/SN) for velocity rotation */ + { + size_t grid_size = (size_t)ny * nx; + char cs_path[MAX_NAME_LEN], sn_path[MAX_NAME_LEN]; + + /* Try AngleCS.data / AngleSN.data first */ + snprintf(cs_path, sizeof(cs_path), "%s/AngleCS.data", directory); + snprintf(sn_path, sizeof(sn_path), "%s/AngleSN.data", directory); + if (!file_exists(cs_path) || !file_exists(sn_path)) { + /* Fall back to CS.data / SN.data */ + snprintf(cs_path, sizeof(cs_path), "%s/CS.data", directory); + snprintf(sn_path, sizeof(sn_path), "%s/SN.data", directory); + } + + if (file_exists(cs_path) && file_exists(sn_path)) { + store->angle_cs = read_binary_floats(cs_path, 0, grid_size); + store->angle_sn = read_binary_floats(sn_path, 0, grid_size); + if (store->angle_cs && store->angle_sn) { + printf("MITgcm: loaded AngleCS/AngleSN for velocity rotation\n"); + } else { + /* Both must succeed */ + free(store->angle_cs); store->angle_cs = NULL; + free(store->angle_sn); store->angle_sn = NULL; + } + } + } + + /* Create USFile */ + USFile *f = calloc(1, sizeof(USFile)); + if (!f) { free(store->depth_values); free(store->hfac); + free(store->angle_cs); free(store->angle_sn); + free(store); return NULL; } + + f->file_type = FILE_TYPE_MITGCM; + strncpy(f->filename, directory, MAX_NAME_LEN - 1); + f->mitgcm_data = store; + + return f; +} + +USMesh *mitgcm_create_mesh(USFile *file) { + if (!file || !file->mitgcm_data) return NULL; + MitgcmStore *store = (MitgcmStore *)file->mitgcm_data; + + int nx = store->nx, ny = store->ny; + size_t n_points = (size_t)ny * nx; + + /* Read XC.data and YC.data */ + char xc_path[MAX_NAME_LEN], yc_path[MAX_NAME_LEN]; + snprintf(xc_path, sizeof(xc_path), "%s/XC.data", store->directory); + snprintf(yc_path, sizeof(yc_path), "%s/YC.data", store->directory); + + float *xc = read_binary_floats(xc_path, 0, n_points); + float *yc = read_binary_floats(yc_path, 0, n_points); + if (!xc || !yc) { + fprintf(stderr, "MITgcm: Failed to read XC/YC grid files\n"); + free(xc); free(yc); + return NULL; + } + + /* Convert float to double for mesh_create */ + double *lon = malloc(n_points * sizeof(double)); + double *lat = malloc(n_points * sizeof(double)); + if (!lon || !lat) { + free(xc); free(yc); free(lon); free(lat); + return NULL; + } + + for (size_t i = 0; i < n_points; i++) { + lon[i] = (double)xc[i]; + lat[i] = (double)yc[i]; + } + free(xc); + free(yc); + + USMesh *mesh = mesh_create(lon, lat, n_points, COORD_TYPE_2D_CURVILINEAR); + if (mesh) { + mesh->orig_nx = nx; + mesh->orig_ny = ny; + printf("MITgcm mesh: %zu points (%dx%d curvilinear)\n", n_points, nx, ny); + } + return mesh; +} + +/* ---- Iteration discovery ---- */ + +/* Compare ints for qsort */ +static int cmp_int(const void *a, const void *b) { + return (*(const int *)a) - (*(const int *)b); +} + +/* + * Find all iteration numbers for a given prefix. + * Scans for files matching .NNNNNNNNNN.meta + */ +static int *find_iterations(const char *directory, const char *prefix, + int *n_out) { + *n_out = 0; + DIR *dir = opendir(directory); + if (!dir) return NULL; + + size_t prefix_len = strlen(prefix); + int capacity = 32; + int *iters = malloc(capacity * sizeof(int)); + if (!iters) { closedir(dir); return NULL; } + + struct dirent *entry; + while ((entry = readdir(dir)) != NULL) { + size_t nlen = strlen(entry->d_name); + /* Match: prefix.DIGITS.meta */ + if (nlen < prefix_len + 7) continue; /* prefix + . + at least 1 digit + .meta */ + if (strncmp(entry->d_name, prefix, prefix_len) != 0) continue; + if (entry->d_name[prefix_len] != '.') continue; + if (nlen < 6 || strcmp(entry->d_name + nlen - 5, ".meta") != 0) continue; + + /* Extract digits between prefix. and .meta */ + const char *digit_start = entry->d_name + prefix_len + 1; + const char *digit_end = entry->d_name + nlen - 5; + if (digit_end <= digit_start) continue; + + /* Check all digits */ + int all_digits = 1; + for (const char *c = digit_start; c < digit_end; c++) { + if (*c < '0' || *c > '9') { all_digits = 0; break; } + } + if (!all_digits) continue; + + char num_buf[32]; + size_t num_len = digit_end - digit_start; + if (num_len >= sizeof(num_buf)) continue; + memcpy(num_buf, digit_start, num_len); + num_buf[num_len] = '\0'; + + int iter = atoi(num_buf); + + if (*n_out >= capacity) { + capacity *= 2; + int *tmp = realloc(iters, capacity * sizeof(int)); + if (!tmp) break; + iters = tmp; + } + iters[(*n_out)++] = iter; + } + closedir(dir); + + if (*n_out > 1) { + qsort(iters, *n_out, sizeof(int), cmp_int); + } + + return iters; +} + +USVar *mitgcm_scan_variables(USFile *f, USMesh *m) { + if (!f || !f->mitgcm_data || !m) return NULL; + MitgcmStore *store = (MitgcmStore *)f->mitgcm_data; + + /* Discover unique prefixes by scanning for *.DIGITS.meta */ + DIR *dir = opendir(store->directory); + if (!dir) return NULL; + + /* Collect unique prefixes */ + char prefixes[64][MAX_NAME_LEN]; + int n_prefixes = 0; + + struct dirent *entry; + while ((entry = readdir(dir)) != NULL) { + size_t nlen = strlen(entry->d_name); + if (nlen < 7 || strcmp(entry->d_name + nlen - 5, ".meta") != 0) continue; + + /* Find first dot that precedes a digit sequence before .meta */ + const char *meta_ext = entry->d_name + nlen - 5; + const char *last_dot = NULL; + for (const char *c = meta_ext - 1; c > entry->d_name; c--) { + if (*c == '.') { last_dot = c; break; } + } + if (!last_dot || last_dot == entry->d_name) continue; + + /* Check that between last_dot+1 and meta_ext are all digits */ + int all_digits = 1; + for (const char *c = last_dot + 1; c < meta_ext; c++) { + if (*c < '0' || *c > '9') { all_digits = 0; break; } + } + if (!all_digits) continue; + + /* Extract prefix */ + size_t plen = last_dot - entry->d_name; + if (plen >= MAX_NAME_LEN) continue; + + char prefix[MAX_NAME_LEN]; + memcpy(prefix, entry->d_name, plen); + prefix[plen] = '\0'; + + /* Check uniqueness */ + int found = 0; + for (int i = 0; i < n_prefixes; i++) { + if (strcmp(prefixes[i], prefix) == 0) { found = 1; break; } + } + if (!found && n_prefixes < 64) { + strncpy(prefixes[n_prefixes++], prefix, MAX_NAME_LEN - 1); + } + } + closedir(dir); + + if (n_prefixes == 0) { + fprintf(stderr, "MITgcm: No diagnostic files found in %s\n", store->directory); + return NULL; + } + + /* For each prefix, parse meta and create variables */ + USVar *head = NULL; + USVar *tail = NULL; + int n_vars = 0; + + for (int pi = 0; pi < n_prefixes; pi++) { + /* Find iterations for this prefix */ + int n_iters = 0; + int *iterations = find_iterations(store->directory, prefixes[pi], &n_iters); + if (!iterations || n_iters == 0) { free(iterations); continue; } + + /* Parse meta from first iteration to get field info */ + char meta_path[MAX_NAME_LEN]; + snprintf(meta_path, sizeof(meta_path), "%s/%s.%010d.meta", + store->directory, prefixes[pi], iterations[0]); + + MetaInfo minfo; + if (parse_meta(meta_path, &minfo) != 0) { + free(iterations); + continue; + } + + int is_3d = (minfo.n_dims >= 3 && minfo.dims[0] > 1); + int var_nx = store->nx; + int var_ny = store->ny; + int var_nz = is_3d ? minfo.dims[0] : 1; + + /* If no fldList, treat as single anonymous field from prefix */ + int actual_fields = minfo.n_fields > 0 ? minfo.n_fields : 1; + + for (int fi = 0; fi < actual_fields; fi++) { + USVar *var = calloc(1, sizeof(USVar)); + if (!var) continue; + + /* Name */ + if (minfo.n_fields > 0 && minfo.fields[fi][0]) { + strncpy(var->name, minfo.fields[fi], MAX_NAME_LEN - 1); + } else { + strncpy(var->name, prefixes[pi], MAX_NAME_LEN - 1); + } + + var->mesh = m; + var->file = f; + var->fill_value = minfo.has_missing ? minfo.missing_value : DEFAULT_FILL_VALUE; + + /* Set up dimensions */ + int dim_idx = 0; + + /* Time dimension (iterations) */ + if (n_iters > 0) { + var->time_dim_id = dim_idx; + strncpy(var->dim_names[dim_idx], "time", MAX_NAME_LEN - 1); + var->dim_sizes[dim_idx] = n_iters; + dim_idx++; + } else { + var->time_dim_id = -1; + } + + /* Depth dimension */ + if (is_3d && var_nz > 1) { + var->depth_dim_id = dim_idx; + strncpy(var->dim_names[dim_idx], "depth", MAX_NAME_LEN - 1); + var->dim_sizes[dim_idx] = var_nz; + dim_idx++; + } else { + var->depth_dim_id = -1; + } + + /* Spatial dimension */ + var->node_dim_id = dim_idx; + strncpy(var->dim_names[dim_idx], "node", MAX_NAME_LEN - 1); + var->dim_sizes[dim_idx] = (size_t)var_ny * var_nx; + dim_idx++; + + var->n_dims = dim_idx; + + /* Private data */ + MitgcmVarData *vd = calloc(1, sizeof(MitgcmVarData)); + if (!vd) { free(var); continue; } + strncpy(vd->prefix, prefixes[pi], MAX_NAME_LEN - 1); + vd->field_index = fi; + vd->iterations = malloc(n_iters * sizeof(int)); + if (vd->iterations) { + memcpy(vd->iterations, iterations, n_iters * sizeof(int)); + } + vd->n_iterations = n_iters; + vd->is_3d = is_3d; + vd->nx = var_nx; + vd->ny = var_ny; + vd->nz = var_nz; + vd->missing_value = minfo.has_missing ? minfo.missing_value : DEFAULT_FILL_VALUE; + vd->has_missing = minfo.has_missing; + var->mitgcm_data = vd; + + /* Append to list */ + var->next = NULL; + if (!head) { + head = var; + tail = var; + } else { + tail->next = var; + tail = var; + } + n_vars++; + } + + free(iterations); + } + + /* Detect velocity components and pair U/V variables */ + if (store->angle_cs) { + /* Mark velocity components by name */ + for (USVar *v = head; v; v = v->next) { + MitgcmVarData *d = (MitgcmVarData *)v->mitgcm_data; + if (v->name[0] == 'U') d->velocity_component = 1; + else if (v->name[0] == 'V') d->velocity_component = 2; + } + + /* Pair U/V variables: same prefix, name differs only in first char */ + for (USVar *u = head; u; u = u->next) { + MitgcmVarData *ud = (MitgcmVarData *)u->mitgcm_data; + if (ud->velocity_component != 1 || ud->paired_velocity) continue; + + for (USVar *v = head; v; v = v->next) { + MitgcmVarData *vvd = (MitgcmVarData *)v->mitgcm_data; + if (vvd->velocity_component != 2 || vvd->paired_velocity) continue; + if (strcmp(ud->prefix, vvd->prefix) != 0) continue; + if (strcmp(u->name + 1, v->name + 1) != 0) continue; + + ud->paired_velocity = v; + vvd->paired_velocity = u; + printf("MITgcm: paired %s / %s for rotation\n", u->name, v->name); + break; + } + } + } + + f->vars = head; + f->n_vars = n_vars; + printf("MITgcm: found %d variables from %d diagnostic groups\n", n_vars, n_prefixes); + + return head; +} + +int mitgcm_read_slice(USVar *var, size_t time_idx, size_t depth_idx, float *data) { + if (!var || !var->mitgcm_data || !var->file || !var->file->mitgcm_data || !data) + return -1; + + MitgcmVarData *vd = (MitgcmVarData *)var->mitgcm_data; + MitgcmStore *store = (MitgcmStore *)var->file->mitgcm_data; + + if ((int)time_idx >= vd->n_iterations) return -1; + + /* Build data file path */ + char data_path[MAX_NAME_LEN]; + snprintf(data_path, sizeof(data_path), "%s/%s.%010d.data", + store->directory, vd->prefix, vd->iterations[time_idx]); + + size_t slab_size = (size_t)vd->ny * vd->nx; + size_t field_size; + + if (vd->is_3d) { + field_size = (size_t)vd->nz * slab_size; + } else { + field_size = slab_size; + } + + /* Compute offset: field_index * field_size + depth_offset */ + size_t offset = (size_t)vd->field_index * field_size; + if (vd->is_3d) { + offset += depth_idx * slab_size; + } + offset *= sizeof(float); + + /* Read the 2D slab */ + FILE *fp = fopen(data_path, "rb"); + if (!fp) { + fprintf(stderr, "MITgcm: Cannot open %s\n", data_path); + return -1; + } + + fseek(fp, (long)offset, SEEK_SET); + size_t nread = fread(data, sizeof(float), slab_size, fp); + fclose(fp); + + if (nread != slab_size) { + fprintf(stderr, "MITgcm: Short read from %s (got %zu, expected %zu)\n", + data_path, nread, slab_size); + return -1; + } + + /* Byte swap */ + if (is_little_endian()) { + swap_endian_float(data, slab_size); + } + + /* Apply CS/SN velocity rotation */ + if (store->angle_cs && vd->velocity_component != 0 && vd->paired_velocity) { + MitgcmVarData *pvd = (MitgcmVarData *)vd->paired_velocity->mitgcm_data; + + /* Compute offset for the paired component in the same data file */ + size_t p_field_size = pvd->is_3d ? (size_t)pvd->nz * slab_size : slab_size; + size_t p_offset = (size_t)pvd->field_index * p_field_size; + if (pvd->is_3d) p_offset += depth_idx * slab_size; + p_offset *= sizeof(float); + + float *other = malloc(slab_size * sizeof(float)); + if (other) { + FILE *fp2 = fopen(data_path, "rb"); + if (fp2) { + fseek(fp2, (long)p_offset, SEEK_SET); + size_t pr = fread(other, sizeof(float), slab_size, fp2); + fclose(fp2); + + if (pr == slab_size) { + if (is_little_endian()) swap_endian_float(other, slab_size); + + if (vd->velocity_component == 1) { + /* This is U, other is V: u_east = CS*U - SN*V */ + for (size_t i = 0; i < slab_size; i++) { + float u = data[i], v = other[i]; + data[i] = store->angle_cs[i] * u - store->angle_sn[i] * v; + } + } else { + /* This is V, other is U: v_north = SN*U + CS*V */ + for (size_t i = 0; i < slab_size; i++) { + float u = other[i], v = data[i]; + data[i] = store->angle_sn[i] * u + store->angle_cs[i] * v; + } + } + } + } + free(other); + } + } + + /* Apply hFacC land masking */ + if (store->hfac && vd->is_3d) { + float *mask = store->hfac + depth_idx * slab_size; + for (size_t i = 0; i < slab_size; i++) { + if (mask[i] == 0.0f) data[i] = var->fill_value; + } + } else if (store->hfac && !vd->is_3d) { + /* Use surface level for 2D vars */ + for (size_t i = 0; i < slab_size; i++) { + if (store->hfac[i] == 0.0f) data[i] = var->fill_value; + } + } + + /* Apply missing value masking */ + if (vd->has_missing) { + for (size_t i = 0; i < slab_size; i++) { + if (data[i] == vd->missing_value || fabsf(data[i]) > INVALID_DATA_THRESHOLD) { + data[i] = var->fill_value; + } + } + } + + return 0; +} + +int mitgcm_estimate_range(USVar *var, float *min_val, float *max_val) { + if (!var || !var->mitgcm_data) return -1; + + MitgcmVarData *vd = (MitgcmVarData *)var->mitgcm_data; + size_t slab_size = (size_t)vd->ny * vd->nx; + + float *data = malloc(slab_size * sizeof(float)); + if (!data) return -1; + + float gmin = 1e30f, gmax = -1e30f; + int found = 0; + + /* Sample first, mid, last timestep at surface (depth 0) */ + int samples[] = {0, vd->n_iterations / 2, vd->n_iterations - 1}; + int n_samples = (vd->n_iterations >= 3) ? 3 : + (vd->n_iterations >= 2) ? 2 : 1; + + for (int s = 0; s < n_samples; s++) { + if (mitgcm_read_slice(var, samples[s], 0, data) != 0) continue; + + for (size_t i = 0; i < slab_size; i++) { + float v = data[i]; + if (v == var->fill_value || fabsf(v) > INVALID_DATA_THRESHOLD) continue; + if (v < gmin) gmin = v; + if (v > gmax) gmax = v; + found = 1; + } + } + + free(data); + + if (!found) { + *min_val = 0.0f; + *max_val = 1.0f; + return -1; + } + + *min_val = gmin; + *max_val = gmax; + return 0; +} + +USDimInfo *mitgcm_get_dim_info(USVar *var, int *n_dims_out) { + if (!var || !var->mitgcm_data || !var->file || !var->file->mitgcm_data) { + *n_dims_out = 0; + return NULL; + } + + MitgcmVarData *vd = (MitgcmVarData *)var->mitgcm_data; + MitgcmStore *store = (MitgcmStore *)var->file->mitgcm_data; + + int n_dims = 0; + if (var->time_dim_id >= 0) n_dims++; + if (var->depth_dim_id >= 0) n_dims++; + + if (n_dims == 0) { + *n_dims_out = 0; + return NULL; + } + + USDimInfo *dims = calloc(n_dims, sizeof(USDimInfo)); + if (!dims) { *n_dims_out = 0; return NULL; } + + int idx = 0; + + /* Time dimension (iteration numbers) */ + if (var->time_dim_id >= 0) { + USDimInfo *di = &dims[idx++]; + strncpy(di->name, "time", MAX_NAME_LEN - 1); + strncpy(di->units, "iteration", MAX_NAME_LEN - 1); + di->size = vd->n_iterations; + di->is_scannable = (vd->n_iterations > 1); + di->values = malloc(vd->n_iterations * sizeof(double)); + if (di->values) { + for (int i = 0; i < vd->n_iterations; i++) { + di->values[i] = (double)vd->iterations[i]; + } + di->min_val = di->values[0]; + di->max_val = di->values[vd->n_iterations - 1]; + } + } + + /* Depth dimension (RC values) */ + if (var->depth_dim_id >= 0) { + USDimInfo *di = &dims[idx++]; + strncpy(di->name, "depth", MAX_NAME_LEN - 1); + strncpy(di->units, "m", MAX_NAME_LEN - 1); + di->size = vd->nz; + di->is_scannable = (vd->nz > 1); + if (store->depth_values) { + di->values = malloc(vd->nz * sizeof(double)); + if (di->values) { + for (int i = 0; i < vd->nz; i++) { + di->values[i] = (double)store->depth_values[i]; + } + di->min_val = di->values[0]; + di->max_val = di->values[vd->nz - 1]; + } + } + } + + *n_dims_out = n_dims; + return dims; +} + +void mitgcm_free_dim_info(USDimInfo *dims, int n_dims) { + if (!dims) return; + for (int i = 0; i < n_dims; i++) { + free(dims[i].values); + } + free(dims); +} + +int mitgcm_read_timeseries(USVar *var, size_t node_idx, size_t depth_idx, + double **times_out, float **values_out, + int **valid_out, size_t *n_out) { + if (!var || !var->mitgcm_data || !var->file || !var->file->mitgcm_data) + return -1; + + MitgcmVarData *vd = (MitgcmVarData *)var->mitgcm_data; + MitgcmStore *store = (MitgcmStore *)var->file->mitgcm_data; + + size_t n_times = vd->n_iterations; + size_t slab_size = (size_t)vd->ny * vd->nx; + + if (node_idx >= slab_size) return -1; + + double *times = malloc(n_times * sizeof(double)); + float *values = malloc(n_times * sizeof(float)); + int *valid = malloc(n_times * sizeof(int)); + if (!times || !values || !valid) { + free(times); free(values); free(valid); + return -1; + } + + size_t field_size; + if (vd->is_3d) { + field_size = (size_t)vd->nz * slab_size; + } else { + field_size = slab_size; + } + + for (size_t t = 0; t < n_times; t++) { + times[t] = (double)vd->iterations[t]; + + /* Build path and read single value */ + char data_path[MAX_NAME_LEN]; + snprintf(data_path, sizeof(data_path), "%s/%s.%010d.data", + store->directory, vd->prefix, vd->iterations[t]); + + size_t offset = (size_t)vd->field_index * field_size; + if (vd->is_3d) { + offset += depth_idx * slab_size; + } + offset += node_idx; + offset *= sizeof(float); + + FILE *fp = fopen(data_path, "rb"); + if (!fp) { + values[t] = var->fill_value; + valid[t] = 0; + continue; + } + + float val; + fseek(fp, (long)offset, SEEK_SET); + if (fread(&val, sizeof(float), 1, fp) == 1) { + if (is_little_endian()) swap_endian_float(&val, 1); + + /* Apply CS/SN velocity rotation */ + if (store->angle_cs && vd->velocity_component != 0 && vd->paired_velocity) { + MitgcmVarData *pvd = (MitgcmVarData *)vd->paired_velocity->mitgcm_data; + size_t p_field_size = pvd->is_3d ? (size_t)pvd->nz * slab_size : slab_size; + size_t p_off = (size_t)pvd->field_index * p_field_size; + if (pvd->is_3d) p_off += depth_idx * slab_size; + p_off += node_idx; + p_off *= sizeof(float); + + float other_val; + fseek(fp, (long)p_off, SEEK_SET); + if (fread(&other_val, sizeof(float), 1, fp) == 1) { + if (is_little_endian()) swap_endian_float(&other_val, 1); + if (vd->velocity_component == 1) { + /* U: u_east = CS*U - SN*V */ + val = store->angle_cs[node_idx] * val + - store->angle_sn[node_idx] * other_val; + } else { + /* V: v_north = SN*U + CS*V */ + float u = other_val, v = val; + val = store->angle_sn[node_idx] * u + + store->angle_cs[node_idx] * v; + } + } + } + + /* Check hFacC mask */ + int masked = 0; + if (store->hfac && vd->is_3d) { + if (store->hfac[depth_idx * slab_size + node_idx] == 0.0f) + masked = 1; + } else if (store->hfac && !vd->is_3d) { + if (store->hfac[node_idx] == 0.0f) + masked = 1; + } + + if (masked || val == vd->missing_value || fabsf(val) > INVALID_DATA_THRESHOLD) { + values[t] = var->fill_value; + valid[t] = 0; + } else { + values[t] = val; + valid[t] = 1; + } + } else { + values[t] = var->fill_value; + valid[t] = 0; + } + fclose(fp); + } + + *times_out = times; + *values_out = values; + *valid_out = valid; + *n_out = n_times; + return 0; +} + +void mitgcm_close(USFile *file) { + if (!file) return; + + if (file->mitgcm_data) { + MitgcmStore *store = (MitgcmStore *)file->mitgcm_data; + free(store->depth_values); + free(store->hfac); + free(store->angle_cs); + free(store->angle_sn); + free(store); + } + + /* Free variables */ + USVar *var = file->vars; + while (var) { + USVar *next = var->next; + if (var->mitgcm_data) { + MitgcmVarData *vd = (MitgcmVarData *)var->mitgcm_data; + free(vd->iterations); + free(vd); + } + free(var); + var = next; + } + + free(file); +} diff --git a/src/file_mitgcm.h b/src/file_mitgcm.h new file mode 100644 index 0000000..14d1bf7 --- /dev/null +++ b/src/file_mitgcm.h @@ -0,0 +1,68 @@ +/* + * file_mitgcm.h - MITgcm MDS binary data reading + * + * Reads MITgcm .data/.meta file pairs (MDS format). + * Always compiled (no external library dependencies). + */ + +#ifndef FILE_MITGCM_H +#define FILE_MITGCM_H + +#include "ushow.defines.h" + +/* + * Check if a path looks like MITgcm data (directory with .meta files, + * or a .data file with companion .meta). + */ +int mitgcm_is_mitgcm(const char *path); + +/* + * Open MITgcm data. Path can be a directory or a .data file. + * Resolves grid dimensions from meta files, loads RC and hFacC if present. + */ +USFile *mitgcm_open(const char *path); + +/* + * Create mesh from XC.data/YC.data grid files. + */ +USMesh *mitgcm_create_mesh(USFile *file); + +/* + * Scan directory for diagnostic variables. + * Returns linked list of USVar. + */ +USVar *mitgcm_scan_variables(USFile *f, USMesh *m); + +/* + * Read a 2D slice of data. + */ +int mitgcm_read_slice(USVar *var, size_t time_idx, size_t depth_idx, float *data); + +/* + * Estimate min/max range by sampling. + */ +int mitgcm_estimate_range(USVar *var, float *min_val, float *max_val); + +/* + * Get dimension info (time iterations, depth levels). + */ +USDimInfo *mitgcm_get_dim_info(USVar *var, int *n_dims_out); + +/* + * Free dimension info. + */ +void mitgcm_free_dim_info(USDimInfo *dims, int n_dims); + +/* + * Read time series at a single spatial node. + */ +int mitgcm_read_timeseries(USVar *var, size_t node_idx, size_t depth_idx, + double **times_out, float **values_out, + int **valid_out, size_t *n_out); + +/* + * Close and free MITgcm file data. + */ +void mitgcm_close(USFile *file); + +#endif /* FILE_MITGCM_H */ diff --git a/src/ushow.c b/src/ushow.c index 130b153..6dac195 100644 --- a/src/ushow.c +++ b/src/ushow.c @@ -11,6 +11,7 @@ #include "regrid_yac.h" #endif #include "file_netcdf.h" +#include "file_mitgcm.h" #ifdef HAVE_ZARR #include "file_zarr.h" #endif @@ -87,6 +88,9 @@ static void on_var_select(int var_index) { /* Set up dimension info panel */ if (current_dim_info) { + if (prev_var && prev_var->file && prev_var->file->file_type == FILE_TYPE_MITGCM) { + mitgcm_free_dim_info(current_dim_info, n_current_dims); + } else #ifdef HAVE_GRIB if (prev_var && prev_var->file && prev_var->file->file_type == FILE_TYPE_GRIB) { grib_free_dim_info(current_dim_info, n_current_dims); @@ -97,6 +101,9 @@ static void on_var_select(int var_index) { } } /* Use fileset dimension info if available (includes virtual time) */ + if (var->file && var->file->file_type == FILE_TYPE_MITGCM) { + current_dim_info = mitgcm_get_dim_info(var, &n_current_dims); + } else #ifdef HAVE_ZARR if (zarr_fileset) { current_dim_info = zarr_get_dim_info_fileset(zarr_fileset, var, &n_current_dims); @@ -305,6 +312,10 @@ static void on_mouse_click(int px, int py) { size_t n_out = 0; int rc; + if (current_var->file && current_var->file->file_type == FILE_TYPE_MITGCM) { + rc = mitgcm_read_timeseries(current_var, node_idx, view->depth_index, + ×, &values, &valid, &n_out); + } else #ifdef HAVE_ZARR if (zarr_fileset) { rc = zarr_read_timeseries_fileset(zarr_fileset, current_var, node_idx, @@ -943,6 +954,14 @@ int main(int argc, char *argv[]) { printf("Opening data file(s)...\n"); if (!use_glob && n_data_files == 1) { + if (mitgcm_is_mitgcm(data_filenames[0])) { + printf("Detected MITgcm data: %s\n", data_filenames[0]); + file = mitgcm_open(data_filenames[0]); + if (!file) { + fprintf(stderr, "Failed to open MITgcm data: %s\n", data_filenames[0]); + return 1; + } + } else #ifdef HAVE_ZARR /* Check if first file is a zarr store */ if (zarr_is_zarr_store(data_filenames[0])) { @@ -1092,6 +1111,14 @@ int main(int argc, char *argv[]) { /* Load mesh */ printf("Loading mesh...\n"); + if (file->file_type == FILE_TYPE_MITGCM) { + mesh = mitgcm_create_mesh(file); + if (!mesh) { + fprintf(stderr, "Failed to load mesh from MITgcm grid files\n"); + mitgcm_close(file); + return 1; + } + } else #ifdef HAVE_ZARR if (file->file_type == FILE_TYPE_ZARR) { /* For zarr, load coordinates from the store itself */ @@ -1169,6 +1196,9 @@ int main(int argc, char *argv[]) { /* Scan for variables */ printf("Scanning for variables...\n"); + if (file->file_type == FILE_TYPE_MITGCM) { + variables = mitgcm_scan_variables(file, mesh); + } else #ifdef HAVE_ZARR if (file->file_type == FILE_TYPE_ZARR) { variables = zarr_scan_variables(file, mesh); @@ -1238,6 +1268,9 @@ int main(int argc, char *argv[]) { int n_init_dims = 0; if (max_var) { /* Use fileset dimension info if available (includes virtual time) */ + if (file->file_type == FILE_TYPE_MITGCM) { + init_dims = mitgcm_get_dim_info(max_var, &n_init_dims); + } else #ifdef HAVE_ZARR if (zarr_fileset) { init_dims = zarr_get_dim_info_fileset(zarr_fileset, max_var, &n_init_dims); diff --git a/src/ushow.defines.h b/src/ushow.defines.h index af92e1e..fd1b7f0 100644 --- a/src/ushow.defines.h +++ b/src/ushow.defines.h @@ -45,6 +45,7 @@ typedef enum { typedef enum { FILE_TYPE_UNKNOWN = 0, FILE_TYPE_NETCDF, + FILE_TYPE_MITGCM, #ifdef HAVE_ZARR FILE_TYPE_ZARR, #endif @@ -143,6 +144,9 @@ struct USVar { USFile *file; int varid; /* NetCDF variable ID */ + /* MITgcm-specific */ + void *mitgcm_data; /* MitgcmVarData* for MITgcm variables */ + #ifdef HAVE_ZARR /* Zarr-specific */ void *zarr_data; /* ZarrArray* for zarr variables */ @@ -166,6 +170,9 @@ struct USFile { /* NetCDF-specific */ int ncid; /* NetCDF file ID */ + /* MITgcm-specific */ + void *mitgcm_data; /* MitgcmStore* for MITgcm files */ + #ifdef HAVE_ZARR /* Zarr-specific */ void *zarr_data; /* ZarrStore* for zarr files */ diff --git a/src/view.c b/src/view.c index 6856135..e6bceb3 100644 --- a/src/view.c +++ b/src/view.c @@ -4,6 +4,7 @@ #include "view.h" #include "file_netcdf.h" +#include "file_mitgcm.h" #ifdef HAVE_ZARR #include "file_zarr.h" #endif @@ -137,6 +138,9 @@ int view_set_variable(USView *view, USVar *var, USMesh *mesh, USRegrid *regrid) /* Estimate data range if not set */ if (!var->range_set) { + if (var->file && var->file->file_type == FILE_TYPE_MITGCM) { + mitgcm_estimate_range(var, &var->global_min, &var->global_max); + } else #ifdef HAVE_ZARR if (var->file && var->file->file_type == FILE_TYPE_ZARR) { zarr_estimate_range(var, &var->global_min, &var->global_max); @@ -432,6 +436,10 @@ int view_update(USView *view) { /* Read data slice - dispatch based on file type */ int read_result; + if (view->variable->file && view->variable->file->file_type == FILE_TYPE_MITGCM) { + read_result = mitgcm_read_slice(view->variable, view->time_index, + view->depth_index, view->raw_data); + } else #ifdef HAVE_ZARR if (view->fileset && view->fileset->files[0]->file_type == FILE_TYPE_ZARR) { /* Zarr multi-file */ diff --git a/tests/Makefile b/tests/Makefile index dd62339..23ff7e8 100644 --- a/tests/Makefile +++ b/tests/Makefile @@ -87,7 +87,7 @@ endif SRCDIR = ../src # Test executables -TEST_TARGETS = test_kdtree test_mesh test_regrid test_colormaps test_file_netcdf test_integration test_term_render_mode test_range_popup test_timeseries test_latband +TEST_TARGETS = test_kdtree test_mesh test_regrid test_colormaps test_file_netcdf test_file_mitgcm test_integration test_term_render_mode test_range_popup test_timeseries test_latband # Add zarr test if enabled ifdef WITH_ZARR @@ -110,6 +110,7 @@ MESH_OBJ = $(SRCDIR)/mesh.c REGRID_OBJ = $(SRCDIR)/regrid.c COLORMAPS_OBJ = $(SRCDIR)/colormaps.c FILE_NETCDF_OBJ = $(SRCDIR)/file_netcdf.c +FILE_MITGCM_OBJ = $(SRCDIR)/file_mitgcm.c FILE_GRIB_OBJ = $(SRCDIR)/file_grib.c # Zarr support files (only when WITH_ZARR=1) @@ -143,6 +144,9 @@ test_colormaps: test_colormaps.c $(COLORMAPS_OBJ) test_file_netcdf: test_file_netcdf.c $(FILE_NETCDF_OBJ) $(MESH_DEPS) $(KDTREE_OBJ) $(CC) $(CFLAGS) -o $@ $^ $(LIBS) +test_file_mitgcm: test_file_mitgcm.c $(FILE_MITGCM_OBJ) $(MESH_DEPS) $(KDTREE_OBJ) + $(CC) $(CFLAGS) -o $@ $^ $(LIBS) + test_integration: test_integration.c $(FILE_NETCDF_OBJ) $(MESH_DEPS) $(KDTREE_OBJ) $(REGRID_OBJ) $(COLORMAPS_OBJ) $(CC) $(CFLAGS) -o $@ $^ $(LIBS) @@ -216,6 +220,9 @@ test-colormaps: test_colormaps test-netcdf: test_file_netcdf ./test_file_netcdf +test-mitgcm: test_file_mitgcm + ./test_file_mitgcm + test-integration: test_integration ./test_integration @@ -245,6 +252,7 @@ clean: rm -f $(TEST_TARGETS) test_file_zarr test_file_grib test_latband test_yac_3d rm -f /tmp/test_ushow_*.nc rm -rf /tmp/test_ushow_zarr_*.zarr + rm -rf /tmp/test_mitgcm_* # Verbose build for debugging verbose: CFLAGS += -v diff --git a/tests/test_file_mitgcm.c b/tests/test_file_mitgcm.c new file mode 100644 index 0000000..fb01205 --- /dev/null +++ b/tests/test_file_mitgcm.c @@ -0,0 +1,973 @@ +/* + * test_file_mitgcm.c - Unit tests for MITgcm MDS binary data reading + * + * Creates synthetic MITgcm .data/.meta file pairs in a temp directory, + * then exercises the file_mitgcm API. + */ + +#include "test_framework.h" +#include "../src/ushow.defines.h" +#include "../src/file_mitgcm.h" +#include "../src/mesh.h" +#include +#include +#include +#include +#include +#include + +/* ---- Helpers to create synthetic MITgcm data ---- */ + +static char test_dir[256]; + +static void swap_to_big_endian(float *data, size_t n) { + /* If the machine is little-endian, swap bytes to produce big-endian output */ + uint32_t check = 1; + if (*(unsigned char *)&check == 0) return; /* already big-endian */ + uint32_t *p = (uint32_t *)data; + for (size_t i = 0; i < n; i++) { + uint32_t v = p[i]; + p[i] = ((v >> 24) & 0xFF) | ((v >> 8) & 0xFF00) | + ((v << 8) & 0xFF0000) | ((v << 24) & 0xFF000000); + } +} + +static int write_binary(const char *path, float *data, size_t n) { + float *buf = malloc(n * sizeof(float)); + if (!buf) return -1; + memcpy(buf, data, n * sizeof(float)); + swap_to_big_endian(buf, n); + FILE *fp = fopen(path, "wb"); + if (!fp) { free(buf); return -1; } + size_t written = fwrite(buf, sizeof(float), n, fp); + fclose(fp); + free(buf); + return (written == n) ? 0 : -1; +} + +static int write_meta(const char *path, const char *content) { + FILE *fp = fopen(path, "w"); + if (!fp) return -1; + fputs(content, fp); + fclose(fp); + return 0; +} + +/* Small grid: 4x3 with 2 depth levels */ +#define TEST_NX 4 +#define TEST_NY 3 +#define TEST_NZ 2 +#define TEST_SLAB (TEST_NX * TEST_NY) + +static void setup_test_dir(void) { + snprintf(test_dir, sizeof(test_dir), "/tmp/test_mitgcm_%d", getpid()); + mkdir(test_dir, 0755); +} + +static void cleanup_test_dir(void) { + /* Remove all files in the test directory */ + char cmd[512]; + snprintf(cmd, sizeof(cmd), "rm -rf %s", test_dir); + (void)system(cmd); +} + +/* Create XC.meta + XC.data */ +static int create_grid_xy(void) { + char path[512]; + + /* XC.meta */ + snprintf(path, sizeof(path), "%s/XC.meta", test_dir); + write_meta(path, + " nDims = [ 2 ];\n" + " dimList = [\n" + " 4, 1, 4,\n" + " 3, 1, 3\n" + " ];\n" + " dataprec = [ 'float32' ];\n" + " nrecords = [ 1 ];\n"); + + /* XC.data: lon values */ + float xc[TEST_SLAB]; + for (int j = 0; j < TEST_NY; j++) + for (int i = 0; i < TEST_NX; i++) + xc[j * TEST_NX + i] = -180.0f + 90.0f * i; /* -180, -90, 0, 90 */ + + snprintf(path, sizeof(path), "%s/XC.data", test_dir); + if (write_binary(path, xc, TEST_SLAB) != 0) return -1; + + /* YC.meta */ + snprintf(path, sizeof(path), "%s/YC.meta", test_dir); + write_meta(path, + " nDims = [ 2 ];\n" + " dimList = [\n" + " 4, 1, 4,\n" + " 3, 1, 3\n" + " ];\n" + " dataprec = [ 'float32' ];\n" + " nrecords = [ 1 ];\n"); + + /* YC.data: lat values */ + float yc[TEST_SLAB]; + for (int j = 0; j < TEST_NY; j++) + for (int i = 0; i < TEST_NX; i++) + yc[j * TEST_NX + i] = -60.0f + 60.0f * j; /* -60, 0, 60 */ + + snprintf(path, sizeof(path), "%s/YC.data", test_dir); + return write_binary(path, yc, TEST_SLAB); +} + +/* Create RC.meta + RC.data (depth levels) */ +static int create_rc(void) { + char path[512]; + + snprintf(path, sizeof(path), "%s/RC.meta", test_dir); + write_meta(path, + " nDims = [ 3 ];\n" + " dimList = [\n" + " 1, 1, 1,\n" + " 1, 1, 1,\n" + " 2, 1, 2\n" + " ];\n" + " dataprec = [ 'float32' ];\n" + " nrecords = [ 1 ];\n"); + + float rc[TEST_NZ] = {-5.0f, -15.0f}; + snprintf(path, sizeof(path), "%s/RC.data", test_dir); + return write_binary(path, rc, TEST_NZ); +} + +/* Create hFacC.meta + hFacC.data (land mask) */ +static int create_hfac(void) { + char path[512]; + + snprintf(path, sizeof(path), "%s/hFacC.meta", test_dir); + write_meta(path, + " nDims = [ 3 ];\n" + " dimList = [\n" + " 4, 1, 4,\n" + " 3, 1, 3,\n" + " 2, 1, 2\n" + " ];\n" + " dataprec = [ 'float32' ];\n" + " nrecords = [ 1 ];\n" + " timeStepNumber = [ 0 ];\n"); + + /* hFacC: all ocean except point [0,0] at both depths */ + size_t total = TEST_NZ * TEST_SLAB; + float *hfac = calloc(total, sizeof(float)); + if (!hfac) return -1; + for (size_t i = 0; i < total; i++) hfac[i] = 1.0f; + hfac[0] = 0.0f; /* depth 0, point (0,0) = land */ + hfac[TEST_SLAB + 0] = 0.0f; /* depth 1, point (0,0) = land */ + + snprintf(path, sizeof(path), "%s/hFacC.data", test_dir); + int rc = write_binary(path, hfac, total); + free(hfac); + return rc; +} + +/* Create a 2D diagnostic file with 2 fields at a given iteration */ +static int create_diag2d(int iteration) { + char path[512]; + char meta_buf[1024]; + + snprintf(path, sizeof(path), "%s/diags2D.%010d.meta", test_dir, iteration); + snprintf(meta_buf, sizeof(meta_buf), + " nDims = [ 2 ];\n" + " dimList = [\n" + " 4, 1, 4,\n" + " 3, 1, 3\n" + " ];\n" + " dataprec = [ 'float32' ];\n" + " nrecords = [ 2 ];\n" + " timeStepNumber = [ %d ];\n" + " missingValue = [ -9.99000000000000E+02 ];\n" + " nFlds = [ 2 ];\n" + " fldList = {\n" + " 'ETAN ' 'MXLDEPTH'\n" + " };\n", iteration); + write_meta(path, meta_buf); + + /* Data: 2 fields * ny * nx */ + size_t total = 2 * TEST_SLAB; + float *data = malloc(total * sizeof(float)); + if (!data) return -1; + + /* ETAN: simple ramp */ + for (int i = 0; i < TEST_SLAB; i++) + data[i] = 0.5f * (float)i + (float)iteration * 0.01f; + + /* MXLDEPTH: constant */ + for (int i = 0; i < TEST_SLAB; i++) + data[TEST_SLAB + i] = 50.0f + (float)iteration * 0.1f; + + snprintf(path, sizeof(path), "%s/diags2D.%010d.data", test_dir, iteration); + int rc = write_binary(path, data, total); + free(data); + return rc; +} + +/* Create a 3D diagnostic file with 1 field at a given iteration */ +static int create_diag3d(int iteration) { + char path[512]; + char meta_buf[1024]; + + snprintf(path, sizeof(path), "%s/diags3D.%010d.meta", test_dir, iteration); + snprintf(meta_buf, sizeof(meta_buf), + " nDims = [ 3 ];\n" + " dimList = [\n" + " 4, 1, 4,\n" + " 3, 1, 3,\n" + " 2, 1, 2\n" + " ];\n" + " dataprec = [ 'float32' ];\n" + " nrecords = [ 1 ];\n" + " timeStepNumber = [ %d ];\n" + " missingValue = [ -9.99000000000000E+02 ];\n" + " nFlds = [ 1 ];\n" + " fldList = {\n" + " 'THETA '\n" + " };\n", iteration); + write_meta(path, meta_buf); + + /* Data: nz * ny * nx */ + size_t total = TEST_NZ * TEST_SLAB; + float *data = malloc(total * sizeof(float)); + if (!data) return -1; + + for (int z = 0; z < TEST_NZ; z++) + for (int i = 0; i < TEST_SLAB; i++) + data[z * TEST_SLAB + i] = 20.0f - (float)z * 5.0f + (float)i * 0.1f; + + snprintf(path, sizeof(path), "%s/diags3D.%010d.data", test_dir, iteration); + int rc = write_binary(path, data, total); + free(data); + return rc; +} + +/* Create a full test dataset */ +static int create_full_dataset(void) { + setup_test_dir(); + if (create_grid_xy() != 0) return -1; + if (create_rc() != 0) return -1; + if (create_hfac() != 0) return -1; + if (create_diag2d(100) != 0) return -1; + if (create_diag2d(200) != 0) return -1; + if (create_diag3d(100) != 0) return -1; + if (create_diag3d(200) != 0) return -1; + return 0; +} + +/* ---- Tests ---- */ + +/* Detection */ + +TEST(mitgcm_is_mitgcm_null) { + ASSERT_FALSE(mitgcm_is_mitgcm(NULL)); + return 1; +} + +TEST(mitgcm_is_mitgcm_nonexistent) { + ASSERT_FALSE(mitgcm_is_mitgcm("/nonexistent/path")); + return 1; +} + +TEST(mitgcm_is_mitgcm_directory) { + ASSERT_EQ(create_full_dataset(), 0); + ASSERT_TRUE(mitgcm_is_mitgcm(test_dir)); + cleanup_test_dir(); + return 1; +} + +TEST(mitgcm_is_mitgcm_data_file) { + ASSERT_EQ(create_full_dataset(), 0); + char path[512]; + snprintf(path, sizeof(path), "%s/diags2D.0000000100.data", test_dir); + ASSERT_TRUE(mitgcm_is_mitgcm(path)); + cleanup_test_dir(); + return 1; +} + +/* Open */ + +TEST(mitgcm_open_null) { + USFile *f = mitgcm_open(NULL); + ASSERT_NULL(f); + return 1; +} + +TEST(mitgcm_open_directory) { + ASSERT_EQ(create_full_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + ASSERT_EQ(f->file_type, FILE_TYPE_MITGCM); + ASSERT_NOT_NULL(f->mitgcm_data); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +TEST(mitgcm_open_from_data_file) { + ASSERT_EQ(create_full_dataset(), 0); + char path[512]; + snprintf(path, sizeof(path), "%s/diags3D.0000000100.data", test_dir); + USFile *f = mitgcm_open(path); + ASSERT_NOT_NULL(f); + ASSERT_EQ(f->file_type, FILE_TYPE_MITGCM); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +/* Mesh */ + +TEST(mitgcm_create_mesh_basic) { + ASSERT_EQ(create_full_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + ASSERT_EQ_SIZET(m->n_points, TEST_SLAB); + ASSERT_EQ(m->coord_type, COORD_TYPE_2D_CURVILINEAR); + ASSERT_EQ_SIZET(m->orig_nx, TEST_NX); + ASSERT_EQ_SIZET(m->orig_ny, TEST_NY); + + /* Check some coordinate values */ + ASSERT_NEAR(m->lon[0], -180.0, 0.1); /* First point lon */ + ASSERT_NEAR(m->lat[0], -60.0, 0.1); /* First point lat */ + + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +/* Variable scanning */ + +TEST(mitgcm_scan_variables_count) { + ASSERT_EQ(create_full_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + /* Should find: ETAN, MXLDEPTH (2D), THETA (3D) = 3 variables */ + int count = 0; + for (USVar *v = vars; v; v = v->next) count++; + ASSERT_EQ_INT(count, 3); + + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +TEST(mitgcm_scan_variables_2d_dims) { + ASSERT_EQ(create_full_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + /* Find ETAN */ + USVar *etan = NULL; + for (USVar *v = vars; v; v = v->next) { + if (strcmp(v->name, "ETAN") == 0) { etan = v; break; } + } + ASSERT_NOT_NULL(etan); + + /* 2D var: time + node */ + ASSERT_TRUE(etan->time_dim_id >= 0); + ASSERT_EQ_INT(etan->depth_dim_id, -1); + ASSERT_EQ_SIZET(etan->dim_sizes[etan->time_dim_id], 2); /* 2 iterations */ + ASSERT_EQ_SIZET(etan->dim_sizes[etan->node_dim_id], TEST_SLAB); + + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +TEST(mitgcm_scan_variables_3d_dims) { + ASSERT_EQ(create_full_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + /* Find THETA */ + USVar *theta = NULL; + for (USVar *v = vars; v; v = v->next) { + if (strcmp(v->name, "THETA") == 0) { theta = v; break; } + } + ASSERT_NOT_NULL(theta); + + /* 3D var: time + depth + node */ + ASSERT_TRUE(theta->time_dim_id >= 0); + ASSERT_TRUE(theta->depth_dim_id >= 0); + ASSERT_EQ_SIZET(theta->dim_sizes[theta->time_dim_id], 2); + ASSERT_EQ_SIZET(theta->dim_sizes[theta->depth_dim_id], TEST_NZ); + ASSERT_EQ_SIZET(theta->dim_sizes[theta->node_dim_id], TEST_SLAB); + + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +/* Read slice */ + +TEST(mitgcm_read_slice_2d) { + ASSERT_EQ(create_full_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + /* Find ETAN */ + USVar *etan = NULL; + for (USVar *v = vars; v; v = v->next) { + if (strcmp(v->name, "ETAN") == 0) { etan = v; break; } + } + ASSERT_NOT_NULL(etan); + + float data[TEST_SLAB]; + int rc = mitgcm_read_slice(etan, 0, 0, data); + ASSERT_EQ_INT(rc, 0); + + /* Point 0 is land (hFacC=0), should be fill value */ + ASSERT_NEAR(data[0], etan->fill_value, 0.1); + + /* Point 1 is ocean, should have real data */ + /* Expected: 0.5 * 1 + 100 * 0.01 = 1.5 */ + ASSERT_NEAR(data[1], 1.5f, 0.01); + + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +TEST(mitgcm_read_slice_3d_surface) { + ASSERT_EQ(create_full_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + /* Find THETA */ + USVar *theta = NULL; + for (USVar *v = vars; v; v = v->next) { + if (strcmp(v->name, "THETA") == 0) { theta = v; break; } + } + ASSERT_NOT_NULL(theta); + + float data[TEST_SLAB]; + int rc = mitgcm_read_slice(theta, 0, 0, data); + ASSERT_EQ_INT(rc, 0); + + /* Point 0 (land) should be fill */ + ASSERT_NEAR(data[0], theta->fill_value, 0.1); + + /* Point 1 (ocean, depth=0): 20.0 - 0*5.0 + 1*0.1 = 20.1 */ + ASSERT_NEAR(data[1], 20.1f, 0.01); + + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +TEST(mitgcm_read_slice_3d_deep) { + ASSERT_EQ(create_full_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + USVar *theta = NULL; + for (USVar *v = vars; v; v = v->next) { + if (strcmp(v->name, "THETA") == 0) { theta = v; break; } + } + ASSERT_NOT_NULL(theta); + + float data[TEST_SLAB]; + int rc = mitgcm_read_slice(theta, 0, 1, data); + ASSERT_EQ_INT(rc, 0); + + /* Point 1 (ocean, depth=1): 20.0 - 1*5.0 + 1*0.1 = 15.1 */ + ASSERT_NEAR(data[1], 15.1f, 0.01); + + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +TEST(mitgcm_read_slice_second_timestep) { + ASSERT_EQ(create_full_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + USVar *etan = NULL; + for (USVar *v = vars; v; v = v->next) { + if (strcmp(v->name, "ETAN") == 0) { etan = v; break; } + } + ASSERT_NOT_NULL(etan); + + float data[TEST_SLAB]; + int rc = mitgcm_read_slice(etan, 1, 0, data); + ASSERT_EQ_INT(rc, 0); + + /* Point 1 (ocean, iteration=200): 0.5 * 1 + 200 * 0.01 = 2.5 */ + ASSERT_NEAR(data[1], 2.5f, 0.01); + + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +/* Range estimation */ + +TEST(mitgcm_estimate_range_basic) { + ASSERT_EQ(create_full_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + USVar *etan = NULL; + for (USVar *v = vars; v; v = v->next) { + if (strcmp(v->name, "ETAN") == 0) { etan = v; break; } + } + ASSERT_NOT_NULL(etan); + + float min_val, max_val; + int rc = mitgcm_estimate_range(etan, &min_val, &max_val); + ASSERT_EQ_INT(rc, 0); + ASSERT_LT(min_val, max_val); + ASSERT_GT(min_val, -100.0f); /* Should not be fill value */ + + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +/* Dimension info */ + +TEST(mitgcm_get_dim_info_2d) { + ASSERT_EQ(create_full_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + USVar *etan = NULL; + for (USVar *v = vars; v; v = v->next) { + if (strcmp(v->name, "ETAN") == 0) { etan = v; break; } + } + ASSERT_NOT_NULL(etan); + + int n_dims = 0; + USDimInfo *dims = mitgcm_get_dim_info(etan, &n_dims); + ASSERT_NOT_NULL(dims); + ASSERT_EQ_INT(n_dims, 1); /* time only, no depth */ + ASSERT_STR_EQ(dims[0].name, "time"); + ASSERT_EQ_SIZET(dims[0].size, 2); + ASSERT_NOT_NULL(dims[0].values); + ASSERT_NEAR(dims[0].values[0], 100.0, 0.01); + ASSERT_NEAR(dims[0].values[1], 200.0, 0.01); + + mitgcm_free_dim_info(dims, n_dims); + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +TEST(mitgcm_get_dim_info_3d) { + ASSERT_EQ(create_full_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + USVar *theta = NULL; + for (USVar *v = vars; v; v = v->next) { + if (strcmp(v->name, "THETA") == 0) { theta = v; break; } + } + ASSERT_NOT_NULL(theta); + + int n_dims = 0; + USDimInfo *dims = mitgcm_get_dim_info(theta, &n_dims); + ASSERT_NOT_NULL(dims); + ASSERT_EQ_INT(n_dims, 2); /* time + depth */ + ASSERT_STR_EQ(dims[0].name, "time"); + ASSERT_STR_EQ(dims[1].name, "depth"); + ASSERT_EQ_SIZET(dims[1].size, TEST_NZ); + ASSERT_NOT_NULL(dims[1].values); + ASSERT_NEAR(dims[1].values[0], -5.0, 0.01); + ASSERT_NEAR(dims[1].values[1], -15.0, 0.01); + + mitgcm_free_dim_info(dims, n_dims); + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +/* Time series */ + +TEST(mitgcm_read_timeseries_basic) { + ASSERT_EQ(create_full_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + USVar *etan = NULL; + for (USVar *v = vars; v; v = v->next) { + if (strcmp(v->name, "ETAN") == 0) { etan = v; break; } + } + ASSERT_NOT_NULL(etan); + + double *times = NULL; + float *values = NULL; + int *valid = NULL; + size_t n_out = 0; + + /* Read timeseries at ocean point (index 1) */ + int rc = mitgcm_read_timeseries(etan, 1, 0, ×, &values, &valid, &n_out); + ASSERT_EQ_INT(rc, 0); + ASSERT_EQ_SIZET(n_out, 2); + ASSERT_NOT_NULL(times); + ASSERT_NOT_NULL(values); + ASSERT_NOT_NULL(valid); + + ASSERT_NEAR(times[0], 100.0, 0.01); + ASSERT_NEAR(times[1], 200.0, 0.01); + ASSERT_EQ_INT(valid[0], 1); + ASSERT_EQ_INT(valid[1], 1); + + /* Values: iteration 100 -> 0.5*1 + 100*0.01 = 1.5, iteration 200 -> 0.5*1 + 200*0.01 = 2.5 */ + ASSERT_NEAR(values[0], 1.5f, 0.01); + ASSERT_NEAR(values[1], 2.5f, 0.01); + + free(times); free(values); free(valid); + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +TEST(mitgcm_read_timeseries_land_point) { + ASSERT_EQ(create_full_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + USVar *etan = NULL; + for (USVar *v = vars; v; v = v->next) { + if (strcmp(v->name, "ETAN") == 0) { etan = v; break; } + } + ASSERT_NOT_NULL(etan); + + double *times = NULL; + float *values = NULL; + int *valid = NULL; + size_t n_out = 0; + + /* Read timeseries at land point (index 0) */ + int rc = mitgcm_read_timeseries(etan, 0, 0, ×, &values, &valid, &n_out); + ASSERT_EQ_INT(rc, 0); + ASSERT_EQ_SIZET(n_out, 2); + ASSERT_EQ_INT(valid[0], 0); /* Land: not valid */ + ASSERT_EQ_INT(valid[1], 0); + + free(times); free(values); free(valid); + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +/* ---- Velocity rotation helpers ---- */ + +/* Create AngleCS.data + AngleSN.data with constant angle */ +static int create_angle_cs_sn(float cs_val, float sn_val) { + char path[512]; + float cs[TEST_SLAB], sn[TEST_SLAB]; + + for (int i = 0; i < TEST_SLAB; i++) { + cs[i] = cs_val; + sn[i] = sn_val; + } + + snprintf(path, sizeof(path), "%s/AngleCS.data", test_dir); + if (write_binary(path, cs, TEST_SLAB) != 0) return -1; + + snprintf(path, sizeof(path), "%s/AngleSN.data", test_dir); + return write_binary(path, sn, TEST_SLAB); +} + +/* Create a 3D velocity diagnostic with UVEL and VVEL fields */ +static int create_vel_diag(int iteration, float uval, float vval) { + char path[512]; + char meta_buf[1024]; + + snprintf(path, sizeof(path), "%s/velDiag.%010d.meta", test_dir, iteration); + snprintf(meta_buf, sizeof(meta_buf), + " nDims = [ 3 ];\n" + " dimList = [\n" + " 4, 1, 4,\n" + " 3, 1, 3,\n" + " 2, 1, 2\n" + " ];\n" + " dataprec = [ 'float32' ];\n" + " nrecords = [ 2 ];\n" + " timeStepNumber = [ %d ];\n" + " nFlds = [ 2 ];\n" + " fldList = {\n" + " 'UVEL ' 'VVEL '\n" + " };\n", iteration); + write_meta(path, meta_buf); + + /* Data: 2 fields * nz * ny * nx */ + size_t field_size = TEST_NZ * TEST_SLAB; + size_t total = 2 * field_size; + float *data = malloc(total * sizeof(float)); + if (!data) return -1; + + /* UVEL: constant value everywhere */ + for (size_t i = 0; i < field_size; i++) + data[i] = uval; + + /* VVEL: constant value everywhere */ + for (size_t i = 0; i < field_size; i++) + data[field_size + i] = vval; + + snprintf(path, sizeof(path), "%s/velDiag.%010d.data", test_dir, iteration); + int rc = write_binary(path, data, total); + free(data); + return rc; +} + +/* Create dataset with angle files and velocity diagnostics */ +static int create_rotation_dataset(void) { + setup_test_dir(); + if (create_grid_xy() != 0) return -1; + if (create_rc() != 0) return -1; + if (create_hfac() != 0) return -1; + /* CS=0.6, SN=0.8 (3-4-5 triangle, CS^2+SN^2=1) */ + if (create_angle_cs_sn(0.6f, 0.8f) != 0) return -1; + /* UVEL=3.0, VVEL=4.0 */ + if (create_vel_diag(100, 3.0f, 4.0f) != 0) return -1; + if (create_vel_diag(200, 6.0f, 8.0f) != 0) return -1; + return 0; +} + +/* Velocity rotation tests */ + +TEST(mitgcm_velocity_pairing) { + ASSERT_EQ(create_rotation_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + /* Find UVEL and VVEL */ + USVar *uvel = NULL, *vvel = NULL; + for (USVar *v = vars; v; v = v->next) { + if (strcmp(v->name, "UVEL") == 0) uvel = v; + if (strcmp(v->name, "VVEL") == 0) vvel = v; + } + ASSERT_NOT_NULL(uvel); + ASSERT_NOT_NULL(vvel); + + /* Check they are paired */ + ASSERT_NOT_NULL(uvel->mitgcm_data); + ASSERT_NOT_NULL(vvel->mitgcm_data); + + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +TEST(mitgcm_velocity_rotation_slice_u) { + ASSERT_EQ(create_rotation_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + USVar *uvel = NULL; + for (USVar *v = vars; v; v = v->next) { + if (strcmp(v->name, "UVEL") == 0) { uvel = v; break; } + } + ASSERT_NOT_NULL(uvel); + + float data[TEST_SLAB]; + int rc = mitgcm_read_slice(uvel, 0, 0, data); + ASSERT_EQ_INT(rc, 0); + + /* Point 0 is land → fill value */ + ASSERT_NEAR(data[0], uvel->fill_value, 0.1); + + /* Ocean points: u_east = CS*U - SN*V = 0.6*3.0 - 0.8*4.0 = 1.8 - 3.2 = -1.4 */ + ASSERT_NEAR(data[1], -1.4f, 0.01); + ASSERT_NEAR(data[2], -1.4f, 0.01); + + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +TEST(mitgcm_velocity_rotation_slice_v) { + ASSERT_EQ(create_rotation_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + USVar *vvel = NULL; + for (USVar *v = vars; v; v = v->next) { + if (strcmp(v->name, "VVEL") == 0) { vvel = v; break; } + } + ASSERT_NOT_NULL(vvel); + + float data[TEST_SLAB]; + int rc = mitgcm_read_slice(vvel, 0, 0, data); + ASSERT_EQ_INT(rc, 0); + + /* Ocean points: v_north = SN*U + CS*V = 0.8*3.0 + 0.6*4.0 = 2.4 + 2.4 = 4.8 */ + ASSERT_NEAR(data[1], 4.8f, 0.01); + ASSERT_NEAR(data[2], 4.8f, 0.01); + + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +TEST(mitgcm_velocity_rotation_timeseries) { + ASSERT_EQ(create_rotation_dataset(), 0); + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + USVar *uvel = NULL; + for (USVar *v = vars; v; v = v->next) { + if (strcmp(v->name, "UVEL") == 0) { uvel = v; break; } + } + ASSERT_NOT_NULL(uvel); + + double *times = NULL; + float *values = NULL; + int *valid = NULL; + size_t n_out = 0; + + /* Read timeseries at ocean point (index 1) */ + int rc = mitgcm_read_timeseries(uvel, 1, 0, ×, &values, &valid, &n_out); + ASSERT_EQ_INT(rc, 0); + ASSERT_EQ_SIZET(n_out, 2); + ASSERT_EQ_INT(valid[0], 1); + ASSERT_EQ_INT(valid[1], 1); + + /* t=100: u_east = 0.6*3.0 - 0.8*4.0 = -1.4 */ + ASSERT_NEAR(values[0], -1.4f, 0.01); + /* t=200: u_east = 0.6*6.0 - 0.8*8.0 = 3.6 - 6.4 = -2.8 */ + ASSERT_NEAR(values[1], -2.8f, 0.01); + + free(times); free(values); free(valid); + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +TEST(mitgcm_no_rotation_without_angles) { + /* Create dataset WITHOUT angle files */ + setup_test_dir(); + ASSERT_EQ(create_grid_xy(), 0); + ASSERT_EQ(create_rc(), 0); + ASSERT_EQ(create_hfac(), 0); + ASSERT_EQ(create_vel_diag(100, 3.0f, 4.0f), 0); + + USFile *f = mitgcm_open(test_dir); + ASSERT_NOT_NULL(f); + USMesh *m = mitgcm_create_mesh(f); + ASSERT_NOT_NULL(m); + USVar *vars = mitgcm_scan_variables(f, m); + ASSERT_NOT_NULL(vars); + + USVar *uvel = NULL; + for (USVar *v = vars; v; v = v->next) { + if (strcmp(v->name, "UVEL") == 0) { uvel = v; break; } + } + ASSERT_NOT_NULL(uvel); + + float data[TEST_SLAB]; + int rc = mitgcm_read_slice(uvel, 0, 0, data); + ASSERT_EQ_INT(rc, 0); + + /* Without angles: raw UVEL=3.0 at ocean points (no rotation) */ + ASSERT_NEAR(data[1], 3.0f, 0.01); + + mesh_free(m); + mitgcm_close(f); + cleanup_test_dir(); + return 1; +} + +/* Close null */ + +TEST(mitgcm_close_null) { + mitgcm_close(NULL); /* Should not crash */ + return 1; +} + +RUN_TESTS("MITgcm File I/O")