From 2380d502f101c2bec51894939134fcae82da960f Mon Sep 17 00:00:00 2001 From: Joaquim Date: Thu, 25 Jun 2026 18:28:41 +0100 Subject: [PATCH 01/29] First cut of the conversion of the Mirone NSWING code to a GMT supplement --- src/CMakeLists.txt | 2 +- src/nswing/CMakeLists.txt | 32 + src/nswing/README.nswing | 23 + src/nswing/nswing.c | 5790 +++++++++++++++++++++++++++++++++++++ 4 files changed, 5846 insertions(+), 1 deletion(-) create mode 100644 src/nswing/CMakeLists.txt create mode 100644 src/nswing/README.nswing create mode 100644 src/nswing/nswing.c diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a209666b74f..9db401cdb23 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -712,7 +712,7 @@ if (BUILD_SUPPLEMENTS) # supplement directories to be built as supplement libraries # SUPPL_EXTRA_DIRS are directories for custom supplements. # See cmake/ConfigUserAdvancedTemplate.cmake for details. - set (GMT_SUPPL_DIRS geodesy gsfml gshhg img mgd77 potential segy seis spotter x2sys windbarbs ${SUPPL_EXTRA_DIRS}) + set (GMT_SUPPL_DIRS geodesy gsfml gshhg img mgd77 potential segy seis spotter x2sys windbarbs nswing ${SUPPL_EXTRA_DIRS}) # process supplement directories and collect information for each supplement libraries set (GMT_SUPPL_LIBRARIES) # empty the library list diff --git a/src/nswing/CMakeLists.txt b/src/nswing/CMakeLists.txt new file mode 100644 index 00000000000..a73ee660072 --- /dev/null +++ b/src/nswing/CMakeLists.txt @@ -0,0 +1,32 @@ +# +# +# Copyright (c) 1991-2020 by the GMT Team (https://www.generic-mapping-tools.org/team.html) +# See LICENSE.TXT file for copying and redistribution conditions. +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation; version 3 or any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# Contact info: www.generic-mapping-tools.org +#------------------------------------------------------------------------------- +# +# CMake settings for supplemental package: +# +# SUPPL_NAME: name of the supplemental package +# SUPPL_HEADERS: header files. Will be installed if BUILD_SUPPLEMENTS is TRUE +# SUPPL_PROGS_SRCS: list of C source codes for supplemental modules +# SUPPL_LIB_SRCS: list of C source codes for supplemental library +# SUPPL_EXAMPLE_FILES: README and other example files +# SUPPL_EXAMPLE_PROGS: Example scripts +# + +set (SUPPL_NAME nswing) +#set (SUPPL_LIB_NAME nswing) +set (SUPPL_PROGS_SRCS nswing.c) +set (SUPPL_LIB_SRCS ${SUPPL_PROGS_SRCS} nswing.c) +set (SUPPL_EXAMPLE_FILES README.nswing) diff --git a/src/nswing/README.nswing b/src/nswing/README.nswing new file mode 100644 index 00000000000..97563070ea5 --- /dev/null +++ b/src/nswing/README.nswing @@ -0,0 +1,23 @@ + +NSWING Supplemental Package +Distributed under the MIT License; see file +LICENSE.TXT in main GMT directory. + +------------------------------------------------------ +Author: Joaquim Luis, Miguel Miranda +Date: June 1, 2026 +Version: 6 (for GMT 6 release). +------------------------------------------------------ + +This directory contains: + + README.nswing This document + nswing.c nswing C program code + + +Documentation via RST files produce HTML and MAN pages that +give full description of the program interfaces. + + +* Acknowledgement + diff --git a/src/nswing/nswing.c b/src/nswing/nswing.c new file mode 100644 index 00000000000..d5bc2bd7ffc --- /dev/null +++ b/src/nswing/nswing.c @@ -0,0 +1,5790 @@ +/*-------------------------------------------------------------------- + * $Id: nswing.c 11410 2019-02-23 17:36:11Z j $ + * + * Copyright (c) 2012-2022 by J. Luis and J. M. Miranda + * + * This program is part of Mirone and is free software; you can redistribute + * it and/or modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * Contact info: w3.ualg.pt/~jluis/mirone + *--------------------------------------------------------------------*/ + + +/* + * Original Fortran version of core hydrodynamic code by J.M. Miranda and COMCOT + * + * + * The nesting grids algorithm + LOOP_OVER_ALL_CYCLES { + mass_L0 + openbd + LOOP_OVER_N_NESTED { + LOOP_TIME_REFINE_1 { + interp_edges_L1 + mass_L1 + LOOP_TIME_REFINE_2 { + interp_edges_L2 + mass_L2 + LOOP_TIME_REFINE_3 { + ... + } + moment_L2 + upscale_L2 + update_L2 + } + moment_L1 + upscale_L1 + update_L1 + } + } + moment_L0 + update_L0 + } + * + * + * To compile, do for example + * cl nswing.c -IC:\programs\compa_libs\netcdf_GIT\compileds\VC12_64\include + * C:\programs\compa_libs\netcdf_GIT\compileds\VC12_64\lib\netcdf.lib /DI_AM_C + * /DHAVE_NETCDF /nologo /D_CRT_SECURE_NO_WARNINGS /fp:precise /Ox + * + * Rewritten in C, added number options, etc... By + * Joaquim Luis - 2013 + * + */ + +#include "gmt_dev.h" /* GMT development API (grids, options, I/O) */ + +#define THIS_MODULE_CLASSIC_NAME "nswing" +#define THIS_MODULE_MODERN_NAME "nswing" +#define THIS_MODULE_LIB "supplements" +#define THIS_MODULE_PURPOSE "A tsunami maker" +#define THIS_MODULE_KEYS " +# define NAN _Nan._Double +# else + static const double _NAN = 20;//(HUGE_VAL-HUGE_VAL); +# define NAN _NAN +# endif +#endif + +/* At least the Intel compiler screws anf the NAN macro returns 0.0 So we need this hack */ +/* http://stackoverflow.com/questions/5714131/nan-literal-in-c */ +union {uint64_t i; double d;} loc_nan = {0x7ff8000000000000}; + +/* err_trap is NOT a netCDF function: it is a local macro wrapping nc_strerror(). It used + to be defined only under #ifdef HAVE_NETCDF, but GMT never defines that symbol, so the + macro vanished and every err_trap(...) call became an undefined external -> LNK2019. + GMT always links netCDF, so just define it unconditionally here. */ +#define err_trap(status) if (status) {fprintf(stderr, "NSWING: error at line %d: %s\n", __LINE__, nc_strerror(status));} + + +#if HAVE_OPENMP +#include +#endif + +#define DEPTH_THRESHOLD -40.0 /* When interpolating nested regions stop at this depth to avoid land contamination */ +#ifndef GMT_CONV_LIMIT +#define GMT_CONV_LIMIT 1.0e-8 /* Fairly tight convergence limit or "close to zero" limit */ +#endif +#define NORMAL_GRAV 9.806199203 /* Moritz's 1980 IGF value for gravity at 45 degrees latitude */ +#define EQ_RAD 6378137.0 /* WGS-84 */ +#define flattening 1.0/298.2572235630 +#define LIMIT_DISCHARGE /* If defined the moment functions will limit the discharge */ +#define ECC2 2 * flattening - flattening * flattening +#define ECC4 ECC2 * ECC2 +#define ECC6 ECC2 * ECC4 +#define EPS12 1e-12 +#define EPS10 1e-10 +#define EPS6 1e-6 +#define EPS5 1e-5 +#define EPS4_ 1e-4 +#define EPS3 1e-3 +#define EPS2 1e-2 +#define EPS1 1e-1 + +static double EPS4 = EPS4_; /* Kinda trick to be able to change EPS4 via a command line option */ + +#define MAXRUNUP -50 /* Do not waste time computing flood above this altitude */ +#define V_LIMIT 20 /* Upper limit of maximum velocity */ + +#define CNULL ((char *)NULL) +#define Loc_copysign(x,y) ((y) < 0.0 ? -fabs(x) : fabs(x)) + +#define ijs(i,j,n) ((i) + (j)*n) +#define ijc(i,j) ((i) + (j)*n_ptmar) +#define ij_grd(col,row,hdr) ((col) + (row)*hdr.nx) + +typedef void (*PFV) (); /* PFV declares a pointer to a function returning void */ + +struct tracers { /* For tracers (oranges) */ + double *x; /* x coordinate */ + double *y; /* y coordinate */ +}; + +struct srf_header { /* Surfer file hdr structure */ + char id[4]; /* ASCII Binary identifier (DSAA/DSBB) */ + short int nx; /* Number of columns */ + short int ny; /* Number of rows */ + double x_min; /* Minimum x coordinate */ + double x_max; /* Maximum x coordinate */ + double y_min; /* Minimum y coordinate */ + double y_max; /* Maximum y coordinate */ + double z_min; /* Minimum z value */ + double z_max; /* Maximum z value */ +}; + +struct grd_header { /* Generic grid hdr structure */ + int nx; /* Number of columns */ + int ny; /* Number of rows */ + unsigned int nm; /* nx * ny */ + double x_inc; + double y_inc; + double x_min; /* Minimum x coordinate */ + double x_max; /* Maximum x coordinate */ + double y_min; /* Minimum y coordinate */ + double y_max; /* Maximum y coordinate */ + double z_min; /* Minimum z value */ + double z_max; /* Maximum z value */ + double lat_min4Coriolis; /* PRECISA SOLUCAO. POR AGORA SERA Cte = 0 */ +}; + +struct nestContainer { /* Container for the nestings */ + int do_upscale; /* If false, do not upscale the parent grid */ + int do_long_beach; /* If true, compute a mask with ones over the "dryed beach" */ + int do_short_beach; /* If true, compute a mask with ones over the "innundated beach" */ + int do_linear; /* If true, use linear approximation */ + int do_max_level; /* If true, inform nestify() on the need to update max level at every inner iteration */ + int do_max_velocity; /* If true, inform nestify() on the need to update max velocity at each inner iteration */ + int do_Coriolis; /* If true, compute the Coriolis effect */ + int out_velocity_x; /* To know if we must compute the vex,vey velocity arrays */ + int out_velocity_y; + int out_momentum; /* To know if save the momentum in the 3D netCDF grid. Mutually exclusive with out_velocity_x|y */ + int isGeog; /* 0 == Cartesian, otherwise Geographic coordinates */ + int n_threads; /* Number of threads for multi-threading */ + int first, last, jupe; /* To hold start and ending rows in the moment conservation functions */ + int writeLevel; /* Store info about which level is (if) to be writen [0] */ + int bnc_pos_nPts; /* Number of points in a external boundary condition file */ + int bnc_var_nTimes; /* Number of time steps in the external boundary condition file */ + int bnc_border[4]; /* Each will be set to true if boundary condition on that border W->0, S->1, E->2, N->3 */ + int level[10]; /* 0 Will mean base level, others the nesting level */ + int LLrow[10], LLcol[10], ULrow[10], ULcol[10], URrow[10], URcol[10], LRrow[10], LRcol[10]; + int incRatio[10]; + short *long_beach[10]; /* Mask arrays for storing the "dry beaches" */ + short *short_beach[10]; /* Mask arrays for storing the "dry beaches" */ + float *work, *wmax; /* Auxiliary pointers (not direcly allocated) to compute max level of nested grids */ + float *vmax; /* Pointer to array storing the max velocity */ + double run_jump_time; /* Time to hold before letting the nested grids start to iterate */ + double lat_min4Coriolis; /* South latitute when computing the Coriolis effect on a cartesian grid */ + double manning_depth; /* Do not use manning if depth is deeper than this value */ + double manning[10]; /* Manning coefficient. Set to zero if no friction */ + double LLx[10], LLy[10], ULx[10], ULy[10], URx[10], URy[10], LRx[10], LRy[10]; + double dt[10]; /* Time step at current level */ + double *bat[10]; /* Bathymetry of current level */ + double *fluxm_a[10], *fluxm_d[10]; /* t-1/2 & t+1/2 fluxes arrays along X */ + double *fluxn_a[10], *fluxn_d[10]; /* t-1/2 & t+1/2 fluxes arrays along Y */ + double *htotal_a[10], *htotal_d[10]; /* t-1/2 & t+1/2 total water depth */ + double *vex[10], *vey[10]; /* X,Y velocity components */ + double *etaa[10], *etad[10]; /* t-1/2 & t+1/2 water height (eta) arrays */ + double *edge_col[10], *edge_colTmp[10]; + double *edge_row[10], *edge_rowTmp[10]; + double *edge_col_P[10], *edge_col_Ptmp[10]; + double *edge_row_P[10], *edge_row_Ptmp[10]; + double *r0[10], *r1m[10], *r1n[10], *r2m[10], *r2n[10], *r3m[10], *r3n[10], *r4m[10], *r4n[10]; + double time_h; + double *bnc_pos_x; + double *bnc_pos_y; + double *bnc_var_t; + double **bnc_var_z; + double *bnc_var_zTmp; + double *bnc_var_z_interp; + struct grd_header hdr[10]; +}; + +/* Argument struct for threading */ +typedef struct { + struct nestContainer *nest; /* Pointer to a nestContainer struct */ + int iThread; /* Thread index */ + int lev; /* Level of nested grid */ + int row_start, row_end; /* Start and end rows for the BY slices in PARALLEL mode */ +} ThreadArg; + +void no_sys_mem(char *where, unsigned int n); +int count_col(char *line); +int read_grd_info_ascii(char *file, struct srf_header *hdr); +int read_header_bin (FILE *fp, struct srf_header *hdr); +int write_grd_bin(char *name, double x_min, double y_min, double x_inc, double y_inc, unsigned int i_start, + unsigned int j_start, unsigned int i_end, unsigned int j_end, unsigned int nX, float *work); +int read_grd_ascii (char *file, struct srf_header *hdr, double *work, int sign); +int read_grd_bin(char *file, struct srf_header *hdr, double *work, int sign); +int read_maregs(struct grd_header hdr, char *file, unsigned int *lcum_p, char *names[]); +int read_tracers(struct grd_header hdr, char *file, struct tracers *oranges); +int count_n_maregs(char *file); +int decode_R(char *item, double *w, double *e, double *s, double *n); +int check_region(double w, double e, double s, double n); +double ddmmss_to_degree (char *text); +void openb(struct grd_header hdr, double *bat, double *fluxm_d, double *fluxn_d, double *etad, struct nestContainer *nest); +void wave_maker(struct nestContainer *nest); +void wall_it(struct nestContainer *nest); +void wall_two(struct nestContainer *nest, int ot1, int ot2, int in1, int in2); +int initialize_nestum(struct nestContainer *nest, int isGeog, int lev); +int intp_lin (double *x, double *y, int n, int m, double *u, double *v); +void inisp(struct nestContainer *nest); +void inicart(struct nestContainer *nest); +void interp_edges(struct nestContainer *nest, double *flux_L1, double *flux_L2, char *what, int lev, int i_time); +void sanitize_nestContainer(struct nestContainer *nest); +void nestify(struct nestContainer *nest, int nNg, int recursionLevel, int isGeog); +void resamplegrid(struct nestContainer *nest, int nNg); +void edge_communication(struct nestContainer *nest, int lev, int i_time); +void mass(struct nestContainer *nest, int lev); +void mass_sp(struct nestContainer *nest, int lev); +void mass_conservation(struct nestContainer *nest, int isGeog, int m); +void moment_conservation(struct nestContainer *nest, int isGeog, int m); +void update(struct nestContainer *nest, int lev); +void upscale(struct nestContainer *nest, double *out, int lev, int i_tsr); +void upscale_(struct nestContainer *nest, double *out, int lev, int i_tsr); +void replicate(struct nestContainer *nest, int lev); +void moment_M(struct nestContainer *nest, int lev); +void moment_N(struct nestContainer *nest, int lev); +void moment_M_slice(int lev, int row_start, int row_end, struct nestContainer *nest); +void moment_N_slice(int lev, int row_start, int row_end, struct nestContainer *nest); +void moment_sp_M(struct nestContainer *nest, int lev); +void moment_sp_N(struct nestContainer *nest, int lev); +void moment_sp_M_slice(int lev, int row_start, int row_end, struct nestContainer *nest); +void moment_sp_N_slice(int lev, int row_start, int row_end, struct nestContainer *nest); +void mass_sp_slice(struct nestContainer *nest, int lev, int row_start, int row_end); +void mass_slice(struct nestContainer *nest, int lev, int row_start, int row_end); +void free_arrays(struct nestContainer *nest, int isGeog, int lev); +int check_paternity(struct nestContainer *nest); +int check_binning(double x0P, double x0D, double dxP, double dxD, double tol, double *suggest); +int read_bnc_file(struct nestContainer *nest, char *file); +int interp_bnc (struct nestContainer *nest, double t); +void total_energy(struct nestContainer *nest, float *work, int lev); +void power(struct nestContainer *nest, float *work, int lev); +void vtm (double lat0, double *t_c1, double *t_c2, double *t_c3, double *t_c4, double *t_e2, double *t_M0); +void deform (struct srf_header hdr, double x_inc, double y_inc, int isGeog, double fault_length, + double fault_width, double th, double dip, double rake, double d, double top_depth, + double xl, double yl, double *z); +void kaba_source(struct srf_header hdr, double x_inc, double y_inc, double x_min, double x_max, + double y_min, double y_max, int type, double *z); +void tm (double lon, double lat, double *x, double *y, double central_meridian, double t_c1, + double t_c2, double t_c3, double t_c4, double t_e2, double t_M0); +double uscal(double x1, double x2, double x3, double c, double cc, double dp); +double udcal(double x1, double x2, double x3, double c, double cc, double dp); +unsigned int gmt_bcr_prep (struct grd_header hdr, double xx, double yy, double wx[], double wy[]); +double GMT_get_bcr_z(double *grd, struct grd_header hdr, double xx, double yy); +void update_max(struct nestContainer *nest); +void update_max_velocity(struct nestContainer *nest); + + +void write_most_slice(struct nestContainer *nest, int *ncid_most, int *ids_most, unsigned int i_start, + unsigned int j_start, unsigned int i_end, unsigned int j_end, float *work, size_t *start, + size_t *count, double *slice_range, int isMost, int lev); +int open_most_nc(struct nestContainer *nest, float *work, char *basename, char *name_var, char hist[], int *ids, + unsigned int nx, unsigned int ny, double xMinOut, double yMinOut, int isMost, int lev); +int open_anuga_sww(struct nestContainer *nest, char *fname_sww, char history[], int *ids, unsigned int i_start, + unsigned int j_start, unsigned int i_end, unsigned int j_end, double xMinOut, double yMinOut, int lev); +void write_anuga_slice(struct nestContainer *nest, int ncid, int z_id, unsigned int i_start, unsigned int j_start, + unsigned int i_end, unsigned int j_end, float *work, size_t *start, size_t *count, + float *slice_range, int idx, int with_land, int lev); +int write_maregs_nc(struct nestContainer *nest, char *fname, float *work, double *t, unsigned int *lcum_p, + char *names[], char hist[], int n_maregs, unsigned int count_time_maregs_timeout, int lev); +int write_greens_nc(struct nestContainer *nest, char *fname, float *work, size_t *start, size_t *count, + double *t, unsigned int *lcum_p, char *names[], char hist[], int *ids, int n_maregs, + unsigned int n_times, int lev); +void err_trap_(int status); + +#if (defined(WIN32) || defined(_WIN32) || defined(_WIN64)) && (defined(DO_MULTI_THREAD) || defined(PARALLEL)) +/* Prototypes for threading related functions */ +unsigned __stdcall MT_cart(void *Arg_p); +unsigned __stdcall MT_sp(void *Arg_p); +unsigned __stdcall MT_moment_sp_M_slice(void *Arg_p); +unsigned __stdcall MT_moment_sp_N_slice(void *Arg_p); +unsigned __stdcall MT_moment_M_slice(void *Arg_p); +unsigned __stdcall MT_moment_N_slice(void *Arg_p); +unsigned __stdcall MT_mass_sp(void *Arg_p); +unsigned __stdcall MT_mass(void *Arg_p); +#endif +int GetLocalNThread(void); + +/* Function pointers to M & N moment components */ +PFV call_moment[2]; +PFV call_moment_sp[2]; +PFV call_moment_sp_M_slice; +PFV call_moment_sp_N_slice; +PFV call_moment_M_slice; +PFV call_moment_N_slice; +PFV call_mass; +PFV call_mass_sp; + +/* ------------------------------------------------------------------------------ */ +/* Read a grid through the GMT API. Works for both on-disk files and in-memory + GMT_GRID objects (virtual files, e.g. handed in from Julia). Fills the + Surfer-style header used internally and returns the GMT_GRID (the caller must + GMT_Destroy_Data() it after copying the data with gmtnswing_copy_grid()). */ +GMT_LOCAL struct GMT_GRID *gmtnswing_get_grid(void *API, char *fname, struct srf_header *h) { + struct GMT_GRID *G = NULL; + if ((G = GMT_Read_Data(API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, fname, NULL)) == NULL) + return NULL; + h->nx = (short int)G->header->n_columns; + h->ny = (short int)G->header->n_rows; + h->x_min = G->header->wesn[XLO]; h->x_max = G->header->wesn[XHI]; + h->y_min = G->header->wesn[YLO]; h->y_max = G->header->wesn[YHI]; + h->z_min = G->header->z_min; h->z_max = G->header->z_max; + return G; +} + +/* Copy a GMT grid (north-up, padded float) into an nswing array (south-up, plain + double, scanline order). 'sign' flips elevation -> depth (bathy uses sign = -1, + while source/momentum grids use sign = +1). */ +GMT_LOCAL void gmtnswing_copy_grid(struct GMT_GRID *G, double *work, int sign) { + uint64_t row, col, ij; + unsigned int nx = G->header->n_columns, ny = G->header->n_rows; + for (row = 0; row < ny; row++) { + for (col = 0; col < nx; col++) { + ij = gmt_M_ijp(G->header, row, col); /* GMT node: row 0 == north */ + work[(ny - 1 - row) * nx + col] = sign * G->data[ij]; + } + } +} + +/* --------------------------------------------------------------------------- */ +GMT_LOCAL int usage(struct GMTAPI_CTRL *API, int level) { + const char *name = gmt_show_name_and_purpose(API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE); + if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR); + GMT_Usage(API, 0, "usage: %s bathy.grd initial.grd [-1] [-2] [-3<...>] [-G|Z[+lev],] " + "[-A] [-B] [-C] [-D] [-E[p][m][,decim]] [-Fdip/strike/rake/slip/length/width/topDepth/x_epic/y_epic] " + "[-Fk[c]] [-H] [-H[,t]] [-J[+run_time_jump]] [-L[name1,name2]] " + "[-M[-|+[]]] [-N] [-O,] [%s] [-S[x|y|n][+m][+s]] [-T,[,]] " + "[-Q] [-X] -t
[%s] [-x] [%s]\n", name, GMT_Rgeo_OPT, GMT_V_OPT, GMT_f_OPT); + + if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS); + + GMT_Message(API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n"); + GMT_Usage(API, 1, "\n are the base level bathymetry and initial-condition (source) grids."); + GMT_Usage(API, 1, "\n-t
Time step for the simulation."); + GMT_Message(API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n"); + GMT_Usage(API, 1, "\n-1 -2 ... -9 Nested bathymetry grids (one per nesting level)."); + GMT_Usage(API, 1, "\n-A Save result as a .SWW ANUGA format file."); + GMT_Usage(API, 1, "\n-n Basename for MOST triplet files (no extension)."); + GMT_Usage(API, 1, "\n-B Name of a BoundaryCondition ASCII file."); + GMT_Usage(API, 1, "\n-C Add Coriolis effect."); + GMT_Usage(API, 1, "\n-D Write grids with the total water depth."); + GMT_Usage(API, -2, "These grids will have wave height on ocean and water thickness on land."); + GMT_Usage(API, 1, "\n-E[p][m][,decim] Write grids with energy or power (-Ep)."); + GMT_Usage(API, -2, "Append 'm' to save only one grid with the max values. This can noticeably slow the run, " + "so optionally append a decimator factor after the comma (causes aliasing visible on shaded illumination). " + "The file name comes from in -G/-Z complemented with a '_max' prefix; saving of multiple grids is disabled. " + "A 3D netCDF file with wave heights is still possible with -Z."); + GMT_Usage(API, 1, "\n-F Okada fault parameters."); + GMT_Usage(API, -2, "Dip, Azimuth, Rake, Slip(m), length, width and depth from sea-bottom; x_epic, y_epic are the " + "X and Y coordinates of the beginning of the fault trace. All dimensions must be in km."); + GMT_Usage(API, 1, "\n-Fk Build a prism source with these limits and height of 1 meter."); + GMT_Usage(API, -2, "-Fkc: alternatively give the prism size as centre x/y and nx/ny half-width cell numbers. " + "-Fk.../RxC: loop over a matrix of size R x C starting at the Lower Left Corner given by w/e/s/n. " + "-Fk.../dx[/dy]: given the w/e/s/n region (pixel registration) loop over the prisms obtained by dividing the " + "region in increments of dx/dy (if not given, dy = dx). Using -Fk sets the output maregraph file to netCDF, unless rows = cols = 1."); + GMT_Usage(API, 1, "\n-G Write grids at the intervals; files are named #.grd."); + GMT_Usage(API, -2, "When doing nested grids, append +lev to save that particular level (only one level is allowed)."); + GMT_Usage(API, 1, "\n-H Write grids with the momentum (velocity times water depth)."); + GMT_Usage(API, -2, "-H[,t]: Hot start using these moment grids. Optional 't' is the hot-start " + "time (also needs the surface displacement corresponding to the time of these grids)."); + GMT_Usage(API, 1, "\n-J Do not write grids or maregraphs for times before time_jump (seconds)."); + GMT_Usage(API, -2, "When doing nested grids, append +
[%s] [-x] [%s]\n", name, GMT_Rgeo_OPT, GMT_V_OPT, GMT_f_OPT); if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS); @@ -418,11 +418,11 @@ GMT_LOCAL int usage(struct GMTAPI_CTRL *API, int level) { GMT_Usage(API, 1, "\n-E[p][m][,decim] Write grids with energy or power (-Ep)."); GMT_Usage(API, -2, "Append 'm' to save only one grid with the max values. This can noticeably slow the run, " "so optionally append a decimator factor after the comma (causes aliasing visible on shaded illumination). " - "The file name comes from in -G/-Z complemented with a '_max' prefix; saving of multiple grids is disabled. " - "A 3D netCDF file with wave heights is still possible with -Z."); - GMT_Usage(API, 1, "\n-F Okada fault parameters."); - GMT_Usage(API, -2, "Dip, Azimuth, Rake, Slip(m), length, width and depth from sea-bottom; x_epic, y_epic are the " - "X and Y coordinates of the beginning of the fault trace. All dimensions must be in km."); + "The file name comes from in -G complemented with a '_max' prefix; saving of multiple grids is disabled. " + "A 3D netCDF file with wave heights is still possible with -G."); + GMT_Usage(API, 1, "\n-F Okada fault parameters."); + GMT_Usage(API, -2, "x_epic, y_epic are the X and Y coordinates of the beginning of the fault trace; Dip, Azimuth, " + "Rake, Slip(m), length, width and depth from sea-bottom follow. All dimensions must be in km."); GMT_Usage(API, 1, "\n-Fk Build a prism source with these limits and height of 1 meter."); GMT_Usage(API, -2, "-Fkc: alternatively give the prism size as centre x/y and nx/ny half-width cell numbers. " "-Fk.../RxC: loop over a matrix of size R x C starting at the Lower Left Corner given by w/e/s/n. " @@ -434,27 +434,32 @@ GMT_LOCAL int usage(struct GMTAPI_CTRL *API, int level) { GMT_Usage(API, 1, "\n-H Write grids with the momentum (velocity times water depth)."); GMT_Usage(API, -2, "-H[,t]: Hot start using these moment grids. Optional 't' is the hot-start " "time (also needs the surface displacement corresponding to the time of these grids)."); - GMT_Usage(API, 1, "\n-J Do not write grids or maregraphs for times before time_jump (seconds)."); - GMT_Usage(API, -2, "When doing nested grids, append +
[%s] [-x] [%s]\n", name, GMT_Rgeo_OPT, GMT_V_OPT, GMT_f_OPT); if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS); @@ -411,7 +411,6 @@ GMT_LOCAL int usage(struct GMTAPI_CTRL *API, int level) { GMT_Usage(API, 1, "\n-1 -2 ... -9 Nested bathymetry grids (one per nesting level)."); GMT_Usage(API, 1, "\n-A Save result as a .SWW ANUGA format file."); //GMT_Usage(API, 1, "\n-n Basename for MOST triplet files (no extension)."); - //GMT_Usage(API, 1, "\n-B Name of a BoundaryCondition ASCII file."); GMT_Usage(API, 1, "\n-C Add Coriolis effect."); GMT_Usage(API, 1, "\n-D Write grids with the total water depth."); GMT_Usage(API, -2, "These grids will have wave height on ocean and water thickness on land."); @@ -434,10 +433,10 @@ GMT_LOCAL int usage(struct GMTAPI_CTRL *API, int level) { GMT_Usage(API, 1, "\n-H Write grids with the momentum (velocity times water depth)."); GMT_Usage(API, -2, "-H[,t]: Hot start using these moment grids. Optional 't' is the hot-start " "time (also needs the surface displacement corresponding to the time of these grids)."); - GMT_Usage(API, 1, "\n-J"); + GMT_Usage(API, 1, "\n-P"); GMT_Usage(API, -2, "Do not write grids or maregraphs for times before (seconds)."); //GMT_Usage(API, -2, "When doing nested grids, append +
[%s] [-x] [%s]\n", name, GMT_Rgeo_OPT, GMT_V_OPT, GMT_f_OPT); if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS); @@ -427,8 +427,9 @@ GMT_LOCAL int usage(struct GMTAPI_CTRL *API, int level) { GMT_Usage(API, -2, "Use x or y to save only one component, or n for no velocity grids (maregs only). Append +m to also write " "velocity (vx,vy) at maregraph locations (needs -T). Append +s to write the max speed (|v|) ('_max_speed' suffix). " "Use the 'n' flag to NOT output the U and V components, e.g. -Sn+s."); - GMT_Usage(API, 1, "\n-T[+o][+t] Save maregraph (virtual tide-gauge) time series."); - GMT_Usage(API, -2, " is the file with the (x y) locations of the virtual maregraphs. Append +o to set the " + GMT_Usage(API, 1, "\n-T|[+o][+t] Save maregraph (virtual tide-gauge) time series."); + GMT_Usage(API, -2, " is the file with the (x y) locations of the virtual maregraphs. For a single maregraph the " + "location may be given directly as instead of a file name. Append +o to set the " "output file name [Default is maregs_out.dat]. A '.dat' extension is added when has none; use a '.nc' " "extension to write the maregraphs as a netCDF file instead. Append +t to save every simulation time " "steps (set by -t) [Default is every time step]."); @@ -453,7 +454,7 @@ struct NSWING_CTRL { bool do_tracers, out_maregs_nc, out_oranges_nc, do_HotStart, write_grids, isGeog; bool maregs_in_input, out_momentum, got_R, deform_only, with_land, saveNested; bool verbose, out_velocity, out_velocity_x, out_velocity_y, out_velocity_r, out_maregs_velocity; - bool cumpt, do_2Dgrids, do_maxs; + bool cumpt, do_2Dgrids, do_maxs, mareg_xy; char *bathy, *fonte, *fname_sww, *basename_most, *bnc_file; char *nesteds[10]; char hcum[256]; @@ -469,6 +470,7 @@ struct NSWING_CTRL { double kaba_xmax, kaba_ymin, kaba_ymax, time_jump, dt, f_dip; double f_azim, f_rake, f_slip, f_length, f_width, f_topDepth; double x_epic, y_epic, dfXmin, dfXmax, dfYmin, dfYmax; + double mareg_x, mareg_y; }; GMT_LOCAL void *New_Ctrl(struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */ @@ -493,7 +495,8 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC int grn = 0, cumint = 0, decimate_max = 1, error = 0; int do_Kaba = 0, n_of_cycles = 1010, KbGridCols = 1, KbGridRows = 1; size_t len; - bool cumpt = false, do_2Dgrids = false, do_maxs = false; + bool cumpt = false, do_2Dgrids = false, do_maxs = false, mareg_xy = false; + double mareg_x = 0, mareg_y = 0; bool out_energy = false, max_energy = false, out_power = false, max_power = false; bool out_sww = false, out_most = false, out_3D = false; bool surf_level = true, max_level = false, max_velocity = false, water_depth = false; @@ -803,7 +806,7 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC dt = atof(opt->arg); nest.dt[0] = dt; break; - case 'T': /* Maregraph xy positions file, with optional +o and +t modifiers */ + case 'T': /* Maregraph xy positions file (or a single x/y pair), with optional +o and +t modifiers */ if (cumpt) { GMT_Report(GMT->parent, GMT_MSG_ERROR, "NSWING: Error, -T option given more than once.\n"); GMT_Report(GMT->parent, GMT_MSG_WARNING, " Ignoring it.\n"); @@ -812,17 +815,19 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC cumint = 1; /* Default: save maregraphs at every time step */ hcum[0] = '\0'; if (opt->arg[0] == '\0') { - GMT_Report(GMT->parent, GMT_MSG_ERROR, "NSWING: Error, -T option, must provide the maregraphs xy file name\n"); + GMT_Report(GMT->parent, GMT_MSG_ERROR, "NSWING: Error, -T option, must provide the maregraphs xy file name (or a x/y pair)\n"); error++; break; } sscanf(opt->arg, "%127s", str_tmp); if ((tok = strtok(str_tmp, "+")) == NULL || tok[0] == '\0') { - GMT_Report(GMT->parent, GMT_MSG_ERROR, "NSWING: Error, -T option, must provide the maregraphs xy file name\n"); + GMT_Report(GMT->parent, GMT_MSG_ERROR, "NSWING: Error, -T option, must provide the maregraphs xy file name (or a x/y pair)\n"); error++; break; } strcpy(maregs, tok); + if (sscanf(maregs, "%lf/%lf%c", &mareg_x, &mareg_y, txt) == 2) + mareg_xy = true; /* A single maregraph given inline as -Tx/y; no positions file */ while ((tok = strtok(NULL, "+")) != NULL) { if (tok[0] == 'o') strcpy(hcum, &tok[1]); @@ -997,6 +1002,9 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC Ctrl->cumpt = cumpt; Ctrl->do_2Dgrids = do_2Dgrids; Ctrl->do_maxs = do_maxs; + Ctrl->mareg_xy = mareg_xy; + Ctrl->mareg_x = mareg_x; + Ctrl->mareg_y = mareg_y; Ctrl->add_const = add_const; Ctrl->time_h = time_h; Ctrl->dxKb = dxKb; @@ -1057,7 +1065,8 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { int i, j, k, n; int grn = 0, cumint = 0, decimate_max = 1, iprc, r_bin_b, r_bin_f, r_bin_mM, r_bin_mN; int error = 0; - bool w_bin = true, cumpt = false, do_2Dgrids = false, do_maxs = false; + bool w_bin = true, cumpt = false, do_2Dgrids = false, do_maxs = false, mareg_xy = false; + double mareg_x = 0, mareg_y = 0; bool out_energy = false, max_energy = false, out_power = false, max_power = false; bool first_anuga_time = true, out_sww = false, out_most = false, out_3D = false; bool surf_level = true, max_level = false, max_velocity = false, water_depth = false; @@ -1212,6 +1221,9 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { cumpt = Ctrl->cumpt; do_2Dgrids = Ctrl->do_2Dgrids; do_maxs = Ctrl->do_maxs; + mareg_xy = Ctrl->mareg_xy; + mareg_x = Ctrl->mareg_x; + mareg_y = Ctrl->mareg_y; add_const = Ctrl->add_const; time_h = Ctrl->time_h; dxKb = Ctrl->dxKb; @@ -1319,8 +1331,9 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { Return(-1); } if (!error && !maregs_in_input) { - n_mareg = count_n_maregs(API, maregs); /* Count maragraphs number */ - if (n_mareg < 0) { + if (mareg_xy) /* A single maregraph given inline via -Tx/y */ + n_mareg = 1; + else if ((n_mareg = count_n_maregs(API, maregs)) < 0) { /* Count maragraphs number */ Return(-1); /* Warning already issued in count_n_maregs() */ } else if (n_mareg == 0) { @@ -1540,7 +1553,21 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { if (cumpt && !maregs_in_input) { lcum_p = (unsigned int *)calloc((size_t)(1024), sizeof(unsigned int)); /* We wont ever use these many */ mareg_names = calloc((size_t)(1024), sizeof(char *)); - if ((n_mareg = read_maregs(API, nest.hdr[writeLevel], maregs, lcum_p, mareg_names)) < 1) { /* Read maregraph locations */ + if (mareg_xy) { /* A single maregraph given inline via -Tx/y */ + if (mareg_x < nest.hdr[writeLevel].wesn[XLO] || mareg_x > nest.hdr[writeLevel].wesn[XHI] || + mareg_y < nest.hdr[writeLevel].wesn[YLO] || mareg_y > nest.hdr[writeLevel].wesn[YHI]) + n_mareg = 0; + else { + i = irint((mareg_x - nest.hdr[writeLevel].wesn[XLO]) / nest.hdr[writeLevel].inc[GMT_X]); + j = irint((mareg_y - nest.hdr[writeLevel].wesn[YLO]) / nest.hdr[writeLevel].inc[GMT_Y]); + lcum_p[0] = j * nest.hdr[writeLevel].n_columns + i; + mareg_names[0] = strdup("NoName"); + n_mareg = 1; + } + } + else + n_mareg = read_maregs(API, nest.hdr[writeLevel], maregs, lcum_p, mareg_names); /* Read maregraph locations */ + if (n_mareg < 1) { GMT_Report(API, GMT_MSG_WARNING, "NSWING - WARNING: No maregraphs inside the (inner?) grid\n"); n_mareg = 0; if (lcum_p) {free(lcum_p); lcum_p = NULL;} From 7af47bc2d9fc664235401ac49f2cb6db50d2e0c0 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Thu, 2 Jul 2026 21:25:57 +0100 Subject: [PATCH 20/29] nswing: moment_N zeroed fluxm_d instead of fluxn_d The #ifndef PARALLEL prologue of the cartesian moment_N() did memset(fluxm_d, ...) - moment_M's output array - instead of its own fluxn_d. Under DO_MULTI_THREAD moment_M and moment_N run concurrently, so the wipe landed at a random point while moment_M was filling fluxm_d: cartesian results were nondeterministic run to run (up to ~0.27 m differences at a gauge). In a serial build the wipe simply zeroed the x-fluxes each step. The three sibling functions (moment_M, moment_sp_M, moment_sp_N) and the PARALLEL-mode hoisted memsets all use the right array, which is why geographic (-fg) runs were always deterministic and unaffected. Verified after the one-word fix: identical cartesian commands are now bit-identical run to run, a serial build matches the threaded one, and geographic output is unchanged. Co-Authored-By: Claude Fable 5 --- src/nswing/nswing.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nswing/nswing.c b/src/nswing/nswing.c index f76a7f8948c..c8bc9b8eaf2 100644 --- a/src/nswing/nswing.c +++ b/src/nswing/nswing.c @@ -4209,7 +4209,7 @@ void moment_N_slice(int lev, int row_start, int row_end, struct nestContainer *n #ifndef PARALLEL row_start = 0; row_end = hdr.n_rows - 1; - memset(fluxm_d, 0, hdr.nm * sizeof(double)); + memset(fluxn_d, 0, hdr.nm * sizeof(double)); /* fluxn_d, NOT fluxm_d: zeroing fluxm_d here raced with moment_M writing it */ #endif dtdx = dt / hdr.inc[GMT_X]; From a2af52a7bf5b71dce522803226b5bbef5186fdb3 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Fri, 3 Jul 2026 00:17:39 +0100 Subject: [PATCH 21/29] nswing: replace home-grown threading with OpenMP Drop both the DO_MULTI_THREAD path (2 Windows threads created and destroyed every time step, moment only) and the never-enabled PARALLEL path (per-step slice threads whose integer-division row splits silently skipped up to n_threads-1 rows). In their place, OpenMP parallel-for with schedule(dynamic, 8) on the six row loops: mass, mass_sp, moment_M, moment_N, moment_sp_M and moment_sp_N. mass is now parallelized too, all cores participate, and threads persist in the OpenMP pool instead of being respawned each step. About 330 lines of scaffolding removed: MT_* wrappers, *_slice function variants and their mid-function #ifdef body-sharing tricks, ThreadArg, the PFV function pointers. -x now maps to omp_set_num_threads() [default all cores]. dynamic scheduling matters on hybrid P/E-core CPUs, where static splits made every per-step barrier wait for the slowest efficiency core and ran slower than the old 2-thread code. Results are bit-identical to the old build and across schedules and thread counts (the loops have no reductions), so runs are finally reproducible regardless of -x. The four moment loops also zero dd/df per cell now: the goto L121/L201 shortcuts could read them stale from the previous cell (a latent NaN source for -S velocities, and a loop-carried dependency the parallel loops must not have). Timings (20-core machine, maregraphs-only runs): 765x476 geographic, 1010 steps: 5.8-7.3 s old -> 3.5 s; 1201x1201 cartesian, 300 steps: 12.0 s old -> 4.4 s. Non-Windows builds gain parallelism for the first time (old threading was Windows-only). Co-Authored-By: Claude Fable 5 --- src/nswing/nswing.c | 465 ++++---------------------------------------- 1 file changed, 41 insertions(+), 424 deletions(-) diff --git a/src/nswing/nswing.c b/src/nswing/nswing.c index c8bc9b8eaf2..777c9ba9446 100644 --- a/src/nswing/nswing.c +++ b/src/nswing/nswing.c @@ -53,18 +53,10 @@ #define THIS_MODULE_NEEDS "" #define THIS_MODULE_OPTIONS "-RVf" -//#define PARALLEL - -#if defined(WIN32) || defined(_WIN32) || defined(_WIN64) -# define DO_MULTI_THREAD /* ISTO TEM DE SER AUTOMATIZADO, OU VIA COMPILA */ -#else +#if !(defined(WIN32) || defined(_WIN32) || defined(_WIN64)) # define strtok_s strtok_r #endif -#ifdef PARALLEL -# undef DO_MULTI_THREAD -#endif - /* This file is now built exclusively as a GMT supplement (module GMT_nswing). The former MEX gateway and stand-alone main() builds were removed. All grid INPUT now goes through the GMT API (GMT_Read_Data), which transparently reads @@ -93,7 +85,7 @@ union {uint64_t i; double d;} loc_nan = {0x7ff8000000000000}; #define err_trap(api, status) if (status) {GMT_Report(api, GMT_MSG_ERROR, "NSWING: error at line %d: %s\n", __LINE__, nc_strerror(status));} -#if HAVE_OPENMP +#ifdef _OPENMP #include #endif @@ -129,8 +121,6 @@ static double EPS4 = EPS4_; /* Kinda trick to be able to change EPS4 via a comm #define ijc(i,j) ((i) + (j)*n_ptmar) #define ij_grd(col,row,hdr) ((col) + (row)*hdr.n_columns) -typedef void (*PFV) (); /* PFV declares a pointer to a function returning void */ - struct tracers { /* For tracers (oranges) */ double *x; /* x coordinate */ double *y; /* y coordinate */ @@ -188,14 +178,6 @@ struct nestContainer { /* Container for the nestings */ struct GMT_GRID_HEADER hdr[10]; }; -/* Argument struct for threading */ -typedef struct { - struct nestContainer *nest; /* Pointer to a nestContainer struct */ - int iThread; /* Thread index */ - int lev; /* Level of nested grid */ - int row_start, row_end; /* Start and end rows for the BY slices in PARALLEL mode */ -} ThreadArg; - void no_sys_mem(void *API, char *where, unsigned int n); int count_col(char *line); int read_maregs(void *API, struct GMT_GRID_HEADER hdr, char *file, unsigned int *lcum_p, char *names[]); @@ -227,14 +209,8 @@ void upscale_(struct nestContainer *nest, double *out, int lev, int i_tsr); void replicate(struct nestContainer *nest, int lev); void moment_M(struct nestContainer *nest, int lev); void moment_N(struct nestContainer *nest, int lev); -void moment_M_slice(int lev, int row_start, int row_end, struct nestContainer *nest); -void moment_N_slice(int lev, int row_start, int row_end, struct nestContainer *nest); void moment_sp_M(struct nestContainer *nest, int lev); void moment_sp_N(struct nestContainer *nest, int lev); -void moment_sp_M_slice(int lev, int row_start, int row_end, struct nestContainer *nest); -void moment_sp_N_slice(int lev, int row_start, int row_end, struct nestContainer *nest); -void mass_sp_slice(struct nestContainer *nest, int lev, int row_start, int row_end); -void mass_slice(struct nestContainer *nest, int lev, int row_start, int row_end); void free_arrays(struct nestContainer *nest, int isGeog, int lev); int check_paternity(void *API, struct nestContainer *nest); int check_binning(double x0P, double x0D, double dxP, double dxD, double tol, double *suggest); @@ -275,29 +251,8 @@ int write_greens_nc(void *API, struct nestContainer *nest, char *fname, float *w unsigned int n_times, int lev); void err_trap_(void *API, int status); -#if (defined(WIN32) || defined(_WIN32) || defined(_WIN64)) && (defined(DO_MULTI_THREAD) || defined(PARALLEL)) -/* Prototypes for threading related functions */ -unsigned __stdcall MT_cart(void *Arg_p); -unsigned __stdcall MT_sp(void *Arg_p); -unsigned __stdcall MT_moment_sp_M_slice(void *Arg_p); -unsigned __stdcall MT_moment_sp_N_slice(void *Arg_p); -unsigned __stdcall MT_moment_M_slice(void *Arg_p); -unsigned __stdcall MT_moment_N_slice(void *Arg_p); -unsigned __stdcall MT_mass_sp(void *Arg_p); -unsigned __stdcall MT_mass(void *Arg_p); -#endif int GetLocalNThread(void); -/* Function pointers to M & N moment components */ -PFV call_moment[2]; -PFV call_moment_sp[2]; -PFV call_moment_sp_M_slice; -PFV call_moment_sp_N_slice; -PFV call_moment_M_slice; -PFV call_moment_N_slice; -PFV call_mass; -PFV call_mass_sp; - /* ------------------------------------------------------------------------------ */ /* Read a grid through the GMT API. Works for both on-disk files and in-memory GMT_GRID objects (virtual files, e.g. handed in from Julia). Fills the @@ -1136,26 +1091,6 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { FILE *fp = NULL, *fp_oranges = NULL; clock_t tic; - call_moment[0] = (PFV)moment_M; - call_moment[1] = (PFV)moment_N; - call_moment_sp[0] = (PFV)moment_sp_M; - call_moment_sp[1] = (PFV)moment_sp_N; - -#ifdef PARALLEL - call_moment_sp_M_slice = (PFV)moment_sp_M_slice; - call_moment_sp_N_slice = (PFV)moment_sp_N_slice; - call_mass_sp = (PFV)mass_sp_slice; - call_mass = (PFV)mass_slice; -#endif - -#ifdef DO_MULTI_THREAD - if ((k = GetLocalNThread()) == 1) { - GMT_Report(API, GMT_MSG_ERROR, "NSWING: This version of the program is build for multi-threading but " - "this machine has only one core. Don't know what will happen here.\n"); - } -#endif - - /*----------------------- Standard module initialization and parsing ----------------------*/ if (API == NULL) return (GMT_NOT_A_SESSION); @@ -1179,6 +1114,10 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { Ctrl = New_Ctrl(GMT); /* Allocate the control structure */ if ((error = parse(GMT, Ctrl, &nest, options)) != 0) Return(error); +#ifdef _OPENMP + if (nest.n_threads > 0) omp_set_num_threads(nest.n_threads); /* Honor -x [Default is all cores] */ +#endif + /* Unpack the parsed values into the local variables used by the main code below */ writeLevel = Ctrl->writeLevel; grn = Ctrl->grn; @@ -1788,8 +1727,8 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { GMT_Report(API, GMT_MSG_INFORMATION, "Computing a grid of prisms with size %d (rows) x %d (cols)\n", KbGridRows, KbGridCols); if (EPS4 != EPS4_) GMT_Report(API, GMT_MSG_INFORMATION, "Using a modified EPS4 const of %g\n", EPS4); -#ifdef PARALLEL - GMT_Report(API, GMT_MSG_INFORMATION, "\nUsing %d cores\n", nest.n_threads); +#ifdef _OPENMP + GMT_Report(API, GMT_MSG_INFORMATION, "\nUsing %d OpenMP threads\n", omp_get_max_threads()); #endif #ifdef LIMIT_DISCHARGE GMT_Report(API, GMT_MSG_INFORMATION, "\nUsing DISCHARGE limit to minimize sources of instability\n"); @@ -3707,13 +3646,8 @@ void err_trap_(void *API, int status) { * * Updates only etad and htotal_d * -------------------------------------------------------------------- */ -#ifdef PARALLEL -void mass_slice(struct nestContainer *nest, int lev, int row_start, int row_end) { -#else void mass(struct nestContainer *nest, int lev) { int row_start, row_end; -#endif - int row, col; int cm1, rm1; /* previous column (cm1 = col -1) and row (rm1 = row - 1) */ unsigned int ij; @@ -3727,10 +3661,11 @@ void mass(struct nestContainer *nest, int lev) { dtdx = nest->dt[lev] / nest->hdr[lev].inc[GMT_X]; dtdy = nest->dt[lev] / nest->hdr[lev].inc[GMT_Y]; -#ifndef PARALLEL row_start = 0; row_end = nest->hdr[lev].n_rows; -#endif +#ifdef _OPENMP +#pragma omp parallel for private(col, ij, rm1, cm1, zzz, dd) schedule(dynamic, 8) +#endif for (row = row_start; row < row_end; row++) { ij = row * nest->hdr[lev].n_columns; rm1 = (row == 0) ? 0 : nest->hdr[lev].n_columns; @@ -4008,34 +3943,14 @@ void moment_M(struct nestContainer *nest, int lev) { row_start = 0; row_end = hdr.n_rows - nest->last; memset(fluxm_d, 0, hdr.nm * sizeof(double)); -#ifdef PARALLEL -} /* Than the moment_M() fun ends here */ -void moment_M_slice(int lev, int row_start, int row_end, struct nestContainer *nest) { - unsigned int ij; - int row, col; - int valid_vel; - int cm1, rm1; /* previous column (cm1 = col - 1) and row (rm1 = row - 1) */ - int cp1, rp1; /* next column (cp1 = col + 1) and row (rp1 = row + 1) */ - int rm2, cp2; - double xp, xqe, xqq, ff = 0, dd, df, f_limit; - double advx, dtdx, dtdy, advy, rlat; - double dpa_ij, dpa_ij_rp1, dpa_ij_rm1, dpa_ij_cm1, dpa_ij_cp1; - double dt, manning, *bat, *htotal_a, *htotal_d, *etad, *fluxm_a, *fluxm_d, *fluxn_a, *fluxn_d, *vex, *r4m; - struct GMT_GRID_HEADER hdr; - - hdr = nest->hdr[lev]; vex = nest->vex[lev]; - dt = nest->dt[lev]; manning = nest->manning[lev]; - bat = nest->bat[lev]; etad = nest->etad[lev]; - htotal_a = nest->htotal_a[lev]; htotal_d = nest->htotal_d[lev]; - fluxm_a = nest->fluxm_a[lev]; fluxm_d = nest->fluxm_d[lev]; - fluxn_a = nest->fluxn_a[lev]; fluxn_d = nest->fluxn_d[lev]; - r4m = nest->r4m[lev]; -#endif - dtdx = dt / hdr.inc[GMT_X]; dtdy = dt / hdr.inc[GMT_Y]; manning *= dt; /* Finish the cte part now that we know 'dt': manning * manning * dt * 4.9 */ +#ifdef _OPENMP +#pragma omp parallel for private(col, ij, valid_vel, cm1, rm1, cp1, rp1, rm2, cp2, xp, xqe, xqq, ff, dd, df, \ + f_limit, advx, advy, dpa_ij, dpa_ij_rp1, dpa_ij_rm1, dpa_ij_cm1, dpa_ij_cp1) schedule(dynamic, 8) +#endif for (row = row_start; row < row_end; row++) { /* main computation cycle fluxm_d */ rp1 = (row < hdr.n_rows - 1) ? hdr.n_columns : 0; rm1 = (row == 0) ? 0 : hdr.n_columns; @@ -4051,7 +3966,7 @@ void moment_M_slice(int lev, int row_start, int row_end, struct nestContainer *n /* Looks weird but it's faster than an IF case (branch prediction?) */ dpa_ij = (dpa_ij = (htotal_d[ij] + htotal_a[ij] + htotal_d[ij+cp1] + htotal_a[ij+cp1]) * 0.25) > EPS5 ? dpa_ij : 0; - xp = 0; + xp = 0; dd = df = 0; /* dd/df must not carry over from the previous cell (goto L121 paths read them) */ valid_vel = true; if (htotal_d[ij] > EPS5 && htotal_d[ij+cp1] > EPS5) { /* case wet-wet */ @@ -4179,14 +4094,8 @@ void moment_M_slice(int lev, int row_start, int row_end, struct nestContainer *n } /* -------------------------------------------------------------------- */ -#ifndef PARALLEL void moment_N(struct nestContainer *nest, int lev) { int row_start, row_end; -#else -void moment_N(struct nestContainer *nest, int lev) {} // UGLY. Just for it to exist -void moment_N_slice(int lev, int row_start, int row_end, struct nestContainer *nest) { -#endif - unsigned int ij; int row, col; int valid_vel; @@ -4207,15 +4116,17 @@ void moment_N_slice(int lev, int row_start, int row_end, struct nestContainer *n fluxn_a = nest->fluxn_a[lev]; fluxn_d = nest->fluxn_d[lev]; r4n = nest->r4n[lev]; -#ifndef PARALLEL row_start = 0; row_end = hdr.n_rows - 1; memset(fluxn_d, 0, hdr.nm * sizeof(double)); /* fluxn_d, NOT fluxm_d: zeroing fluxm_d here raced with moment_M writing it */ -#endif dtdx = dt / hdr.inc[GMT_X]; dtdy = dt / hdr.inc[GMT_Y]; manning *= dt; /* Finish the cte part now that we know 'dt': manning * manning * dt * 4.9 */ +#ifdef _OPENMP +#pragma omp parallel for private(col, ij, valid_vel, cm1, rm1, cp1, rp1, cm2, rp2, xq, xpe, xpp, ff, dd, df, \ + f_limit, advx, advy, dqa_ij, dqa_ij_rp1, dqa_ij_rm1, dqa_ij_cm1, dqa_ij_cp1) schedule(dynamic, 8) +#endif for (row = row_start; row < row_end; row++) { /* main computation cycle fluxn_d */ rp1 = hdr.n_columns; rp2 = (row < hdr.n_rows - 2) ? 2*hdr.n_columns : hdr.n_columns; @@ -4230,7 +4141,7 @@ void moment_N_slice(int lev, int row_start, int row_end, struct nestContainer *n /* Looks weird but it's faster than an IF case (branch prediction?) */ dqa_ij = (dqa_ij = (htotal_d[ij] + htotal_a[ij] + htotal_d[ij+rp1] + htotal_a[ij+rp1]) * 0.25) > EPS5 ? dqa_ij : 0; - xq = 0; + xq = 0; dd = df = 0; /* dd/df must not carry over from the previous cell (goto L201 paths read them) */ /* moving boundary - Imamura algorithm following cho 2009 */ valid_vel = true; @@ -4423,16 +4334,9 @@ void mass_sp(struct nestContainer *nest, int lev) { row_start = 0; row_end = nest->hdr[lev].n_rows; -#ifdef PARALLEL -} /* Than the mass_sp() fun end here */ -void mass_sp_slice(struct nestContainer *nest, int lev, int row_start, int row_end) { - unsigned int ij; - int row, col; - int cm1, rm1, rowm1; /* previous column (cm1 = col -1) and row (rm1 = row - 1) */ - double etan, dd; - double r2m_r, r2n_r, r1n_r, r1n_r1; +#ifdef _OPENMP +#pragma omp parallel for private(col, ij, rm1, cm1, rowm1, etan, dd, r2m_r, r2n_r, r1n_r, r1n_r1) schedule(dynamic, 8) #endif - for (row = row_start; row < row_end; row++) { ij = row * nest->hdr[lev].n_columns; rm1 = ((row == 0) ? 0 : 1) * nest->hdr[lev].n_columns; @@ -4510,39 +4414,11 @@ void moment_sp_M(struct nestContainer *nest, int lev) { memset(fluxm_d, 0, hdr.nm * sizeof(double)); -#ifdef PARALLEL -} /* Than the moment_sp_M() fun end here */ -void moment_sp_M_slice(int lev, int row_start, int row_end, struct nestContainer *nest) { - unsigned int ij; - int valid_vel; - int row, col; - int cm1, rm1; /* previous column (cm1 = col - 1) and row (rm1 = row - 1) */ - int cp1, rp1; /* next column (cp1 = col + 1) and row (rp1 = row + 1) */ - int rm2, cp2; - double ff = 0; - double dd, df, xp, xqe, xqq, advx, advy, f_limit; - double dpa_ij, dpa_ij_rp1, dpa_ij_rm1, dpa_ij_cm1, dpa_ij_cp1; - double dt, manning, *htotal_a, *htotal_d, *bat, *etad, *fluxm_a, *fluxn_a, *fluxm_d, *fluxn_d, *vex; - double *r0, *r2m, *r3m, *r4m, r2m_r; - double bat__ij; - double htotal_d__ij; - double htotal_d__ij_cp1; - double etad__ij, etad__ij_cp1; - double fluxm_a__ij; - struct GMT_GRID_HEADER hdr; - - hdr = nest->hdr[lev]; vex = nest->vex[lev]; - dt = nest->dt[lev]; manning = nest->manning[lev]; - bat = nest->bat[lev]; etad = nest->etad[lev]; - htotal_a = nest->htotal_a[lev]; htotal_d = nest->htotal_d[lev]; - fluxm_a = nest->fluxm_a[lev]; fluxm_d = nest->fluxm_d[lev]; - fluxn_a = nest->fluxn_a[lev]; fluxn_d = nest->fluxn_d[lev]; - r0 = nest->r0[lev]; r2m = nest->r2m[lev]; - r3m = nest->r3m[lev]; r4m = nest->r4m[lev]; - - manning *= dt; /* Finish the cte part now that we know 'dt': manning * manning * dt * 4.9 */ +#ifdef _OPENMP +#pragma omp parallel for private(col, ij, valid_vel, cm1, rm1, cp1, rp1, rm2, cp2, ff, dd, df, xp, xqe, xqq, \ + advx, advy, f_limit, dpa_ij, dpa_ij_rp1, dpa_ij_rm1, dpa_ij_cm1, dpa_ij_cp1, \ + r2m_r, bat__ij, htotal_d__ij, htotal_d__ij_cp1, etad__ij, etad__ij_cp1, fluxm_a__ij) schedule(dynamic, 8) #endif - for (row = row_start; row < row_end; row++) { /* - main computation cycle fluxm_d */ rp1 = (row < hdr.n_rows - 1) ? hdr.n_columns : 0; rm1 = (row == 0) ? 0 : hdr.n_columns; @@ -4564,7 +4440,7 @@ void moment_sp_M_slice(int lev, int row_start, int row_end, struct nestContainer etad__ij = etad[ij]; etad__ij_cp1 = etad[ij+cp1]; fluxm_a__ij = fluxm_a[ij]; - xp = 0; + xp = 0; dd = df = 0; /* dd/df must not carry over from the previous cell (goto L121 paths read them) */ /* Looks weird but it's faster than an IF case (branch prediction?) */ dpa_ij = (dpa_ij = (htotal_d__ij + htotal_a[ij] + htotal_d__ij_cp1 + htotal_a[ij+cp1]) * 0.25) > EPS5 ? dpa_ij : 0; @@ -4719,39 +4595,11 @@ void moment_sp_N(struct nestContainer *nest, int lev) { memset(fluxn_d, 0, hdr.nm * sizeof(double)); -#ifdef PARALLEL -} /* Than the moment_sp_N() fun ends here */ -void moment_sp_N_slice(int lev, int row_start, int row_end, struct nestContainer *nest) { - unsigned int ij; - int valid_vel; - int row, col; - int cm1, rm1; /* previous column (cm1 = col - 1) and row (rm1 = row - 1) */ - int cp1, rp1; /* next column (cp1 = col + 1) and row (rp1 = row + 1) */ - int cm2, rp2; - double ff = 0; - double dd, df, xq, xpe, xpp, advx, advy, f_limit; - double dqa_ij, dqa_ij_rp1, dqa_ij_rm1, dqa_ij_cm1, dqa_ij_cp1; - double dt, manning, *htotal_a, *htotal_d, *bat, *etad, *fluxm_a, *fluxn_a, *fluxm_d, *fluxn_d, *vey; - double *r0, *r2n, *r3n, *r4n, r2n_r; - double bat__ij; - double htotal_d__ij; - double htotal_a__ij_rp1, htotal_d__ij_rp1; - double etad__ij, etad__ij_rp1; - double fluxn_a__ij; - struct GMT_GRID_HEADER hdr; - - hdr = nest->hdr[lev]; vey = nest->vey[lev]; - dt = nest->dt[lev]; manning = nest->manning[lev]; - bat = nest->bat[lev]; etad = nest->etad[lev]; - htotal_a = nest->htotal_a[lev]; htotal_d = nest->htotal_d[lev]; - fluxm_a = nest->fluxm_a[lev]; fluxm_d = nest->fluxm_d[lev]; - fluxn_a = nest->fluxn_a[lev]; fluxn_d = nest->fluxn_d[lev]; - r0 = nest->r0[lev]; r2n = nest->r2n[lev]; - r3n = nest->r3n[lev]; r4n = nest->r4n[lev]; - - manning *= dt; /* Finish the cte part now that we know 'dt': manning * manning * dt * 4.9 */ +#ifdef _OPENMP +#pragma omp parallel for private(col, ij, valid_vel, cm1, rm1, cp1, rp1, cm2, rp2, ff, dd, df, xq, xpe, xpp, \ + advx, advy, f_limit, dqa_ij, dqa_ij_rp1, dqa_ij_rm1, dqa_ij_cm1, dqa_ij_cp1, \ + r2n_r, bat__ij, htotal_d__ij, htotal_d__ij_rp1, htotal_a__ij_rp1, etad__ij, etad__ij_rp1, fluxn_a__ij) schedule(dynamic, 8) #endif - for (row = row_start; row < row_end; row++) { /* - main computation cycle fluxn_d */ rp1 = hdr.n_columns; rp2 = (row < hdr.n_rows - 2) ? 2*hdr.n_columns : hdr.n_columns; @@ -4775,7 +4623,7 @@ void moment_sp_N_slice(int lev, int row_start, int row_end, struct nestContainer htotal_a__ij_rp1 = htotal_a[ij+rp1]; etad__ij_rp1 = etad[ij+rp1]; fluxn_a__ij = fluxn_a[ij]; - xq = 0; + xq = 0; dd = df = 0; /* dd/df must not carry over from the previous cell (goto L201 paths read them) */ /* Looks weird but it's faster than an IF case (branch prediction?) */ dqa_ij = (dqa_ij = (htotal_d__ij + htotal_a[ij] + htotal_d__ij_rp1 + htotal_a__ij_rp1) * 0.25) > EPS5 ? dqa_ij : 0; @@ -5276,275 +5124,44 @@ void edge_communication(void *API, struct nestContainer *nest, int lev, int i_ti /* ------------------------------------------------------------------------------ */ void mass_conservation(struct nestContainer *nest, int isGeog, int m) { /* m is the level of nesting which starts counting at one for FIRST nesting level */ - int i, row_start, row_end, n_rows_block; - -#ifdef PARALLEL - HANDLE ThreadList[64]; /* Handles to the worker threads */ - ThreadArg Arg_List[64], *Arg_p; - - n_rows_block = (int)ceil(nest->hdr[m].n_rows / nest->n_threads); - if (isGeog) { - for (i = 0; i < nest->n_threads; i++) { - row_start = i * n_rows_block; - row_end = MIN(row_start + n_rows_block, nest->hdr[m].n_rows); - Arg_p = &Arg_List[i]; - Arg_p->nest = nest; - Arg_p->lev = m; - Arg_p->row_start= row_start; - Arg_p->row_end = row_end; - ThreadList[i] = (HANDLE)_beginthreadex(NULL, 0, MT_mass_sp, Arg_p, 0, NULL); - } - } - else { - for (i = 0; i < nest->n_threads; i++) { - row_start = i * n_rows_block; - row_end = MIN(row_start + n_rows_block, nest->hdr[m].n_rows); - Arg_p = &Arg_List[i]; - Arg_p->nest = nest; - Arg_p->lev = m; - Arg_p->row_start= row_start; - Arg_p->row_end = row_end; - ThreadList[i] = (HANDLE)_beginthreadex(NULL, 0, MT_mass, Arg_p, 0, NULL); - } - } - /* Wait until all threads are ready and close the handles */ - WaitForMultipleObjects(nest->n_threads, ThreadList, true, INFINITE); - for (i = 0; i < nest->n_threads; i++) - CloseHandle(ThreadList[i]); -#else if (isGeog) mass_sp(nest, m); else mass(nest, m); -#endif } /* ------------------------------------------------------------------------------ */ void moment_conservation(struct nestContainer *nest, int isGeog, int m) { /* m is the level of nesting which starts counting at one for FIRST nesting level */ - int i; - -#ifdef DO_MULTI_THREAD - HANDLE ThreadList[2]; /* Handles to the worker threads */ - ThreadArg Arg_List[2], *Arg_p; - - if (m > 0) {nest->jupe = 0; nest->first = 1; nest->last = 0;} - else {nest->jupe = 10; nest->first = 0; nest->last = 1;} - if (nest->do_linear) nest->jupe = 1e6; /* A tricky way of imposing linearity */ - - if (isGeog == 0) { - for (i = 0; i < 2; i++) { - Arg_p = &Arg_List[i]; - Arg_p->nest = nest; - Arg_p->iThread = i; - Arg_p->lev = m; - ThreadList[i] = (HANDLE)_beginthreadex(NULL, 0, MT_cart, Arg_p, 0, NULL); - } - } - else { - for (i = 0; i < 2; i++) { - Arg_p = &Arg_List[i]; - Arg_p->nest = nest; - Arg_p->iThread = i; - Arg_p->lev = m; - ThreadList[i] = (HANDLE)_beginthreadex(NULL, 0, MT_sp, Arg_p, 0, NULL); - } - } - - /* Wait until all threads are ready and close the handles */ - WaitForMultipleObjects(2, ThreadList, true, INFINITE); - for (i = 0; i < 2; i++) - CloseHandle(ThreadList[i]); - -#elif defined(PARALLEL) - HANDLE ThreadList[64]; /* Handles to the worker threads */ - ThreadArg Arg_List[64], *Arg_p; - int row_start, row_end, last, n_rows_block; - /* - fixes the width of the lateral buffer for linear aproximation */ /* - if jupe>nnx/2 and jupe>nny/2 linear model will be applied */ if (m > 0) {nest->jupe = 0; nest->first = 1; nest->last = 0;} else {nest->jupe = 10; nest->first = 0; nest->last = 1;} if (nest->do_linear) nest->jupe = 1e6; /* A tricky way of imposing linearity */ - - last = (m > 0) ? 0 : 1; - if (isGeog) { - memset(nest->fluxm_d[m], 0, nest->hdr[m].nm * sizeof(double)); - n_rows_block = (int)ceil((nest->hdr[m].n_rows - last) / nest->n_threads); - for (i = 0; i < nest->n_threads; i++) { - row_start = i * n_rows_block; - row_end = MIN(row_start + n_rows_block, nest->hdr[m].n_rows - last); - Arg_p = &Arg_List[i]; - Arg_p->nest = nest; - Arg_p->iThread = i; - Arg_p->lev = m; - Arg_p->row_start= row_start; - Arg_p->row_end = row_end; - ThreadList[i] = (HANDLE)_beginthreadex(NULL, 0, MT_moment_sp_M_slice, Arg_p, 0, NULL); - } - /* Wait until all threads are ready and close the handles */ - WaitForMultipleObjects(nest->n_threads, ThreadList, true, INFINITE); - for (i = 0; i < nest->n_threads; i++) - CloseHandle(ThreadList[i]); - - /* Now the N component */ - memset(nest->fluxn_d[m], 0, nest->hdr[m].nm * sizeof(double)); - n_rows_block = (int)ceil((nest->hdr[m].n_rows - 1) / nest->n_threads); - for (i = 0; i < nest->n_threads; i++) { - row_start = i * n_rows_block; - row_end = MIN(row_start + n_rows_block, nest->hdr[m].n_rows - 1); - Arg_p = &Arg_List[i]; - Arg_p->nest = nest; - Arg_p->iThread = i; - Arg_p->lev = m; - Arg_p->row_start= row_start; - Arg_p->row_end = row_end; - ThreadList[i] = (HANDLE)_beginthreadex(NULL, 0, MT_moment_sp_N_slice, Arg_p, 0, NULL); - } - /* Wait until all threads are ready and close the handles */ - WaitForMultipleObjects(nest->n_threads, ThreadList, true, INFINITE); - for (i = 0; i < nest->n_threads; i++) - CloseHandle(ThreadList[i]); + /* The M and N components are internally parallelized with OpenMP over their row loops */ + if (isGeog == 0) { + moment_M(nest, m); + moment_N(nest, m); } else { - memset(nest->fluxm_d[m], 0, nest->hdr[m].nm * sizeof(double)); - n_rows_block = (int)ceil((nest->hdr[m].n_rows - last) / nest->n_threads); - for (i = 0; i < nest->n_threads; i++) { - row_start = i * n_rows_block; - row_end = MIN(row_start + n_rows_block, nest->hdr[m].n_rows - last); - Arg_p = &Arg_List[i]; - Arg_p->nest = nest; - Arg_p->iThread = i; - Arg_p->lev = m; - Arg_p->row_start= row_start; - Arg_p->row_end = row_end; - ThreadList[i] = (HANDLE)_beginthreadex(NULL, 0, MT_moment_M_slice, Arg_p, 0, NULL); - } - /* Wait until all threads are ready and close the handles */ - WaitForMultipleObjects(nest->n_threads, ThreadList, true, INFINITE); - for (i = 0; i < nest->n_threads; i++) - CloseHandle(ThreadList[i]); - - /* Now the N component */ - memset(nest->fluxn_d[m], 0, nest->hdr[m].nm * sizeof(double)); - n_rows_block = (int)ceil((nest->hdr[m].n_rows - 1) / nest->n_threads); - for (i = 0; i < nest->n_threads; i++) { - row_start = i * n_rows_block; - row_end = MIN(row_start + n_rows_block, nest->hdr[m].n_rows - 1); - Arg_p = &Arg_List[i]; - Arg_p->nest = nest; - Arg_p->iThread = i; - Arg_p->lev = m; - Arg_p->row_start= row_start; - Arg_p->row_end = row_end; - ThreadList[i] = (HANDLE)_beginthreadex(NULL, 0, MT_moment_N_slice, Arg_p, 0, NULL); - } - /* Wait until all threads are ready and close the handles */ - WaitForMultipleObjects(nest->n_threads, ThreadList, true, INFINITE); - for (i = 0; i < nest->n_threads; i++) - CloseHandle(ThreadList[i]); - } - -#else - if (m > 0) {nest->jupe = 0; nest->first = 1; nest->last = 0;} - else {nest->jupe = 10; nest->first = 0; nest->last = 1;} - if (nest->do_linear) nest->jupe = 1e6; /* A tricky way of imposing linearity */ - - if (isGeog == 0) { - for (i = 0; i < 2; i++) - call_moment[i](nest, m); - } - else { - for (i = 0; i < 2; i++) - call_moment_sp[i](nest, m); + moment_sp_M(nest, m); + moment_sp_N(nest, m); } -#endif } -#ifdef DO_MULTI_THREAD -/* ------------------------------------------------------------------------------ */ -unsigned __stdcall MT_cart(void *Arg_p) { - /* Convert input from (void *) to (ThreadArg *), call the moment and stop the thread. */ - ThreadArg *Arg = (ThreadArg *)Arg_p; - call_moment[Arg->iThread](Arg->nest, Arg->lev); - _endthreadex(0); - return (0); -} -/* ------------------------------------------------------------------------------ */ -unsigned __stdcall MT_sp(void *Arg_p) { - /* Convert input from (void *) to (ThreadArg *), call the moment and stop the thread. */ - ThreadArg *Arg = (ThreadArg *)Arg_p; - call_moment_sp[Arg->iThread](Arg->nest, Arg->lev); - _endthreadex(0); - return (0); -} -#endif - -#if defined(DO_MULTI_THREAD) || defined(PARALLEL) /* ------------------------------------------------------------------------------ */ int GetLocalNThread(void) { /* Get number of processors from the environment variable NUMBER_OF_PROCESSORS. */ - char *pStr; int localNThread = 1; /* Default */ - + if ((pStr = getenv("NUMBER_OF_PROCESSORS")) != NULL) sscanf(pStr, "%d", &localNThread); - + return (localNThread); } -#else -int GetLocalNThread(void) { - return 1; -} -#endif -#ifdef PARALLEL -/* ------------------------------------------------------------------------------ */ -unsigned __stdcall MT_moment_sp_M_slice(void *Arg_p) { - /* Convert input from (void *) to (ThreadArg *), call the moment and stop the thread. */ - ThreadArg *Arg = (ThreadArg *)Arg_p; - call_moment_sp_M_slice(Arg->lev, Arg->row_start, Arg->row_end, Arg->nest); - _endthreadex(0); - return (0); -} -/* ------------------------------------------------------------------------------ */ -unsigned __stdcall MT_moment_sp_N_slice(void *Arg_p) { - ThreadArg *Arg = (ThreadArg *)Arg_p; - call_moment_sp_N_slice(Arg->lev, Arg->row_start, Arg->row_end, Arg->nest); - _endthreadex(0); - return (0); -} -/* ------------------------------------------------------------------------------ */ -unsigned __stdcall MT_moment_M_slice(void *Arg_p) { - ThreadArg *Arg = (ThreadArg *)Arg_p; - call_moment_M_slice(Arg->lev, Arg->row_start, Arg->row_end, Arg->nest); - _endthreadex(0); - return (0); -} -/* ------------------------------------------------------------------------------ */ -unsigned __stdcall MT_moment_N_slice(void *Arg_p) { - ThreadArg *Arg = (ThreadArg *)Arg_p; - call_moment_N_slice(Arg->lev, Arg->row_start, Arg->row_end, Arg->nest); - _endthreadex(0); - return (0); -} -/* ------------------------------------------------------------------------------ */ -unsigned __stdcall MT_mass_sp(void *Arg_p) { - ThreadArg *Arg = (ThreadArg *)Arg_p; - call_mass_sp(Arg->nest, Arg->lev, Arg->row_start, Arg->row_end); - _endthreadex(0); - return (0); -} -/* ------------------------------------------------------------------------------ */ -unsigned __stdcall MT_mass(void *Arg_p) { - ThreadArg *Arg = (ThreadArg *)Arg_p; - call_mass(Arg->nest, Arg->lev, Arg->row_start, Arg->row_end); - _endthreadex(0); - return (0); -} -#endif /* ---------------------------------------------------------------------------------------- */ void kaba_source(struct GMT_GRID_HEADER hdr, double x_inc, double y_inc, double x_min, double x_max, From 7ed1cc8965b1833ce84d3ed9231c59452b43d4b6 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Fri, 3 Jul 2026 00:41:30 +0100 Subject: [PATCH 22/29] nswing: auto-detect geographic coordinates when -fg is not given Two-stage detection, applied in both the simulation path and the deformation-only mode: first trust GMT's own grid metadata (gmt_M_is_geographic, set by the netCDF reader when the grid carries degree units), then fall back to the old -180/360, -90/90 range heuristic for unit-less grids (e.g. Mirone-written ones). Either way a warning is issued and isGeog is set. The deformation-only mode previously had NO detection at all - it returns before the bathymetry-reading branch where the heuristic lived, so 'nswing -R -F... -G' on a lon/lat region silently ran the Okada deformation in Cartesian, producing wrong results with no warning. Verified against an explicit -fg run: bit-identical now. Co-Authored-By: Claude Fable 5 --- src/nswing/nswing.c | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/src/nswing/nswing.c b/src/nswing/nswing.c index 777c9ba9446..36503cf6b28 100644 --- a/src/nswing/nswing.c +++ b/src/nswing/nswing.c @@ -1226,6 +1226,17 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { hdr_b.z_min = hdr_b.z_max = 0; nm_def = (uint64_t)hdr_b.n_columns * (uint64_t)hdr_b.n_rows; + if (!isGeog) { /* No -fg given. Trust GMT's -R metadata detection first, fall back to a range heuristic */ + if (gmt_M_is_geographic(GMT, GMT_IN)) { + GMT_Report(API, GMT_MSG_WARNING, "NSWING: -R region is in geographic coordinates. Assuming -fg.\n"); + isGeog = true; + } + else if (dfXmin >= -180 && dfXmax <= 360 && dfYmin > -90 && dfYmax < 90) { + GMT_Report(API, GMT_MSG_WARNING, "NSWING: Warning, -R region seams to be in Geographical coords but -f was not set.\n"); + isGeog = true; + } + } + if ((def = (double *)calloc((size_t)nm_def, sizeof(double))) == NULL) {no_sys_mem(API, "(deform)", (unsigned int)nm_def); Return(-1);} @@ -1310,9 +1321,15 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { Return(-1); } - if (!isGeog && hdr_b.wesn[XLO] >= -180 && hdr_b.wesn[XHI] <= 360 && hdr_b.wesn[YLO] > -90 && hdr_b.wesn[YHI] < 90) { - GMT_Report(API, GMT_MSG_WARNING, "NSWING: Warning, the bathymetry grid seams to be in Geographical coords but -f was not set.\n"); - isGeog = true; + if (!isGeog) { /* No -fg given. Trust GMT's grid-metadata detection first, fall back to a range heuristic */ + if (gmt_M_is_geographic(GMT, GMT_IN)) { + GMT_Report(API, GMT_MSG_WARNING, "NSWING: Bathymetry grid metadata says geographic coordinates. Assuming -fg.\n"); + isGeog = true; + } + else if (hdr_b.wesn[XLO] >= -180 && hdr_b.wesn[XHI] <= 360 && hdr_b.wesn[YLO] > -90 && hdr_b.wesn[YHI] < 90) { + GMT_Report(API, GMT_MSG_WARNING, "NSWING: Warning, the bathymetry grid seams to be in Geographical coords but -f was not set.\n"); + isGeog = true; + } } if (!do_Okada && !do_Kaba) { /* Otherwise we will compute initial condition later down after arrays are allocated */ From af333fc1b0e99b1804d662c739c956603cba8252 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Fri, 3 Jul 2026 01:44:18 +0100 Subject: [PATCH 23/29] Update the manual --- doc/rst/source/supplements/nswing/nswing.rst | 77 ++++++++++++-------- 1 file changed, 47 insertions(+), 30 deletions(-) diff --git a/doc/rst/source/supplements/nswing/nswing.rst b/doc/rst/source/supplements/nswing/nswing.rst index 30a9dd01390..9875003f102 100644 --- a/doc/rst/source/supplements/nswing/nswing.rst +++ b/doc/rst/source/supplements/nswing/nswing.rst @@ -24,7 +24,7 @@ Synopsis [ |-E|\ [**p**][**m**][,\ *decim*] ] [ |-F|\ *x_epic/y_epic/dip/strike/rake/slip/length/width/topDepth* ] [ |-F|\ **k**\ [**c**]\ *w/e/s/n* ] -[ |-G|\ *stem*\ [**+l**\ *ev*],\ *int* ] +[ |-G|\ *name*\ [**+m**],\ *int* ] [ |-H|\ [*momentM,momentN*\ [,\ *t*]] ] [ |-P|\ *time_jump*\ [**+t**\ *run_time_jump*] ] [ |-L|\ [*name1,name2*] ] @@ -34,13 +34,11 @@ Synopsis [ |-Q|\ *z_offset* ] [ |SYN_OPT-R| ] [ |-S|\ [**x**\|\ **y**\|\ **n**][**+m**][**+s**] ] -[ |-T|\ *mareg*\ [**+o**\ *outmaregs*][**+t**\ *int*] ] +[ |-T|\ *mareg*\|\ *x/y*\ [**+o**\ *outmaregs*][**+t**\ *int*] ] [ |-X|\ *manning0*\ [,\ *...*] ] -[ |-Z|\ *name*\ [**+l**\ *ev*],\ *int* ] -[ **-i** ] [ |SYN_OPT-V| ] [ |SYN_OPT-f| ] -[ **-x**\ [[-]\ *n*] ] +[ **-x**\ *n* ] [ |SYN_OPT--| ] |No-spaces| @@ -78,9 +76,8 @@ Optional Arguments **-1**\ *bat_lev1* **-2**\ *bat_lev2* ... **-9**\ *bat_lev9* Nested bathymetry grids, one per nesting level. Each level refines the solution inside the area covered by its grid. Warning, the grids must be aligned and have a cell size - that is an integer fraction of the previous level. The first level (level 1) is the base - bathymetry grid, and the last level (level 9) is the finest grid. The number of levels - is determined by the number of **-n** options given. These grids are not trivial to create, + that is an integer fraction of the previous level. The number of levels is determined + by the number of **-1** ... **-9** options given. These grids are not trivial to create, the best way is to use the Mirone TINTOL tool. .. _-A: @@ -109,9 +106,9 @@ Optional Arguments Append **m** to save only a single grid holding the max values. This can noticeably slow the run, so optionally append a *decim* decimator factor after the comma (causes aliasing visible under shaded illumination). The - file name comes from *name* in **-G**/**-Z** complemented with a *_max* + file name comes from *name* in **-G** complemented with a *_max* prefix; saving of multiple grids is then disabled. A 3D netCDF file with - wave heights is still possible with **-Z**. + wave heights is still possible with **-G**. .. _-F: @@ -141,10 +138,12 @@ Optional Arguments .. _-G: -**-G**\ *stem*\ [**+l**\ *ev*],\ *int* - Write grids at the *int* time intervals; files are named *stem#*\ **.grd**. - When doing nested grids, append **+l**\ *ev* to save that particular level - (only one level is allowed). +**-G**\ *name*\ [**+m**],\ *int* + Save the water level every *int* time steps in a single 3D netCDF file + called *name*\ **.nc** (the extension is appended when *name* has none). + Append **+m** to instead write each saved step as a separate grid; files + are then named *name#*\ **.grd**. When doing nested grids the finest level + is the one saved. .. _-H: @@ -176,12 +175,12 @@ Optional Arguments .. _-M: **-M**\ [**-**\|\ **+**\ [*maskname*]] - Write a grid with the max water level (name from *name* in **-Z**, *_max* + Write a grid with the max water level (name from *name* in **-G**, *_max* prefix). Append **-** to instead compute the maximum water retreat, written to a mask file (default *long_beach.grd*; append a name after **-** to change it, e.g. **-M-**\ *beach_long.grd*). Append **+** for a mask with the Run In extent (behaves like **-M-**). **-M** may be repeated, e.g. - **-M -M- -M+** computes all three. With **-Z** the *long* and *short* beach + **-M -M- -M+** computes all three. With **-G** the *long* and *short* beach arrays are also saved in the *.nc* file. .. _-N: @@ -218,14 +217,16 @@ Output grids only in the sub-region enclosed by *west/east/south/north*. .. _-T: -**-T**\ *mareg*\ [**+o**\ *outmaregs*][**+t**\ *int*] +**-T**\ *mareg*\|\ *x/y*\ [**+o**\ *outmaregs*][**+t**\ *int*] Save time series (maregraphs) at virtual tide-gauge locations. *mareg* is - the file with the (x y) locations of the virtual maregraphs. Append - **+o**\ *outmaregs* to set the output file name [Default is + the file with the (x y) locations of the virtual maregraphs. For a single + maregraph the location may be given directly as *x/y* instead of a file + name. Append **+o**\ *outmaregs* to set the output file name [Default is *maregs_out.dat*]. A *.dat* extension is added when *outmaregs* has none; use a *.nc* extension to write the maregraphs as a netCDF file instead. Append **+t**\ *int* to save every *int* simulation time steps (set by - **-t**) [Default is every time step]. + **-t**) [Default is every time step]. **-T** alone (without **-G**) is + allowed and runs a simulation that only outputs the maregraph series. .. _-X: @@ -235,15 +236,6 @@ Output grids only in the sub-region enclosed by *west/east/south/north*. **+**\ *depth* to apply Manning only at depths shallower than *depth* (positive up). -.. _-Z: - -**-Z**\ *name*\ [**+l**\ *ev*],\ *int* - Same as **-G** but saves the result in a 3D netCDF file. - -**-i** - Do not interpolate the initial surface of nested grids, at time zero, from - the mother grids. - .. |Add_-V| replace:: |Add_-V_links| .. include:: ../../explain_-V.rst_ :start-after: **Syntax** @@ -252,6 +244,15 @@ Output grids only in the sub-region enclosed by *west/east/south/north*. .. |Add_-f| unicode:: 0x20 .. just an invisible code .. include:: ../../explain_-f.rst_ + When **-fg** is not given, **nswing** first checks the grid's own metadata + (netCDF degree units) and, failing that, falls back to a coordinate-range + heuristic; if either indicates geographic coordinates a warning is issued + and **-fg** is set implicitly. + +**-x**\ *n* + Number of OpenMP threads to use [Default is all available cores]. Results + are identical for any number of threads. + .. include:: ../../explain_core.rst_ .. include:: ../../explain_help.rst_ @@ -269,7 +270,23 @@ To run the same simulation but generate the source on the fly from Okada fault parameters and sample the wave field at the virtual maregraphs listed in *gauges.dat* every 5 time steps, try:: - gmt nswing bathy.grd dummy.grd -t5 -F-8/37/12/90/90/3/100/50/10 -Tgauges.dat+t5 + gmt nswing bathy.grd -t5 -F-8/37/12/90/90/3/100/50/10 -Tgauges.dat+t5 + +To record only the time series at a single virtual tide gauge (no grids at +all), give its location directly to **-T**:: + + gmt nswing bathy.grd source.grd -t5 -T-10.7/37.3 + +To compute just the Okada co-seismic deformation over the geometry of the grid +*bathy.grd* (no simulation), try:: + + gmt nswing -Rbathy.grd -F94.3/2.8/25/330/90/10/250/65/10 -Gdeform.grd + +A less hypothetical example that generates only a few layers in the tsu.nc cube:: + + gmt grdcut @earth_relief_02m_g -R-15/-7.5/34/39.5 -Gbat.grd + gmt nswing -Rbat.grd -F-12.49593345/35.93634937/25/58.2/90/10/215/53.75/10 -Ginit.grd + gmt nswing bat.grd init.grd -t3 -Gtsu,10 -N100 -V See Also -------- From dec53d7ed7520ab5d6e3f03e7d93bc9356b715dd Mon Sep 17 00:00:00 2001 From: Joaquim Date: Fri, 3 Jul 2026 01:47:57 +0100 Subject: [PATCH 24/29] Another man update --- doc/rst/source/supplements/nswing/nswing.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/rst/source/supplements/nswing/nswing.rst b/doc/rst/source/supplements/nswing/nswing.rst index 9875003f102..529f668f4bb 100644 --- a/doc/rst/source/supplements/nswing/nswing.rst +++ b/doc/rst/source/supplements/nswing/nswing.rst @@ -292,4 +292,4 @@ See Also -------- :doc:`gmt `, -:doc:`grdseamount ` +:doc:`grdinterpolate ` From 2e4efe30d2f95fa3f3221f0d1d6b9bac3bfc0d2a Mon Sep 17 00:00:00 2001 From: Joaquim Date: Tue, 7 Jul 2026 02:02:10 +0100 Subject: [PATCH 25/29] Update the KEYS. It now accepts GMTgrids in input --- src/gmt_api.c | 4 +- src/nswing/nswing.c | 109 +++++++++++++++++++++++++------------------- 2 files changed, 63 insertions(+), 50 deletions(-) diff --git a/src/gmt_api.c b/src/gmt_api.c index 359fb696377..f59f30a9b92 100644 --- a/src/gmt_api.c +++ b/src/gmt_api.c @@ -1629,7 +1629,7 @@ GMT_LOCAL char ** gmtapi_process_keys (void *V_API, const char *string, char typ size_t len, k, kk, n; int o_id = GMT_NOTSET, family = GMT_NOTSET, geometry = GMT_NOTSET; bool change_type = false; - char **s = NULL, *next = NULL, *tmp = NULL, magic = 0, revised[GMT_LEN64] = {""}; + char **s = NULL, *next = NULL, *tmp = NULL, magic = 0, revised[GMT_LEN128] = {""}; struct GMT_OPTION *opt = NULL; struct GMTAPI_CTRL *API = gmtapi_get_api_ptr (V_API); @@ -1795,7 +1795,7 @@ GMT_LOCAL char ** gmtapi_process_keys (void *V_API, const char *string, char typ n = kk; /* May have lost some NULLs. Make a revised string for debug output */ for (k = 0; k < n; k++) { if (k) strcat (revised, ","); - strncat (revised, s[k], GMT_LEN64-1); + strncat (revised, s[k], GMT_LEN128-1); } if (revised[0]) GMT_Report (API, GMT_MSG_DEBUG, "gmtapi_process_keys: Revised keys string is %s\n", revised); *n_items = (unsigned int)n; /* Total number of remaining keys for this module */ diff --git a/src/nswing/nswing.c b/src/nswing/nswing.c index 36503cf6b28..fbd05da7d46 100644 --- a/src/nswing/nswing.c +++ b/src/nswing/nswing.c @@ -49,7 +49,9 @@ #define THIS_MODULE_MODERN_NAME "nswing" #define THIS_MODULE_LIB "supplements" #define THIS_MODULE_PURPOSE "A tsunami maker" -#define THIS_MODULE_KEYS "u, 2->d, 3->r, 4->q, 5->c, 6->s, 7->e, 8->o + so that we can still use -1, -2, etc. in command line and via Julia. */ +#define THIS_MODULE_KEYS "next) { - if (opt->option == 'R') /* Already consumed by GMT_Parse_Common into GMT->common.R; not one of nswing's own options */ - continue; switch (opt->option) { case GMT_OPT_INFILE: /* bathy / source grid (positional, no leading '-') */ if (bathy == NULL) @@ -502,9 +502,6 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC error = true; } break; - case 'c': - add_const = atof(opt->arg); - break; case 'f': /* */ isGeog = true; break; @@ -527,10 +524,7 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC if (opt->arg[0] == 'l') /* Output land nodes in SWW file */ with_land = true; break; - case 'O': /* File with the boundary condition (experimental) */ - bnc_file = opt->arg; - break; - case 'C': /* Output momentum grids */ + case 'C': /* Coriolis */ nest.do_Coriolis = true; if (opt->arg[0]) nest.lat_min4Coriolis = atof(opt->arg); @@ -626,6 +620,7 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC } if ((pch = strstr(stem,"+m")) != NULL) { /* Write grids at grn intervals */ write_grids = true; + pch[0] = '\0'; /* Strip the "+m" part */ } else { out_3D = true; @@ -712,6 +707,9 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC case 'N': /* Number of cycles to compute */ n_of_cycles = atoi(opt->arg); break; + case 'O': /* File with the boundary condition (experimental) */ + bnc_file = opt->arg; + break; case 'Q': /* Vertical offset (simulate tide) */ if (opt->arg[0]) sscanf(opt->arg, "%lf", &z_offset); @@ -829,12 +827,29 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC case 'V': verbose = true; break; - case '1': + case '1': case 'u': nesteds[0] = opt->arg; break; - case '2': case '3': case '4': case '5': - case '6': case '7': case '8': case '9': - nesteds[opt->option - '0' - 1] = opt->arg; + case '2': case 'd': + nesteds[1] = opt->arg; + break; + case '3': case 'r': + nesteds[2] = opt->arg; + break; + case '4': case 'q': + nesteds[3] = opt->arg; + break; + case '5': case 'c': + nesteds[4] = opt->arg; + break; + case '6': case 's': + nesteds[5] = opt->arg; + break; + case '7': case 'e': + nesteds[6] = opt->arg; + break; + case '8': case 'o': + nesteds[7] = opt->arg; break; default: GMT_Report(GMT->parent, GMT_MSG_ERROR, "NSWING: Unknown option -%c%s\n", opt->option, opt->arg); @@ -1036,7 +1051,7 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { int ncid_3D[3], ids_z[10], ids_3D[3], ncid_Mar, ids_Mar[8]; int n_of_cycles = 1010; /* Default number of cycles to compute */ int num_of_nestGrids = 0; /* Number of nesting grids */ - bool bat_in_input = false, source_in_input = false, write_grids = false, isGeog = false; + bool write_grids = false, isGeog = false; bool maregs_in_input = false, out_momentum = false, got_R = false, deform_only = false; bool with_land = false, do_nestum = false, saveNested = false, verbose = false; bool out_velocity = false, out_velocity_x = false, out_velocity_y = false, out_velocity_r = false; @@ -1304,45 +1319,43 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { if (out_momentum && (out_sww || out_most)) out_momentum = false; if ((out_velocity || out_velocity_x || out_velocity_y || out_velocity_r) && (out_sww || out_most)) out_velocity = false; - if (!bat_in_input && !source_in_input) { /* If bathymetry & source where not given as arguments, load them */ - if (!bathy || (!fonte && !bnc_file && !do_Okada && !do_Kaba)) { - GMT_Report(API, GMT_MSG_ERROR, "NSWING: error, bathymetry and/or source grids were not provided.\n"); - Return(-1); - } + if (!bathy || (!fonte && !bnc_file && !do_Okada && !do_Kaba)) { + GMT_Report(API, GMT_MSG_ERROR, "NSWING: error, bathymetry and/or source grids were not provided.\n"); + Return(-1); + } - /* Read base bathymetry through the GMT API (real file or in-memory grid) */ - if ((Gb = gmtnswing_get_grid(API, bathy, &hdr_b)) == NULL) { - GMT_Report(API, GMT_MSG_ERROR, "NSWING: %s Invalid bathymetry grid.\n", bathy); - Return(-1); + /* Read base bathymetry through the GMT API (real file or in-memory grid) */ + if ((Gb = gmtnswing_get_grid(API, bathy, &hdr_b)) == NULL) { + GMT_Report(API, GMT_MSG_ERROR, "NSWING: %s Invalid bathymetry grid.\n", bathy); + Return(-1); + } + if (hdr_b.n_columns < 2 || hdr_b.n_rows < 2) { + GMT_Report(API, GMT_MSG_ERROR, "NSWING: bathymetry grid must have at least 2 rows and 2 columns.\n"); + GMT_Destroy_Data(API, &Gb); + Return(-1); + } + + if (!isGeog) { /* No -fg given. Trust GMT's grid-metadata detection first, fall back to a range heuristic */ + if (gmt_M_is_geographic(GMT, GMT_IN)) { + GMT_Report(API, GMT_MSG_WARNING, "NSWING: Bathymetry grid metadata says geographic coordinates. Assuming -fg.\n"); + isGeog = true; } - if (hdr_b.n_columns < 2 || hdr_b.n_rows < 2) { - GMT_Report(API, GMT_MSG_ERROR, "NSWING: bathymetry grid must have at least 2 rows and 2 columns.\n"); - GMT_Destroy_Data(API, &Gb); - Return(-1); + else if (hdr_b.wesn[XLO] >= -180 && hdr_b.wesn[XHI] <= 360 && hdr_b.wesn[YLO] > -90 && hdr_b.wesn[YHI] < 90) { + GMT_Report(API, GMT_MSG_WARNING, "NSWING: Warning, the bathymetry grid seams to be in Geographical coords but -f was not set.\n"); + isGeog = true; } + } - if (!isGeog) { /* No -fg given. Trust GMT's grid-metadata detection first, fall back to a range heuristic */ - if (gmt_M_is_geographic(GMT, GMT_IN)) { - GMT_Report(API, GMT_MSG_WARNING, "NSWING: Bathymetry grid metadata says geographic coordinates. Assuming -fg.\n"); - isGeog = true; - } - else if (hdr_b.wesn[XLO] >= -180 && hdr_b.wesn[XHI] <= 360 && hdr_b.wesn[YLO] > -90 && hdr_b.wesn[YHI] < 90) { - GMT_Report(API, GMT_MSG_WARNING, "NSWING: Warning, the bathymetry grid seams to be in Geographical coords but -f was not set.\n"); - isGeog = true; - } + if (!do_Okada && !do_Kaba) { /* Otherwise we will compute initial condition later down after arrays are allocated */ + if (!bnc_file && (Gf = gmtnswing_get_grid(API, fonte, &hdr_f)) == NULL) { /* Read source grid */ + GMT_Report(API, GMT_MSG_ERROR, "NSWING: %s Invalid source grid.\n", fonte); + Return(-1); } - if (!do_Okada && !do_Kaba) { /* Otherwise we will compute initial condition later down after arrays are allocated */ - if (!bnc_file && (Gf = gmtnswing_get_grid(API, fonte, &hdr_f)) == NULL) { /* Read source grid */ - GMT_Report(API, GMT_MSG_ERROR, "NSWING: %s Invalid source grid.\n", fonte); - Return(-1); - } - - if (!bnc_file && (hdr_f.n_columns != hdr_b.n_columns || hdr_f.n_rows != hdr_b.n_rows)) { - GMT_Report(API, GMT_MSG_ERROR, "Bathymetry and source grids have different rows/columns\n"); - GMT_Report(API, GMT_MSG_ERROR, "%d %d %d %d\n", hdr_b.n_rows, hdr_f.n_rows, hdr_b.n_columns, hdr_f.n_columns); - error++; - } + if (!bnc_file && (hdr_f.n_columns != hdr_b.n_columns || hdr_f.n_rows != hdr_b.n_rows)) { + GMT_Report(API, GMT_MSG_ERROR, "Bathymetry and source grids have different rows/columns\n"); + GMT_Report(API, GMT_MSG_ERROR, "%d %d %d %d\n", hdr_b.n_rows, hdr_f.n_rows, hdr_b.n_columns, hdr_f.n_columns); + error++; } } From 7021eacc1deae968893e0223e7ded01a83c875f6 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Tue, 7 Jul 2026 10:38:48 +0100 Subject: [PATCH 26/29] Fix option -T such that we can pass it a matrix via externals --- src/nswing/nswing.c | 63 ++++++++++++++++++++++++++++++--------------- 1 file changed, 42 insertions(+), 21 deletions(-) diff --git a/src/nswing/nswing.c b/src/nswing/nswing.c index fbd05da7d46..61d296c8af0 100644 --- a/src/nswing/nswing.c +++ b/src/nswing/nswing.c @@ -51,7 +51,7 @@ #define THIS_MODULE_PURPOSE "A tsunami maker" /* We can't use 1G(, etc because -1 is taken leterally, so julia changes: 1->u, 2->d, 3->r, 4->q, 5->c, 6->s, 7->e, 8->o so that we can still use -1, -2, etc. in command line and via Julia. */ -#define THIS_MODULE_KEYS "n_records; - GMT_Destroy_Data(API, &D); + if (D_out) + *D_out = D; + else + GMT_Destroy_Data(API, &D); return (n); } /* -------------------------------------------------------------------- */ -int read_maregs(void *API, struct GMT_GRID_HEADER hdr, char *file, unsigned int *lcum_p, char *names[]) { - /* Read maregraph positions (real file or virtual dataset) and convert them to vector linear indices */ +int read_maregs(void *API, struct GMT_GRID_HEADER hdr, char *file, unsigned int *lcum_p, char *names[], struct GMT_DATASET *D_in) { + /* Read maregraph positions (real file or virtual dataset) and convert them to vector linear indices. + * If D_in is not NULL it is a dataset already read by count_n_maregs() and is reused here instead of + * re-reading the file (a GMT virtual-file reference can only be read once). */ int i = 0, ix, jy; uint64_t tbl, seg, row; double x, y; - struct GMT_DATASET *D = NULL; + struct GMT_DATASET *D = D_in; struct GMT_DATASEGMENT *S = NULL; - if ((D = GMT_Read_Data(API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_READ_NORMAL, NULL, file, NULL)) == NULL) { + if (D == NULL && (D = GMT_Read_Data(API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_READ_NORMAL, NULL, file, NULL)) == NULL) { GMT_Report(API, GMT_MSG_ERROR, "NSWING: Unable to open file %s - exiting\n", file); return (-1); } @@ -2692,8 +2708,11 @@ int read_maregs(void *API, struct GMT_GRID_HEADER hdr, char *file, unsigned int for (row = 0; row < S->n_rows; row++) { x = S->data[GMT_X][row]; y = S->data[GMT_Y][row]; - if (x < hdr.wesn[XLO] || x > hdr.wesn[XHI] || y < hdr.wesn[YLO] || y > hdr.wesn[YHI]) + if (x < hdr.wesn[XLO] || x > hdr.wesn[XHI] || y < hdr.wesn[YLO] || y > hdr.wesn[YHI]) { + GMT_Report(API, GMT_MSG_WARNING, "NSWING: maregraph position %g/%g in %s falls outside grid W/E/S/N = %g/%g/%g/%g - skipped\n", + x, y, file, hdr.wesn[XLO], hdr.wesn[XHI], hdr.wesn[YLO], hdr.wesn[YHI]); continue; + } ix = irint((x - hdr.wesn[XLO]) / hdr.inc[GMT_X]); jy = irint((y - hdr.wesn[YLO]) / hdr.inc[GMT_Y]); lcum_p[i] = jy * hdr.n_columns + ix; @@ -2707,15 +2726,17 @@ int read_maregs(void *API, struct GMT_GRID_HEADER hdr, char *file, unsigned int } /* -------------------------------------------------------------------- */ -int read_tracers(void *API, struct GMT_GRID_HEADER hdr, char *file, struct tracers *oranges) { - /* Read tracers positions (real file or virtual dataset) */ +int read_tracers(void *API, struct GMT_GRID_HEADER hdr, char *file, struct tracers *oranges, struct GMT_DATASET *D_in) { + /* Read tracers positions (real file or virtual dataset). If D_in is not NULL it is a dataset + * already read by count_n_maregs() and is reused here instead of re-reading the file (a GMT + * virtual-file reference can only be read once). */ int i = 0; uint64_t tbl, seg, row; double x, y; - struct GMT_DATASET *D = NULL; + struct GMT_DATASET *D = D_in; struct GMT_DATASEGMENT *S = NULL; - if ((D = GMT_Read_Data(API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_READ_NORMAL, NULL, file, NULL)) == NULL) { + if (D == NULL && (D = GMT_Read_Data(API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_READ_NORMAL, NULL, file, NULL)) == NULL) { GMT_Report(API, GMT_MSG_ERROR, "NSWING: Unable to open file %s - exiting\n", file); return (-1); } From 90f3477e3ba821932a497504b95de4688087789d Mon Sep 17 00:00:00 2001 From: Joaquim Date: Tue, 7 Jul 2026 19:21:20 +0100 Subject: [PATCH 27/29] nswing: add +a modifier to -E/-S so energy/velocity can append to the 3D cube Previously -E's Energy/Power silently overwrote the z (sea-surface) variable in the -G 3D netCDF cube, mislabeled with z's own name/units ("Sea surface [meters]"), while -S always wrote velocity alongside z with no way to opt out. Add +a to both: absent, the option's quantity is the cube's only variable (correctly named/labeled); with +a, it's written as an extra channel next to z, matching -S's old default behavior. Note: -S without +a now omits the z variable (previously always included). Co-Authored-By: Claude Sonnet 5 --- doc/rst/source/supplements/nswing/nswing.rst | 20 ++- src/nswing/nswing.c | 172 ++++++++++++++----- 2 files changed, 140 insertions(+), 52 deletions(-) diff --git a/doc/rst/source/supplements/nswing/nswing.rst b/doc/rst/source/supplements/nswing/nswing.rst index 529f668f4bb..38946f334bb 100644 --- a/doc/rst/source/supplements/nswing/nswing.rst +++ b/doc/rst/source/supplements/nswing/nswing.rst @@ -21,7 +21,7 @@ Synopsis [ |-A|\ *fname.sww* ] [ |-C| ] [ |-D| ] -[ |-E|\ [**p**][**m**][,\ *decim*] ] +[ |-E|\ [**p**][**m**][**+a**][,\ *decim*] ] [ |-F|\ *x_epic/y_epic/dip/strike/rake/slip/length/width/topDepth* ] [ |-F|\ **k**\ [**c**]\ *w/e/s/n* ] [ |-G|\ *name*\ [**+m**],\ *int* ] @@ -33,7 +33,7 @@ Synopsis [ |-O|\ *BCfile* ] [ |-Q|\ *z_offset* ] [ |SYN_OPT-R| ] -[ |-S|\ [**x**\|\ **y**\|\ **n**][**+m**][**+s**] ] +[ |-S|\ [**x**\|\ **y**\|\ **n**][**+m**][**+s**][**+a**] ] [ |-T|\ *mareg*\|\ *x/y*\ [**+o**\ *outmaregs*][**+t**\ *int*] ] [ |-X|\ *manning0*\ [,\ *...*] ] [ |SYN_OPT-V| ] @@ -101,14 +101,18 @@ Optional Arguments .. _-E: -**-E**\ [**p**][**m**][,\ *decim*] +**-E**\ [**p**][**m**][**+a**][,\ *decim*] Write grids with energy, or with power if **p** is appended (**-Ep**). Append **m** to save only a single grid holding the max values. This can noticeably slow the run, so optionally append a *decim* decimator factor after the comma (causes aliasing visible under shaded illumination). The file name comes from *name* in **-G** complemented with a *_max* - prefix; saving of multiple grids is then disabled. A 3D netCDF file with - wave heights is still possible with **-G**. + prefix; saving of multiple grids is then disabled. + + With **-G**'s 3D netCDF cube, Energy/Power replaces the sea-surface + (*z*) variable unless **+a** is appended, in which case Energy/Power is + written as an extra variable alongside *z* (same idea as **-S**'s + **+a**). .. _-F: @@ -207,7 +211,7 @@ Output grids only in the sub-region enclosed by *west/east/south/north*. .. _-S: -**-S**\ [**x**\|\ **y**\|\ **n**][**+m**][**+s**] +**-S**\ [**x**\|\ **y**\|\ **n**][**+m**][**+s**][**+a**] Write grids with the velocity (names get *_U* and *_V* suffixes). Use **x** or **y** to save only one component, or **n** for no velocity grids (maregraphs only). Append **+m** to also write velocity (vx,vy) at maregraph @@ -215,6 +219,10 @@ Output grids only in the sub-region enclosed by *west/east/south/north*. (|v|) (*_max_speed* suffix). Use the **n** flag to NOT output the U and V components, e.g. **-Sn+s**. + With **-G**'s 3D netCDF cube, velocity is written as the only variable + unless **+a** is appended, in which case the sea-surface (*z*) variable + is also written (*z* + velocity). + .. _-T: **-T**\ *mareg*\|\ *x/y*\ [**+o**\ *outmaregs*][**+t**\ *int*] diff --git a/src/nswing/nswing.c b/src/nswing/nswing.c index 61d296c8af0..ff5201ad5ac 100644 --- a/src/nswing/nswing.c +++ b/src/nswing/nswing.c @@ -51,7 +51,7 @@ #define THIS_MODULE_PURPOSE "A tsunami maker" /* We can't use 1G(, etc because -1 is taken leterally, so julia changes: 1->u, 2->d, 3->r, 4->q, 5->c, 6->s, 7->e, 8->o so that we can still use -1, -2, etc. in command line and via Julia. */ -#define THIS_MODULE_KEYS "] [-2] [-3<...>] -G[+m], " - "[-A] [-C] [-D] [-E[p][m][,decim]] [-Fx_epic/y_epic/dip/strike/rake/slip/length/width/topDepth] " + "[-A] [-C] [-D] [-E[p][m][+a][,decim]] [-Fx_epic/y_epic/dip/strike/rake/slip/length/width/topDepth] " "[-Fk[c]] [-H] [-H[,t]] [-P] [-L[name1,name2]] " - "[-M[-|+[]]] [-N] [-O] [%s] [-S[x|y|n][+m][+s]] [-T|[+o][+t]] " + "[-M[-|+[]]] [-N] [-O] [%s] [-S[x|y|n][+m][+s][+a]] [-T|[+o][+t]] " "[-Q] [-X] -t
[%s] [-x] [%s]\n", name, GMT_Rgeo_OPT, GMT_V_OPT, GMT_f_OPT); if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS); @@ -341,14 +345,15 @@ GMT_LOCAL int usage(struct GMTAPI_CTRL *API, int level) { GMT_Usage(API, 1, "\n-C Add Coriolis effect."); GMT_Usage(API, 1, "\n-D Write grids with the total water depth."); GMT_Usage(API, -2, "These grids will have wave height on ocean and water thickness on land."); - GMT_Usage(API, 1, "\n-E[p][m][,decim] Write grids with energy or power (-Ep)."); - GMT_Usage(API, -2, "Append 'm' to save only one grid with the max values. This can noticeably slow the run, " - "so optionally append a decimator factor after the comma (causes aliasing visible on shaded illumination). " - "The file name comes from in -G complemented with a '_max' prefix; saving of multiple grids is disabled. " - "A 3D netCDF file with wave heights is still possible with -G."); - GMT_Usage(API, 1, "\n-F Okada fault parameters."); - GMT_Usage(API, -2, "x_epic, y_epic are the X and Y coordinates of the beginning of the fault trace; Dip, Azimuth, " - "Rake, Slip(m), length, width and depth from sea-bottom follow. All dimensions must be in km. " + GMT_Usage(API, 1, "\n-E[p][m][+a][,decim]"); + GMT_Usage(API, -2, "Write grids with energy or power (-Ep). Append 'm' to save only one grid with the max values. " + "This can noticeably slow the run, so optionally append a decimator factor after the comma (causes aliasing " + "visible on shaded illumination). The file name comes from in -G complemented with a '_max' prefix; " + "saving of multiple grids is disabled. With -G's 3D netCDF file, Energy/Power replaces the sea-surface " + "variable unless +a is appended, in which case it is written as an extra variable alongside it (like -S does)."); + GMT_Usage(API, 1, "\n-F"); + GMT_Usage(API, -2, "Okada fault parameters. x_epic, y_epic are the X and Y coordinates of the beginning of the fault trace;" + " Dip, Azimuth, Rake, Slip(m), length, width and depth from sea-bottom follow. All dimensions must be in km. " "If no bathymetry grid is given, -F (or -Fk) together with -R and -G computes the deformation over -R's " "grid geometry (e.g. -R) and writes it straight to -G's name: no simulation is run, -t and -G's " "saving interval are not needed."); @@ -372,18 +377,22 @@ GMT_LOCAL int usage(struct GMTAPI_CTRL *API, int level) { GMT_Usage(API, -2, "-L,"); GMT_Usage(API, -2, "Do Lagrangian tracers, where is the tracers initial-position file " "and the file to hold the results."); - GMT_Usage(API, 1, "\n-M Write a grid with the max water level (name from in -G, '_max' prefix)."); - GMT_Usage(API, -2, "Append '-' to compute instead the maximum water retreat, written to a mask file (default 'long_beach.grd'; " + GMT_Usage(API, 1, "\n-M"); + GMT_Usage(API, -2, " Write a grid with the max water level (name from in -G, '_max' prefix). " + "Append '-' to compute instead the maximum water retreat, written to a mask file (default 'long_beach.grd'; " "append a name after '-' to change it, e.g. -M-beach_long.grd). Append '+' for a mask with the Run In extent (behaves like -M-). " "-M may be repeated, e.g. -M -M- -M+ computes all three. With -G the 'long' and 'short' beach arrays are also saved in the .nc file."); - GMT_Usage(API, 1, "\n-N Number of cycles [Default 1010]."); + GMT_Usage(API, 1, "\n-N"); + GMT_Usage(API, -2, "Number of ccycles [Default 1010]."); GMT_Usage(API, 1, "\n-O Name of a BoundaryCondition ASCII file (experimental)."); GMT_Option(API, "R"); GMT_Usage(API, -2, "Output grids only in the sub-region enclosed by ."); - GMT_Usage(API, 1, "\n-S Write grids with the velocity (names get _U and _V suffixes)."); - GMT_Usage(API, -2, "Use x or y to save only one component, or n for no velocity grids (maregs only). Append +m to also write " + GMT_Usage(API, 1, "\n-S[x|y|n][+m][+s][+a]"); + GMT_Usage(API, -2, " Write grids with the velocity (names get _U and _V suffixes). Use x or y to save " + "only one component, or n for no velocity grids (maregs only). Append +m to also write " "velocity (vx,vy) at maregraph locations (needs -T). Append +s to write the max speed (|v|) ('_max_speed' suffix). " - "Use the 'n' flag to NOT output the U and V components, e.g. -Sn+s."); + "Use the 'n' flag to NOT output the U and V components, e.g. -Sn+s. With -G's 3D netCDF file, velocity is written " + "as the only variable unless +a is appended, in which case the sea-surface variable is also written (z + velocity)."); GMT_Usage(API, 1, "\n-T|[+o][+t] Save maregraph (virtual tide-gauge) time series."); GMT_Usage(API, -2, " is the file with the (x y) locations of the virtual maregraphs. For a single maregraph the " "location may be given directly as instead of a file name. Append +o to set the " @@ -392,10 +401,12 @@ GMT_LOCAL int usage(struct GMTAPI_CTRL *API, int level) { "steps (set by -t) [Default is every time step]."); GMT_Usage(API, 1, "\n-Q"); GMT_Usage(API, -2, "Apply a vertical offset to ALL bathymetry grids (e.g. to simulate tide)."); - GMT_Usage(API, 1, "\n-X]> Manning friction coefficients."); - GMT_Usage(API, -2, "If only one is provided, use it for all nesting levels; otherwise specify one per level, comma separated. " - "Append + to apply Manning only at depths shallower than (positive up)."); - GMT_Usage(API, 1, "\n-x Number of cores to use in a parallel run [Default is the max in the machine]."); + GMT_Usage(API, 1, "\n-X]>"); + GMT_Usage(API, -2, "Manning friction coefficients. If only one is provided, use it for all nesting levels; " + "otherwise specify one per level, comma separated. Append + to apply Manning only at depths " + "shallower than (positive up)."); + GMT_Usage(API, 1, "\n-x"); + GMT_Usage(API, -2, "Number of cores to use in a parallel run [Default is the max in the machine]."); GMT_Option(API, "V,f,."); return (GMT_MODULE_USAGE); @@ -412,6 +423,7 @@ struct NSWING_CTRL { bool maregs_in_input, out_momentum, got_R, deform_only, with_land, saveNested; bool verbose, out_velocity, out_velocity_x, out_velocity_y, out_velocity_r, out_maregs_velocity; bool cumpt, do_2Dgrids, do_maxs, mareg_xy; + bool append_z; char *bathy, *fonte, *fname_sww, *basename_most, *bnc_file; char *nesteds[10]; char hcum[256]; @@ -463,6 +475,7 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC bool with_land = false, saveNested = false, verbose = false; bool out_velocity = false, out_velocity_x = false, out_velocity_y = false, out_velocity_r = false; bool out_maregs_velocity = false; + bool E_append = false, S_append = false; char *bathy = NULL; char hcum[256] = ""; char maregs[256] = ""; @@ -538,6 +551,8 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC EPS4 = atof(&opt->arg[4]); break; } + if (strstr(opt->arg, "+a")) /* Append to (rather than replace) the -G 3D cube's sea-surface variable */ + E_append = true; if (opt->arg[0] == 'p') { if (opt->arg[1] == 'm') max_power = true; @@ -717,6 +732,8 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC break; case 'S': /* Output velocity grids */ strncpy(str_tmp, opt->arg, sizeof(str_tmp)-1); str_tmp[sizeof(str_tmp)-1] = '\0'; + if (strstr(str_tmp, "+a")) /* Append to (rather than be the only variable in) the -G 3D cube */ + S_append = true; if ((pch = strstr(str_tmp,"+m")) != NULL) { /* Velocity at maregraphs */ out_maregs_velocity = true; out_velocity_x = out_velocity_y = true; @@ -874,7 +891,7 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC do_maxs = (max_level || max_energy || max_power); do_2Dgrids = (write_grids || out_velocity || out_velocity_x || out_velocity_y || out_velocity_r || out_momentum - || max_level || max_velocity || max_energy || out_power || max_power || nest.do_long_beach + || max_level || max_velocity || max_energy || out_energy || out_power || max_power || nest.do_long_beach || nest.do_short_beach); if (!deform_only && !(do_2Dgrids || out_sww || out_most || out_3D || cumpt)) { @@ -942,6 +959,7 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC Ctrl->out_energy = out_energy; Ctrl->max_energy = max_energy; Ctrl->out_power = out_power; + Ctrl->append_z = E_append || S_append; Ctrl->max_power = max_power; Ctrl->out_sww = out_sww; Ctrl->out_most = out_most; @@ -1049,7 +1067,7 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { bool do_HotStart = false; /* For when doing a Hot Start */ int n_arg_no_char = 0; int ncid, ncid_most[3], z_id = -1, ids[13], ids_ha[6], ids_ua[6], ids_va[6], ids_most[3]; - int ncid_3D[3], ids_z[10], ids_3D[3], ncid_Mar, ids_Mar[8]; + int ncid_3D[3], ids_z[10], ids_3D[4], ncid_Mar, ids_Mar[8]; int n_of_cycles = 1010; /* Default number of cycles to compute */ int num_of_nestGrids = 0; /* Number of nesting grids */ bool write_grids = false, isGeog = false; @@ -1057,6 +1075,7 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { bool with_land = false, do_nestum = false, saveNested = false, verbose = false; bool out_velocity = false, out_velocity_x = false, out_velocity_y = false, out_velocity_r = false; bool out_maregs_velocity = false; + bool append_z = false, wrote_z = true; int KbGridCols = 1, KbGridRows = 1; /* Number of rows & columns IF computing a grid of 'Kabas' */ int cntKabas = 0; /* Counter of the number of Kabas (prisms) already processed */ int n_mareg, n_ptmar, n_oranges, n_oranges_alloc = 0; @@ -1098,7 +1117,7 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { double z_offset = 0; /* To apply to bathymetry to simulate a tide */ double manning[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; /* Manning coefficients */ - double actual_range[6] = {1e30, -1e30, 1e30, -1e30, 1e30, -1e30}; + double actual_range[8] = {1e30, -1e30, 1e30, -1e30, 1e30, -1e30, 1e30, -1e30}; float stage_range[2], xmom_range[2], ymom_range[2], *tmp_slice; struct GMT_GRID_HEADER hdr_b, hdr_f, hdr_mM, hdr_mN; struct GMT_GRID_HEADER hdr; @@ -1146,6 +1165,7 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { out_energy = Ctrl->out_energy; max_energy = Ctrl->max_energy; out_power = Ctrl->out_power; + append_z = Ctrl->append_z; max_power = Ctrl->max_power; out_sww = Ctrl->out_sww; out_most = Ctrl->out_most; @@ -1462,6 +1482,10 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { nest.out_velocity_x = out_velocity_x; nest.out_velocity_y = out_velocity_y; nest.out_momentum = out_momentum; + nest.out_energy = out_energy; + nest.out_power = out_power; + nest.append_z = append_z; + wrote_z = !((out_velocity_x || out_velocity_y) && !append_z); /* -S without +a omits the sea-surface variable */ nest.isGeog = isGeog; nest.writeLevel = writeLevel; if (initialize_nestum(API, &nest, isGeog, 0)) {free_arrays(&nest, isGeog, num_of_nestGrids); Return(-1);} @@ -1681,16 +1705,20 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { tmp_slice = (float *)malloc(sizeof(float) * (nx * ny)); /* To use inside slice writing */ } else if (out_3D) { + char *var3D = "z"; + if (!append_z && out_energy) var3D = "Energy"; /* -E without +a: Energy/Power replaces the z variable */ + else if (!append_z && out_power) var3D = "Power"; nx = nest.hdr[writeLevel].n_columns; ny = nest.hdr[writeLevel].n_rows; - ncid_3D[0] = open_most_nc(API, &nest, work, fname3D, "z", history, ids_z, nx, ny, xMinOut, yMinOut, false, writeLevel); + ncid_3D[0] = open_most_nc(API, &nest, work, fname3D, var3D, history, ids_z, nx, ny, xMinOut, yMinOut, false, writeLevel); if (ncid_3D[0] == -1) { GMT_Report(API, GMT_MSG_ERROR, "NSWING: failure to create netCDF file\n"); Return(-1); } - ids_3D[0] = ids_z[3]; /* ID of z vriable */ + ids_3D[0] = ids_z[3]; /* ID of z (or Energy/Power, if it replaced z) vriable */ ids_3D[1] = ids_z[5]; /* ID of Vx vriable (only used when it exists) */ ids_3D[2] = ids_z[6]; /* ID of Vy vriable (only used when it exists) */ + ids_3D[3] = ids_z[9]; /* ID of the appended Energy/Power vriable (only used when it exists) */ /* To be used when writing the data slices */ count1_M[0] = 1; count1_M[1] = ny; count1_M[2] = nx; @@ -2026,10 +2054,12 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { } } - if (out_energy) - total_energy(&nest, work, writeLevel); - else if (out_power) - power(&nest, work, writeLevel); + if (!(append_z && out_3D)) { /* With -E+a on the 3D cube, Energy/Power gets its own channel instead (write_most_slice) */ + if (out_energy) + total_energy(&nest, work, writeLevel); + else if (out_power) + power(&nest, work, writeLevel); + } if (write_grids) { sprintf(prenome, "%s%05d.grd", stem, irint(time_h)); @@ -2136,11 +2166,14 @@ EXTERN_MSC int GMT_nswing(void *V_API, int mode, void *args) { err_trap(API, nc_close(ncid_most[2])); } else if (out_3D) { /* Uppdate range values and close 3D file */ - err_trap(API, nc_put_att_double(ncid_3D[0], ids_z[3], "actual_range", NC_DOUBLE, 2U, actual_range)); + if (wrote_z) + err_trap(API, nc_put_att_double(ncid_3D[0], ids_z[3], "actual_range", NC_DOUBLE, 2U, actual_range)); if (out_velocity_x) err_trap(API, nc_put_att_double(ncid_3D[0], ids_z[5], "actual_range", NC_DOUBLE, 2U, &actual_range[2])); if (out_velocity_y) err_trap(API, nc_put_att_double(ncid_3D[0], ids_z[6], "actual_range", NC_DOUBLE, 2U, &actual_range[4])); + if (append_z && (out_energy || out_power)) + err_trap(API, nc_put_att_double(ncid_3D[0], ids_z[9], "actual_range", NC_DOUBLE, 2U, &actual_range[6])); /* It's ugly to have this here but I have no means to tell write_most_slice() when it's the last call */ if (nest.do_long_beach || nest.do_short_beach) { /* Write the mask(s) */ @@ -2991,6 +3024,11 @@ int open_most_nc(void *API, struct nestContainer *nest, float *work, char *base, unsigned int m, n, ij; float dummy = -1e34f; double *x, *y; + /* -S without +a omits the sea-surface variable entirely (velocity/momentum keep their own dedicated + channels regardless); in every other case the primary variable (named z, or Energy/Power when -E + replaced it) is written. */ + int write_z = !((nest->out_velocity_x || nest->out_velocity_y) && !nest->append_z); + int write_energy_chan = (nest->append_z && (nest->out_energy || nest->out_power)); basename = (char *)malloc(strlen(base) + 8); /* +8: room for the NUL and a "_ha.nc" type suffix */ strcpy(basename, base); @@ -3013,6 +3051,14 @@ int open_most_nc(void *API, struct nestContainer *nest, float *work, char *base, long_name = "Sea surface"; units = "meters"; } + else if (!strcmp(name_var,"Energy")) { + long_name = "Total Wave Energy"; + units = "Joules/m^2"; + } + else if (!strcmp(name_var,"Power")) { + long_name = "Wave Power"; + units = "Watts/m"; + } if ((status = nc_create(basename, NC_NETCDF4, &ncid)) != NC_NOERR) { GMT_Report(API, GMT_MSG_ERROR, "NSWING: Unable to create file %s - exiting\n", basename); @@ -3037,7 +3083,8 @@ int open_most_nc(void *API, struct nestContainer *nest, float *work, char *base, } else { err_trap(API, nc_def_var(ncid, "time", NC_DOUBLE,1, &dim0[2], &ids[2])); - err_trap(API, nc_def_var(ncid, name_var, NC_FLOAT, 3, dim3, &ids[3])); + if (write_z) + err_trap(API, nc_def_var(ncid, name_var, NC_FLOAT, 3, dim3, &ids[3])); if (nest->out_momentum) { err_trap(API, nc_def_var(ncid, "Mlon", NC_FLOAT,3, dim3, &ids[5])); err_trap(API, nc_def_var(ncid, "Mlat", NC_FLOAT,3, dim3, &ids[6])); @@ -3046,6 +3093,8 @@ int open_most_nc(void *API, struct nestContainer *nest, float *work, char *base, err_trap(API, nc_def_var(ncid, "Vlon", NC_FLOAT,3, dim3, &ids[5])); if (nest->out_velocity_y) err_trap(API, nc_def_var(ncid, "Vlat", NC_FLOAT,3, dim3, &ids[6])); + if (write_energy_chan) + err_trap(API, nc_def_var(ncid, (nest->out_energy ? "Energy" : "Power"), NC_FLOAT,3, dim3, &ids[9])); dim3[0] = dim0[1]; dim3[1] = dim0[0]; /* Bathym array is rank 2 */ err_trap(API, nc_def_var(ncid, "bathymetry",NC_FLOAT,2, dim3, &ids[4])); } @@ -3066,7 +3115,8 @@ int open_most_nc(void *API, struct nestContainer *nest, float *work, char *base, } else { err_trap(API, nc_def_var(ncid, "time", NC_DOUBLE,1, &dim0[2], &ids[2])); - err_trap(API, nc_def_var(ncid, name_var,NC_FLOAT,3, dim3, &ids[3])); + if (write_z) + err_trap(API, nc_def_var(ncid, name_var,NC_FLOAT,3, dim3, &ids[3])); if (nest->out_momentum) { err_trap(API, nc_def_var(ncid, "Mx",NC_FLOAT,3, dim3, &ids[5])); err_trap(API, nc_def_var(ncid, "My",NC_FLOAT,3, dim3, &ids[6])); @@ -3075,6 +3125,8 @@ int open_most_nc(void *API, struct nestContainer *nest, float *work, char *base, err_trap(API, nc_def_var(ncid, "Vx",NC_FLOAT,3, dim3, &ids[5])); if (nest->out_velocity_y) err_trap(API, nc_def_var(ncid, "Vy",NC_FLOAT,3, dim3, &ids[6])); + if (write_energy_chan) + err_trap(API, nc_def_var(ncid, (nest->out_energy ? "Energy" : "Power"), NC_FLOAT,3, dim3, &ids[9])); dim3[0] = dim0[1]; dim3[1] = dim0[0]; /* Bathym array is rank 2 */ err_trap(API, nc_def_var(ncid, "bathymetry",NC_FLOAT,2, dim3, &ids[4])); } @@ -3087,8 +3139,11 @@ int open_most_nc(void *API, struct nestContainer *nest, float *work, char *base, err_trap(API, nc_def_var(ncid, "ShortBeach", NC_UBYTE, 2, dim3, &ids[8])); } - /* Set a deflation level of 5 (4 zero based) and shuffle for z variable */ - id = (isMost) ? 5 : 3; + /* Set a deflation level of 5 (4 zero based) and shuffle for the primary variable */ + if (isMost) id = 5; + else if (write_z) id = 3; + else if (nest->out_velocity_x) id = 5; + else id = 6; err_trap(API, nc_def_var_deflate(ncid, ids[id], 1, 1, 4)); /* ---- Variables Attributes --------- */ @@ -3128,11 +3183,13 @@ int open_most_nc(void *API, struct nestContainer *nest, float *work, char *base, err_trap(API, nc_put_att_text(ncid, ids[1], "units", 6, "meters")); } err_trap(API, nc_put_att_text (ncid, ids[2], "units", 7, "Seconds")); - err_trap(API, nc_put_att_text (ncid, ids[3], "long_name", strlen(long_name), long_name)); - err_trap(API, nc_put_att_text (ncid, ids[3], "units", strlen(units), units)); - err_trap(API, nc_put_att_float (ncid, ids[3], "missing_value", NC_FLOAT, 1, &nan)); - err_trap(API, nc_put_att_float (ncid, ids[3], "_FillValue", NC_FLOAT, 1, &nan)); - err_trap(API, nc_put_att_double(ncid, ids[3], "actual_range", NC_DOUBLE, 2U, dummy)); + if (write_z) { + err_trap(API, nc_put_att_text (ncid, ids[3], "long_name", strlen(long_name), long_name)); + err_trap(API, nc_put_att_text (ncid, ids[3], "units", strlen(units), units)); + err_trap(API, nc_put_att_float (ncid, ids[3], "missing_value", NC_FLOAT, 1, &nan)); + err_trap(API, nc_put_att_float (ncid, ids[3], "_FillValue", NC_FLOAT, 1, &nan)); + err_trap(API, nc_put_att_double(ncid, ids[3], "actual_range", NC_DOUBLE, 2U, dummy)); + } err_trap(API, nc_put_att_text (ncid, ids[4], "long_name", 10, "bathymetry")); err_trap(API, nc_put_att_text (ncid, ids[4], "units", 6, "meters")); @@ -3177,6 +3234,15 @@ int open_most_nc(void *API, struct nestContainer *nest, float *work, char *base, err_trap(API, nc_put_att_float (ncid, ids[6], "_FillValue", NC_FLOAT, 1, &nan)); err_trap(API, nc_put_att_double(ncid, ids[6], "actual_range", NC_DOUBLE, 2U, dummy)); } + if (write_energy_chan) { /* -E +a: Energy/Power appended alongside z, 3D case */ + long_name = nest->out_energy ? "Total Wave Energy" : "Wave Power"; + units = nest->out_energy ? "Joules/m^2" : "Watts/m"; + err_trap(API, nc_put_att_text (ncid, ids[9], "long_name", strlen(long_name), long_name)); + err_trap(API, nc_put_att_text (ncid, ids[9], "units", strlen(units), units)); + err_trap(API, nc_put_att_float (ncid, ids[9], "missing_value", NC_FLOAT, 1, &nan)); + err_trap(API, nc_put_att_float (ncid, ids[9], "_FillValue", NC_FLOAT, 1, &nan)); + err_trap(API, nc_put_att_double(ncid, ids[9], "actual_range", NC_DOUBLE, 2U, dummy)); + } if (nest->do_long_beach) { float act_range[2] = {0, 0}; @@ -3229,15 +3295,18 @@ void write_most_slice(void *API, struct nestContainer *nest, int *ncid, int *ids double *slice_range, int isMost, int lev) { /* Write a slice of _ha.nc, _va.nc & _ua.nc MOST netCDF files */ unsigned int col, row, n, ij, k; + int write_z = !((nest->out_velocity_x || nest->out_velocity_y) && !nest->append_z); if (!isMost) { - /* ----------- Find the min/max of this slice --------- */ - for (ij = 0; ij < nest->hdr[lev].nm; ij++) { - slice_range[0] = MIN(work[ij], slice_range[0]); - slice_range[1] = MAX(work[ij], slice_range[1]); - } + if (write_z) { + /* ----------- Find the min/max of this slice --------- */ + for (ij = 0; ij < nest->hdr[lev].nm; ij++) { + slice_range[0] = MIN(work[ij], slice_range[0]); + slice_range[1] = MAX(work[ij], slice_range[1]); + } - err_trap(API, nc_put_vara_float(ncid[0], ids[0], start, count, work)); + err_trap(API, nc_put_vara_float(ncid[0], ids[0], start, count, work)); + } /* Conditionally write the Vx & Vy velocity components */ if (nest->out_velocity_x) { @@ -3276,6 +3345,17 @@ void write_most_slice(void *API, struct nestContainer *nest, int *ncid, int *ids } err_trap(API, nc_put_vara_float (ncid[0], ids[2], start, count, work)); } + if (nest->append_z && (nest->out_energy || nest->out_power)) { /* -E +a: Energy/Power as an extra channel */ + if (nest->out_energy) + total_energy(nest, work, nest->writeLevel); + else + power(nest, work, nest->writeLevel); + for (ij = 0; ij < nest->hdr[nest->writeLevel].nm; ij++) { + slice_range[6] = MIN(work[ij], slice_range[6]); + slice_range[7] = MAX(work[ij], slice_range[7]); + } + err_trap(API, nc_put_vara_float (ncid[0], ids[3], start, count, work)); + } } else { for (n = 0; n < 3; n++) { /* Loop over, Amplitude, Xmomentum & Ymomentum */ From ae431890e36313c63fc36233b1106381c8801b7d Mon Sep 17 00:00:00 2001 From: Joaquim Date: Tue, 7 Jul 2026 22:55:47 +0100 Subject: [PATCH 28/29] Restore parsing of option -R --- src/nswing/nswing.c | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/nswing/nswing.c b/src/nswing/nswing.c index ff5201ad5ac..905a31e9f81 100644 --- a/src/nswing/nswing.c +++ b/src/nswing/nswing.c @@ -868,9 +868,8 @@ GMT_LOCAL int parse(struct GMT_CTRL *GMT, struct NSWING_CTRL *Ctrl, struct nestC case '8': case 'o': nesteds[7] = opt->arg; break; - default: - GMT_Report(GMT->parent, GMT_MSG_ERROR, "NSWING: Unknown option -%c%s\n", opt->option, opt->arg); - error = true; + default: /* Let GMT catch common options (like -R) already consumed by GMT_Parse_Common */ + error += gmt_default_option_error(GMT, opt); break; } } From 0fd77f567f291dc7c667d058735ed276ab023407 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Wed, 8 Jul 2026 17:49:02 +0100 Subject: [PATCH 29/29] Fix the access to layers in a nc cube file that has more that one variable. Before, this cube file wiyj two variables: 'z' and 'bathymetry' always returned the 'bathymetry' layer. grdinterpolate tsu.nc?z -T570 -Glay570.grd Now it works, (note the 'tsu.nc' file was created with the new nswing supplement). Assisted-by: Sonnet 5.0 --- src/grdinterpolate.c | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/grdinterpolate.c b/src/grdinterpolate.c index 03d23dc54be..de8d9786c6c 100644 --- a/src/grdinterpolate.c +++ b/src/grdinterpolate.c @@ -772,8 +772,12 @@ EXTERN_MSC int GMT_grdinterpolate (void *V_API, int mode, void *args) { } else if (C[GMT_IN] == NULL) { /* Read the cube */ gmt_M_free (GMT, level); /* Free this one now since it will be re-read by GMT_Read_Data */ - if ((C[GMT_IN] = GMT_Read_Data (API, GMT_IS_CUBE, GMT_IS_FILE, GMT_IS_VOLUME, GMT_CONTAINER_AND_DATA, wesn, Ctrl->In.file[0], NULL)) == NULL) - Return (GMT_DATA_READ_ERROR); + if (cube_layer[0]) /* Restore the requested variable name that was chopped off earlier */ + sprintf(file, "%s?%s", Ctrl->In.file[0], cube_layer); + else + strncpy(file, Ctrl->In.file[0], PATH_MAX-1); + if ((C[GMT_IN] = GMT_Read_Data(API, GMT_IS_CUBE, GMT_IS_FILE, GMT_IS_VOLUME, GMT_CONTAINER_AND_DATA, wesn, file, NULL)) == NULL) + Return(GMT_DATA_READ_ERROR); } if (convert_to_cube) { /* Just want to build cube from input stack */