Disclaimer and Free Software Statement

The computer codes that I post to the JParkCodes blog are provided as free software.

You can redistribute the codes and/or modify them under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,

but WITHOUT ANY WARRANTY; without even the implied warranty of

MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

GNU General Public License for more details.

You may obtain a copy of the GNU General Public License

along with these program at the link above. If not active, see http://www.gnu.org/licenses/.

grn_propzon.f -- code to compute greensfunctions for zonally-symmetric earth models.

click title to see Google Drive directory JParkCodes

c program grn_propzon -- three-component gfcts
c for zonally-symmetric earth models
c f77 -o /Users/jjpark/bin/grn_propzon grn_propzon.f /Users/jjpark/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/grn_propzon grn_propzon.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/jlib.a
c the code reads a special compressed hybrid eigenvector file
c produced by program propzon.
c This routine does a coupling calculation in frequency intervals, storing all
c hybrid singlets in a central interval of the coupling bandwidth.
c the trick is to read the coupling partners for each block, calculate
c receiver and source scalars, then read the eigenvectors one by one.

As with the spherical-earth greensfunction code, this code creates waveforms for a source-receiver geometry stored in the SAC header of a common seismic data file. Unlike the sph-earth code, the output is tied to geographic coordinates (east/north), not radial/transverse. The code syndat.f takes the greensfunctions and sums with a source mechanism to form synthetic seismograms.

sample input:

0 --> # of point to compute, 0=same number as in the data file
ffc.lhz --> data file for source-receiver-station info
gzB0.05_ffc --> filename root for greenfct files, code will append ".lhz" ".lhe" ".lhn"
../Vary/ymodel10_B0.05.zon20 --> input file of hybrid singlets, computed with zonchn.f
/Users/jjpark/Mineos/Models/premiss --> earth model (isotropic PREM in this case)
/Users/jjpark/Mineos/sphnl --> N,L file of spheroidal mode parameters
/Users/jjpark/Mineos/sphe --> direct-access Fortran file of spheroidal modal eigenfunctions
/Users/jjpark/Mineos/tornl --> N,L file of toroidal mode parameters
/Users/jjpark/Mineos/tore --> direct-access Fortran file of toroidal modal eigenfunctions
stop

zonsubs.f -- the subroutines needed to make zonchn.f work

click title to see Google Drive directory JParkCodes

This set of subroutines serves zonchn.f, see the earlier post. The first routine yread reads the parameters of the zonally symmetric model example file here

1 1 ! useless parameters for option I never implemented
0.0 0.05 0.0 0.0 0.0 A,B,C,D,E relative perts to the Backus anisotropy parameters for azim anisotropy with an axis of symmetry.
1 152 181 ! lithosphere in model premiss
10 ! k=10 wavenumber on the sphere

Some relevant passages of code are the following:

c angint is original theta integral for cos{(l+1/2)\theta} structure
c angint2 uses 3-j tricks to speed up the calculation
c angint2X uses 3-j tricks, but for either spherical cap or equatorial band
c angint2XX is s=20 expansion of equatorial band model, no discontinuites
c call angint2(l1,l2,ls(kii),angi)
call angint2X(l1,l2,ls(kii),angi)
c call angint2XX(l1,l2,ls(kii),angi)


c combine isotropic and anisotropic parts, use b,c,e & a,d
c **************************** note that code can incorporate only
c ***************************** one aspherical wavenumber, despite
c ***************************** the do loops
c NOTE: some lines have alternate versions for phi and theta as axis of sym
if(ick.eq.3) then ! tor-tor
do j=1,nom
coup(j)=coup(j)
x +angi(3,j)*(c*ysc(1,1)+e*ysc(2,1)) ! anisotropic
c x -angi(5,j)*e*ysc(3,1)-angi(6,j)*c*ysc(1,1) ! anisotropic-phi
x +angi(5,j)*e*ysc(3,1)-angi(6,j)*c*ysc(1,1) ! anisotropic-theta
x +angi(2,j)*d*ysc(3,1)+angi(3,j)*d*ysc(2,1) ! isotropic
end do

