USERS’ MANUAL

Program outline and flow chart

The code was written using Fortran 90 with the c preprocessor (cpp) statements for separation of the source code. Arrays are dynamically allocated at runtime. Precision is selected using the selected_real_kind Fortran intrinsic function defined in the makefile. The default precision is single.

The present version of NearCoM-TVD includes a number of options including (1) choice of serial or parallel code (2) Cartesian or curvilinear coordinate, (3) samples.

The flow chart is shown in Fig. 4.

_images/chart.jpg

Fig. 4 Flow chart of the main program.

Permanent variables associated with coupling

Depth(): still water depth h at element point

DepthX(): still water depth h at x-interface

DepthY(): still water depth h at y-interface

Eta(): surface elevation, for dry point, Eta() = MinDepth - Depth(), MinDepth is specified in input.txt.

Eta0(): \(\eta\) at previous time level

MASK(): 1 - wet, 0 - dry

MASK_STRUC(): 0 - permanent dry point

U(): depth-averaged u

V(): depth-averaged v

Uc(): contravariant component of depth-averaged velocity in \(\xi_1\) direction

Vc(): contravariant component of depth-averaged velocity in \(\xi_2\) direction

HU(): \((h+\eta)u\) at element

HV(): \((h+\eta)v\) at element

P(): \((h+\eta)U\) at x- interface for Cartesian, and \((h+\eta)Uc\) at \(\xi_1\) -interface curvilinear

Q(): \((h+\eta)V\) at y- interface for Cartesian, and \((h+\eta)Vc\) at \(\xi_2\) -interface curvilinear

Fx(): numerical flux F at x-interface

Fy(): numerical flux F at y-interface

Gx(): numerical flux G at x-interface

Gy(): numerical flux G at y-interface

Ubar(): \(HU\) for Cartesian and \(JHU\) for curvilinear

Vbar(): \(HV\) for Cartesian and \(JHV\) for curvilinear

EtaRxL(): \(\eta\) Left value at x-interface

EtaRxR(): \(\eta\) Right value at x-interface

EtaRyL(): \(\eta\) Left value at y-interface

EtaRyR(): \(\eta\) Right value at y-interface

HxL(): total depth Left value at x-interface

HxR(): total depth Right value at x-interface

HyL(): total depth Left value at y-interface

HyR(): total depth Right value at y-interface

HUxL(): \((h+\eta)u\) Left value at x-interface

HUxR(): \((h+\eta)u\) Right value at x-interface

HVyL(): \((h+\eta)v\) Left value at y-interface

HVyR(): \((h+\eta)v\) Right value at y-interface

PL(): Left P value at x-interface

PR(): Right P value at x-interface

QL(): Left Q value at y-interface

QR(): Right Q value at y-interface

Installation and compilation

NearCoM-TVD is distributed in a compressed fie. To install the programs, first, uncompress the package. Then use

\(>\) tar xvf \(*\).tar

to extract files from the uncompressed package. The exacted files will be distributed in two new directories: /CIRC_SWAN and /work.

To compile the program, go to /CIRC_SWAN and modify Makefile if needed. There are several necessary flags in Makefile needed to specify below.

-DDOUBLE_PRECISION: use double precision, default is single precision.

-DPARALLEL: use parallel mode, default is serial mode.

-DSAMPLES: include all samples, default is no sample included.

-DCURVILINEAR: curvilinear version, otherwise Cartesian.

\(\bf NOTE:\) setting curvilinear is a must for SWAN and SHORECIRC coupled model.

-DSEDIMENT: include sediment and seabed modules.

-DINTEL: INTEL compiler.

-DRESIDUAL: include tidal residual calculation.

-DSTATIONARY: stationary mode for SHORECIRC

CPP: path to CPP directory.

FC: Fortran compiler.

Then execute

\(>\) make clean

\(>\) make

The executable file ‘nearcom’ will be generated and copied from /CIRC_SWAN to /work/. Note: use ‘make clean’ after any modification of Makefile.

To run the model, go to /work. Modify INPUT if needed and run.

Input

Following are descriptions of parameters in input.txt (\(\bf NOTE:\) all parameter names are capital sensitive).

  • SWAN INPUT: refer to SWAN manual. Model run time is set in SWAN model. For example,

    COMPUTE NONSTAT 20081114.160000 1 MI 20081114.230000

    The above setting means model run start from 2008 11 14 16:00 to 2008 11 14 23:00. The model call swan at \(DT_{\mbox{swan}}\) = 1 minute. The loop number for SHORECIRC and SEDIMENT is estimated by \(DT_{\mbox{swan}}\) and the time step of SHORECIRC (time varying).

    IMPORTANT SETTING IN SWAN:

    1. in SET, always set CARTESIAN in order to make a grid orientation consistent with SHORECIRC

    2. in SET, always set [inrhog] as 1 to get a true wave energy dissipation.

    3. in COMPUTE, always set NONSTAT mode.

  • WAVE CURRENT INTERACTION

