recfunk.f -- original recfunk version without moveout corrections

recfunk.f -- original recfunk version without moveout corrections

Click the title of the post to access the code recfunk.f on Google Drive directory JParkCodes.

c  program recfunk

c  5/25/99 JJP updated for GFORTRAN compiler 07/27/15

c

c  original code written for the paper 

c  Park, J., and V. Levin, Receiver functions from multiple-taper 

c  spectral correlation estimates, Bull. Seismo. Soc Amer., v90, 

c  pp1507-1520, 2000.

c

c  this version of the RF code reads datafile names from standard input

c  then the time intervals in the filename file

c  the data must be binary SAC format

c  horizontals must be rotated to radial and transverse

c

c  for start times in the file:

c

c  1997.070.19.33.bh?  <-- code="" div="" r="" replaces="" t="" with="" z="">

c  57 52       <-- analysis="" div="" duration="" of="" sec="" start="" time="" window="">

c  1997.076.08.15.bh?

c  62 62 

c  ...

c  ...

c  ...

c  stop                <-- 399="" code="" data="" div="" events="" finished="" is="" max="" tells="" that="">

c

c  has kluge to cheat the pre-event noise for synthetic records  3/12/00 JJP

c  check to see if the kluge is commented out

c

c  multitaper program to invert for receiver functions

c  SINGLE STATION WITH SPECTRAL WEIGHTING

c  

c  uses spectral estimates of pre-event noise 

c  to weight the inversion in the frequency domain, thereby avoiding the 

c  usual problem with microseism noise in 0.1-0.5 Hz

c

c  xf77 -o /park/backup/bin/recfunk recfunk.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a

c xf77 -o /Users/jjpark/bin/recfunk recfunk.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a

rfmigrate.f -- recfunk version with migration

rfmigrate.f -- recfunk version with migration

Click the title of the post to access the code rfmigrate.f on Google Drive directory JParkCodes.

c  program rfmigrate

c  9/20/00 JJP -- adapted from recfunk  --- UPDATED 08/19/07 JJP

c   updated again 07/22/15 JJP

c  migrates MTC receiver functions in the frequency domain.

c  requires a stacking model in the anirec format

c  such a model may have anisotropy parameters, 

c  but migration code only uses the isotropic velocities. 

c

c  has kluge to cheat the pre-event noise for synthetic records  3/12/00 JJP

c  check to see if the kluge is commented out, because anirec_synth was updated 

c  in 2015 to have a proper-length pre-event time window.

c

c  this version of the RF code reads a file of data filenames 

c  you have two choices: either read the time intervals in the filename file

c  or read them in the sac header

c  the data must be binary SAC format

c  horizontals must be rotated to radial and transverse

c

c  for start times in the file:

c  the file defaults to "in_recfunk" and has lines of the form:

c

c  1997.070.19.33.bh?  <-- code="" div="" r="" replaces="" t="" with="" z="">

c  57 52       <-- analysis="" div="" duration="" of="" sec="" start="" time="" window="">

c  1997.076.08.15.bh?

c  62 62 

c  ...

c  ...

c  ...

c  stop                <-- 799="" code="" data="" div="" events="" finished="" is="" max="" tells="" that="">

c

c

c  for start times in the SAC header, the default file is "in_recpick"

c  reads seismic record start times from the sac header

c  will search the SAC header for specific markers of P phases

c  T1 - P, Pdiff    ahead(12)

c  T2 - PKP,PKIKP   ahead(13)

c  T3 - PP          ahead(14)

c  T1=T2=T3=0 ==> use original A-marker ahead(9)

c

c  legacy data has timing in the

c   A-marker (ahead(9)) which doesnt indicate which phase was marked, 

c  For the A-marker, the phase is identified as P/Pdiff for DELTA<120 div="" nbsp="">

c  PKP/PKIKP for DELTA>120

c

c  code does NOT combine data with different sample rates

c  it is possible to spline interpolate the data files for a uniform sample rate

c  see the code in rfmig_mcboot.f, which should be spliceable into this code.

c  data files are limited to 99K pnts. To increase, see common block /datastuff/

c

c   many intermediate quantities are plotted with PLOTIT as the code proceeds.

c  other intermediate quantities can be plotted by uncommenting calls to PLOTIT

c

c  the code writes the BAZ- and EPICEN-dependent RFs to files

c  in a format easily digested by GMT (traces are separated by ">" lines)

c

c  filenames: out[rt]_baz.grid oum[rt]_epi.grid oum[rt]_baz.grid

c

c  out --> no moveout correction;  oum --> moveout correction

c

c  these files are overwritten everytime you run the program

c  so rename them to save them

c

c  xf77 -o /park/backup/bin/rfmigrate rfmigrate.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a

c  xf77 -o /Users/jjpark/bin/rfmigrate rfmigrate.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a

c

anirec_synth_circle.f -- make P-coda synthetics at evenly-spaced back azimuth and constant epicentral distance

