diff --git a/doc/rst/source/modules-classic.rst b/doc/rst/source/modules-classic.rst index c79a2e45269..89903753327 100644 --- a/doc/rst/source/modules-classic.rst +++ b/doc/rst/source/modules-classic.rst @@ -145,6 +145,7 @@ All modules are requested via a call to the :doc:`gmt` program. supplements/mgd77/mgd77path supplements/mgd77/mgd77sniffer supplements/mgd77/mgd77track-classic + supplements/nswing/nswing supplements/potential/gmtflexure supplements/potential/gmtgravmag3d supplements/potential/gravfft @@ -315,6 +316,7 @@ Supplemental Modules - :doc:`/supplements/mgd77/mgd77path` - :doc:`/supplements/mgd77/mgd77sniffer` - :doc:`/supplements/mgd77/mgd77track-classic` + - :doc:`nswing ` - :doc:`gmtflexure ` - :doc:`gmtgravmag3d ` - :doc:`/supplements/potential/gravfft` @@ -675,6 +677,13 @@ MGD77 | :doc:`/supplements/mgd77/mgd77track-classic` | |mgd77track_purpose| | +-----------------------------------------------+------------------------+ +nswing +------ + ++-----------------------------------------------+------------------------+ +| :doc:`nswing ` | |nswing_purpose| | ++-----------------------------------------------+------------------------+ + potential --------- diff --git a/doc/rst/source/modules.rst b/doc/rst/source/modules.rst index 69123d14b37..224001eac40 100644 --- a/doc/rst/source/modules.rst +++ b/doc/rst/source/modules.rst @@ -143,6 +143,7 @@ All modules are requested via a call to the :doc:`gmt` program. supplements/mgd77/mgd77path supplements/mgd77/mgd77sniffer supplements/mgd77/mgd77track + supplements/nswing/nswing supplements/potential/gmtflexure supplements/potential/gmtgravmag3d supplements/potential/gravfft @@ -323,6 +324,7 @@ Supplemental Modules - :doc:`/supplements/mgd77/mgd77path` - :doc:`/supplements/mgd77/mgd77sniffer` - :doc:`/supplements/mgd77/mgd77track` + - :doc:`/supplements/nswing/nswing` - :doc:`/supplements/potential/gmtflexure` - :doc:`/supplements/potential/gravfft` - :doc:`/supplements/potential/gmtgravmag3d` @@ -703,6 +705,13 @@ MGD77 | :doc:`/supplements/mgd77/mgd77track` | |mgd77track_purpose| | +------------------------------------------+------------------------+ +nswing +------ + ++------------------------------------------+------------------------+ +| :doc:`/supplements/nswing/nswing` | |nswing_purpose| | ++------------------------------------------+------------------------+ + potential --------- diff --git a/doc/rst/source/supplements/module_supplements_purpose.rst_ b/doc/rst/source/supplements/module_supplements_purpose.rst_ index bf5333ea0bd..a88a1fabf8f 100644 --- a/doc/rst/source/supplements/module_supplements_purpose.rst_ +++ b/doc/rst/source/supplements/module_supplements_purpose.rst_ @@ -106,6 +106,8 @@ .. |rotsmoother_purpose| replace:: Get mean rotations and covariance matrices from set of finite rotations +.. |nswing_purpose| replace:: A tsunami maker + .. |grdbarb_purpose| replace:: Plot wind barb field from two component grids .. |barb_purpose| replace:: Plot wind barbs in 2-D and 3-D diff --git a/doc/rst/source/supplements/nswing/nswing.rst b/doc/rst/source/supplements/nswing/nswing.rst new file mode 100644 index 00000000000..3f873b27a43 --- /dev/null +++ b/doc/rst/source/supplements/nswing/nswing.rst @@ -0,0 +1,308 @@ +.. index:: ! nswing +.. include:: ../module_supplements_purpose.rst_ + +****** +nswing +****** + +.. only:: not man + + |nswing_purpose| + +Synopsis +-------- + +.. include:: ../../common_SYN_OPTs.rst_ + +**nswing** +*bathy.grd* *initial.grd* +**-t**\ *dt* +[ **-1**\ *bat_lev1* ] [ **-2**\ *bat_lev2* ] [ **-3**\ *...* ] +[ |-A|\ *fname.sww* ] +[ |-C| ] +[ |-D| ] +[ |-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* ] +[ |-H|\ [*momentM,momentN*\ [,\ *t*]] ] +[ |-P|\ *time_jump*\ [**+t**\ *run_time_jump*] ] +[ |-L|\ [*name1,name2*] ] +[ |-M|\ [**-**\|\ **+**\ [*maskname*]] ] +[ |-N|\ *n_cycles* ] +[ |-O|\ *BCfile* ] +[ |-Q|\ *z_offset* ] +[ |SYN_OPT-R| ] +[ |-S|\ [**x**\|\ **y**\|\ **n**][**+m**][**+s**][**+a**] ] +[ |-T|\ *mareg*\|\ *x/y*\ [**+o**\ *outmaregs*][**+t**\ *int*] ] +[ |-X|\ *manning0*\ [,\ *...*] ] +[ |SYN_OPT-V| ] +[ |SYN_OPT-f| ] +[ **-x**\ *n* ] +[ |SYN_OPT--| ] + +|No-spaces| + +Description +----------- + +**nswing** is a Non-linear Shallow Water model that propagates a tsunami over a +bathymetry grid. Starting from a base level bathymetry and an initial-condition +(source) grid it integrates the shallow water equations in time and writes the +resulting wave field (water level, velocity, momentum, energy, ...) either as a +series of GMT grids, a 3D netCDF file, or as time series sampled at virtual +maregraph (tide-gauge) locations. Nested grids of increasing resolution may be +used to refine the solution near the coast, and a tsunami source can be generated +on the fly from Okada fault parameters. + +Required Arguments +------------------ + +*bathy.grd* + The base level bathymetry grid (positive up; negative = below sea level). + (See :ref:`Grid File Formats `). +*initial.grd* + The initial-condition (source) grid holding the sea surface displacement at + time zero. (See :ref:`Grid File Formats `). + +.. _-t: + +**-t**\ *dt* + Time step (in seconds) for the simulation. + +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 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: + +**-A**\ *fname.sww* + Save the result as a *.sww* ANUGA-format file. + +**-n**\ *base* + Basename for MOST triplet files (no extension). + +.. _-C: + +**-C** + Add the Coriolis effect. + +.. _-D: + +**-D** + Write grids with the total water depth. These grids will have wave height + on the ocean and water thickness on land. + +.. _-E: + +**-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. + + 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: + +**-F**\ *x_epic/y_epic/dip/strike/rake/slip/length/width/topDepth* + Okada fault parameters used to build the source. *x_epic*, *y_epic* are the + x and y coordinates of the beginning of the fault trace; *dip*, *strike* + (azimuth), *rake*, *slip* (m), *length*, *width* and *topDepth* (depth from + the 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**\ *grid*) and writes it straight to **-G**'s name: no simulation is + run, so **-t** and **-G**'s saving interval are not needed. + +**-F**\ **k**\ *west/east/south/north* + Build a prism source with these limits and a height of 1 meter. + + - **-Fkc**\ *x/y/nx/ny* - 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* by *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 *dy* is not given, *dy* = *dx*). + + Using **-Fk** sets the output maregraph file to netCDF, unless rows = cols = 1. + +.. _-G: + +**-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: + +**-H** + Write grids with the momentum (velocity times water depth). + + - **-H**\ *fname_momentM,fname_momentN*\ [,\ *t*] - hot start using these + moment grids. The optional *t* is the hot-start time (also needs the + surface displacement corresponding to the time of these grids). + +.. _-P: + +**-P**\ *time_jump*\ [**+t**\ *run_time_jump*] + Do not write grids or maregraphs for times before *time_jump* (seconds). + When doing nested grids, append **+**\ *time* to NOT start nested-grid + computations before this time has elapsed. Allowed forms: **-P**\ *t1*, + **-P+**\ *t2*, **-P**\ *t1*\ **+**\ *t2* or **-P**\ *t1* **-P+**\ *t2*. + +.. _-L: + +**-L** + Use the linear approximation in the moment conservation equations (faster + but less accurate). + + - **-L**\ *in_fname,out_fname* - do Lagrangian tracers, where *in_fname* is + the tracers initial-position file and *out_fname* the file to hold the + results. + +.. _-M: + +**-M**\ [**-**\|\ **+**\ [*maskname*]] + 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 **-G** the *long* and *short* beach + arrays are also saved in the *.nc* file. + +.. _-N: + +**-N**\ *n_cycles* + Number of cycles in the simulation [Default 1010]. Total simulation time is *n_cycles* times the time step *dt*. + +.. _-O: + +**-O**\ *BCfile* + Name of a Boundary Condition ASCII file (experimental). + +.. _-Q: + +**-Q**\ *z_offset* + Apply a vertical offset to ALL bathymetry grids (e.g. to simulate tide). + +.. _-R: + +.. |Add_-Rgeo| replace:: |Add_-R_auto_table| +.. include:: ../../explain_-Rgeo.rst_ + +Output grids only in the sub-region enclosed by *west/east/south/north*. + +.. _-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 + 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**. + + 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*] + Save time series (maregraphs) at virtual tide-gauge locations. *mareg* 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** alone (without **-G**) is + allowed and runs a simulation that only outputs the maregraph series. + +.. _-X: + +**-X**\ *manning0*\ [,\ *manning1*\ [,\ *...*]][**+**\ *depth*] + Manning friction coefficients. If only one is provided, use it for all + nesting levels; otherwise specify one per level, comma separated. Append + **+**\ *depth* to apply Manning only at depths shallower than *depth* + (positive up). + +.. |Add_-V| replace:: |Add_-V_links| +.. include:: ../../explain_-V.rst_ + :start-after: **Syntax** + :end-before: **Description** + +.. |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_ + +Examples +-------- + +To propagate a tsunami over the bathymetry *bathy.grd* given the source +*source.grd*, using a 5 second time step and saving the water level every 10 +time steps to grids named *wave#.grd*, try:: + + gmt nswing bathy.grd source.grd -t5 -Gwave,10 + +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 -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 +-------- + +:doc:`gmt `, +:doc:`grdinterpolate ` + +Reference +--------- + +[Validation of NSWING, a multi-core finite difference code for tsunami propagation and run-up]( https://www.researchgate.net/publication/275349940_Validation_of_NSWING_a_multi-core_finite_difference_code_for_tsunami_propagation_and_run-up) 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/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/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..8897fb078eb --- /dev/null +++ b/src/nswing/README.nswing @@ -0,0 +1,19 @@ + +NSWING Supplemental Package +Distributed under the MIT License + +------------------------------------------------------ +Authors: Joaquim Luis, Miguel Miranda +Date: July 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. + diff --git a/src/nswing/nswing.c b/src/nswing/nswing.c new file mode 100644 index 00000000000..905a31e9f81 --- /dev/null +++ b/src/nswing/nswing.c @@ -0,0 +1,5637 @@ +/*-------------------------------------------------------------------- + * 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 MIT License + * + * Contact info: www.generic-mapping-tools.org + *--------------------------------------------------------------------*/ + +/* + * 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 + } + * + * 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" +/* 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 " +# 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(api, status) if (status) {GMT_Report(api, GMT_MSG_ERROR, "NSWING: error at line %d: %s\n", __LINE__, nc_strerror(status));} + + +#ifdef _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.n_columns) + +struct tracers { /* For tracers (oranges) */ + double *x; /* x coordinate */ + double *y; /* y coordinate */ +}; + +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 out_energy; /* To know if -E (or -Ep) was selected: write Energy or Power in the 3D netCDF grid */ + int out_power; + int append_z; /* -E|-S "+a": true = write z (sea surface) together with the extra quantity; false = the extra + quantity replaces z (-E) or z is omitted altogether (-S) */ + 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 GMT_GRID_HEADER hdr[10]; +}; + +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[], struct GMT_DATASET *D_in); +int read_tracers(void *API, struct GMT_GRID_HEADER hdr, char *file, struct tracers *oranges, struct GMT_DATASET *D_in); +int count_n_maregs(void *API, char *file, struct GMT_DATASET **D_out); +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 GMT_GRID_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(void *API, struct nestContainer *nest, int isGeog, int lev); +int intp_lin (void *API, 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(void *API, struct nestContainer *nest, double *flux_L1, double *flux_L2, char *what, int lev, int i_time); +void sanitize_nestContainer(struct nestContainer *nest); +void nestify(void *API, struct nestContainer *nest, int nNg, int recursionLevel, int isGeog); +void resamplegrid(struct nestContainer *nest, int nNg); +void edge_communication(void *API, 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_sp_M(struct nestContainer *nest, int lev); +void moment_sp_N(struct nestContainer *nest, int lev); +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); +int read_bnc_file(void *API, struct nestContainer *nest, char *file); +int interp_bnc (void *API, 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 GMT_GRID_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 GMT_GRID_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 GMT_GRID_HEADER hdr, double xx, double yy, double wx[], double wy[]); +double GMT_get_bcr_z(double *grd, double *bat, struct GMT_GRID_HEADER hdr, double xx, double yy); +void update_max(struct nestContainer *nest); +void update_max_velocity(struct nestContainer *nest); + + +void write_most_slice(void *API, 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(void *API, 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(void *API, 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(void *API, 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(void *API, 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(void *API, 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_(void *API, int status); + +int GetLocalNThread(void); + +/* ------------------------------------------------------------------------------ */ +/* 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 + internal header copy 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 GMT_GRID_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 = *G->header; /* Struct copy; we only ever read the plain public fields */ + h->hidden = NULL; /* The hidden part stays with (and dies with) G */ + 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]; + } + } +} + +/* Write an nswing array (south-up, scanline order) to disk as a GMT grid through + the GMT API (replaces the old Surfer-binary writer). Writes the sub-region + [i_start:i_end) x [j_start:j_end) of a grid that has nX columns per row. */ +GMT_LOCAL int gmtnswing_write_grid(void *API, 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) { + uint64_t row, ij; + unsigned int nx = i_end - i_start, ny = j_end - j_start, i, j, col; + double wesn[4], inc[2]; + struct GMT_GRID *G = NULL; + + wesn[XLO] = x_min; wesn[XHI] = x_min + (nx - 1) * x_inc; + wesn[YLO] = y_min; wesn[YHI] = y_min + (ny - 1) * y_inc; + inc[GMT_X] = x_inc; inc[GMT_Y] = y_inc; + + if ((G = GMT_Create_Data(API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, wesn, inc, + GMT_GRID_NODE_REG, GMT_NOTSET, NULL)) == NULL) + return (-1); + + for (j = j_start; j < j_end; j++) { /* nswing row j == y_min + j*y_inc (south-up) */ + row = ny - 1 - (j - j_start); /* GMT row 0 == north */ + for (i = i_start; i < i_end; i++) { + col = i - i_start; + ij = gmt_M_ijp(G->header, row, col); + G->data[ij] = work[ijs(i,j,nX)]; + } + } + + if (GMT_Write_Data(API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, name, G) != GMT_NOERROR) { + GMT_Destroy_Data(API, &G); + return (-1); + } + GMT_Destroy_Data(API, &G); + return (0); +} + +/* --------------------------------------------------------------------------- */ +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[+m], " + "[-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][+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); + + 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-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][+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."); + 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[+m],int"); + GMT_Usage(API, -2, "Output name of a 3D netCDF file is .nc. Save steps at the intervals. Optionally " + "append '+m' to save each time step as a separate file. Files are then named #.grd."); + 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-P"); + GMT_Usage(API, -2, "Do not write grids or maregraphs for times before (seconds)."); + //GMT_Usage(API, -2, "When doing nested grids, append +