SWAN_RUN: logical parameter to run SWAN

SHORECIRC_RUN: logical parameter to run SHORECIRC

WC_BOUND_WEST: west bound region (number of grid point) in which wave-current is inactive.

WC_BOUND_EAST : east bound region (number of grid point) in which wave-current is inactive.

WC_BOUND_SOUTH : south bound region (number of grid point) in which wave-current is inactive.

WC_BOUND_NORTH: north bound region (number of grid point) in which wave-current is inactive.

WC_LAG : time delay for wave-current interaction

  • TITLE:

title for SHORECIRC log file

  • SPECIFICATION OF MULTI-PROCESSORS

PX: processor numbers in X

PY: processor numbers in Y

\(\bf NOTE:\) PX and PY must be consistency with number of processors defined in mpirun command, e.g., mpirun -np n (where n = px \(\times\) py).

  • SPECIFICATION OF WATER DEPTH

DEPTH_TYPE: depth input type.

DEPTH_TYPE=DATA: from a depth file.
The program includes several simple bathymetry configurations such as
DEPTH_TYPE=FLAT: flat bottom, need DEPTH_FLAT
DEPTH_TYPE=SLOPE: plane beach along \(x\) direction. It needs three parameters: slope,SLP, slope starting point, Xslp and flat part of depth, DEPTH_FLAT
DEPTH_FILE: bathymetry file if DEPTH_TYPE=DATA, file dimension should be Mglob x Nglob

with the first point as the south-west corner. The read format in the code is shown below.

DO J=1,Nglob
 READ(1,*)(Depth(I,J),I=1,Mglob)
ENDDO
DEPTH_FLAT: water depth of flat bottom if DEPTH_TYPE=FLAT or DEPTH_TYPE=SLOPE

(flat part of a plane beach).

SLP: slope if DEPTH_TYPE=SLOPE

Xslp: starting \(x\) (m) of a slope, if DEPTH_TYPE=SLOPE

  • SPECIFICATION OF RESULT FOLDER

RESULT_FOLDER: result folder name, e.g., RESULT_FOLDER = /Users/fengyanshi/tmp/

  • SPECIFICATION OF DIMENSION

Mglob: global dimension in \(x\) direction.

Nglob: global dimension in \(y\) direction.

\(\bf NOTE:\) For parallel runs, Mglob and Nglob can be divided by PX and PY, respectively. MAX(Mglob,Nglob) can be divided by PX \(\times\) PY.

  • SPECIFICATION OF STATIONARY MODE

N_ITERATION: the iteration number for stationary mode of SHORECIRC

(set -DSTATIONARY in Makefile).

WATER_LEVEL_FILE: the file name of water level file containing time and water level, for

stationary mode. The following example shows the format.

water levels for stationary mode
5 - number of water level data
0.0 0.0 ! Time (s), Level (m)
3600.0 0.5
7200.0 0.8660
10800.0 1.0
14400.0 0.866
18000.0 0.5
  • SPECIFICATION OF TIME

PLOT_INTV: output interval in seconds (Note, output time is not exact because adaptive dt is used.)

SCREEN_INTV: time interval (s) of screen print.

PLOT_INTV_STATION: time interval (s) of gauge output

  • SPECIFICATION OF GRID

DX: grid size(m) in \(x\) direction, for Cartesian mode

DY: grid size(m) in \(y\) direction, for Cartesian mode

X_FILE: name of file to store x for curvilinear mode

Y_FILE: name of file to store y for curvilinear mode

\(\bf NOTE:\) data format is the same as the depth data shown above.

CORI_CONSTANT: logical parameter for constant Coriolis parameter

LATITUDE: latitude if constant Coriolis parameter is used

LATITUDE_FILE: name of file to store latitude at every grid point if not constant Coriolis

\(\bf NOTE:\) data format is the same as the depth data shown above.

  • BOUNDARY CONDITIONS

ETA_CLAMPED: logical parameter for surface elevation clamped condition

V_CLAMPED: logical parameter for velocity clamped condition

FLUX_CLAMPED: logical parameter for flux clamped condition

TIDE_FILE: name of file to store tidal constituents