anirec_synth_circle.f -- make P-coda synthetics at evenly-spaced back azimuth and constant epicentral distance

click title to see Google Drive directory JParkCodes

c  anirec_synth_circle - anirec to compute synthetics for a specified model and
c  a specified station location from one seismic record read from in_recfunk
c  the code will take the station location and generate a ring of 120 events
c  that are 90 degrees away, one for every 3 degrees back-azimuth.
c  To do this we generate a great circle path
c  with the station at the pole of the great circle
c
c  To minimize aliasing in the time domain, we create synthetics that last 200 seconds
c  The inverse FFT of the frequency-domain reflectivity response is time shifted to
c  place the P wave just beyond the 100-sec mark, using a wraparound of the inverse FFT.
c  The code can use several pulse shapes, according to the flag ipulse
c   ipulse=1 -- one-sided pulse (squared sine wave) with user-chosen period
c   ipulse=2 -- two-sided sine pulse with user-chosen period
c   ipulse=3 -- two-sided dispersive pulse, 160-points
c   ipulse=4 -- two-sided 10-second dispersive pulse with random-walk phase
c   ipulse=5 -- random-walk pulse (integration of RAND) -- emergent 10 sec pulse with tapered end
c
c   This code is hardwired for incident P waves (ity=1), but the reflectivity calculation has info
c   to construct incident SV (ity=2) and SH (ity=3) synthetics as well.  The exact timing of these
c   S-coda are different from the P-coda, so you would need to check the wraparound timing to get it right.
c
c  code accepts as input SACfiles with start time info in the header
c  legacy data has timing in the
c   A-marker (ahead(9)) which doesnt indicate which phase was marked,
c  For the A-marker, the phase is identified as P/Pdiff for DELTA<120 nbsp="" p="">c  PKP/PKIKP for DELTA>120
c  the value of A-marker is not important.
c  In the output files of synthetics, the A-marker and/or the T1,T2,T3-marker
c  in the header are set to 95.0 seconds to allow for some pre-event noise
c  that is generated by intrinsic function RAND
c  for A.ne.0 by recfunk codes that use the header information for timing.
c  code divides between T1 (P) or T2 (PKP) or T3 (PP) markers based on epicentral distance,
c  subroutine call_ttimes reads files that give slowness values for 0c  Preactically, these wont be used for all DELTA, but consistency!
c  the synthetic files are s_filename.bh[zrt], where filename is the name of the datafile.
c  the code writes in_recpick_synth, a list of the filenames for recfunk_pick
c  and other recfunk codes that expect the start time in the SAC header
c
c  written originally 12/3/03 JJP
c  3/3/04 -- Added Kluge -- I have added code that forces the phase velocity to
c  fluctuate with BZ to mimic the effect of a dipping layer,
c  changing its incidence angle with a cos(\phi) dependence.
c  This kluge mimics the P-SV conversion variation caused by a dipping interface
c  but not the P-SH conversion behavior.
c  This kluge was used for the last check on RF interpretation in
c  Park et al 2004 in JGR, the paper on the Cascadia subduction zone.
c  Look for the word "kluge" in the code and make certain that the code
c  between begin-kluge and end-kluge comment-lines is commented
c
c  xf77 -o /Users/jjpark/bin/anirec_synth_circle anirec_synth_circle.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c
c  anirec - program to calculate receiver-function response of
c  a stack of anisotropic layers to a plane wave incident from below
c  cannibalized from aniprop.f  11/18/95
c  for hexagonally symmetric media
c  reads fast axis orientation, constants A,B,C,D,E from file animodel
c  calculate quadratic eigenvalue problem based on the Christoffel matrix
c  see appendix of P. Shearer's thesis
c
c  read model, phase velocity of incident wave, P, SV, or SH
c
c  calc the eigenvector decomps for the layers
c  loop over frequency, calc reflection/transmission matrices
c  calc 3-comp transfer fct response at surface
c  find distortion of reference wavelet
c

anirec_synth.f --- a code to compute synthetic P coda from a collection of SAC-format data

anirec_synth.f --- a code to compute synthetic P coda from a collection of SAC-format data

click title to see Google Drive directory JParkCodes

