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

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.