DATA FORMAT: please refer to mk_tide.f90. The formula of surface elevation at a tidal boundary can be expressed by

(67)\[\eta_0 (t) = \sum_{n=1}^Na_{0}({\bf x}, n) f_c (n) \cos \left(\frac{2\pi}{T(n)} t - \phi({\bf x}, n) + (V_0 +u_0)(n) \right)\]

where \(a_0\) and \(\phi\) represent amplitude and phase lag, respectively, for a harmonic constituent at location \(\bf x\). \(T\) is tidal period. \(f_c\) and \((V_0+u_0)\) are the lunar node factor and the equilibrium argument, respectively, for a constituent.

The following is an example of M2 + O1.
tidal boundary conditions
150 — number of days from Jan 1, to simulation date
2 — number of constituents
1.000 0.000 — \(f_c\) and \((V_0+u_0)\) for M2
0.980 0.000 — \(f_c\) and \((V_0+u_0)\) for O1
80 — number of tidal boundary points
1 , 1 — (i,j) grid location of tidal boundary
12.420 1.200 21.000 — \(T\), amplitude \(a_0\) and phase lag \(\phi\) for M2
24.000 0.3 30.100 – \(T\), amplitude \(a_0\) and phase lag \(\phi\) for O1
2 , 1 — (i,j) grid location of tidal boundary
12.420 1.200 21.000 — \(T\), amplitude \(a_0\) and phase lag \(\phi\) for M2
24.000 0.3 30.100 – \(T\), amplitude \(a_0\) and phase lag \(\phi\) for O1
3 , 1
…

FLUX_FILE: name of file to store time series of flux (e.g., unit width river flux)

DATA FORMAT:
title
Number of data, Number of flux point
I, J, River orientation
Time, Flux, Angle in Cartesian
…
where (I,J) represent grid points of river location. River orientation represents the direction which a river flows from in the IMAGE domain (for curvilinear coordinates). Use W,E,S and N for the orientation. For example, ‘W’ represents a river flowing into the domain from the west boundary (in IMAGE domain for curvilinear coordinates).

Please refer to mk_river.f90. The following is an example.
river flux boundary condition
5 2 ! NumTimeData, NumFluxPoint
1 38 W ! I, J, River_Orientation
0.000 0.200 0.000
360000.000 0.200 0.000
720000.000 0.200 0.000
1080000.000 0.200 0.000
1440000.000 0.200 0.000
1 39 W ! I, J, River_Orientation
0.000 0.200 0.000
360000.000 0.200 0.000
720000.000 0.200 0.000
1080000.000 0.200 0.000
1440000.000 0.200 0.000
end of file
  • WIND CONDITION

    Spatially uniform wind field is assumed in this version.

WindForce: logical parameter for wind condition, T or F.

WIND_FILE: name of file for a time series of wind speed.

DATA FORMAT: the following is an example of wind data.

wind data
100 - number of data
0.0 , -10.0 0.0 — time(s), wu, wv (m/s)
2000.0, -10.0, 0.0
8000.0, -10.0, 0.0
…

Cdw: wind stress coefficient for the quadratic formula.

  • SPECIFICATION OF INITIAL CONDITION

INT_UVZ : logical parameter for initial condition, default is FALSE

ETA_FILE: name of file for initial \(\eta\), e.g., ETA_FILE= /Users/fengyanshi/work/input/CVV_H.grd, data format is the same as depth data.

U_FILE: name of file for initial \(u\), e.g.,U_FILE= /Users/fengyanshi/work/input/CVV_U.grd, data format is the same as depth data.

V_FILE: name of file for initial \(v\), e.g., V_FILE= /Users/fengyanshi/work/input/CVV_V.grd, data format is the same as depth data.

  • SPECIFICATION OF WAVEMAKER

There is no wavemaker implemented in SHORECIRC.

  • SPECIFICATION OF PERIODIC BOUNDARY CONDITION

    (Note: only south-north periodic condition was implemented)

PERIODIC_X: logical parameter for periodic boundary condition in x direction, T - periodic, F - wall boundary condition.

PERIODIC_Y: logical parameter for periodic boundary condition in x direction.

Num_Transit: grid numbers needed to make periodic condition for SWAN. The reason to set this parameter is that SWAN doesn’t have an option for periodic boundary condition. In this implementation, a periodic boundary condition is implemented by making a transition from a left array ( count to Num_Transit from left boundary) to a right array.

  • SPECIFICATION OF SPONGE LAYER

SPONGE_ON: logical parameter, T - sponge layer, F - no sponge layer.

