TsunAWI output files¶

Netcdf output¶

As the most common format in oceanography, netcdf (https://www.unidata.ucar.edu/software/netcdf/) is also used for TsunAWI. The output file contains meta data on the model run, input data, and of course the results. The most important parameters and fields are

nodes_2D

Number of mesh nodes.

elements_2D

Number of mesh elements.

iteration

Number of written time steps (determined by T_out).

time(iteration)

Model time of the output time steps, in seconds since the start of the simulation.

surface_locations(nodes_2D,2)

Node coordinates, either as geographical degrees or in meters.

triangles(elements_2D,3)

Triangulation as in elem_2D.out.

tri_area(elements_2D)

Areas of the triangles [m:math:^2].

initial_topography(nodes_2D)

Bathymetry [m] at the nodes, before smoothing and before adding the uplift, subtracting the subsidence due to the earthquake.

topography(nodes_2D)

Bathymetry [m] at the nodes, after smoothing and adding/substracting the uplift/subsidence.

ssh(iteration,nodes_2D)

If enable_snapshots=.true.: Sea surface elevation $\zeta_h$ at the nodes for each output time step T_out [s].

This forms usually the largest part of the netcdf output. To improve compression of the resulting file, $\zeta_h$ can be written with reduced accuracy (skipping too detailed information) specified as nc_snapshot_accuracy_ssh [m] in the namelist. Reasonable values are 0.01 (1cm) or 0.001 (1mm).

first_arrival(nodes_2D)

Arrival time [s], defined as the first time to reach the threshold of 0.01m.

max_wave_height(nodes_2D)

For each node, determine the maximum ssh [m] over all timesteps.

max_abs_vel(nodes_2D)

Determine the maximum absolute velocity in each edge, interpolated to the nodes. To be taken with care in the inundation zone, when the velocity is extrapolated.

max_flux(nodes_2D)

In each node, multply the absolute velocity and the water depth and determine the maximum value over all timesteps.

Further fields can be written on request:

velocity_in_nodes(iteration,nodes_2D,2): If write_velocity=.true., the horizontal velocity $\vec{v}_h$ is interpolated from the edges to the nodes and written.
first_arrival_ttt(nodes_2D): The estimated travel time of the wave purely determined by bathymetry with velocity $\sqrt{g h}$ ) is computed and written to the output before the main computation starts. Compared to first_arrival derived from the simulation, this simple estimate usually results in faster arrival and in a smooth field by definition. This makes it handy to use as a warning product: The threat is overestimated, and the isolines look nice.
epot(iteration), ekin(iteration): In each timestep, intergrate the potential and kinetic energy [J].
cfl(iteration, nodes_2D): If write_cfl=.true., the Courant–Friedrichs–Lewy condition can be examined: $|\vec{v}_h|/h$ is computed for each edge, $h$ =edgelength, and interpolated to the nodes.

Netcdf output of inundated area¶

Switch on with write_nc_inundation=.true.. For faster postprocessing, the information can be reduced to the inundated area, i.e., the triangulation is reduced to elements that were fully on land before the earthquake elevated or subsided the topography, and at least one node of each triangle was inundated with.

The following parameters and fields are available:

nodes_2D, elements_2D: The number of nodes and elements of the reduced triangulation
longitude(nodes_2D), latitude(nodes_2D): Coordinates of these nodes
triangles(elements_2D, 3), tri_area(elements_2D): The reduced triangulation and the area of each element, which might be handy in post processing
initial_topography(nodes_2D), topography(nodes_2D): The topography before and after the rupture(s) and the smoothing step
max_wave_height(nodes_2D): As we only regard nodes on land: the maximum depth of the inundation at each node
max_abs_vel(nodes_2D), max_flux(nodes_2D): Maximum velocity and flux at each node.

Mareogram at tide gauge location¶

Mareograms at tide gauge locations can be written into ascii files. In the namelist, set write_tidegauge_data=.true. and the interval in T_out_tidegauge, if it differs from the snapshot interval T_out. The tide gauge locations are read from the file TG_list.txt in the working directory or any other file specified by TideGaugeFile in the namelist.

Each row defines one tide gauge with the location’s longitude and latitude. The tide gauge name may be the urn from the GITEWS project, a short identifier, or some own definition without any space and possibly a short identifier separted by a “:”.

urn:gfz:gitews:def:procedure:tg:station:como  43.2481  -11.7035
urn:gfz:gitews:def:procedure:tg:station:meul  96.2167    4.3167
musc                                          23.6333   58.5667
Northsea-Germany-Cuxhaven:cux                 53.8678    8.7175