comment or uncomment the phi and theta terms as you wish, but dont comment both or uncomment both.

zonchn.f -- code to couple free oscillations for a zonally symmetric earth model

click title to see Google Drive directory JParkCodes

zonchn.f is a code to couple together free oscillations of Earth for a set of highly restricted aspherical earth models. The models are zonally symmetric, dependent only on co-latitude. The code has subroutines that integrate the modes against two types of models (1) a model that depends on a single wavenumber k, and (2) a model that is an equatorial band or a symmetric spherical cap. The first model type was used to examine the wavenumber behavior of surface-wave coupling, and the second type simulates the interaction of surface waves with a sharp boundary. The paper that describes the code's workings is

Park, J., The sensitivity of seismic free oscillations to upper mantle anisotropy I. Zonal symmetry, J. Geophys. Res., v98, 19933-19949, 1993.

Zonal symmetry in the model allows the coupling to vanish for all free-oscillation singlets with differing azimuthal numbers m, m'. In this manner, if there are two modes nSl and n'Sl' with a total of 2l'+1+2l+1 = 2(l+l'+1) singlets, the coupling problem splits apart into (for l'>l) 2l+1 2x2 coupling matrices, and l'-l uncoupled singlets. As a result of this symmetry, zonchn can readily accept a large number of modes in a narrow frequency band without breaking the computer memory or the CPU. The general scheme is to march up in frequency in steps of df and Df: All modes in a half-bandwidth Df about some fiducial frequency fo are coupled, but only those hybrid modes in a narrower bandwidth df about fo are saved. The code increases fo by 2*df and iterates up to some cutoff frequency. The user specifies how many overtones n are included. The uploaded version of zonchn.f has been tested with 0.2.le.f.le.20 mHz and df=0.1 mHz, Df=0.5 mHz and n-max=5.

The code depends on a set of subroutines in the file zonsubs.f, which are posted next in this blog.

The code reads from an input file "phzonchn" that contains

print *,'open(15,file=phzonchn,form=formatted)'
open(15,file='phzonchn',form='formatted')
call yread
print *,'output chain file '
read(15,102)namchn
print *,namchn
102 format(80a)
print *,'initial open of output file in append mode (0=no) '
read(15,*) iappend
print *,'iappend= ',iappend
print *,'enter inner freq halfwidth, outer freq halfwidth (mHz)'
c outer band is modes that couple, inner band contains modes that are saved
read(15,*) hfin,hfout
print *,hfin,hfout
print *,'enter freq band for hybrid modes (mHz)'
read(15,*) fb1,fb2
print *,fb1,fb2
fb1=dmax1(fb1,0.25d0)
print *,'enter max overtone branch to consider'
read(15,*) nbrmax

The option of writing the output file in "append" mode is an artifact of slow computers of the 1990s, in which one might need to stop and restart a long computation. I doubt that the user will need this option in 2009, though my test calculation required the greater part of an hour to complete.

example for phzonchn:

ymodel10_B0.05 --> aspherical model file
ymodel10_B0.05.zon20 --> output file of hybrid modes (linear comb of sph-earth modes, new frequencies)
0 --> dont open output file in append mode (usual case)
0.1 0.5 --> inner band contains modes that are saved, outer band is modes that couple (mHz)
0.2 20 --> frequency band to compute (mHz)
5 --> max number of overtones
/Users/jjpark/Mineos/Models/premiss --> earth model (isotropic PREM in this case)
/Users/jjpark/Mineos/sphnl --> N,L file of spheroidal mode parameters
/Users/jjpark/Mineos/sphe --> direct-access Fortran file of spheroidal modal eigenfunctions
/Users/jjpark/Mineos/tornl --> N,L file of toroidal mode parameters
/Users/jjpark/Mineos/tore --> direct-access Fortran file of toroidal modal eigenfunctions


