Elcirc 4.01 User Manual
Input files
np: # of nodes in the horizontal grid;
ne: # of elements in the horizontal grid;
ns: # of sides in the horizontal grid;
nvrt: # of levels in the vertical grid
Horizontal grid (hgrid.gr3)
In gredit 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)
2 3 3 4 1
3 3
3 5 4
4 3 1 4 6
5 3 4 7 6
...........................................
(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 boudnaries.
Vertical grid (vgrid.in)
43 4825.1 : total # of vertical levels, z-coordinate
of MSL (mean sea level).
1 3667.00 3667.00 : level #, layer
thickness, z-coordinate of level line
2 1000.00 4667.00
3 100.00
4767.00
4 10.00 4777.00
5 9.40 4786.40
6 5.00 4791.40
7 5.00
4796.40
8 3.00 4799.40
9 2.00 4801.40
10 2.00 4803.40
11 1.50
4804.90
12 1.00 4805.90
13 0.80 4806.70
14 0.80 4807.50
15 0.80
4808.30
16 0.80 4809.10
17 0.80 4809.90
18 0.80 4810.70
19 0.80
4811.50
20 0.70 4812.20
21 0.70 4812.90
22 0.70 4813.60
23 0.70
4814.30
24 0.70 4815.00
25 0.70 4815.70
26 0.70 4816.40
27 0.70
4817.10
28 0.70 4817.80
29 0.70 4818.50
30 0.70 4819.20
31 0.70
4819.90
32 0.70 4820.60
....................................................
Notes:
Explanation of each line:
48-character string
description of the version.
48-character start time
info string.
ipre: pre-processing flag.
=1: output centers.bp,
sidecenters.bp and obe.out
(centers build point, sidcenters build point, and list of open boundary
elements), and stop. =0: skip this part.
nscreen = screen output
on/off switch (0: off; 1: on). In either case, mirror output messages will be
directed into mirror.out.
ihot =
hot start flag. If ihot=0, cold start; if ihot/=0, hot start from hotstart.in.
ics = coordinate frame
flag. If ics=1, Cartesian coordinates are used; if ics=2, degrees
latitude/longitude are used.
slam0, sfea0 = centers of
projection used to convert lat/long to Cartesian coordinates. These are only
used for ics=2, or a variable Coriolis parameter is employed (ncor=2).
theta0 = implicitness
parameter (between 0.5 and 1).
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 ramp-up period in days (not used if nramp=0).
dt = external time step(in
sec) for momentum and continuity equations.
nsubfl = 0, 1, or 2: flag
to determine how the number of subdivisions in backtracking is computed. For
nsubfl=0, a constant value, specified below, is used throughout the domain;
for nsubfl=1, the number of subdivisions is read in from ndelt.gr3;
for nsubfl=2, the number of subdivisions is automatically calculated in the
code based on the local velocity gradient, subject to the max. and min.
specified below.
If nsubfl=0, the next line
is: ndelt = constant number of subdivisions.
If nsubfl=2, the next line
is: ndeltmin, ndeltmax: minimum and maximum number of subdivisions
allowable.
nadv = advection on/off
switch. If nadv=1, advection is on for the whole domain (and this is the
default). If nadv=0, advection is selectively turned off based on the input
file adv.gr3.
h0 = minimum depth (in m)
(recommended value: 1cm). When the total depth is less than h0, the
corresponding sides/elements are considered as dry; it also controls the
minimum layer thickness. This should always be positive to prevent
underflow.
ntau = bottom friction
option. If ntau=0, a constant drag coefficient is used for bottom friction
parameterization. If ntau=1, a logarithmic law is applied to the
parameterization. If ntau=2, the drag coefficients are read in from drag.bp.
If ntau=0, the next line
is: Cd0 = constant drag coefficient.
ncor = Coriolis option. If
ncor=0, a constant Coriolis parameter is used. If ncor=2, a variable Coriolis
parameter, based on a beta-plane approximation, is used, with the lat/long.
coordinates of the grid 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,
spatially and temporally variable wind is applied and the input consists of a
number of hdf files customized for CORIE.
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 customized heat budget model is invoked for CORIE domain. The latter entails
a number of hdf files for radiation flux input.
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=1, a step
function model is used. If itur=2, the zero-equation Pacanowski and Philander
closure is used. If itur=3, then the two-equation Mellor-Yamada-Galperin
closure scheme is used (default scheme). If itur=-1,
horizontally homogeneous but vertically varying diffusivities are used, which
are read in from vvd.dat.
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: vdiff_max, vdiff_min, tp, tdiff_min. Eddy viscosity is computed as: vdiff=vdiff_max/(1+rich)^2+vdiff_min, and diffusivity tdiff=vdiff_max/tp/(1+rich)^2+tdiff_min, where rich is a Richardson number.
If itur=3, the next line
is: bgdiff, diffmax, hest_my, hcont_my, xlmin_est, xlmin_sea (background
diffusivity, maximum diffusivity, cut depths for estuary and open sea and
associatedmin. mixing lengths). The min. mixing length is computed as: L_min=
xlmin_est for h<hest_my; L_min= xlmin_sea for h>hcont_my; L_min=
xlmin_est + (xlmin_sea - xlmin_est) *
(h-hest_my)/(hcont_my-hest_my) for hest_my <= h<= hcont_my (h is
the local depth).
ihorcon, horcon:
horizontal diffusion option for momentum equation. If ihorcon=0, a constant
diffusion coefficient (horcon) is used. If ihorcon=1, Smagorinsky
paramertization is used and horcon is the dimensionless constant used in the
scheme.
Next 2 lines are not active at the moment. **********
******************************************
ictemp, icsalt = options
for specifying initial temperature and salinity field. If ictemp (or
icsalt)=1, a vertically homogeneous but horizontally varying initial
temperature (or salinity) field is contained in temp.ic
(or salt.ic).
If ictemp (or icsalt)=2, a horizontally homogeneous but vertically
varying initial temperature (or salinity) field is contained in temp.ic
(or salt.ic).
ispon_e, ispon_v,
ispon_st: sponge layer options for elevation, velocity, and S,T (experiment
stage). 0: off.
ntip = total # of tidal
potential forcing frequencies.
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 ocean boundary;
For j=1, nope
netaelem(j), iettype(j), ifltype(j), itetype(j), isatype(j) = # of elements on the open boundary segment j, 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, netaelem(j)
emo(ietaelem(j,i),k),
efa(ietaelem(j,i),k) !amplitude
and phase for each element on this open boundary;
end for
end for;
else
no input in this file; elevations are computed as average of surrounding eevations.
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
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.
tthconst = imposed temperature for inflow
endif
Salintiy b.c. closely follows 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 if the flag nhstar
is turned on
below). Therefore the output is named as [1,2,3,...]_salt.63 etc.
next 20 lines are global
output options. They share the same
structure, and thus 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 n_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 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 test variable.
nhstar= hot start output control parameter. If nhstar=0, no hot start output (it_hotstart) is generated. If nhstar=1, hot start output is spooled to it_hotstart every ihfskip time steps, where it is the corresponding time step stamp. If a run needs to be hot started from step it, the user can create a symbolic link of hotstart.in to it_hotstart as the code expects the hot start input file to be the former.
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 for itmax1=1000, iremove=0, zeta=5.e-6, tol=1.e-13.
iflux, ihcheck = parameter
for checking volume, heat and salt budget balances. If turned on (=1), the
conservation will be checked in a region specified, respectively, by fluxflag.gr3
and hcheck.gr3.
iwmode: flag to decide
either the continuity eq. or the vertical momentum eq. is used to solve for
w.
iwmode=1: use
continuity eq;
iwmode=2: use
vertical momentum eq;
iwmode=3:
use one of the 2 eqs. based on the flags (1 or 2) specified in wswitch.gr3.
nsplit, dvel00, mnosm = mode splitting factor (i.e., the ratio between internal and external time steps), minimum vel. gradient for smoothing, max. # of iteration in smoothing (in general, no smoothing should be used unless absolutely necessary, and so one should set mnosm=0 and the value of dvel00 is not important in this case).
Depending on the values of icsalt and ictemp (see parameter input file):
If ntau=2 in param.in, this input is needed. It takes the form of a build point file:
0.0025 to 0.0045 from TP to Woody !file decription
27918
!total # of nodes
1 386738.500000 285939.060000 0.004500 !node #, x, y, drag
coefficient
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
5.
8.660254
5. 8.660254
.......
Note that the speed varies linearly in time, and the time interval 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
...........
This file is generated with the pre-processing flag in param.in, and info contained here is needed in param.in (e.g., # of open boundary elements) or in calculating tidal amplitudes and phases.
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
........
This is needed when nsubfl=1. It is the same as hgrid.gr3 except the depth indicates the # of subdivisions used in backtracking.
If nadv=0, the advection on/off flags are the "bathymetry" (0: off; 1: on) in this grid file, which is otherwise similar to hgrid.gr3.
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 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.
When iwmode=3, the "depth" (either 1 or 2) of this grid file indicates the method used to solve for w. If depth=1, the continuity eq. is used; if depth=2, the vertical momentum eq. is used.
Conservation check files (fluxflag.gr3 and hcheck.gr3)
vvd.dat, hvd.mom, and hvd.tran
Amplitudes and phases of boundary forcings
To generate amplitudes and phases for each element on a particular open boundary , follow these steps (scripts and sample input can be found in amb24:~yinglong/ElcircScripts/TIDES/):
Nodal factor and equilibrium arguments
Scripts (tid_e) and sample input (tide_colu.com) can be found in amb24:~yinglong/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.
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). Depending on the value of the turbulence closure flag itur, it assumes 2 different forms:
If itur=3:
Total record length
ihot_len=nbyte*(3+4*ne+2*ne*(nvrt+1)+4*ns*(nvrt+1)+4*ns*nvrt+3*np+7*np*(nvrt+1)+8*np*nvrt+1)+12
The variables in order are:
time,iths,(eta1(i),eta2(i), (we(i,j),j=0,nvrt),i=1,ne),
((vn2(i,j),vt2(i,j),j=0,nvrt),(tsd(i,j),ssd(i,j),j=1,nvrt),i=1,ns)
,(peta(i),ibad(i),(uu1(i,j),vv1(i,j),ww1(i,j),nosm(i,j),j=0,nvrt),
(tnd(i,j),snd(i,j),tem0(i,j),sal0(i,j),j=1,nvrt),i=1,np),((q2(i,j),xl(i,j),j=0,nvrt),
i=1,ns),ifile,ifile_char
where (eta1(i),eta2(i), (we(i,j),j=0,nvrt),i=1,ne) is equivalent to:
do i=1,ne
eta1(i)
eta2(i)
do j=0,nvrt
we(i,j)
enddo
enddo
etc.
and ifile_char is a 12-character string corresponding to ifile.
If itur/=3:
Delete q2 and xl in the list, i.e., ((q2(i,j),xl(i,j),j=0,nvrt), i=1,ns). The total record length ihot_len is reduced by ns*(nvrt+1)*16.
Output files
There are 4 types of output in Elcirc4.01, which correspond to the following 4 types of suffixes:
Every output variable is defined on hgrid.gr3, i.e. nodes. The output can be
in either ASCII or binary format, only the latter being described here. The
header part contains grid and other useful info:
Vertical grid part:
Horizontal grid part:
Time iteration part:
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
..............................................