Tuesday, July 24, 2012


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.




Tuesday, June 14, 2011

GMT script to generate Figure 4 of Park and Royer (2011)

#!/bin/csh


gmtset ANNOT_FONT_SIZE 12p LABEL_FONT_SIZE 15p HEADER_FONT_SIZE 15p

gmtset ANNOT_FONT Times-Roman LABEL_FONT Times-Roman HEADER_FONT Times-Roman


awk '{print $1, $2, $3/1000}' out_ppdf_control.dat >! blob.dat

surface blob.dat -R0.2/10.6/0.025/0.135 -I0.4/0.01 -Gblobg -T0.25

grdimage blobg -JX3/2 -Cppdf.cpt -R -Ba2f1:'CO@-2@- Sensitivity (degC)':/0.02:'ACT'::'.Fraction of ModelRuns within 3-@~s@~':/SWen -X1.0 -Y1.0 -K >! aaplot


awk '{print $1, $5, $6/1000}' out_ppdf_control.dat >! blob.dat

surface blob.dat -R0.2/10.6/0.15/0.85 -I0.4/0.05 -Gblobg -T0.25

grdimage blobg -JX3/2 -Cppdf.cpt -R -Ba2f1:'CO@-2@- Sensitivity (degC)':/0.1:'FERT'::'.Fraction of ModelRuns within 3-@~s@~':/SWen -X0.0 -Y3.5 -K -O >> aaplot


awk '{print $1, $8, $9/1000}' out_ppdf_control.dat >! blob.dat

surface blob.dat -R0.2/10.6/0.05/0.55 -I0.4/0.025 -Gblobg -T0.25

grdimage blobg -JX3/2 -Cppdf.cpt -R -Ba2f1:'CO@-2@- Sensitivity (degC)':/0.1:'LIFE'::'.Fraction of ModelRuns within 3-@~s@~':/SWen -X4.0 -Y0.0 -K -O >> aaplot


awk '{print $1, $11, $12/1000}' out_ppdf_control.dat >! blob.dat

surface blob.dat -R0.2/10.6/0.45/1.25 -I0.4/0.05 -Gblobg -T0.25

grdimage blobg -JX3/2 -Cppdf.cpt -R -Ba2f1:'CO@-2@- Sensitivity (degC)':/0.2:'GYM'::'.Fraction of ModelRuns within 3-@~s@~':/SWen -X0.0 -Y-3.5 -K -O >> aaplot


psscale -D3.5/1/2/0.25c -Cppdf.cpt -Ba0.2f0.1:'.Fraction <3-@~s@~': -O >> aaplot


gv aaplot &


CO2 proxy data set for Park and Royer (2011)


click title to see Google Drive directory JParkCodes

carbon dioxide proxy data set (co2proxy_2010.dat) used by

Park, J., and D. L. Royer, Geologic constraints on the glacial amplification of Phanerozoic climate sensitivity, American J. Science, 311, 1-26, doi:10.2475/01.2011.01, 2011.

The sources of data are enumerated in the reference list. The dataset is an updated version of the data compilation reported by

Royer DL. 2006. CO2-forced climate thresholds during the Phanerozoic. Geochimica et Cosmochimica Acta, 70: 5665-5675.



CO2 Proxy data set from Royer et al 2007

click title to see Google Drive directory JParkCodes

Carbon dioxide proxy data set for the Phanerozoic (co2proxy_2007.dat) used by

Royer, D. L., R. A. Berner, and J. Park, Climate sensitivity constrained by carbon dioxide concentrations over the last 420 million years, Nature, v446, 530-532, 2007.

data compilation described in


Royer DL. 2006. CO2-forced climate thresholds during the Phanerozoic. Geochimica et Cosmochimica Acta, 70: 5665-5675.  

post-processing of GEOCARBSULFvolc simulation results

click title to see Google Drive directory JParkCodes


c program gcsv_ppdf.f
c 5/30/10 JJP
c f77 -o /Users/jjpark/bin/gcsv_ppdf gcsv_ppdf.f
c program to read variance reduction matrices for gcsfdata_export/search
c reads relative variance residual for CO2-proxy datafit to ndat*10000
c geocarb model runs

GEOCARBSULFvolc carbon-cycle simulation code

Google Drive JParkCodes folder gcsv10_export.f
click title to see Google Drive directory JParkCodes

