SELFE

SELFE 1.4a User Manual

Input files:
Output files:

Input files

Important symbols

"!" is used to add comments after actual input (as customary in FORTRAN 90);

np: # of nodes in the horizontal grid;

ne: # of elements in the horizontal grid;

ns: # of sides in the horizontal grid;

nvrt: total # of levels in the vertical grid

go back

Horizontal grid (hgrid.gr3): same as ELCIRC

In xmgredit grid format.

Below is a sample:

fort_12062002.14 : alphanumeric description

55880 30001 : # of elements and nodes in the horizontal grid

(Coordinate part):

1 386380.409604 286208.187634 5.122 : node #, x,y, depth
2 386460.352736 285995.038877 9.167
3 386687.720000 286213.590000 1.000
4 386460.076848 286367.779818 2.209
5 386678.380000 286483.440000 1.614
6 386180.219063 286405.956765 4.627
7 386409.007263 286563.660632 2.629
8 386186.575437 286680.225393 4.195
9 385958.392423 286604.196847 4.177

.............................................

(Connectivity part)

1 3 1 2 3 : element #, element type (triangles only), node 1, node 2, node 3
2 3 3 4 5
3 3 3 5 6
4 3 1 4 6
5 3 4 7 9

...........................................

(Boundary condition part)

3 : Number of open boundaries
95 : Total number of open boundary nodes
85 : Number of nodes for open boundary 1
15185 : first node

........................................

25 = number of land boundaries
4079 = Total number of land boundary nodes
1452 0 = Number of nodes for land boundary 1
19947 : first node
19945
19943

.......................................

Note: the boundary condition (b.c.) part can be generated with xmgredit5 --> GridDEM --> Create open/land boundaries.

go back

Vertical grid (vgrid.in)

This version uses hybrid S-Z coordinates in the vertical, with S on top of Z.

