forked from geodynamics/hc
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathREADME.TXT
387 lines (243 loc) · 11.8 KB
/
README.TXT
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
README FOR HC
Version 1.0.15 - 2022
Thorsten Becker - [email protected]
HC is a global mantle circulation solver following Hager & O'Connell
(1981) which can compute velocities, tractions, and the geoid for
arbirtrary density distributions and radial viscosity variations.
This particular implementation illustrates one possible way to combine
the HC solver routines. There is a forward tool example, as included
in the graphical user interface SEATREE, and a parameter space
exploration tool to scan viscosity distributions in terms of their
predictions of the geoid.
This code is based on Fortran routine by Brad Hager, Richard
O'Connell, and most recently Bernhard Steinberger, who made several
modifications. This version is by Thorsten Becker with contributions
to the plate velocity inversion tool by Craig O'Neill (the latter
remains incomplete).
The Solid Earth Teaching and Research Environment (SEATREE,
https://github.com/thwbecker/seatree/), as preinstalled with
software as a VirtualBox available from the Unified Geodynamics Earth
Science Computing Environment (UGESCE,
http://www-udc.ig.utexas.edu/external/becker/ugesce.html) provides a
convenient graphical user interface for HC to select seismic
tomography models, edit viscosity structures, etc.
INSTALLATION
The code should compile on any basic UNIX/Linux system with the GMT
tools (BY DEFAULT, VERSION > 4.2.1, AND ONLY VERSION 4.X.X IS
SUPPORTED!) and libraries installed. See the Makefile for comments and
what to modify.
What is described below is the hc Hager & O'Connell (1981) forward
flow and geoid computation. For plate velocity inversions following
the Ricard & Vigny (1989) method, see the hcplates subdirectory and
the README there. The plate inversion code is not quite working yet.
The Makefile assumes that the following environment variables are
predefined:
1) GMTHOME: needs to point to the base directory of the GMT installation,
e.g.
/usr/local/src/GMTdev/GMT4.5.18/
if you installed GMT yourself, or
/usr/lib/gmt
if you installed GMT with a package manager.
NOTE: ONLY GMT VERSION 4 IS SUPPORTED.
2) NETCDFHOME: for the netcdf libraries. This could be e.g.
/usr/local/src/netcdf-3.5.0/
or
/usr
3) HC_HOME: By default, the object files and binaries will be
installed in the "objects/" and "bin/" directories in the current
directory. If HC_HOME is set (e.g. /usr/local/), then they will
be installed in $HC_HOME/objects and $HC_HOME/bin.
4) ARCH: If you like, you can put the object files and binaries in
different directories depending on your architecture. This may be
useful if your directories are NFS mounted on different machines.
One way is to use the output of uname
setenv ARCH `uname -m | gawk '{print(tolower($1))}'` # sh/tcsh
export ARCH=`uname -m | gawk '{print(tolower($1))}'` # bash
On a 32 bit Intel machine, this will put the binaries in bin/i686.
Also, the Makefile uses the commonly defined compiler variables CC,
F90, CFLAGS, LD, and LDFLAGS. So to make static executables, set
LDFLAGS="-static"
By default, the Makefile is set up for the new syntax of GMT version
4.1.2 and higher. The alternative is to use GMT3, which can be done by
defining -DUSE_GMT3 and modifying the two corresponding lines in the
Makefile.include.
With all things set up, you should be able to type
make all
to compile the programs.
IF YOU DO NOT WANT TO COMPILE USING THE GTM4 LIBRARIES, use
make all_no_gmt
this will compile the main HC program, but not the sh_exp type of
routines, as those rely on the netcdf/grd I/O capabilities via GMT.
HC CODE usage
Described in the help page that is displayed for "hc -h" as below.
Also see SEATREE (https://github.com/thwbecker/seatree/) for a
graphical user interface, and example plotting scripts as provided
below.
Example input data is provided in subdirectory example_data/
>>>
hc - perform Hager & O'Connell flow computation
This code can compute velocities, tractions, and geoid for simple
density distributions and plate velocities using the semi-analytical
approach of Hager & O'Connell (1981). This particular implementation
illustrates *one possible way* to combine the HC solver routines,
hc_visc_scan (see below) another.
Based on code by Brad Hager, Richard O'Connell, and Bernhard
Steinberger. This version by Thorsten Becker and Craig O'Neill
usage example:
bin/hc -vvv
Compute mantle flow solution using the default input files:
viscosity profile visc.dat
density profile dens.sh.dat
earth model prem/prem.dat
and provide lots of output. Default setting is quiet operation.
See README.TXT in the installation directory for example for how to plot output, and
http://geosys.usc.edu/projects/seatree/ for a graphical user interface.
http://www-udc.ig.utexas.edu/external/becker/sdata.html for a VirtualBox install.
density anomaly options:
-dens name use name as a SH density anomaly model (dens.sh.dat)
All density anomalies are in units of 1% of PREM, all SH coefficients
in Dahlen & Tromp convention.
-dshs use the short, Becker & Boschi (2002) format for the SH density model (OFF)
-ds density scaling factor (0.2)
-dnp do not scale density anomalies with PREM but rather mean density (OFF)
-dsf file read depth dependent density scaling from file
(overrides -ds, OFF), use pdens.py to edit
Earth model options:
-prem name set Earth model to name (prem/prem.dat)
-vf name viscosity structure filename (visc.dat), use pvisc.py to edit
This file is in non_dim_radius viscosity[Pas] format
boundary condition options:
-fs perform free slip computation (ON)
-ns perform no slip computation (OFF)
-pvel name set prescribed surface velocities from file name (OFF)
The file (e.g. pvel.sh.dat) is based on a DT expansion of cm/yr velocity fields.
-vshs use the short format (only lmax in header) for the plate velocities (OFF)
-vdir velocities are given in files name/vel.1.ab to vel.140.ab for different times,
-140 to -1 Ma before present, where name is from -pvel
-vtime time use this particular time step of the plate velocities (-1)
solution procedure and I/O options:
-cbckl val will modify CMB boundary condition for all l > val with solver kludge (2147483647)
-ng do not compute and print the geoid (1)
-ag compute geoid at all layer depths, as opposed to the surface only
-rg name compute correlation of surface geoid with that in file "name",
this will not print out the geoid file, but only correlations (OFF)
-pptsol print pol[6] and tor[2] solution vectors (OFF)
-px print the spatial solution to file (OFF)
-rtrac compute srr,srt,srp tractions [MPa] instead of velocities [cm/yr] (default: vel)
-htrac compute stt,stp,spp tractions [MPa] instead of velocities [cm/yr] (default: vel)
-v -vv -vvv: verbosity levels (0)
<<<
OTHER BINARIES
hc_visc_scan Illustrates how to do a parameter space scan for
viscosities and compute geoid correlations on the
fly. drive_visc_scan runs a bunch of tests, and
plot_visc_scan visualizes the output.
KNOWN LIMITATIONS
The propagator matrix approach can be unstable for moderately high
maximum degrees. This behavior can be addressed partially by compiling
HC in quadruple precision, and kinematic boundary can be addressed
with a kludge, see -cbckl. If you compile with this trick and
quadruple precision (the default within SEATREE), you should be good
to up to L=127.
SPHERICAL HARMONICS FORMAT
(A) Regular (long) format, which allows for both scalar and vector
harmonics
All single layer spherical harmonics are in the fully normalized,
physical convention, e.g. Dahlen & Tromp (1998).
All spherical harmonics files have an SH_HEADER
lmax layer_i zlabel_i nlayer nrset type
lmax: maximum l of expansion
layer_i: layer number, from 0....nlayer
zlabel_i: depth label of this layer, in km (positive, from 0: surface
to 2871 for CMB)
nlayer: number of layers
nrset: number of spherical harmonic coefficient sets
type: 0 for Rick, 1 for HealPix etc. (will determine only the
internal representation, external all is physical
convention)
1) Scalar fields with layers (e.g.: hc_assign_density, which calls:
sh_read_parameters, sh_init_expansion, sh_read_coefficients)
SH_HEADER (see above), with nrset == 1
a00 b00
a10 b10
a11 b11
a20 b20
a21 b21
....
Unformatted file.
2) Surface velocity fields
SH_HEADER, e.g: 127 0 0 1 2 0 (nrset == 2 for poloidal and toroidal fields)
a00p b00p a00t b00t
a10p b10p a10t b10t
...
(B) Short format
As above, but will only expect a single integer, lmax, as the header
for a spherical harmonic expansion.
OTHER INPUT FILE FORMATS
1) Viscosity structure (hc_assign_viscosity)
r_i e_i
Unformatted list of radii (radius of layer/Earth radius) and viscosity
(in Pas) values, reads until end of file. Values determine each layer
viscosity upward until the next entry. Use the graphical tool
pvisc.py to edit such files.
2) Depth dependent density scaling file
r_i d_i
Format as for the viscosity file, but d_i are the depth-dependent
scaling factors (this overridings -ds). Use the graphical tool
pdens.py to edit such files.
OUTPUT FILES
After a regular run, file sol.bin will have the velocity solution
(cm/yr) in binary format. This file can be extracted using
hc_extract_sh_layer (see below).
File geoid.ab will have the geoid in meters in a spherical harmonic
expansion.
USING THE OUTPUT
1) Extracting spherical harmonics solutions.
a) extract SH expansion of radial velocity at layer 2
hc_extract_sh_layer vel.sol.bin 2 1
b) extract SH expansion of radial velocity at layer 2 and convert
to spatial expansion
hc_extract_sh_layer vel.sol.bin 2 1 0 | sh_syn
c) extract SH expansion of poloidal and toroidal velocity at layer 5 and convert
to spatial expansion of v_theta v_phi
hc_extract_sh_layer vel.sol.bin 5 2 0 | sh_syn
Also see the script calc_vel_and_plot for some suggestions on
how to convert the output.
d) convert velocity output and density anomalies into a binary VTK file
for paraview
hc_extract_spatial vel.sol.bin -2 6 dscaled.sol.bin > vel.vtk
Load a view into paraview (pvsm file is included here)
paraview --state=hc_world.pvsm
to get the view in hc_world.png.
USING THE SPHERICAL HARMONICS TOOLS AS STANDALONE
1) Convert scalar values from a GMT/Netcdf grd file into a spherical
harmonics expansion of degree 31
a) obtain scalar data at the Gaussian intergration points
sh_ana -127 | grdtrack -Lx -G$datadir/etopo5/etopo5.1.grd > etopo5.dat
b) expand
cat etopo5.dat | sh_ana 127 > etopo5.ab
c) Alternative: use the grd file directly
sh_ana 127 $datadir/etopo5/etopo5.1.grd
2) convert spherical harmonics to spatial expansion
cat etopo5.ab | sh_syn > etopo5.127.dat
Note that sh_syn and sh_ana are only example implementations of the
subroutines, there's very limited actual functionality. For a more
useful spherical harmonics package, see shansyn at
http://www-udc.ig.utexas.edu/external/becker/sdata.html
That being said, also note helper programs sh_corr and sh_power.
SEATREE
HC is a module of the Solid Earth Teaching and Research Environment
(SEATREE) which provides a graphical user interface to flow
computations and plotting.
https://github.com/thwbecker/seatree
UGESCE
The Unified Geodynamics Earth Science Computing Environment (UGESCE,
http://www-udc.ig.utexas.edu/external/becker/ugesce.html) provides a VirtualBox
Linux install that includes SEATREE, HC, and a range of other Earth
Science data and software, all in one (big) package, ready to go.
COPYRIGHT
Versions of this software might include Numerical Recipes code (HC
does not rely on it, only some related tools do) - copyright is with
these authors, do not distribute without permission.
For all of the HC code, copyright by Thorsten Becker,
[email protected], under GPL of 1991.