Sponge_west_width: width (m) of sponge layer at west boundary.

Sponge_east_width: width (m) of sponge layer at east boundary.

Sponge_south_width: width (m) of sponge layer at south boundary.

Sponge_north_width width (m) of sponge layer at north boundary

R_sponge: decay rate in sponge layer. Its values are between 0.85 \(\sim\) 0.95.

A_sponge: maximum damping magnitude. The value is \(\sim\) 5.0.

  • SPECIFICATION OF OBSTACLES

OBSTACLE_FILE: name of obstacle file. 1 - water point, 0 - permanent dry point. Data dimension is (Mglob \(\times\) Nglob). Data format is the same as the depth data.

  • SPECIFICATION OF PHYSICS

Cd: quadratic bottom friction coefficient

nu_bkgd : background eddy viscosity parameter.

  • SPECIFICATION OF NUMERICS

Time_Scheme: stepping option, Runge_Kutta or Predictor_Corrector (not suggested for this version).

HIGH_ORDER: spatial scheme option, FOURTH for the fourth-order, THIRD for the third-order, and SECOND for the second-order (not suggested for Boussinesq modeling).

CONSTRUCTION: construction method, HLL for HLL scheme, otherwise for averaging scheme.

CFL: CFL number, CFL \(\sim\) 0.5.

FroudeCap: cap for Froude number in velocity calculation for efficiency. The value could be 5 \(\sim\) 10.0.

MinDepth: minimum water depth (m) for wetting and drying scheme. Suggestion: MinDepth = 0.001 for lab scale and 0.01 for field scale.

MinDepthFrc: minimum water depth (m) to limit bottom friction value. Suggestion: MinDepthFrc = 0.01 for lab scale and 0.1 for field scale.

  • SPECIFICATION OF TIDAL RESIDUAL

T_INTV_mean: time-averaging interval for Eulerian mean current and elevation. Note: use -DRESIDUAL in Makefile to make this option active.

  • SPECIFICATION OF SEDIMENT CALCULATION

    Note: set -DSEDIMENT in Makefile to make sediment module active

T_INTV_sed: time interval to call sediment module

Factor_Morpho: morphology factor.

D_50 : \(D_{50}\)

D_90 : \(D_{90}\)

por: sediment porosity

RHO: water density

nu_water: water eddy viscosity

S_sed: specific gravity

SOULSBY: logical parameter for Soulsby (1997) total load formula, T = true, F = false

z0: \(z_0\), bed roughness length.

KOBAYASHI: logical parameter for KOBAYASHI’s formula, T = true, F = false

angle_x_beach: coordinate rotation angle defined in Figure 1.

eB: \(e_B\), suspension efficiency for energy dissipation rate due to wave breaking

ef: \(e_f\), suspension efficiency for energy dissipation rate due to bottom friction

a_k: \(a\), empirical suspended load parameter.

b_k: \(b\), empirical bedload parameter.

TanPhi: \(\tan \phi\), where \(\phi\) is the angle of internal friction of the sediment.

Gm: \(G_m\) for slope function (\(G_m=10\)).

frc: friction coefficient in Kobayashi.

Si_c: a coefficient in calculating \(P_b\).

  • SPECIFICATION OF OUTPUT VARIABLES

NumberStations: number of station for output. If NumberStations \(> 0\), need input i,j in STATION_FILE

DEPTH_OUT: logical parameter for output depth. T or F.

U: logical parameter for output \(u\). T or F.

V: logical parameter for output \(v\). T or F.

ETA: logical parameter for output \(\eta\). T or F.

HS: logical parameter for output of significant wave height \(H_s\). T or F.

WFC: logical parameter for output of wave force. T or F.

WDIR: logical parameter for output of peak wave direction. T or F.

WBV: logical parameter for output of wave orbital velocity. T or F.

MASK: logical parameter for output wetting-drying MASK. T or F.

SourceX: logical parameter for output source terms in \(x\) direction. T or F.

SourceY: logical parameter for output source terms in \(y\) direction. T or F.

UV3D: logical parameter for output 3D structure. T or F.

Qstk: logical parameter for output Stokes mass flux. T or F.

DepDt: logical parameter for output depth variation rate. T or F.

Qsed: logical parameter for output sediment transport rate. T or F.

Output

The output files are saved in the result directory defined by RESULT_FOLDER in INPUT. For outputs in ASCII, a file name is a combination of variable name and an output series number such eta_0001, eta_0002, …. The format and read/write algorithm are consistent with a depth file. Output for stations is a series of numbered files such as sta_0001, sta_0002 ….