c program gcsv10_export.f
c fortran version of bob berner's GEOCARBSULF BASIC code
c gcsv10 stands for GeoCarbSulfVolc algorithm for 2010
c THIS CODE CAN READ EITHER THE 2006/7 proxy-CO2 data file
c OR THE 2010 DATA FILE.
c 8/17/2010 -- confirmed the VNV parameter as 4, though Berner asked for 5,
c as motivated by Aaron Taylors thesis -- but other Taylor source is vnv=3
c this code developed in MacOS 10.4, using gcc compiler,
c takes 90minutes to run on a 2006-vintage Macbook Pro with 2 GHz Intel Core Duo processor
c takes 45minutes to run on a 2009-vintage Macbook Pro with 2.53 GHz Intel Core 2 Duo processor
c -- looped over range of input parameters, set up to fit proxy data for CO2 ppm
c 11/22/06 JJP
c proxy data has uncertainties that are sometimes large. To make the data-fitting robust,
c compute chi-squared using an uncertainty in the log-domain.
c xf77 -o /Users/jjpark/bin/gcsv10_export gcsv10_export.f /Users/jjpark/Plotxy/plotlib.a
c
c code uses plotit, a standalone fortran/C plotting subroutine in Xwindows
c this code can be found on JPark's website http://earth.geology.yale.edu/~jjpark/Code/plotit.tar
c instructions for PLOTIT at http://earth.geology.yale.edu/~jjpark/Code/plotit.html
c
c code updated from the code used for the paper
c
c Royer, D. L., Berner, R. A., and Park, J., 2007, Climate sensitivity constrained by carbon
c dioxide concentrations over the last 420 million years: Nature, v.446, p.530-532.
c
c and used for this paper
c
c Park, J., and D. L. Royer, Geologic constraints on the glacial amplification of Phanerozoic
c climate sensitivity, American J. Science, V.311, 2011, P. 1-26, DOI 10.2475/01.2011.01
c
c there are many test computations in Park and Royer (2011) that were effected
c with modified versions of this gcsv10_export.f code
c
c Output of this code is written to two files
c outoutout.dat is a simple table of CO2proxy data values averaged in 10My intervals
c out_gcsv10.dat is a large ASCII file (>50 Mb) that contains the GEOCARB carbon-cycle simulations
c for all parameter choices. This allows the user to examine the statistical behavior of the
c carbon-cycle simulations offline this code (e.g. in gcsv_ppdf.f)
c
c code converted from BASIC to run grid search over parameters
c Newton-Raphson damping applied here, to avoid NaN.
c divergent transition over 380-350My is interpolated better,
c correcting earlier convergence error
c the code reads a file proxydata.txt that contains proxy CO2 PPM data
c the GEOCARBSULF carbon-cycle model is computed at time steps of 1-My from 570Ma to present
c the proxy data ranges only from 420Ma to the present
c
c Five nested loops comprise the heart of the code, from outermost to innermost
c DeltaT = temp increase with CO2 doubling
c ACT = dimensionless activation-energy for silicate weathering
c FERT = plant CO2-fertilization coefficient
c LIFE = liverwort-based weathering factor -- less than vascular plants
c GYM = gymnosperm-based weathering factor (1=angiosperm weathering parameter)
c gymnosperms (e.g. conifers) might be better at weathering than angiosperms
c
c for each choice of parameters, we compute the carbon flux balance between terms
c that depend on silicate weathering (fBBS) and carbon burial/degassing (fB),
c each normalized with common factors derived by Bob Berner many years ago
c The carbon-cycle balance depends on CO2 level,
c because the weathering of silicates accelerates with warmer temps induced by greenhouse gases
c
c the code solves for the CO2 level that achieves the balance of carbon fluxes
c The CO2-dependent parameter fB is calculated at each time step from carbon and
c carbon isotope mass balance and values of all other parameters, based on
c geological and biological data, that affect the carbon cycle. Then the value of
c RCO2 (CO2 concentration normalized to the mean value for the past 1 million
c years = 250 ppm) is calculated by inversion from a complex expression for fB based on the
c greenhouse effect, CO2 fertilization of plant weathering, solar evolution, and
c changes in land temperature).
c
c for detailed exposition of the theoretical GEOCARB models, see
c
c Berner, R. A., The Phanerozoic Carbon Cycle: CO2 and O2, 150pp. Oxford University Press, 2004.
c
c and a collection of GEOCARB and GEOCARBSULF papers published over the past decade by
c Berner and Colleagues.
c
c Berner, R. A., 2006a, GEOCARBSULF: A combined model for Phanerozoic atmospheric O2 and CO2:
c Geochimica et Cosmochimica Acta, v.70, p.5653-5664.
c
c Berner, R. A., 2006b, Inclusion of the weathering of volcanic rocks in the GEOCARBSULF model:
c American Journal of Science, v.306, p.295-302, doi:10.2475/052006.01.
c
c Berner, R. A., 2008, Addendum to Inclusion of the weathering of volcanic rocks in the
c GEOCARBSULF model (Berner, R. A., 2006, v.306, p.295-302) American Journal of Science,
c v.308, p.100-103, doi:10.2475/01.2008.04.
c
c Berner, R. A., 2009, Phanerozoic atmospheric oxygen: New results using the GEOCARBSULF Model:
c American Journal of Science, v.309, p.603-606, doi:10.2475/07.2009.03.
c
c The code computes the aggregate data fit of GEOCARBSULF for each choice
c of the combined parameters, expressed as the fractional variance residual when expressed
c in the log domain without uncertainty weighting.
c description of the proxy data can be found in
c
c Royer et al, CO2 as a primary driver of Phanerozoic climate, GSA Today, March 2004.
c

Reminder for the Fortran plotting subroutine PLOTIT

click title to see Google Drive directory JParkCodes

plotit.tar

information README at

plotit.html


 
Link