SELFE 1.5k7 User Manual
Input files
"!" 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;
mne_ke: max. # of elements used in local Kriging;
nope: total # of open boundary segments;
ntracers: total # of passive tracers.
Horizontal grid (hgrid.gr3)
In xmgredit grid format.
Here is an example; below are explanations of format with this grid:
hgrid.gr3 ! alphanumeric description 60356 31082 ! of elements and nodes in the horizontal grid !Following is coordinate
part: 1 402672.000000 282928.000000 2.0000000e+01
! node #, x,y, depth ............................................. 31082 331118.598253 112401.547031 2.3000000e-01 !last node! Following is connectivity part: ........................................... |
!Following is boundary condition part (needed for hgrid.gr3 only):
3 : Number of open boundaries
95 ! Total number of open boundary nodes
3 ! Number of nodes for open boundary 1
29835 ! first node on this segment
29834 ! 2nd node on this segment
29836 ! 3rd node on this segment
........................................
2 ! last node on this open segment
16 ! number of land boundaries
1743 ! Total number of land boundary nodes (including islands)
753 0 ! Number of nodes for land boundary 1 ('0' means the exterior boundary)
30381 ! first node on this segment
30380
.......................................
.......................................
10 1 ! Number of nodes for island boundary 1 ('1' means island)
29448 ! first node on this island
29451
29525
.......................................
.......................................
Note: (1) the boundary condition (b.c.) part can be generated with xmgredit5 --> GridDEM --> Create open/land boundaries;
(2) if you have no open boundary, you can create two land boundary segments that are linked to each other;
(3) although not required, we recommend you follow the following convention when genrating the boundary segments. For the exterior boundary (open+land), go in counter-clockwise direction; for interior boundaries (islands), go in clockwise direction;
(4) The use of a Shapiro filter places some constraints on the boundary sides. In particular, the center of any internal sides must be inside the quad formed by the centers of its 4 adjacent sides (see Fig. 1). If not, the code will try to enlarge the stencil, but if the side is near the boundary, fatal error will occur. To find out all violating boundary elements, just prepare hgrid.gr3 (note that the open boundary info needs no be correct at this stage), vgrid.in and param.in up to ihorcon, and run the code with ipre=1 in param.in. You'll find a list of all such sides in fort.11 (the two node numbers of a side will be shown). Method to eliminate this problem includes: (1) move node, (2) swap the side for a pair of elements, and (3) refine or coarsen. Note that in most grid editors, the first 2 methods won't change the node numbering and so you'd try them first before Method (3), to save time. After the pre-processing run is successful (with a screen message indicating so), you can then proceed to prepare other inputs.
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. !
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,
20 -0.944444
.......
54 0. !last
Notes:
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).
ihorcon: flag to
use non-zero horizontal viscosity. If
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:
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).
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.
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
netcdf
files in the directory sflux/.
If
nws>0, the next two line are: (1) nrampwind, drampwind = ramp option and period
(in days) for wind.
ihconsv,
isconsv = heat budget and salt conservation models flags. If ihconsv=0, the heat budget model is not used. If
ihconsv=1, a heat budget model is invoked, and a number of
netcdf files for radiation flux
input are read in from he directory sflux/.
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=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 xlsc.gr3.
If itur=4, GOTM turbulence model is invoked (recommended). But the user need
to compile the GOTM libraries first (see FAQ or README inside GOTM/ for
instructions), and turn on
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.
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
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;
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) ! nudge to 3D time series in temp3D.th
tobc: nudging factor (between 0 and 1).
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
25+ 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 precipitation rate (*prcp.61).
output options for wind speed (*wind.62).
output options for wind stresses (*wist.62).
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 eddy diffusivity (*tdff.63).
output options for eddy viscosity (*vdff.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 horizontal velocity (*hvel.64).
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
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,
int_mom: linear (lq=1) or quadratic (lq=2) interpolation for T, S in backtracking,
and linear (
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
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.
inactive line.
mmm:
order of vertical integration used in calculating the baroclinic
pressure gradient for
idrag: bottom drag option. idrag=1: linear drag formulation; idrag=2:
quadratic drag
inactive line.
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,
vis_coe1,
rmaxvel: maximum velocity. This is needed mainly for the air-water exchange as the model may blow up if the water velocity is above 20m^2/s.
ntr: # of passive tracers as defined in the global module. If ntr>0, additional lines are needed that specify the transport method (upwind or TVD) and horizontal boundary conditions etc. Consult the source code for details. See FAQ for how to interface your own code to SELFE.
Depending on the values of icst (see parameter input file):
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
..........
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.
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
...........
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
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
........
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.
The depth specifies the Kriging option for each node: 0 means no Kriging; 1 means applying Kriging there. The order of the generalized covariance function is specified
in param.in.
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.
The depth specifies the surface mixing length used when itur=3. It specifies the portion of surface layer thickness; e.g., 0.5 mean 1/2 of the layer is used as surface mixing length. Recommend value: 0.5.
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.
Lat/long coordinates (hgrid.ll)
This file is identical to hgrid.gr3 except the x,y coordinates are replaced by lattitudes and longitudes.
This consists of a suite of input for wind and radiation fluxes found in a sub-directory sflux/. 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 netcdf files for various periods have been pre-computed by Mike Zulauf and deposited in a data base.
Conservation check files (fluxflag.gr3)
vvd.dat, hvd.mom, and hvd.tran
Amplitudes and phases of boundary forcings
To generate amplitudes and phases for each node on a particular open boundary, see SELFE Utilites.
Nodal factor and equilibrium arguments
See SELFE Utilites.
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+4*nvrt*ntracers+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), (trel0(j,i,l),trel(j,i,l),l=1,ntracers),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
enddo
etc.
and ifile_char is a 12-character string corresponding to ifile. See the source code for the meaning of each internal variable.
Bed deformation input (bdef.gr3)
This file is needed if imm=1 (tsunami option) and is in a grid format:
"alphanumeric description";
ne, np: same as in hgrid.gr3;
do i=1,np
i,x(i),y(i),bdef(i) !bdef(i) is the local defomration (positive for uplift)
enddo
Output files
There are 4 types of output in SELFE, which correspond to the following 4 types of suffixes:
All output variables are defined on hgrid.gr3, i.e. @ nodes and in binary
format. Please consult the script read_output*.f90 for a complete description of
the format. The header
part contains grid and other useful info:
Vertical grid part:
sigma(k): sigma-coordinates of each S-level;
enddo
Horizontal grid part:
enddo
The header is followed by time iteration part:
do it=1,ntime
enddo !it
These structures can also be seen in the simple I/O utility code read_output*.f90 included in the package.
Warning message (fort.12) contains non-fatal warnings, while fatal message file (fort.11) is useful for debugging.
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
..............................................