The model file is based on the A,B,C,D,E Backus parameterization of hexagonally-symmetric anisotropy. A and D are isotropic P and S wavespeed perturbations (peak2peak variation: A=0.05 implies 5%). B is the 2\phi azimuthal Vp term, C is the 4\phi azimuthal Vp term, and E is the 2\phi azimuthal Vs term. The input, as explained in the post for zonsubs.f iare the values of the ABCDE and the wavenumber k of a corrugation in colatitude. The code does not allow the user to choose the equatorial band model on input. You must change a subroutine call and recompile the code. I may fix this soon, as it is cumbersome. Similarly, there is a slight alteration in the coupling formulas that switches the cases for an anisotropic symmetry axis w-hat aligned with theta-hat, or with phi-hat, that is, with colatitude and longitude, respectively. Due to symmetry restrictions, only w-hat parallel to theta-hat or phi-hat is allowed by this code. No tilted axes of symmetry, no general azimuthal angle.

syndat.f -- sum greens functions with moment tensor and centroid time correction

click title to see Google Drive directory JParkCodes

NOTE:  This code was updated on 7/24/2012.  The title of the post points to the newer code.  The older code can be found at Google Drive JParkCodes/syndat_2009.f

This code is a version of the old sac_syndat.f that incorporates the IRIS-suppied EVRESP subroutine to compute and apply instrument response correction to convert displacement to instrument counts.

c f77 -o /Users/jjpark/bin/syndat syndat.f /Users/jjpark/Ritz/libevresp.a /Users/jjpark/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/syndat syndat.f /Users/jjpark/Ritz/libevresp.a /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/jlib.a
c program to make synthetics using green fcts
c the greenfcts are raw displacement in meters. This code applies the instrument
c response in the frequency domain with a FFT --> multiply by H(f) --> IFFT process.

program asks for the centroid time delay, which is applied with a phase-ramp in the frequency domain
The centroid time delay occurs because the start time of a large event is typically not the best time to place the stepfunction onset of the long-period course. If moment release occurs over T seconds, then a better fit to data is found by delaying the synthetic onset by T/2 seconds. (Assuming a symmetric source time function, natch).

print *,'enter time delay e.g. centroid-time correction'
read(5,*) tconst
40 print *,'do you want to use strike,dip and rake (0=no) '
read(5,*) ick
if(ick.ne.0) then
call fault(sol)
else
print *,'enter f(x6) : '
read(5,*) (sol(i),i=1,6)
endif
print *,'routine expects the gfct filename to start with g'
print *,'the code substitutes s for g in the output filename'
30 print *,'green fns file name(=stop) : '

The program will read a succession of greensfct files, one per station-component, until it reads 'stop' e.g.

0 --> zero centroid time delay
1 --> read strike,dip,slip
329 8 110 --> strike dip slip
400 --> seismic moment (units of 10**27 dyne-cm/10**20 nt-m)
/Users/jjpark/Synth/gsph_ffc.lhz --> greensfct file
/Users/jjpark/Synth/gsph_ffc.lhr
/Users/jjpark/Synth/gsph_ffc.lht
stop

There are PLOTIT calls within the code (commented out for now) that allow the user to view his/her handiwork, testing the causality of the instrument response, etc.

The EVRESP routine is found in the libevresp.a library that you can create from the IRIS EVALRESP software distribution. For reasons only the DMC programmers know, the latest versions of the EVALRESP distribution hide the libevresp.a library in a subdirectory ".lib" that wont show up in a "ls" command unless you are looking for it. Just another convenience feature, I guess.

Evalresp follies