54 18 100. !nvrt; kz (# of Z-levels); h_s (transition depth between S and Z)
Z levels !Z-levels first
1 -5000. !level index, z-coordinates
2 -2300.
3 -1800.
4 -1400.
5 -1000.
6 -770.
7 -570.
8 -470.
9 -390.
10 -340.
11 -290.
12 -240.
13 -190.
14 -140.
15 -120.
16 -110.
17 -105.
18 -100. !z-coordinate of the last Z-level must be -h_s
S levels !S-levels below
30. 0.7 10. ! constants used in S-coordinates: h_c, theta_b, theta_f (see notes below)
18 -1. !first S-level (S-coordinate must be -1)
19 -0.972222 !levels index, S-coordinate
20 -0.944444
.......
54 0. !last S-coordinate must be 0

Notes:

Origin of the vertical axis is at MSL;
h_c, theta_b, theta_f: constants used in Song and Haidvogel's (1994) S-coordinate system. h_c controls surface/bottom boundary layer thickness that requires fine resolution; theta_f-->0 leads to traditional sigma-coordinates, while theta_f >> 1 skews resolution towards surface and/or bottom; theta_b=1 leads to both bottom and surface being resolved, while theta_b=0 resolves only surface.
Global output is from the bottom (variable in space) to surface (level nvrt) at each node;
If a "pure S" model is desired, use only 1 Z-level and set h_s to a very large number (e.g., 1.e6).

go back

Parameter input (param.in):

Explanation of each line:

48-character string description of the version.
48-character start time info string, e.g., 04/23/2002 00:00:00 PST (only used for visualization with xmvis)
ipre: pre-processing flag. ipre=1: code will output centers.bp, sidecenters.bp, obe.out (centers build point, sidcenters build point, and list of open boundary elements), and mirror.out and stop. This is useful also for checking z-levels at various depths (in mirror.out) for any given choice of vgrid. ipre=0: normal run.
nscreen = screen output on/off switch (0: off; 1: on). In either case, mirror output messages will be directed into mirror.out.

iwrite: writing destination option for COIRE system. Default: 0.

imm: tsunami option. Default: 0 (no bed deformation); 1: with bed deformation (needs bdef.bp).

If imm=1, this line is ibdef: total # of deformation steps (i.e., the bed will change from the initial position to the position specified in bdef.bp in ibdef steps).

ihot = hot start flag. If ihot=0, cold start; if ihot/=0, hot start from hotstart.in. If ihot=1, the time and time step are reset to zero, and outputs start from T=0 accordingly. If ihot=2, the run (and output) will continue from the time specified in hotstart.in.

ics = coordinate frame flag. If ics=1, Cartesian coordinates are used; if ics=2, degrees latitude/longitude are used (but the output will still be in Cartesian coordinates).

slam0, sfea0 = centers of projection used to convert lat/long to Cartesian coordinates. These are used for ics=2, or a variable Coriolis parameter is employed (ncor=1), or the heat exchange sub-model is invoked (ihconsv=1).

theta0 = implicitness parameter (between 0.5 and 1). Recommended value: 0.6.

ibcc, ibtp = barotropic/baroclinic flags. If ibcc=0, a baroclinic model is used and regardless of the value for ibtp, the transport equation is solved. If ibcc=1, a barotropic model is used, and the transport equation may (when ibtp=1) or may not (when ibtp=0) be solved; in the former case, S and T are treated as passive tracers.

If ibcc=0, the next line is: nrampbc, drampbc: ramp option flag and ramp-up period (in days). If nrampbc=0, drampbc is not used. A hyperbolic tangent function is used for ramp-up.

rnday = total # of run days.

nramp, dramp = ramp option for the tides and some boundary conditions, and ramp-up period in days (not used if nramp=0).

dt = time step (in sec).

Inactive (will be removed eventually).

nadv = advection on/off switch. If nadv=0, advection is selectively turned off based on the input file adv.gr3. If nadv=1 or 2, advection is on for the whole domain, and backtracking is done using either Euler or 5th-order Runge-Kutta (expensive) scheme.

dtb1, dtb2: sub-steps used in btrack. For nadv /= 0, only dtb1 is used, which is the minimum sub-step used in Euler scheme (nadv=1) or maximum sub-step used in 5th-order Runge-Kutta scheme (nadv=2). For nadv=0, dtb1 is the minimum sub-step used in Euler scheme (when local depth=1 in adv.gr3) and dtb2 is the maximum sub-step used in 5th-order Runge-Kutta scheme (when local depth=2 in adv.gr3).

h0 = minimum depth (in m) for wetting and drying (recommended value: 1cm). When the total depth is less than h0, the corresponding nodes/sides/elements are considered dry. It should always be positive to prevent underflow.

nchi = bottom friction option. If nchi=0, spatially varying drag coefficients are read in from drag.gr3 (as depth info). For nchi=1, bottom roughnesses (in meters) are read in from rough.gr3. In all cases, the logarithmic law is assumed.

If nchi=1, the next line is: Cdmax = max. drag coefficient (to prevent exaggeration of drag coefficient in shallow areas).

ncor = Coriolis option. If ncor=0 or -1, a constant Coriolis parameter is used (see next line). If ncor=1, a variable Coriolis parameter, based on a beta-plane approximation, is used, with the lat/long. coordinates read in from hgrid.ll. In this case, the center of CPP projection must be correctly specified (see above).

If ncor=0, the next line is: cori = constant Coriolis parameter. If ncor=-1, the next line is the reference lattitude in degrees.

nws, wtiminc = wind forcing options and the interval (in seconds) with which the wind input is read in. If nws=0, no wind is applied (and wtiminc becomes immaterial). If nws=1, constant wind is applied to the whole domain at any given time, and the time history of wind is input from wind.th. If nws=2 or 3, spatially and temporally variable wind is applied and the input consists of a number of hdf files in the directory hdf/. The option nws=3 is only for checking heat conservation and needs sflux.th.

If nws>0, the next line is: nrampwind, drampwind = ramp option and period (in days) for wind.

ihconsv, isconsv = heat budget and salt conservation models flags (the latter is inactive at the moment). If ihconsv=0, the heat budget model is not used. If ihconsv=1, a heat budget model is invoked, and a number of hdf files for radiation flux input are read in from he directory hdf/. If ELM option is used for heat transport (see iupwind_t below), this is only used to invoke Zeng's bulk aerodynamic model, and the heat budget model is bypassed in the code.

itur = turbulence closure model selection. If itur=0, constant diffusivities are used for momentum and transport (and the values are specified in the next line). If itur=-2, vertically homogeneous but horizontally varying diffusivities are used, which are read in from hvd.mom.and hvd.tran. If itur=-1, horizontally homogeneous but vertically varying diffusivities are used, which are read in from vvd.dat. If itur=2, the zero-equation Pacanowski and Philander closure is used. If itur=3, then the two-equation closure schemes (Mellor-Yamada-Galperin, K-epsilon, Umlauf and Burchard etc.) are used.

If itur=0, the next line is: vdiff, tdiff = constant diffusivities for momentum and transport.

If itur=2, the next line is: hestu_pp, vdmax1, vdmin1, tdmin1, hcont_pp, vdmax2, vdmin2, tdmin2. Eddy viscosity is computed as: vdiff=vdiff_max/(1+rich)^2+vdiff_min, and diffusivity tdiff=vdiff_max/(1+rich)^2+tdiff_min, where rich is a Richardson number. The limits (vdiff_max, vdiff_min and tdiff_min) vary linearly with depth between depths hestu_pp and hcont_pp.

If itur=3, the next line is:

mid, stab: choice of model description ("MY"-Mellor & Yamada, "KL"-GLS as k-kl, "KE"-GLS as k-epsilon, "KW"-GLS as k-omega, or "UB"-Umlauf & Burchard's optimal), and stability function ("GA"-Galperin's, or "KC"-Kantha & Clayson's for GLS models). In this case, the minimum and maximum viscosity/diffusivity are specified in diffmin.gr3 and diffmax.gr3, and the surface mixing length is specified in xlfs.gr3.

icst = options for specifying initial temperature and salinity field for cold start. If icst=1, a vertically homogeneous but horizontally varying initial temperature and salinity field is contained in temp.ic and salt.ic. If icst=2, a horizontally homogeneous but vertically varying initial temperature and salinity field, prescribed in a series of z-levels, is contained in ts.ic. For general 3D initial S,T fields, use the hot start option.

ntip, tip_dp: # of constituents used in earth tidal potential; cut-off depth for applying tidal potential (i.e., it is not calculated when depth < tip_dp).

For k=1, ntip

talpha(k) = tidal constituent name;

jspc(k), tamp(k), tfreq(k), tnf(k), tear(k) = tidal species # (0: declinational; 1: diurnal; 2: semi-diurnal), amplitude constants, frequency, nodal factor, earth equilibrium argument (in degrees);

end for;

nbfr = total # of tidal boundary forcing frequencies.

For k=1, nbfr

alpha(k) = tidal constituent name;

amig(k), ff(k), face(k) = forcing frequency, nodal factor, earth equilibrium argument (in degrees) for constituents forced on the open boundary;

end for;

nope: # of open boundary segments;

For j=1, nope

neta(j), iettype(j), ifltype(j), itetype(j), isatype(j) = # of nodes on the open boundary segment j (from hgrid.gr3), b.c. flags for elevation, normal velocity, temperature, and salinity;

if (iettype(j) == 1) !time history of elevation on this boundary

no input in this file; time history of elevation is read in from elev.th;

else if (iettype(j) == 2) !this boundary is forced by a constant elevation

ethconst: constant elevation

else if (iettype(j) == 3) !this boundary is forced by tides

for k=1, nbfr

alpha(k) = tidal constituent name;

for i=1, neta(j)

emo(ietaelem(j,i),k), efa(ietaelem(j,i),k) !amplitude and phase for each node on this open boundary;

end for

end for;

else if (iettype(j) == 4) !space- and time-varying input

no input in this file; time history of elevation is read in from elev3D.th;

else if (iettype(j) == 0)

elevations are not specified for this boundary (in this case the velocity must be specified).

endif

if (ifltype(j) == 0) !nornal vel. not specified

no input needed

else if (ifltype(j) == 1) !time history of discharge on this boundary

no input in this file; time history of discharge is read in from flux.th;

else if (ifltype(j) == 2) !this boundary is forced by a constant discharge

vthconst: constant discharge (note that a negative number means inflow)

else if (ifltype(j) == 3) !vel. is forced in frequency domain

for k=1, nbfr

vmo(j,k),vfa(j,k) !uniform amplitude and phase along each boundary segment

end for;

else if (ifltype(j) == 4) !3D input

no input in this file; time history of discharge is read in from uv3D.th;

else if (ifltype(j) == -1) !Flanther type radiation b.c. (iettype must be 0 in this case)

eta_m0,qthcon(j): mean elevation and discharge on the jth boundary

endif

if (itetype(j) == 0) !temperature not specified

no input needed

else if (itetype(j) == 1) !time history of temperature on this boundary

no input in this file; time history of temperature is read in from temp.th;

else if (itetype(j) == 2) !this boundary is forced by a constant temperature

tthconst = constant temperature

else if (itetype(j) == 3) !keep initial temperature profile

no input is needed

else if(itetype(j) == -1) !open b.c.; nudge to initial condition

tobc: nudging factor (between 0 and 1).

else if(itetype(j) == 4) !3D input

no input in this file; time history of temperature is read in from temp3D.th;

endif

Salintiy b.c. is similar to temperature:

if (isatype(j) == 0) !salinity not specified

.........

endif

nspool, ihfskip: Global output skips. Output is done every nspool steps, and a new output file is opened every ihfskip steps (and in addition, a hotstart file is output at the same time step if the flag nhstar is turned on below). Therefore the outputs are named as [1,2,3,...]_salt.63 etc.

next 22 lines are global output (in machine-dependent binary) options. The outputs share the same structure. Only the first line is detailed here.

noutge = global elevation output control. If noutge=0, no global elevation is recorded. If noutge= 1, global elevation for each node in the grid is recorded in [1,2,3...]_elev.61 in binary format. The output is either starting from scratch or appended to existing ones depending on ihot.
output options for atmospheric pressure (*pres.61).
output options for air temperature (*airt.61).
output options for specific humidity (*shum.61).
output options for solar radiation (*srad.61).
output options for short wave radiation (*flsu.61).
output options for long wave radiation (*fllu.61).
output options for upward heat flux (*radu.61).
output options for downward flux (*radd.61).
output options for total flux (*flux.61).
output options for wind speed (*wind.62).
output options for wind stresses (*wist.62).
output options for horizontal velocity (*hvel.64).
output options for vertical velocity (*vert.63).
output options for temperature (*temp.63).
output options for salinity (*salt.63).
output options for density (*conc.63).
output options for diffusivity (*tdff.63).
output options for turbulent kinetic energy (*kine.63).
output options for macroscale mixing length (*mixl.63).
output options for z coordinates at each node (*zcor.63). This output will always be there even if the flag is turned off.
output options for depth-averaged velocity (*dahv.62).
output options for test variable (test.60). The user may choose any internal variable by modifying the source code

nhstar= hot start output control parameter. If nhstar=0, no hot start output is generated. If nhstar=1, hot start output is spooled to it_hotstart every ihfskip time steps, where it is the corresponding time iteration number. If a run needs to be hot started from step it, the user can copy it_hotstart to hotstart.in, as the code expects the hot start input file to be hotstart.in.

isolver, itmax1, iremove, zeta, tol = ITPACK solver control parameters.

· If isolver=1, the Jacobian Conjugate Gradient Method is used (recommended);

· If isolver=2, the Jacobian Semi-Iteration Method is used;

· If isolver=3, the Successive Over-relaxation Conjugate Gradient Method is used;

· If isolver=4, the Successive Over-relaxation Semi-Iteration Method is used;

Recommended values: isolver=1, itmax1=1000, iremove=0, zeta=5.e-6, tol=1.e-13.

iflux= parameter for checking volume and salt conservation. If turned on (=1), the conservation will be checked in regions specified by fluxflag.gr3.

lq_t, lq_s, int_mom: linear (1) or quadratic (2) interpolation for T, S in backtracking, and linear (int_mom=0) or linear with element splitting (int_mom=1) for backtracking velocity. If lq_t (or lq_s)=0, the option is read in from lq_t.gr3 (or lq_s.gr3) as the depth info. Note that int_mom=1 is currently not recommended.

h_bcc1: depth thresholds used in calculation of bariclonic force. The force will be evaluated using the pressure Jacobian method when h<= h_bcc1 (otherwise the Z-method is used where interpolation back to Z-plane is done). We usually recommend the Jacobian method, and in this case, a very large h_bcc1 (e.g., 1.e6m) can be used.

islip: option for land b.c. islip=0: free slip; =1: no slip.

inu_st, step_nu, vnh1,vnf1,vnh2,vnf2: nudging flag for S,T, nudging step, parameters for vertical nudging. When inu_st=0, no nudging is done. When inu_st=1, nudge to initial conditions. When inu_st=2, nudge to values specified in temp_nu.in and salt_nu.in, given at an interval of step_nu. For inu_st/=0, the horizontal nudging factors are given in t_nudge.gr3 and s_nudge.gr3 (as depths info), and the vertical nudging factors vary linearly along the depth as: min(vnf2,max(vnf1,vnf)), where vnf=vnf1+(vnf2-vnf1)*(h-vnh1)/(vnh2-vnh1). The nudging factor is the sum of the two.

xlmax00: maximum mixing length (e.g., 10m).

mmm: order of vertical integration used in calculating the baroclinic pressure gradient for pressure Jacobian method. Must be between 0 and 3. Use mmm=0 for best efficiency, or higher value for accuracy with significant computational cost.

idrag: bottom drag option. idrag=1: linear drag formulation; idrag=2: quadratic drag formulation (default).

ielad: used only in upwind option (see below). If ielad=0, the time stepping continues even if the ELAD correction for over- and -undershoots does not converge. If ielad=1, the time stepping will stop.

ihhat: wet/dry option. If ihhat=1, the friction-reduced depth will be kept non-negative to ensure robustness (at the expense of accuracy); if ihhat=0, the depth is unrestricted.

iupwind_t, iupwind_s: upwind option for T,S. A value of "0" corresponds to Eulerian-Lagrangian transport option, and "1" for the mass-conservative upwind option.

If the upwind option is invoked for at least one of the variables, then this line is mbig: the number of sub-divisions in time step. The users need to estimate the maximum Courant number and choose the smallest mbig possible for efficiency. The maximum Courant number in each step is output in mirror.out. Upwind option also requires the horizontal grid to be free of very skinny elements (otherwise a large mbig must be used).

go back

Initial condition for S,T

Depending on the values of icst (see parameter input file):

If icst = 1, salt.ic and temp.ic take the following build point format:

33 to 0 from outh to TP. !file description
30001 !# of nodes;
1 386380.409604 286208.187634 0.000000 !node #, x, y, initial salinity (or temperature)
2 386460.352736 285995.038877 0.000000
3 386687.720000 286213.590000 0.000000
4 386460.076848 286367.779818 0.000000
5 386678.380000 286483.440000 0.000000
..............
If icst = 2, ts.ic takes the following format:

43 !total # of vertical levels;
1 -2000. 4. 34. !level #, z-coordinates, temperature, salinity
2 -1000. 5. 34.
3 -500.   6. 33.
4 -200.   7. 32.
5 -100.   8. 31.
...........

go back

Bottom drag (drag.gr3 or rough.gr3)

grid !file decription
40000 27918 !# of elements, # of nodes
1 386738.500000 285939.060000 0.004500 !node #, x, y, drag coefficient Cd (for nchi=0) or roughness (in meters; for nchi=1)
2 386687.720000 286213.590000 0.004500
3 386421.090000 286172.160000 0.004500
4 386471.720000 286376.030000 0.004500
5 386678.380000 286483.440000 0.004500
6 386140.190000 286439.220000 0.004500
7 386387.280000 286557.310000 0.004500
8 386209.840000 286676.470000 0.004500
..........

go back

wind (wind.th)

If nws=1 in param.in, a time history of wind speed must be specified in this file:

5. 8.660254 ! x and y components of wind speed @ 0*wtiminc
5. 8.660254
5. 8.660254

.......

Note that the speed varies linearly in time, and the time interval (wtiminc) is specified in param.in.

go back

Time history input:

This includes elev.th, flux.th, temp.th, salt.th, which share same structure. Below is a sample flux.th:

300. -1613.05005 -6186.60156 !time (in sec), discharge at the 1st boundary with ifltype=1, discharge at the 2nd boundary with ifltype=1
600. -1611.37854 -6208.62549
900. -1609.39612 -6232.22314
1200. -1607.42651 -6254.24707
1500. -1605.45703 -6276.27148
1800. -1603.48743 -6298.2959
2100. -1601.3772 -6321.89307
2400. -1599.40772 -6343.91748
2700. -1597.43811 -6365.94141
3000. -1595.46863 -6387.96582
3300. -1593.49902 -6409.99023
3600. -1591.38879 -6433.5874
3900. -1589.41931 -6452.94287
4200. -1587.2959 -6472.29834

...........

go back

Space- and time-varying time history inputs:

These include elev3D.th, uv3D.th, temp3D.th, salt3D.th, which share similar structure. For example, uv3D.th:

for it=1,nt !all time steps

time stamp (in sec);

for i=1,nope !all open boundary segments that have this type of b.c.

for j=1,nond(k) !all nodes on this boundary

node #, (uth(i,j,k),vth(i,j,k),k=1,nvrt);

end for !j

end for !i

end for !it

go back

obe.out

This file is generated with the pre-processing flag in param.in for debugging purpose only.

3 # of open bnd
Element list:
251 bnd # 1
1 31587
2 31588
3 31589
4 31590
5 31592
6 31595
7 31601
8 31603
9 31605
10 31606

........

4 bnd # 2
1 31583
2 31584
3 31585
4 31586
........

go back

adv.gr3

If nadv=0, the advection on/off flags are the "depths" (0: off; 1: Euler; 2: 5th order Runge-Kutta) in this grid file, which is otherwise similar to hgrid.gr3.

go back

hvis.gr3

The depth specifies the horizontal viscosity used at each node. We recommend using 0 for most cases. If a non-zero viscosity is needed, it must be between -0.5 and 0.5 for stability reason. It usually should be on the order of 0.01 or so.

go back

Min. and max. diffusivity

The depth specifies the minimum and maximum diffusivity imposed at each node. This is needed to further constraint outputs from the GLS model. We generally recommend a constant value of 1.e-6 m^2/s for minimum diffusivity, and 1.e-2 m^2/s for maximum diffusivity inside estuaries and 10 m^2/s otherwise. The minimum diffusivity may also be changed locally, e.g., to create a mixing pool near the end of a river.

go back

Surface mixing length

The depth specifies the surface mixing length used in GLS model in meters. We recommend using a constant value equal to the maximum surface layer thickness in the domain, but the users

may try smaller values as well.

go back

Interpolation mode

The depth specifies the way the vertical interpolation is done locally, i.e., along Z or S plane. If the depth=1, it is done along Z-plane; if the depth=2, along S-plane. We recommend the value of 2 in "pure S" zone, and 1 in SZ zone. You may not use "2" in the SZ zone.

go back

Lat/long coordinates (hgrid.ll)

This file is identical to hgrid.gr3 except the x,y coordinates are replaced by lattitudes and longitudes.

go back

Heat exchange

This consists of a suite of input for wind and radiation fluxes found in a sub-directory hdf/. When nws=2, the wind speed and atmospheric pressure are read in from this directory; when ihconsv=1, various fluxes are read in from it as well. The hdf files for various periods have been pre-computed by Mike Zulauf and deposited in a data base. Running his script inside hdf/ generates links to this data base for specified period:

/home/workspace/ccalmr/mazulauf/scripts/make_wind_links.csh yyyy d1 d2

where yyyy is the 4-digit year, d1 and d2 are the starting and ending Julian dates.

Note that Mike has recently added an alternative route for preparing the inputs (sflux_alt8d.f90; see Makefile). The code is a template the user can use to read in the input files in other formats.

go back

Conservation check files (fluxflag.gr3)

Mass conservation (fluxflag.gr3)

The "depths" of this grid file divide the whole domain into 3 regions: 0, 1 and 2, where regions "1" and "2" must be adjacent to each other. The total volume and various fluxes in regions "1" and "2" are computed.

go back

vvd.dat, hvd.mom, and hvd.tran

vvd.dat (vgrid.in format):

43 !total # of vertical levels;
1 1.e-2 1.e-4 !level #, viscosity, diffusivity
2 1.e-2 1.e-4
3 1.e-2 1.e-4
.........................
hvd.mom & hvd.tran (build point format)

hvd.mom !file decription
27918 !total # of nodes
1 386738.500000 285939.060000 0.0045 !node #, x, y, viscosity/diffusivity
2 386687.720000 286213.590000 0.004
..............................................

go back

Amplitudes and phases of boundary forcings

To generate amplitudes and phases for each node on a particular open boundary , follow README in ambcs01:~yinglong/AMB24Stuff/ElcircScripts/TIDES2/.

go back

Nodal factor and equilibrium arguments

Scripts (tid_e) and sample input (tide_colu.com) can be found in ambcs01:~yinglong/AMB24Stuff/ElcircScripts/TIDES/NodalFactor/. There is also a README there.

To generate nodal factor and equilibrium argument at a particular time, first change the time in the first line of tide_colu.com; e.g., for April 30, 2001:

0830040120 !Hour (PST), Day, Month, Year(1995 -> 9519) ; note the 8 hrs difference between 00GMT and 00PST

Then run the script on amb24 as: tid_e < tide_colu.com > out, and the file out contains nodal factors and arguments for all constituents.

go back

Hotstart input

This file is always in direct-access binary format, and all integers (i.e., those beginning with i-n) occupy nbyte=4 bytes, and all real variables are in double precision (8 bytes).

Total record length ihot_len=nbyte*(3+(6*nvrt+1)*ne+(8*nvrt+1)*ns+3*np+20*np*nvrt+1)+12.

The variables in order are (beware the order of i,j in each variable):

time,iths, (idry_e(i), (we(j,i), tsel(j,i,1), tsel(j,i,2), j=1,nvrt), i=1,ne), (idry_s(i), (su2(j,i),sv2(j,i),tsd(j,i),ssd(j,i),j=1,nvrt), i=1,ns), (eta2(i),idry(i), (tnd(j,i),snd(j,i),tem0(j,i),sal0(j,i),q2(i,j),xl(i,j),dfv(i,j),dfh(i,j),dfq1(i,j),dfq2(i,j), j=1,nvrt), i=1,np), ifile, ifile_char

where (eta1(i),eta2(i), (we(j,i),j=1,nvrt),i=1,ne) is equivalent to:

do i=1,ne

eta1(i)

eta2(i)

do j=1,nvrt

we(j,i)

enddo

etc.

and ifile_char is a 12-character string corresponding to ifile. See the source code for the meaning of each internal variable.

go back

Bed deformation input (bdef.bp)

This file is needed if imm=1 (tsunami option) and is in a build-point format:

"alphanumeric description";

np: # of build points (same as # of nodes in hgrid.gr3);

do i=1,np

i,x(i),y(i),bdef(i) !bdef(i) is the local defomration (positive for uplift)

enddo

go back

Output files

Global output

There are 4 types of output in SELFE, which correspond to the following 4 types of suffixes:

*.61: 2D scalar - no vertical structure (elevation and 9 other variables used in the heat exchange model);
*.62: 2D vector - no vertical structure: wind speed (u,v) and stress (tauxz,tauyz);
*.63: 3D scalar - has vertical structure (vertical vel., temperature, salinity, density, z-coordinates, diffusivity, turbulent kinetic energy and mixing length);
*.64: 3D vector - has vertical structure: horizontal vel.

All output variables are defined on hgrid.gr3, i.e. @ nodes and in binary format. The header part contains grid and other useful info:

Data format description (char*48): e.g., 'DataFormat v5.0'
version (char*48): version of SELFE;
start_time (char*48): start time of the run;
variable_nm (char*48): variable description;
variable_dim (char*48): '2D(3D) scalar(vector)'
# of output time steps (int), output time step (real), skip (int), ivs (=1 or 2 for scalar or vector), i23d (=2 or 3 for
2D or 3D);

Vertical grid part:

nvrt: toatal # of vertical levels;
kz: # of Z-levels;
h0, h_s, hc, theta_b, theta_f: constants used in defining S-Z hybrid grid;
do k=1,kz-1
ztot(k): z-coordinates of each Z-level;
enddo
do k=1,nvrt-kz+1

sigma(k): sigma-coordinates of each S-level;

enddo
Horizontal grid part:

np: # of nodes;
ne: # of elements;
do m=1,np
x(m)
y(m)
h(m): depth
kbp00(m): initial bottom indices
enddo
do m=1,ne
i34: element type (currently must be 3)
do mm=1,i34
nm(m,mm): node # (connectivity table)
enddo

enddo

The header is followed by time iteration part:

do it=1,nt

time (real);
it: iteration #;
do i=1,np
eta(i): surface elevation of node i;
enddo
do i=1,np
if(i23d=2) then !2D output
     if(ivs.eq.1) out1
     if(ivs.eq.2) out1,out2
else !i23d=3 !3D output
     do k=max(1,kbp00(i)),nvrt
           if(ivs.eq.1) out1
           if(ivs.eq.2) out1,out2
     enddo !k
endif
enddo !i

enddo !it

These structures can also be seen in the simple I/O utility code read_output*.f90 included in the package.

go back

Warning and fatal messages

Warning message (fort.12) contains non-fatal warnings, while fatal message file (fort.11) is useful for debugging.

go back

Run info output (mirror.out)

This is a mirror image of screen output and is particularly useful when the latter is suppressed with nscreen=0. Below is a sample:

There are 85902 sides in the grid...
done computing geometry...
done classifying boundaries...
You are using baroclinic model
Check slam0 and sfea0 as variable Coriolis is used
Warning: you have chosen a heat conservation model
which assumes start time at 0:00 PST!
Last parameter in param.in is mnosm= 0
done reading grids...
done initializing outputs
done initializing cold start
hot start at time= 0.00000000000000D+000 0

calculating grid weightings for wind_file_1

calculating grid weightings for wind_file_2

wind file starting Julian date: 127.000000000000
wind file assumed UTC starting time: 8.00000000000000
done initializing variables...
time stepping begins... 1 2016
done computing initial levels...
Total # of faces= 1914122
done computing initial nodal vel...
done computing initial density...

calculating grid weightings for rad fluxes

rad fluxes file starting Julian date: 127.000000000000
rad fluxes file assumed UTC starting time: 8.00000000000000

..............................................

go back