The tide gauge identifier or the last part of the urn or similar definition (extracted as the characters between the last “:” and the space), becomes part of the file name, e.g., for the first station in the list above TG_como_scen-id.out. In the following example output, the tide gauge was vertically moved by the virtual earthquake. The water level is simply extraced from the nearest grid node in water (no interpolation).

%  urn:gfz:gitews:def:procedure:tg:station:gano
% TG longitude, latitude        102.277800   -5.346100
% Index of nearest grid node with depth>1m      573645
% Grid node longitude latitude  102.277243   -5.345681
% Distance [m] to grid node     77.4244
% Water depth [m] at grid node      8.1311
% time since rupture [min]   sea surface elevation [m]
         0.00     -0.0005
         1.00     -0.0005
         2.00     -0.0005
         3.00     -0.0005
         4.00     -0.0005
         5.00     -0.0005
         6.00     -0.0005
         7.00     -0.0005
         8.00     -0.0004
         9.00     -0.0003
        10.00     -0.0002

If you have other needs, e.g., another input format, modify output.F90, subroutines ascii_out_tidegauge_init and ascii_out_tidegauge.

Raster output¶

This experimental TsunAWI feature interpolates to a raster of given resolution. We chose Surfer 6 Binary Grid File Format (see http://grapherhelp.goldensoftware.com/subsys/surfer_6_grid_file_format.htm), because it is a simple binary format that can be written w/o additional libraries. Even if compressed by gzip, if can be read by QGIS and translated to e.g. GeoTiff with GDAL.

Switch on with enable_raster_out=.true., see &output_raster for further settings in the namelist.

Unfortunately, many choices are still hard coded and you have to edit output_raster.F90 e.g. to switch between values only in the ocean, only on land, or both.

Summing up, the raster output produces (replace ID by the ID chosen for the simulation):

ID_mwh.grd: Maximum ssh over all timesteps, linear interpolation.
ID_mwh_maxLand.grd: Maximum inundation, with a more careful maximum value approach for each raster pixel.
ID_ssh_tttttt.grd: Sea surface height, one file each raster_ssh_T_out seconds, here tttttt seconds after the tsunami origin, padded with zeros to 6 digits..
ID_ttt.grd: If enable_write_ttt and enable_raster_out=.true., interpolate the estimated tsunami trave time to raster, too.

We recommend to compress the grd-files with gzip, because usually, it saves a lot of disk space, in particular, if the accuracy was reduced.

The interpolation is a simple linear interpolation from the triangular element which includes the middle of the raster pixel. Clearly, this is only a valid choice for a qualitative visualisation. If the raster is coarser than the computational mesh, and many features of the topography and thus the tsunami are not resolved, it is basically a lottery which value is chosen for a pixel.

Therefore, a more careful approach is used for TSID_mwh_maxLand.grd. Here, the linear interpolation is replaced by the maximum value of the linear interpolation and of all values at nodes inside the pixel. The aim is to not underestimate the inundation by an unfortunate interplay of computational mesh, linear interpolation, and raster resolution.

ASCII output¶

The ASCII output was added for users without a working netcdf implementation. It is enabled by enable_ascii_out=.true. in the namelist, and provides the following output files for variables at the nodes during and at the end of the run

(path=OutputPath/scenarioID_prefix)

Sea surface elevation: Each $n\times$ T_out model seconds, the sea surface height $\zeta_h$ is written to path_ssh_00000n.out ( $n$ is padded with zeros to 6 digits).
Velocity: If write_velocity=.true., in the output step also the horizontal velocity $\vec{v}_h=(u_n,v_h)$ interpolated to the nodes are written to path_vel_00000n.out
Maximum amplitude, arrival time: At the end of the run, the maximum amplitude at each node is written to path_mwh.out and the arrival time (threshold of 0.01m) to path_arrivaltime.out.
Estimated travel time: If write_ttt=.true., this simple estimate (based on a phase speed of $\sqrt{g/h}$ ) is written to path_ttt.out.

Binary output¶

A routine for binary output is implemented (enable_bin_out=.true.), but has not been used or maintained for some while. Please refer to the source code output.F90, subroutine binout and adjust to your needs. As the ASCII output was more recently added, it might be an idea to start from the routines ascii_out_init, ascii_out, ascii_finalize and modify those for binary output.

Restart files¶

Time consuming computations can be divided into chunks, e.g., to meet the time limits of a batch system. If write_restart=.true, restart files are written after a model time of T_chunk seconds and TsunAWI stops. The following model run automatically detects restart files for the current namelist configuration (output_prefix, type of run) and proceeds the computation until 2 $\times$ T_chunk seconds are reached and so on.