c  anirec_synth - anirec to compute synthetics for a specified model and
c  a specified collection of seismic records
c  (read from in_recpick-format file of filenames)
c the goal is to create P-coda synthetics that match the event distribution of your dataset.
c
c  To minimize aliasing in the time domain, we create synthetics that last 200 seconds
c  The inverse FFT of the frequency-domain reflectivity response is time shifted to
c  place the P wave just beyond the 100-sec mark, using a wraparound of the inverse FFT.
c  The code can use several pulse shapes, according to the flag ipulse
c   ipulse=1 -- one-sided pulse (squared sine wave) with user-chosen period
c   ipulse=2 -- two-sided sine pulse with user-chosen period
c   ipulse=3 -- two-sided dispersive pulse, 160-points
c   ipulse=4 -- two-sided 10-second dispersive pulse with random-walk phase
c   ipulse=5 -- random-walk pulse (integration of RAND) -- emergent 10 sec pulse with tapered end
c
c   This code is hardwired for incident P waves (ity=1), but the reflectivity calculation has info
c   to construct incident SV (ity=2) and SH (ity=3) synthetics as well.  The exact timing of these
c   S-coda are different from the P-coda, so you would need to check the wraparound timing to get it right.
c
c  code accepts as input SACfiles with start time info in the header
c  legacy data has timing in the
c   A-marker (ahead(9)) which doesnt indicate which phase was marked,
c  For the A-marker, the phase is identified as P/Pdiff for DELTA<120 font="" nbsp="">
c  PKP/PKIKP for DELTA>120
c  the value of A-marker is not important.
c  In the output files of synthetics, the A-marker and/or the T1,T2,T3-marker
c  in the header are set to 95.0 seconds to allow for some pre-event noise
c  that is generated by intrinsic function RAND
c  for A.ne.0 by recfunk codes that use the header information for timing.
c  code divides between T1 (P) or T2 (PKP) or T3 (PP) markers based on epicentral distance,
c  subroutine call_ttimes reads files that give slowness values for 0
c  Preactically, these wont be used for all DELTA, but consistency!
c  the synthetic files are s_filename.bh[zrt], where filename is the name of the datafile.
c  the code writes in_recpick_synth, a list of the filenames for recfunk_pick
c  and other recfunk codes that expect the start time in the SAC header
c
c  modified 11/30/07 JJP
c  5/12/01 JJP, 7/23/2015
c  xf77 -o /Users/jjpark/bin/anirec_synth anirec_synth.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c
c  anirec - program to calculate receiver-function response of
c  a stack of anisotropic layers to a plane wave incident from below
c  cannibalized from aniprop.f  11/18/95
c  for hexagonally symmetric media
c  reads fast axis orientation, constants A,B,C,D,E from file animodel
c  calculate quadratic eigenvalue problem based on the Christoffel matrix
c  see appendix of P. Shearer's thesis
cawk '{print "s_" $1}' in_recpick >! in_recpick_synth

c  read model, phase velocity of incident wave, P, SV, or SH
c
c  calc the eigenvector decomps for the layers
c  loop over frequency, calc reflection/transmission matrices
c  calc 3-comp transfer fct response at surface
c  find distortion of reference wavelet
c

Updated Fortran Programs for Receiver-Function Processing

Sportsfans,

After upgrading my Macs to MacOS X 10.10.4 Yosemite, I loaded the gcc5 compiler from Fink.org and started modifying my codes to compile properly with the gfortran compiler that comes with Rev5 of gcc.  Although the codes are still written mostly in the oldest flavors of Fortran, they will indeed compile and run in the new computing environment.  I will post descriptions of the updated codes over the next days, and load them into my Google Drive directory.

One update is that the RF codes that implement a moveout correction while stacking have been altered to use P-wave slowness files for the P, PKP-PKIKP and PP cases.  These files are newly computed for the iaspei91 model from the TauP codes that superseded the venerable ttimes program in the IRIS software library.  You will see that the subroutine read_ttimes opens three files in a subdirectory of an Applications folder where TauP sits.  I will upload my own versions of these slowness files to Google Drive. These files average over triplications and allow Pdiff to reach DELTA=180 and follow the PKiKP phase to its zero-slowness value at DELTA=0.

look for

taup_P.dat
taup_PP.dat
taup_PKP.dat

The Geology Department Webserver was taken down, so all of the links in this code blog have been lost.  I have the material saved, and will be putting the codes back into blogspot itself.

DONE.  The codes have been placed in the cloud on Google Drive.  A link in most posts takes you to the Google Drive directory JParkCodes.

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

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

click title to see Google Drive directory JParkCodes

This code is an updated version of a code first posted on 10/30/2009.  The differences are mainly cosmetic, but have practical value.  The revision allows the code to circumvent an annoying convention in horizontal components and spherical-earth mode-sum synthetic seismograms.  Because the code grn_sph.f writes greensfunctions for radial and transverse horizontal components of motion, it needs to apply the instrument response for one of the actual components of the station.  The code was originally written to use the E-component response file.  A snag arises when there is no east component for the station because the horizontals are deflected too far from the compass orientations.  If the RESP-file for the E-component does not exist, syndat looks for the response file of the "1" component instead.  (A station may collect data on the LH1,LH2 components, instead of the LHE,LHN components.)   This workaround seems to work in test examples.

ALSO, the code is using the software library from evalresp3.3.3, obtained from IRIS.  The C-routine evresp.c in that package spits out a comment that "velocity (m/s)' is being converted to instrument counts.  This comment is not correct -- the syndat code tells the evresp() function that its input is displacement, and the application of instrument response looks right.  (I also tested the false specification of "VEL" input, but that choice misfit seismic data greatly.

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.