HC is a global mantle circulation solver following Hager & O'Connell (1981) which can compute velocities, tractions, and geoid for simple density distributions and plate velocities.
Switch branches/tags
Nothing to show
Clone or download
Permalink
Type Name Latest commit message Commit time
Failed to load latest commit information.
anisotropic_viscosity Added basic VTK output functionality and hc_extract_spatial Jul 16, 2011
example_data Added Gauss lat program, moved some files around. Sep 14, 2010
hcplates Internal modifications to structure naming definitions and such. Feb 17, 2008
prem Minor radius compatibility fix May 4, 2010
test Modified treatment of L=1 solution, added testing directory. Nov 27, 2008
.dat No commit message May 5, 2006
COPYING added license file Apr 8, 2009
Makefile Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
Makefile.craig Made GMT 4 the default choice, now USE_GMT3 needs to be defined for o… Apr 8, 2009
Makefile.include Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
README.TXT Change geoid default to EGM 2008 Aug 8, 2017
calc_dyn_topo No commit message Oct 14, 2013
calc_rms No commit message May 5, 2006
calc_vel_and_plot Typo fixed. Feb 21, 2010
change.log No commit message May 5, 2006
cig_bench.tgz Readd benchmark file Oct 14, 2013
compare_solution_methods No commit message May 5, 2006
compile_simple_test No commit message May 5, 2006
compute_and_extract_tractions Fixed minor typos, added another example script Apr 29, 2008
convert_bernhard_dens.c Updated readme. Jul 19, 2011
cpol.gpl No commit message May 5, 2006
drive_visc_scan Some new interpolation features Aug 14, 2017
dtopo_ref.ab No commit message Jun 21, 2013
egm2008-hc-geoid.chambat.31.ab Change geoid default to EGM 2008 Aug 8, 2017
egm2008-hc-geoid.nakiboglu.ab Change geoid default to EGM 2008 Aug 8, 2017
expansion.par - Fixed a few minor bugs. Aug 27, 2007
gaussp.c Added capabilities for hc_visc_scan, attempts for kernel output (not … Aug 3, 2017
geoid_ref.ab Change geoid default to EGM 2008 Aug 8, 2017
ggrd_base.h Removed an unfinished function. Sep 3, 2017
ggrd_grdtrack_util.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
ggrd_grdtrack_util.h Made radius of planet an argument rather than fixed. Oct 13, 2018
ggrd_readgrds.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
ggrd_struc.h Removed an unfinished function. Sep 3, 2017
ggrd_test.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
ggrd_velinterpol.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
grdinttester.c Added support for nearneighbor "interpolation" Dec 29, 2010
hc.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
hc.h Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
hc2gmt Added Gauss lat program, moved some files around. Sep 14, 2010
hc_auto_proto.h Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
hc_constants.h Removed an unfinished function. Sep 3, 2017
hc_extract_sh_layer.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
hc_extract_spatial.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
hc_filenames.h Fixed boolean init problem Aug 7, 2017
hc_ggrd.h Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
hc_ggrd_auto_proto.h Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
hc_init.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
hc_init.old No commit message May 5, 2006
hc_input.c No commit message Jun 21, 2013
hc_invert_dtopo.c Added capabilities for hc_visc_scan, attempts for kernel output (not … Aug 3, 2017
hc_matrix.c Some more format fixes that got away. Jul 16, 2012
hc_misc.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
hc_output.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
hc_polsol.c Added capabilities for hc_visc_scan, attempts for kernel output (not … Aug 3, 2017
hc_propagator.c Improved geographic region precisions for sh_syn Aug 5, 2012
hc_solve.c Added capabilities for hc_visc_scan, attempts for kernel output (not … Aug 3, 2017
hc_torsol.c Updated contact info. Jul 15, 2017
hc_visc_scan.c Fixed boolean init problem Aug 7, 2017
hc_world.png Updated readme. Jul 19, 2011
hc_world.pvsm Updated readme. Jul 19, 2011
healpix.ab No commit message May 5, 2006
itg-hc-geoid.ab No commit message Jun 21, 2013
make_binary_package Add visc.dat to binary package Apr 15, 2009
make_tar Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
manipulate_hc.py Modified the way that the surface boundary conditions are handled on … Jan 2, 2008
orig.ab No commit message May 5, 2006
pdens.py Modified the way that the surface boundary conditions are handled on … Jan 2, 2008
plot_geoid Added tentative computation of geoid at all depths. Nov 18, 2008
plot_spatial_model No commit message May 5, 2006
plot_spectral_model Added Gauss lat program, moved some files around. Sep 14, 2010
plot_visc_scan updated README Aug 7, 2017
pow.sh_ana Minor updates in the hc part, some modifications to the hcplates comp… Jul 12, 2007
pow.shana Minor updates in the hc part, some modifications to the hcplates comp… Jul 12, 2007
prem.h Removed an Sep 3, 2017
prem2dsm.c Got rid of one gcc compiler warning. Oct 23, 2009
prem2r.c Can't remmber. Oct 13, 2018
prem2rp.c Removed an Sep 3, 2017
prem_util.c Removed an Sep 3, 2017
print_gauss_lat.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
pvisc.py Modified the way that the surface boundary conditions are handled on … Jan 2, 2008
rick.ab No commit message May 5, 2006
rick.v.ab No commit message May 5, 2006
rick_fft.f90 No commit message May 5, 2006
rick_fft_c.c Allow for quadruple precision. Jul 16, 2012
rick_sh_c.c Fixed boolean init problem Aug 7, 2017
rms.dat No commit message May 5, 2006
rotvec2vel.c Minute updated for supp program. Jun 20, 2012
run_bench No commit message May 5, 2006
sh.h Allow for quadruple precision. Jul 16, 2012
sh.reexpanded.pt No commit message May 5, 2006
sh.reexpanded.r No commit message May 5, 2006
sh_ana.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
sh_corr.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
sh_exp.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
sh_exp_ggrd.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
sh_extract_layer.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
sh_model.c No commit message Jun 21, 2013
sh_power.c Updated contact info. Jul 15, 2017
sh_rick.h Allow for quadruple precision. Jul 16, 2012
sh_rick_ftrn.h No commit message May 5, 2006
sh_shana.h Minor updates in the hc part, some modifications to the hcplates comp… Jul 12, 2007
sh_spherepack.h Updated readme. Jul 19, 2011
sh_syn.c Removed GMT/netcdf compile to allow obtaining hc without those libaries. Nov 7, 2018
sh_test.c Updated contact info. Jul 15, 2017
shana.ab No commit message May 5, 2006
shana.v.ab No commit message May 5, 2006
shana_sh.c Minor updates in the hc part, some modifications to the hcplates comp… Jul 12, 2007
simple_test.c Allow for quadruple precision. Jul 16, 2012
spherepack_sh.c Updated readme. Jul 19, 2011
state1.pvsm Added binary release. Sep 22, 2011
state2.pvsm Added binary release. Sep 22, 2011
test.ab No commit message May 5, 2006
test_codes No commit message May 5, 2006
test_fft.c No commit message May 5, 2006
test_sh_routines Minor updates in the hc part, some modifications to the hcplates comp… Jul 12, 2007
vdepth.dat Added flag to not scale input density anomalies with the value of PRE… Jun 27, 2008
vel.cpt No commit message May 5, 2006
visc.dat No commit message Jun 21, 2013

README.TXT

README FOR HC

Version 1.0.7 - 2017

Thorsten Becker - twb@ig.utexas.edu


HC is a global mantle circulation solver following Hager & O'Connell
(1981) which can compute velocities, tractions, and geoid for simple
density distributions and plate velocities.

This particular implementation illustrates one possible way to combine
the HC solver routines.

Based on code by Brad Hager, Richard O'Connell, and Bernhard
Steinberger, this version is by Thorsten Becker and Craig O'Neill.

Note also that the Solid Earth Teaching and Research Environment
(SEATREE, http://geosys.usc.edu/projects/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 user interface for HC.


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, due to Craig O'Neill isn't
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.6/

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.



HC CODE usage

Described in the help page that is displayed for "hc -h" as below.
Also see SEATREE (http://geosys.usc.edu/projects/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://geosys.usc.edu/projects/seatree/wiki/ 
     
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.