Skip to content

Commit

Permalink
Harper et al. PNAS resubmission release updates
Browse files Browse the repository at this point in the history
  • Loading branch information
dustintharper committed May 27, 2024
1 parent f9c11a0 commit 4695234
Show file tree
Hide file tree
Showing 51 changed files with 102 additions and 5,167 deletions.
11 changes: 5 additions & 6 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -3,14 +3,13 @@
.Rapp.history

# Session Data files
Harperetal_subm/out


# User-specific files
Harperetal_subm/out/inversion_sum.csv
Harperetal_subm/RevisionApril2024/out_senstest
Harperetal_subm/RevisionApril2024/out_primary/LPEE_prim.rda
Harperetal_subm/RevisionApril2024/.DS_Store
Harperetal_subm/RevisionApril2024/data/.DS_Store
Harperetal_resubm/out_senstest
Harperetal_resubm/out_primary/LPEE_prim.rda
Harperetal_resubm/R.DS_Store
Harperetal_resubm/data/.DS_Store

# Example code in package build process
*-Ex.R
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ parms <- c("sal", "tempC", "press", "xca", "xmg", "xso4", "d11Bsw", "d18Osw", "p
"m.1", "m.2", "c.1", "c.2", "alpha", "d11Bf.1", "d11Bf.2", "d18Of", "mgcaf")

# Read in proxy time series data
prox.in <- readRDS(file = "Harperetal_subm/RevisionApril2024/data/ShatskyLPEE_data.rds")
prox.in <- readRDS(file = "Harperetal_resubm/data/ShatskyLPEE_data.rds")

# Setup age range and bins
ages.prox <- unique(round(prox.in$age))
Expand Down Expand Up @@ -155,7 +155,7 @@ pH.u = 7.75
############################################################################################
# Read in D[CO3=] from file and linearly interpolate for ages associated with each time step