The IRIS-supplied instrument response software was not straightforward to incorporate into our codes. The key C-subroutine EVRESP is built for incorporating into Fortran codes, but seems to have developed a bug over its past revisions, and no one may have complained of it. The subroutine generates the filename of a RESP file, nominally output by rdseed. In previous versions of the subroutine I have used (circa 2005) the routine opens the file, reads the contents and generates a frequency domain instrument response for convolving synthetics or deconvolving data. The 2008-9 versions of EVRESP cant seem to find the file whose name they create (seemingly correctly) and the program fails. Diddling with the SEEDRESP environmental variable does not seem to help. The fix that I found was to feed EVRESP the RESP.* filename explicitly, by generating the filename myself from the same group of character-string variables that EVRESP uses. To wit:

file='RESP.'//net//'.'//sta//'.'//locid//'.'//cha
print *,file
c remove blank spaces
do i=30,7,-1
if(file(i:i).eq.' ') file=file(1:i-1)//file(i+1:80)//' '
end do
print *,file
print *,'entering evresp'
print *,sta,cha,net,locid,datime,unts,nscan,npad,nha
iflag=evresp(sta,cha,net,locid,datime,unts,file,freq,nha,resp,
x rtyp,vbs,start_stage,stop_stage,stdio_flag)
print *,'done evresp'

The character string "file" is passed to the routine and the code looks for it successfully, if the file is there!

sacread.f -- subroutines to read and write SAC-format files -- adjusted to write 6-component greens functions

sacread.f

subroutines to read and write SAC-format files -- adjusted to write 6-component greens functions. The routines provided by SAC are fine, of course, but I use a hydrid file format for greens-function files. These files contain six displacement traces rather than one, one for each independent component of the seismic moment tensor. There is a flag in the SAC header that I took over for indicating that six traces need to be read instead of one.

if(chead(145).ne.' GRN'.and.chead(145).ne.' grn') then
nword=nscan+158
nrec=(nword-1)/128+1
last=nword-128*(nrec-1)
else
print *,'reading greens functions'
nword=6*nscan+158
nrec=(nword-1)/128+1
last=nword-128*(nrec-1)
endif

The subroutines are sacread, sacout, sacread_swap and sacout_swap. The last two routines swap big-endian and little-endian data. You will be happier if you never need to swap bytes, but they are there for when you are painted into that corner.

csolve.f -- solve a system of equations with complex coefficients

click title to see Google Drive directory JParkCodes

csolve.f is a subroutine for doing Gaussian elimination with partial pivoting to solve a system of linear equations A.x=y with complex matrix A and vector y. Useful in several routines.

Snow Leopard transition: Fortran

Just an informational item regarding Snow Leopard's treatment of Fortran. The first comment is that programmers need to install XCode from the Snow Leapard install DVD. I bulled ahead without checking this and ran into a roadblock in trying to compile a Fortran program -- the computer could not find the unix command "as" which assembles object files. Once I loaded the XCodes package from the Install disk, this problem went away, and others cropped up.

The innovation of Snow Leopard (and Leopard, I think, Im not exactly sure when the transition occurred relative to Tiger (MacOS 10.4)) is the full utilization of 64-bit chip architecture. The Snow Leopard gcc compiler creates object files that are 64bit, but the Fortran compiler g77 does not. If you try to assemble Fortran and C routines into one archive or a single program, the computer balks. You can tell the nature of the object file with the command "file *.o". I sampled the internets for hints on solving this problem. The simplest solution, and the one Im adopting for now, is to force the gcc compiler to generate a 32-bit object file, rather than a 64bit version, with the "-m32" runtime option. Less powerful than 64bit, and Im still getting warnings from the compiler when I create large codes, but the executables seem to work OK.

The "gfortran" option, rather than "g77" for the Fortran compiler, might work with the 64bit option, but the internets are not unanimous on its success. There are compilers available for download from the HPC site (http://hpc.sourceforge.net/) that claim to vectorize gcc (version 4.4) in Snow Leopard. I read comments from programmers that full-optimization of the compiler has bugs, so buyer beware. I might try it later, though.

NOTE: In 2012-2013 I switched to compile codes with gfortran compiler.

grn_sph.f for computing greens function for synthetic seismograms using summed free oscillations for spherical earth models

click title to see Google Drive directory JParkCodes

grn_sph.f

The task of computing synthetic seismograms from summed free oscillations for an aspherical earth model is typically accomplished in two or three steps. The first step computes the hybrid free-oscillations for a chosen earth model, comprised of the "simple" spherical-earth modes whose particle motions follow single vector spherical harmonics. The second step takes those hybrid modes and computes greens functions for particle motion at a given source location (allowing six different time series, one for each component of the seismic moment tensor) and station location (for a given sensor and three time series for three components: vertical/horizontal). The last step is the computation of a synthetic seismogram by summing the six time series for each motion component by weighting with the components of the moment tensor.

Historically, I have used the greens function routine as the location for adding the instrument response, but I have reconsidered this choice. Computers are fast enough for the FFTs an instrument response to be computed on the fly in the last step. This also avoids the problem of updating the instrument response codes (notoriously unstable, owing to the various formats for instrument responses) in multiple greensfunction codes. My final code (sac_syndat.f) is used for all, and should be the place where instrument response either is (or is not) applied. The output of this routine has units of displacement in meters.

So, for making these codes more public, lets start with the simple code for synthetic seismograms on a spherical earth model.
The code is grn_sph.f, which has existed in one form or another for over 20 years (this explains the dead code inside it).

print *,'typical input for this program:'
print *,'0 20 --> freq interval to sum (mHz)'
print *,'/home/jjpark/Mineos/Models/premiss --> 1-D earth model'
print *,'/home/jjpark/Mineos/premissnl --> spheroidal info'
print *,'/home/jjpark/Mineos/premiss_eig --> spheroidal modes'
print *,'/home/jjpark/Mineos/premistnl --> toroidal info'
print *,'/home/jjpark/Mineos/premist_eig --> toroidal modes'
print *,'1000 --> # of overtone branches (here: all of them)'
print *,'0 --> # of data points to synthesize (0-> #=data: NPTS)'
print *,'/park/nobackup/Sumatra/VH/PET.VHZ --> data file'
print *,'.....'
print *,'. . . more data files'
print *,'.....'
print *,'stop'
print *,' '
print *,'The format of the greens fct files is impure SAC-format'
print *,'readable only by sac_syndat and related JPark codes'
print *,'It has a SAC header with one ASCII variable set to GRN'
print *,'and six*NPTS data values that follow the header'
print *,'NPTS for each component of the source moment tensor'
print *,'The greens-function file names are created from your input file'
print *,'by prepending gsph_ to the input filename'
print *,'and replacing the last character with z,r,t'

This code reads files full of free-oscillation displacement functions for a standard earth model, computes the strain of the mode at the desired source depth, and sums 6 traces for surface displacement at a specified distance from the source. Because the earth model is a spherical-earth model, only the source depth and source-receiver distance affect the waveforms that represent oscillatory displacement. The code, however, takes a data file to specify the source and receiver locations and creates synthetics for the geometry. A subsequent processing step (syndat.f) will take the greens function file, sum with the components of a moment tensor, and apply the instrument responses appropriate to the station's sensors and channel.

The start time for the earthquake determines where the decaying sinusoids begin. Older versions of green-fct codes used an explicit time that the user had to supply. This version assumes that the information is in the omarker time in the SAC header. There is also a bmarker time that is used in case the start time of the record does not correspond to the reference time in the data header.
c the reference time of the file could be before or after the first data point
c the ymdhmsec parameters define the reference time
c the origin quake of the quake and start of record are wrt reference time
origin=ahead(8) ! event start is the omarker in sac header
recstart=ahead(6) ! record start is the bmarker in sac header