Dco3_in <- as.data.frame(readRDS(file = "Harperetal_subm/RevisionApril2024/data/Dco3_bw.rds"))
Dco3_in <- as.data.frame(readRDS(file = "Harperetal_resubm/data/Dco3_bw.rds"))
Dco3.interp <- approx(Dco3_in$age, Dco3_in$Dco3, xout=ages.prox, method="linear")
Dco3_2s.interp <- approx(Dco3_in$age, Dco3_in$Dco3_2s, xout=ages.prox, method="linear")
Dco3.avg.pri <- Dco3.interp[["y"]]
Expand All @@ -175,15 +175,15 @@ for (i in 1:length(Dco3.avg.pri)){
############################################################################################
# Read in DIC from LOSCAR simulation output: mean of 2 sims for PETM and 2 sims for ETM-2 with ZT19 long-term carb chem

dic_LOSCAR <- readRDS(file = "Harperetal_subm/RevisionApril2024/data/dic_LOSCAR.rds")
dic_LOSCAR <- readRDS(file = "Harperetal_resubm/data/dic_LOSCAR.rds")
# Linearly interpolate DIC for ages associated with each time step using input DIC time series
dic.interp.LOSCAR <- approx(dic_LOSCAR$dic.LOSCAR.x, dic_LOSCAR$dic.LOSCAR.y, xout=ages.prox, method="linear")
dic.interp.LOSCAR.err <- approx(dic_LOSCAR$dic.LOSCAR.x, dic_LOSCAR$dic.LOSCAR.2s, xout=ages.prox, method="linear")
dic.LOSCAR <- dic.interp.LOSCAR[["y"]]
dic.LOSCAR.err <- dic.interp.LOSCAR.err[["y"]] / 2

# Read in DIC time series from Haynes and Hönisch (2020)
dic_HH <- readRDS("Harperetal_subm/RevisionApril2024/data/dic_HH.rds")
dic_HH <- readRDS("Harperetal_resubm/data/dic_HH.rds")
# Linearly interpolate DIC for ages associated with each time step using input DIC time series
dic.interp.HH <- approx(dic_HH$dic.HH.x, dic_HH$dic.HH.y, xout=ages.prox, method="linear")
dic.HH.meanerr <- rowMeans(cbind(dic_HH$dic.HH.2sp, dic_HH$dic.HH.2sn))
Expand Down Expand Up @@ -289,7 +289,7 @@ data <- list("d11Bf.data1" = clean.d11B1$d11B,
############################################################################################
# Run the inversion

system.time({jout = jags.parallel(model.file = "Harperetal_subm/RevisionApril2024/code/_ForamPSMLPEE.R",
system.time({jout = jags.parallel(model.file = "Harperetal_resubm/code/_ForamPSMLPEE.R",
parameters.to.save = parms, data = data, inits = NULL,
n.chains = 9, n.iter = 800000, n.burnin = 500000, n.thin = 200)})
# 500k burn in, 9 chains, 800k iterations takes 10 hours
Expand All @@ -299,10 +299,10 @@ system.time({jout = jags.parallel(model.file = "Harperetal_subm/RevisionApril202
# Display summary statistics and save summary as .csv

#View(jout$BUGSoutput$summary)
write.csv(jout$BUGSoutput$summary, "Harperetal_subm/RevisionApril2024/out_primary/inversion_sum.csv")
write.csv(jout$BUGSoutput$summary, "Harperetal_resubm/out_primary/inversion_sum.csv")

save(ages.prox, file = "Harperetal_subm/RevisionApril2024/out_primary/ages.prox.rda")
#save(jout, file = "Harperetal_subm/RevisionApril2024/out_primary/LPEE_prim.rda")
save(ages.prox, file = "Harperetal_resubm/out_primary/ages.prox.rda")
#save(jout, file = "Harperetal_resubm/out_primary/LPEE_prim.rda")

############################################################################################

Expand Down
Original file line number Diff line number Diff line change
@@ -1,8 +1,10 @@
#
# This contains the R code to plot Figures 2, 3, and S5 in Harper et al., submission to PNAS
#
# Note that full inversion output files used here ('jout') are too big to include in the repository.
# Thus inversions would need to be re-run using "...Driver..." .R files.

jout <- jout_prim
load(file = "Harperetal_resubm/out_primary/ages.prox.rda")

###########################################################################################
# Temp, pCO2, and pH time series draw plots
Expand Down Expand Up @@ -117,7 +119,7 @@ print(c(mean(-LP.Dpco2), quantile(-LP.Dpco2, 0.025), quantile(-LP.Dpco2, 0.975))
###########################################################################################
# LOSCAR Supplemental Figure - 5800 Gt release of -35 per mille C
###########################################################################################
LOSCAR <- readRDS("Harperetal_subm/RevisionApril2024/data/LOSCAR_PETM5800.rds")
LOSCAR <- readRDS("Harperetal_resubm/data/LOSCAR_PETM5800.rds")

# read in and format PSM output
jout_summ <- jout$BUGSoutput$summary
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ parms <- c("sal", "tempC", "press", "xca", "xmg", "xso4", "d11Bsw", "d18Osw", "p
"m.1", "m.2", "c.1", "c.2", "alpha", "d11Bf.1", "d11Bf.2", "d18Of", "mgcaf")

# Read in proxy time series data
prox.in <- readRDS(file = "Harperetal_subm/RevisionApril2024/data/ShatskyLPEE_data.rds")
prox.in <- readRDS(file = "Harperetal_resubm/data/ShatskyLPEE_data.rds")

# Setup age range and bins
ages.prox <- unique(round(prox.in$age))
Expand Down Expand Up @@ -155,7 +155,7 @@ pH.u = 7.75
############################################################################################
# Read in D[CO3=] from file and linearly interpolate for ages associated with each time step

Dco3_in <- as.data.frame(readRDS(file = "Harperetal_subm/RevisionApril2024/data/Dco3_bw.rds"))
Dco3_in <- as.data.frame(readRDS(file = "Harperetal_resubm/data/Dco3_bw.rds"))
Dco3.interp <- approx(Dco3_in$age, Dco3_in$Dco3, xout=ages.prox, method="linear")
Dco3_2s.interp <- approx(Dco3_in$age, Dco3_in$Dco3_2s, xout=ages.prox, method="linear")
Dco3.avg.pri <- Dco3.interp[["y"]]
Expand All @@ -175,15 +175,15 @@ for (i in 1:length(Dco3.avg.pri)){
############################################################################################
# Read in DIC from LOSCAR simulation output: mean of 2 sims for PETM and 2 sims for ETM-2 with ZT19 long-term carb chem

dic_LOSCAR <- readRDS(file = "Harperetal_subm/RevisionApril2024/data/dic_LOSCAR.rds")
dic_LOSCAR <- readRDS(file = "Harperetal_resubm/data/dic_LOSCAR.rds")
# Linearly interpolate DIC for ages associated with each time step using input DIC time series
dic.interp.LOSCAR <- approx(dic_LOSCAR$dic.LOSCAR.x, dic_LOSCAR$dic.LOSCAR.y, xout=ages.prox, method="linear")
dic.interp.LOSCAR.err <- approx(dic_LOSCAR$dic.LOSCAR.x, dic_LOSCAR$dic.LOSCAR.2s, xout=ages.prox, method="linear")
dic.LOSCAR <- dic.interp.LOSCAR[["y"]]
dic.LOSCAR.err <- dic.interp.LOSCAR.err[["y"]] / 2

# Read in DIC time series from Haynes and Hönisch (2020)
dic_HH <- readRDS("Harperetal_subm/RevisionApril2024/data/dic_HH.rds")
dic_HH <- readRDS("Harperetal_resubm/data/dic_HH.rds")
# Linearly interpolate DIC for ages associated with each time step using input DIC time series
dic.interp.HH <- approx(dic_HH$dic.HH.x, dic_HH$dic.HH.y, xout=ages.prox, method="linear")
dic.HH.meanerr <- rowMeans(cbind(dic_HH$dic.HH.2sp, dic_HH$dic.HH.2sn))
Expand Down Expand Up @@ -289,7 +289,7 @@ data <- list("d11Bf.data1" = clean.d11B1$d11B,
############################################################################################
# Run the inversion

system.time({jout = jags.parallel(model.file = "Harperetal_subm/RevisionApril2024/code/_ForamPSMLPEE.R",
system.time({jout = jags.parallel(model.file = "Harperetal_resubm/code/_ForamPSMLPEE.R",
parameters.to.save = parms, data = data, inits = NULL,
n.chains = 9, n.iter = 800000, n.burnin = 500000, n.thin = 200)})
# 500k burn in, 9 chains, 800k iterations takes 10 hours
Expand All @@ -299,8 +299,8 @@ system.time({jout = jags.parallel(model.file = "Harperetal_subm/RevisionApril202
# Display summary statistics and save summary as .csv

#View(jout$BUGSoutput$summary)
#write.csv(jout$BUGSoutput$summary, "Harperetal_subm/RevisionApril2024/out/inversion_sum.csv")
#write.csv(jout$BUGSoutput$summary, "Harperetal_resubm/out_senstest/inversion_sum.csv")

############################################################################################
save(jout, file = "Harperetal_subm/RevisionApril2024/out/LPEE_pHMgCa3.rda")
save(jout, file = "Harperetal_resubm/out_senstest/LPEE_pHMgCa3.rda")

Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ parms <- c("sal", "tempC", "press", "xca", "xmg", "xso4", "d11Bsw", "d18Osw", "p
"m.1", "m.2", "c.1", "c.2", "alpha", "d11Bf.1", "d11Bf.2", "d18Of", "mgcaf")

# Read in proxy time series data
prox.in <- readRDS(file = "Harperetal_subm/RevisionApril2024/data/ShatskyLPEE_data.rds")
prox.in <- readRDS(file = "Harperetal_resubm/data/ShatskyLPEE_data.rds")

# Setup age range and bins
ages.prox <- unique(round(prox.in$age))
Expand Down Expand Up @@ -155,7 +155,7 @@ pH.u = 7.75
############################################################################################
# Read in D[CO3=] from file and linearly interpolate for ages associated with each time step

Dco3_in <- as.data.frame(readRDS(file = "Harperetal_subm/RevisionApril2024/data/Dco3_bw.rds"))
Dco3_in <- as.data.frame(readRDS(file = "Harperetal_resubm/data/Dco3_bw.rds"))
Dco3.interp <- approx(Dco3_in$age, Dco3_in$Dco3, xout=ages.prox, method="linear")
Dco3_2s.interp <- approx(Dco3_in$age, Dco3_in$Dco3_2s, xout=ages.prox, method="linear")
Dco3.avg.pri <- Dco3.interp[["y"]]
Expand All @@ -175,15 +175,15 @@ for (i in 1:length(Dco3.avg.pri)){
############################################################################################
# Read in DIC from LOSCAR simulation output: mean of 2 sims for PETM and 2 sims for ETM-2 with ZT19 long-term carb chem

dic_LOSCAR <- readRDS(file = "Harperetal_subm/RevisionApril2024/data/dic_LOSCAR.rds")
dic_LOSCAR <- readRDS(file = "Harperetal_resubm/data/dic_LOSCAR.rds")
# Linearly interpolate DIC for ages associated with each time step using input DIC time series
dic.interp.LOSCAR <- approx(dic_LOSCAR$dic.LOSCAR.x, dic_LOSCAR$dic.LOSCAR.y, xout=ages.prox, method="linear")
dic.interp.LOSCAR.err <- approx(dic_LOSCAR$dic.LOSCAR.x, dic_LOSCAR$dic.LOSCAR.2s, xout=ages.prox, method="linear")
dic.LOSCAR <- dic.interp.LOSCAR[["y"]]
dic.LOSCAR.err <- dic.interp.LOSCAR.err[["y"]] / 2

# Read in DIC time series from Haynes and Hönisch (2020)
dic_HH <- readRDS("Harperetal_subm/RevisionApril2024/data/dic_HH.rds")
dic_HH <- readRDS("Harperetal_resubm/data/dic_HH.rds")
# Linearly interpolate DIC for ages associated with each time step using input DIC time series
dic.interp.HH <- approx(dic_HH$dic.HH.x, dic_HH$dic.HH.y, xout=ages.prox, method="linear")
dic.HH.meanerr <- rowMeans(cbind(dic_HH$dic.HH.2sp, dic_HH$dic.HH.2sn))
Expand Down Expand Up @@ -289,7 +289,7 @@ data <- list("d11Bf.data1" = clean.d11B1$d11B,
############################################################################################
# Run the inversion

system.time({jout = jags.parallel(model.file = "Harperetal_subm/RevisionApril2024/code/_ForamPSMLPEE.R",
system.time({jout = jags.parallel(model.file = "Harperetal_resubm/code/_ForamPSMLPEE.R",
parameters.to.save = parms, data = data, inits = NULL,
n.chains = 9, n.iter = 800000, n.burnin = 500000, n.thin = 200)})
# 500k burn in, 9 chains, 800k iterations takes 10 hours
Expand All @@ -299,8 +299,8 @@ system.time({jout = jags.parallel(model.file = "Harperetal_subm/RevisionApril202
# Display summary statistics and save summary as .csv

#View(jout$BUGSoutput$summary)
#write.csv(jout$BUGSoutput$summary, "Harperetal_subm/RevisionApril2024/out/inversion_sum.csv")
#write.csv(jout$BUGSoutput$summary, "Harperetal_resubm/out_senstest/inversion_sum.csv")

############################################################################################
save(jout, file = "Harperetal_subm/RevisionApril2024/out/LPEE_pHMgCa7.rda")
save(jout, file = "Harperetal_resubm/out_senstest/LPEE_pHMgCa7.rda")

Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ parms <- c("sal", "tempC", "press", "xca", "xmg", "xso4", "d11Bsw", "d18Osw", "p
"m.1", "m.2", "c.1", "c.2", "alpha", "d11Bf.1", "d11Bf.2", "d18Of", "mgcaf")

# Read in proxy time series data
prox.in <- readRDS(file = "Harperetal_subm/RevisionApril2024/data/ShatskyLPEE_data.rds")
prox.in <- readRDS(file = "Harperetal_resubm/data/ShatskyLPEE_data.rds")

# Setup age range and bins
ages.prox <- unique(round(prox.in$age))
Expand Down Expand Up @@ -155,7 +155,7 @@ pH.u = 7.75
############################################################################################
# Read in D[CO3=] from file and linearly interpolate for ages associated with each time step

Dco3_in <- as.data.frame(readRDS(file = "Harperetal_subm/RevisionApril2024/data/Dco3_bw.rds"))
Dco3_in <- as.data.frame(readRDS(file = "Harperetal_resubm/data/Dco3_bw.rds"))
Dco3.interp <- approx(Dco3_in$age, Dco3_in$Dco3, xout=ages.prox, method="linear")
Dco3_2s.interp <- approx(Dco3_in$age, Dco3_in$Dco3_2s, xout=ages.prox, method="linear")
Dco3.avg.pri <- Dco3.interp[["y"]]
Expand All @@ -175,15 +175,15 @@ for (i in 1:length(Dco3.avg.pri)){
############################################################################################
# Read in DIC from LOSCAR simulation output: mean of 2 sims for PETM and 2 sims for ETM-2 with ZT19 long-term carb chem

dic_LOSCAR <- readRDS(file = "Harperetal_subm/RevisionApril2024/data/dic_LOSCAR.rds")
dic_LOSCAR <- readRDS(file = "Harperetal_resubm/data/dic_LOSCAR.rds")
# Linearly interpolate DIC for ages associated with each time step using input DIC time series
dic.interp.LOSCAR <- approx(dic_LOSCAR$dic.LOSCAR.x, dic_LOSCAR$dic.LOSCAR.y, xout=ages.prox, method="linear")
dic.interp.LOSCAR.err <- approx(dic_LOSCAR$dic.LOSCAR.x, dic_LOSCAR$dic.LOSCAR.2s, xout=ages.prox, method="linear")
dic.LOSCAR <- dic.interp.LOSCAR[["y"]]
dic.LOSCAR.err <- dic.interp.LOSCAR.err[["y"]] / 2

# Read in DIC time series from Haynes and Hönisch (2020)
dic_HH <- readRDS("Harperetal_subm/RevisionApril2024/data/dic_HH.rds")
dic_HH <- readRDS("Harperetal_resubm/data/dic_HH.rds")
# Linearly interpolate DIC for ages associated with each time step using input DIC time series
dic.interp.HH <- approx(dic_HH$dic.HH.x, dic_HH$dic.HH.y, xout=ages.prox, method="linear")
dic.HH.meanerr <- rowMeans(cbind(dic_HH$dic.HH.2sp, dic_HH$dic.HH.2sn))
Expand Down Expand Up @@ -289,7 +289,7 @@ data <- list("d11Bf.data1" = clean.d11B1$d11B,
############################################################################################
# Run the inversion

system.time({jout = jags.parallel(model.file = "Harperetal_subm/RevisionApril2024/code/_ForamPSMLPEE.R",
system.time({jout = jags.parallel(model.file = "Harperetal_resubm/code/_ForamPSMLPEE.R",
parameters.to.save = parms, data = data, inits = NULL,
n.chains = 9, n.iter = 800000, n.burnin = 500000, n.thin = 200)})
# 500k burn in, 9 chains, 800k iterations takes 10 hours
Expand All @@ -299,8 +299,8 @@ system.time({jout = jags.parallel(model.file = "Harperetal_subm/RevisionApril202
# Display summary statistics and save summary as .csv

#View(jout$BUGSoutput$summary)
#write.csv(jout$BUGSoutput$summary, "Harperetal_subm/RevisionApril2024/out/inversion_sum.csv")
#write.csv(jout$BUGSoutput$summary, "Harperetal_resubm/out_senstest/inversion_sum.csv")

############################################################################################
save(jout, file = "Harperetal_subm/RevisionApril2024/out/LPEE_nopHd18O.rda")
save(jout, file = "Harperetal_resubm/out_senstest/LPEE_nopHd18O.rda")

Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ parms <- c("sal", "tempC", "press", "xca", "xmg", "xso4", "d11Bsw", "d18Osw", "p
"m.1", "m.2", "c.1", "c.2", "alpha", "d11Bf.1", "d11Bf.2", "d18Of", "mgcaf")

# Read in proxy time series data
prox.in <- readRDS(file = "Harperetal_subm/RevisionApril2024/data/ShatskyLPEE_data.rds")
prox.in <- readRDS(file = "Harperetal_resubm/data/ShatskyLPEE_data.rds")

# Setup age range and bins
ages.prox <- unique(round(prox.in$age))
Expand Down Expand Up @@ -155,7 +155,7 @@ pH.u = 7.75
############################################################################################
# Read in D[CO3=] from file and linearly interpolate for ages associated with each time step

Dco3_in <- as.data.frame(readRDS(file = "Harperetal_subm/RevisionApril2024/data/Dco3_bw.rds"))
Dco3_in <- as.data.frame(readRDS(file = "Harperetal_resubm/data/Dco3_bw.rds"))
Dco3.interp <- approx(Dco3_in$age, Dco3_in$Dco3, xout=ages.prox, method="linear")
Dco3_2s.interp <- approx(Dco3_in$age, Dco3_in$Dco3_2s, xout=ages.prox, method="linear")
Dco3.avg.pri <- Dco3.interp[["y"]]
Expand All @@ -175,15 +175,15 @@ for (i in 1:length(Dco3.avg.pri)){
############################################################################################
# Read in DIC from LOSCAR simulation output: mean of 2 sims for PETM and 2 sims for ETM-2 with ZT19 long-term carb chem

dic_LOSCAR <- readRDS(file = "Harperetal_subm/RevisionApril2024/data/dic_LOSCAR.rds")
dic_LOSCAR <- readRDS(file = "Harperetal_resubm/data/dic_LOSCAR.rds")
# Linearly interpolate DIC for ages associated with each time step using input DIC time series
dic.interp.LOSCAR <- approx(dic_LOSCAR$dic.LOSCAR.x, dic_LOSCAR$dic.LOSCAR.y, xout=ages.prox, method="linear")
dic.interp.LOSCAR.err <- approx(dic_LOSCAR$dic.LOSCAR.x, dic_LOSCAR$dic.LOSCAR.2s, xout=ages.prox, method="linear")
dic.LOSCAR <- dic.interp.LOSCAR[["y"]]
dic.LOSCAR.err <- dic.interp.LOSCAR.err[["y"]] / 2

# Read in DIC time series from Haynes and Hönisch (2020)
dic_HH <- readRDS("Harperetal_subm/RevisionApril2024/data/dic_HH.rds")
dic_HH <- readRDS("Harperetal_resubm/data/dic_HH.rds")
# Linearly interpolate DIC for ages associated with each time step using input DIC time series
dic.interp.HH <- approx(dic_HH$dic.HH.x, dic_HH$dic.HH.y, xout=ages.prox, method="linear")
dic.HH.meanerr <- rowMeans(cbind(dic_HH$dic.HH.2sp, dic_HH$dic.HH.2sn))
Expand Down Expand Up @@ -289,7 +289,7 @@ data <- list("d11Bf.data1" = clean.d11B1$d11B,
############################################################################################
# Run the inversion

system.time({jout = jags.parallel(model.file = "Harperetal_subm/RevisionApril2024/code/_ForamPSMLPEE_MgCasw.R",
system.time({jout = jags.parallel(model.file = "Harperetal_resubm/code/_ForamPSMLPEE_MgCasw.R",
parameters.to.save = parms, data = data, inits = NULL,
n.chains = 9, n.iter = 800000, n.burnin = 500000, n.thin = 200)})
# 500k burn in, 9 chains, 800k iterations takes 10 hours
Expand All @@ -299,8 +299,8 @@ system.time({jout = jags.parallel(model.file = "Harperetal_subm/RevisionApril202
# Display summary statistics and save summary as .csv

#View(jout$BUGSoutput$summary)
#write.csv(jout$BUGSoutput$summary, "Harperetal_subm/RevisionApril2024/out/inversion_sum.csv")
#write.csv(jout$BUGSoutput$summary, "Harperetal_resubm/out_senstest/inversion_sum.csv")

############################################################################################
save(jout, file = "Harperetal_subm/RevisionApril2024/out/LPEE_MgCasw.rda")
save(jout, file = "Harperetal_resubm/out_senstest/LPEE_MgCasw.rda")

Loading

0 comments on commit 4695234

Please sign in to comment.