# This program does the calculation presented in Rohling et al. (2021 Science Advances)
# HOWEVER: some changes have been made to Antarctic ice volume description (VseqA), 
# to ensure that the sea-level budget is properly closed. Hence the results are not
# completely identical to those in the paper (but are so by close approximation).
library(ggplot2)
library(dplyr)
library(hexbin)
library(gcookbook)
library(zoo)
library(patchwork)
library(phenocamr)
library(pspline)
library(Hmisc)
library(splines)
library(GauPro)
library(stats)
library(forecast)
library(scales)
# <------ SET working directory
setwd("~/Dropbox/aaaa-HIVE/WORK/ARTICLES/deconvolution_in_R/SciAdvmodels in R")
# <------ SET file to be used
Input <- read.delim("WH-pos_t-SL-dOadjtoLR.txt", header=FALSE)
# Input time is set to count FORward in time, starting at 0, and the records need to have been interpolated with replicates removed already 
names(Input)<-c("t","SL","dOc") # assigns column names
t=Input$t
SLtot=Input$SL
dOc=Input$dOc
# 
t0<--max(t) # This tells you when the record starts in negative values for past (i.e., what age does time 0 represent)
realage<-t-max(t) # This gives the real ages in negative values for past
# 
Voc=1.335*10^18   # m^3
Aoc=361.9*10^12   # m^2
Hoc<-Voc/Aoc      # m
Loc=3700          # Pre-set & rounded approximation for  Hoc
eps=2*10^-3       # planoconvex ice sheet aspect ratio
VAnow=57.8        # modern Antarctic volume in mseq
VGnow=7.3         # modern Greenland volume in mseq
VxLmax=70         # LGM Laurentide (LIS) excess max volume in mseq
VxFmax=35         # LGM Fennoscandian (EIS) excess max volume in mseq
VxAgmax=15        # LGM Antarctic (AIS) excess max volume in mseq
VxGgmax=5         # LGM Greenland (GrIS) excess max volume in mseq
Vgtot<-VxLmax+VxFmax+VxAgmax+VxGgmax  # total LGM ice volume excess in mseq
# at any time, total global ice volume (in mseq) is negative of SL. Here we identify the range limits for ice volume in our record relative to present ice volume (in mseq)
Vmax<-max(-SLtot)
Vmin<-min(-SLtot)
## 
SLmax<-max(-SLtot)
## LGM-based scaling of the various ice sheets
VLmax<-VxLmax*(Vmax/Vgtot)
VFmax<-VxFmax*(Vmax/Vgtot)
VAgmax<-VxAgmax*(Vmax/Vgtot)
VGmax<-VxGgmax*(Vmax/Vgtot)
# Here we set base values for the snow precip.-isotope calculations
P=0.9
g=2
eta=10
lamb=63
delta0=-15
d0<-delta0
# For modern ice volume of roughly 65 mseq, gross accum is 0.8m/cy or 8m/kyr. We set "agrmod" with this value.
agrmod=8
Armod<-pi*(((1.5*65*Aoc)/(0.9*pi*eps))^(1/3))^2 # gross accum over total area of modern ice volume of roughly 65 mseq 
Armod
# Here we define current (timestep j in the calculations) and previous (j-1) values needed for increment-sensitive calculations
N<-NROW(SLtot)
N
Q<-N-3
tcurr<-t[2:N]
# tprev<-t[1:(N-1)]
SLcurr<-SLtot[2:N]
# SLprev<-SLtot[1:(N-1)]
SL<-SLtot
# Here we define a heuristic equation used to allow growth rate gradient between LIS and EIS
# exp(1) is Euler's number, the natural logarithm base
Gra<-1-(0.5*exp(1))^(20*SL/Vmax) # equation for growth rate gradient between LIS and EIS
# uses SLtot(j-1) in the Mathcad code because it depends on the preceding ice volume,
# but that is the same as SLtot here (except for the very youngest value).
# For computational reasons (vector lengths), it's easiest to use SLtot here (renamed to SL for consistency with Mathcad code)
# 
############# LAURENTIDE
calcL<-(Gra*(-SL*VLmax/Vmax))-((VAgmax*SL/SLmax)+VAgmax*(SL/SLmax)^2)
VL<-ifelse(-SL>0,calcL,0)
VseqL<-ifelse(VL<0,0, ifelse(VL>VLmax,VLmax,VL))
VseqLcurr<-VseqL[2:N]
VseqLprev<-VseqL[1:(N-1)]
ArL<-pi*(((1.5*VseqL*Aoc)/(0.9*pi*eps))^(1/3))^2 # Area (uses REAL volume; so compensate with 0.9 for density diff relative to mseq) 
ArLcurr<-ArL[2:N]
ArLprev<-ArL[1:(N-1)]
anetLcurr<-VseqLcurr-VseqLprev # net accumulation rate at timestep j
agrLcurr<-agrmod*ArLcurr/Armod
lossLcurr<-agrLcurr-anetLcurr
agrL<-append(0,agrLcurr)  # add a zero gross growth value for t=0 (i.e., between t=-1 and t=0)
agrLprev<-agrL[1:(N-1)]
LocLcurr<-Loc-cumsum(anetLcurr) # imposed change in ocean level (this had a PLUS in Mathcad.....!)
LocL<-append(Loc,LocLcurr) # add a imposed change in ocean level
LocLprev<-LocL[1:(N-1)]
delLcurr<-(0.55*(exp((((P*VseqLcurr)-eta)/lamb)^g))*(100-VseqLcurr))-(56.4+15)
dmeanL<-delLcurr # seed values 
calcL<-delLcurr # seed values
dxL<-delLcurr # seed values
# 
############# EURASIAN (using F for FENNOSCANDIAN)
calcLcomp<-(Gra*(-SL*VLmax/Vmax))
VLcomp<-ifelse(-SL>0,calcLcomp,0)
VF<--SL*((VLmax+VFmax)/Vmax)-VLcomp
VseqF<-ifelse(VF<0,0, ifelse(VF>VFmax,VFmax,VF))
VseqFcurr<-VseqF[2:N]
VseqFprev<-VseqF[1:(N-1)]
ArF<-pi*(((1.5*VseqF*Aoc)/(0.9*pi*eps))^(1/3))^2 # Area (uses REAL volume; so compensate with 0.9 for density diff relative to mseq) 
ArFcurr<-ArF[2:N]
ArFprev<-ArF[1:(N-1)]
anetFcurr<-VseqFcurr-VseqFprev # net accumulation rate at timestep j
agrFcurr<-agrmod*ArFcurr/Armod
lossFcurr<-agrFcurr-anetFcurr
agrF<-append(0,agrFcurr)  # add a zero gross growth value for t=0 (i.e., between t=-1 and t=0)
agrFprev<-agrF[1:(N-1)]
LocFcurr<-Loc-cumsum(anetFcurr) # imposed change in ocean level (this had a PLUS in Mathcad.....!)
LocF<-append(Loc,LocFcurr) # add a imposed change in ocean level
LocFprev<-LocF[1:(N-1)]
delFcurr<-(0.55*(exp((((P*VseqFcurr)-eta)/lamb)^g))*(100-VseqFcurr))-(56.4+15)
dmeanF<-delFcurr # seed values 
calcF<-delFcurr # seed values 
dxF<-delFcurr # seed values 
# 
############# ANTARCTIC
calcA<-VAnow+(VAgmax*(SL/SLmax)^2)
VseqA<-ifelse(SL>(VAnow+VGnow),0, ifelse(SL<=(VAnow+VGnow) & SL>2*VGnow,(VAnow+VGnow-SL), #factor 2 in VGnow because we compare with equal AIS+GrIS when GrIS is maximised
           ifelse(VGnow+(-SL/2)<=VGnow & VGnow+(-SL/2)>0,VAnow+(-SL/2),calcA)))
VseqAcurr<-VseqA[2:N]
VseqAprev<-VseqA[1:(N-1)]
ArA<-pi*(((1.5*VseqA*Aoc)/(0.9*pi*eps))^(1/3))^2 # Area (uses REAL volume; so compensate with 0.9 for density diff relative to mseq)
ArAcurr<-ArA[2:N]
ArAprev<-ArA[1:(N-1)]
anetAcurr<-VseqAcurr-VseqAprev # net accumulation rate at timestep j
agrAcurr<-agrmod*ArAcurr/Armod
lossAcurr<-agrAcurr-anetAcurr
agrA<-append(0,agrAcurr)  # add a zero gross growth value for t=0 (i.e., between t=-1 and t=0)
agrAprev<-agrA[1:(N-1)]
LocAcurr<-Loc-cumsum(anetAcurr) # imposed change in ocean level (this had a PLUS in Mathcad.....!)
LocA<-append(Loc,LocAcurr) # add a imposed change in ocean level
LocAprev<-LocA[1:(N-1)]
delAcurr<-(1*(exp((((P*VseqAcurr)-eta)/lamb)^g))*(100-VseqAcurr))-(102.55+15)
dmeanA<-delAcurr # seed values 
calcA<-delAcurr # seed values 
dxA<-delAcurr # seed values 
# 
############# GREENLAND
calcG<-VGnow+(-SL*VGmax/SLmax)
VseqG<-ifelse(VGnow+(-SL/2)<=VGnow & VGnow+(-SL/2)>0,VGnow+(-SL/2),
         ifelse(VGnow+(-SL/2)<0,0,calcG))
VseqGcurr<-VseqG[2:N]
VseqGprev<-VseqG[1:(N-1)]
ArG<-pi*(((1.5*VseqG*Aoc)/(0.9*pi*eps))^(1/3))^2 # Area (uses REAL volume; so compensate with 0.9 for density diff relative to mseq)
ArGcurr<-ArG[2:N]
ArGprev<-ArG[1:(N-1)]
anetGcurr<-VseqGcurr-VseqGprev # net accumulation rate at timestep j
agrGcurr<-agrmod*ArGcurr/Armod
lossGcurr<-agrGcurr-anetGcurr
agrG<-append(0,agrGcurr)  # add a zero gross growth value for t=0 (i.e., between t=-1 and t=0)
agrGprev<-agrG[1:(N-1)]
LocGcurr<-Loc-cumsum(anetGcurr) # imposed change in ocean level (this had a PLUS in Mathcad.....!)
LocG<-append(Loc,LocGcurr) # add a imposed change in ocean level
LocGprev<-LocG[1:(N-1)]
delGcurr<-(1*(exp((((P*VseqGcurr)-eta)/lamb)^g))*(100-VseqGcurr))-(102.55+15)
dmeanG<-delGcurr # seed values 
calcG<-delGcurr # seed values 
dxG<-delGcurr # seed values 
# # 
# Northern and Southern Hemisphere ice volumes
VNH<-VseqLcurr+VseqFcurr+VseqGcurr
VSH<-VseqAcurr
# 
# Now we do the mean ice-sheet d18O calculations, which are iterative (dependent on their own previous value)
for (i in 2:length(anetLcurr)) {
  calcL[i]<-(((agrLcurr[i]*delLcurr[i])-((lossLcurr[i])*dmeanL[i-1])+(VseqLprev[i]*dmeanL[i-1]))/(VseqLcurr[i]))
  calcF[i]<-(((agrFcurr[i]*delFcurr[i])-((lossFcurr[i])*dmeanF[i-1])+(VseqFprev[i]*dmeanF[i-1]))/(VseqFcurr[i]))
  calcA[i]<-(((agrAcurr[i]*delAcurr[i])-((lossAcurr[i])*dmeanA[i-1])+(VseqAprev[i]*dmeanA[i-1]))/(VseqAcurr[i]))
  calcG[i]<-(((agrGcurr[i]*delGcurr[i])-((lossGcurr[i])*dmeanG[i-1])+(VseqGprev[i]*dmeanG[i-1]))/(VseqGcurr[i]))
  dxL[i]<-if (VseqLcurr[i]<=0) {d0} else if (calcL[i] > d0) {d0} else if (VseqLcurr[i]<=VseqLprev[i]) {
    dmeanL[i-1]} else {calcL[i]}
  dxF[i]<-if (VseqFcurr[i]<=0) {d0} else if (calcF[i] > d0) {d0} else if (VseqFcurr[i]<=VseqFprev[i]) {
    dmeanF[i-1]} else {calcF[i]}
  dxA[i]<-if (VseqAcurr[i]<=0) {d0} else if (calcA[i] > d0) {d0} else if (VseqAcurr[i]<=VseqAprev[i]) {
    dmeanA[i-1]} else {calcA[i]}
  dxG[i]<-if (VseqGcurr[i]<=0) {d0} else if (calcG[i] > d0) {d0} else if (VseqGcurr[i]<=VseqGprev[i]) {
    dmeanG[i-1]} else {calcG[i]}
  dmeanL<-dxL
  dmeanF<-dxF
  dmeanA<-dxA
  dmeanG<-dxG
}
# Then we calculate each ice sheet's impact on global d18Ow, normalised to present
# note the minus before dmeanX, which applies because ice sheet d18O inversely affects ocean d18O
DdocLcurr<--dmeanL*VseqLcurr/LocLcurr
DdocLnorm<-DdocLcurr-DdocLcurr[length(DdocLcurr)]
DdocFcurr<--dmeanF*VseqFcurr/LocFcurr
DdocFnorm<-DdocFcurr-DdocFcurr[length(DdocFcurr)]
DdocAcurr<--dmeanA*VseqAcurr/LocAcurr
DdocAnorm<-DdocAcurr-DdocAcurr[length(DdocAcurr)]
DdocGcurr<--dmeanG*VseqGcurr/LocGcurr
DdocGnorm<-DdocGcurr-DdocGcurr[length(DdocGcurr)]
DdocAll<-DdocLnorm+DdocFnorm+DdocAnorm+DdocGnorm
# 
# Finally, we calculate deep-sea T from the dOc residuals
dOccurr<-dOc[2:N]
dOcnorm<-dOccurr-dOccurr[length(dOccurr)]
dOcT<-(dOcnorm-DdocAll)/-0.25
# 
# Now we need to make a good (intuitive) time axis for plotting
age<-tcurr-tcurr[max(tcurr)]
# 
# And now we plot things up.
plot(age,VseqLcurr,type="l",pch="",col="red", xlab="",ylab="", ylim=c(-5,75),xlim=rev(range(c(min(age),max(age)))))
  lines(age,VseqFcurr,col="blue")
  lines(age,VseqAcurr,col="magenta")
  lines(age,VseqGcurr,col="green")
  minor.tick(nx=4,ny=4,tick.ratio=0.7)
  grid(nx=NULL,ny=NULL,col=alpha("lightgrey",0.5), lty = 1,
       lwd = 1)
  mtext("Ice volume (mseq)", side=2,line=2.5,cex=1)
  mtext("Age (ka)", side=1,line=3,cex=1)
  legend("topright", legend=c("LIS", "EIS", "AIS", "GrIS"),
         col=c("red", "blue", "magenta", "green"), lty=1, cex=0.8, text.font=4)
plot(age,dmeanL,type="l",pch="", col="red", xlab="",ylab="", ylim=c(-60,-10),xlim=rev(range(c(min(age),max(age)))))
  lines(age,dmeanF,col="blue")
  lines(age,dmeanA,col="magenta")
  lines(age,dmeanG,col="green")
  minor.tick(nx=4,ny=5,tick.ratio=0.7)
  grid(nx=NULL,ny=NULL,col=alpha("lightgrey",0.5), lty = 1,
       lwd = 1)
  mtext("Mean ice delta18O (permille)", side=2,line=2.5,cex=1)
  mtext("Age (ka)", side=1,line=3,cex=1)
  legend("bottomright", legend=c("LIS", "EIS", "AIS", "GrIS"),
         col=c("red", "blue", "magenta", "green"), lty=1, cex=0.8, text.font=4)
plot(VseqLcurr,dmeanL,type="l",pch="*",col="red", xlab="",ylab="", xlim=c(0,70), ylim=c(-65,-15))
  lines(VseqFcurr,dmeanF,col="blue")
  lines(VseqAcurr,dmeanA,col="magenta")
  lines(VseqGcurr,dmeanG,col="green")
  minor.tick(nx=5,ny=5,tick.ratio=0.7)
  grid(nx=NULL,ny=NULL,col=alpha("lightgrey",0.5), lty = 1,
       lwd = 1)
  mtext("Mean ice delta18O (permille)", side=2,line=2.5,cex=1)
  mtext("Ice volume (mseq)", side=1,line=3,cex=1)
  legend("topright", legend=c("LIS", "EIS", "AIS", "GrIS"),
         col=c("red", "blue", "magenta", "green"), lty=1, cex=0.8, text.font=4)
plot(age,DdocLnorm,type="l",pch="", col="red", xlab="",ylab="", ylim=c(-1,1.5),xlim=rev(range(c(min(age),max(age)))))
  lines(age,DdocFnorm,col="blue")
  lines(age,DdocAnorm,col="magenta")
  lines(age,DdocGnorm,col="green")
  lines(age,DdocAll,col=alpha("black",0.5))
  minor.tick(nx=4,ny=5,tick.ratio=0.7)
  grid(nx=NULL,ny=NULL,col=alpha("lightgrey",0.5), lty = 1,
       lwd = 1)
  mtext("Global delta18Ow impacts (permille)", side=2,line=2.5,cex=1)
  mtext("Age (ka)", side=1,line=3,cex=1)
  legend("topright", legend=c("LIS", "EIS", "AIS", "GrIS","All"),
         col=c("red", "blue", "magenta", "green","grey"), lty=1, cex=0.8, text.font=4)
plot(age,dOcnorm,type="l",pch="", col="grey", xlab="",ylab="", ylim=c(-3,2),xlim=rev(range(c(min(age),max(age)))))
  lines(age,DdocAll, col="pink")
  minor.tick(nx=4,ny=5,tick.ratio=0.7)
  grid(nx=NULL,ny=NULL,col=alpha("lightgrey",0.5), lty = 1,
       lwd = 1)
  mtext("Global delta18Oc and delta18Ow (permille)", side=2,line=2.5,cex=1)
  mtext("Age (ka)", side=1,line=3,cex=1)
  legend("topright", legend=c("delta18Oc", "delta18Ow"),
         col=c("grey", "pink"), lty=1, cex=0.8, text.font=4)
plot(age,dOcT,type="l",pch="", col="black", xlab="",ylab="", ylim=c(-3,8),xlim=rev(range(c(min(age),max(age)))))
  minor.tick(nx=4,ny=4,tick.ratio=0.7)
  grid(nx=NULL,ny=NULL,col=alpha("lightgrey",0.5), lty = 1,
       lwd = 1)
  mtext("Deep-sea temperature (deg.C)", side=2,line=2.5,cex=1)
  mtext("Age (ka)", side=1,line=3,cex=1)
  legend("bottomright", legend=c("Tw"),
         col=c("black"), lty=1, cex=0.8, text.font=4)
plot(age,VNH,type="l",pch="", col="red", xlab="",ylab="", ylim=c(0,115),xlim=rev(range(c(min(age),max(age)))))
  lines(age,VSH, col="blue")  
  minor.tick(nx=4,ny=4,tick.ratio=0.7)
  grid(nx=NULL,ny=NULL,col=alpha("lightgrey",0.5), lty = 1,
       lwd = 1)
  mtext("Ice volume (mseq)", side=2,line=2.5,cex=1)
  mtext("Age (ka)", side=1,line=3,cex=1)
  legend("topright", legend=c("NH","SH"),
         col=c("red","blue"), lty=1, cex=0.8, text.font=4)
# # END
# 
# ADDITIONAL: sea-level budget closure test
plot(SL,VseqA+VseqG+VseqL+VseqF, pch=".", col="black")
  points(SL,VseqG, pch=".", col="green")
  points(SL,VseqL, pch=".", col="red")
  points(SL,VseqF, pch=".", col="blue")
  points(SL,VseqA, pch=".", col="magenta")
  abline(v=2*VGnow, lty=2)
  abline(v=0, lty=2)
  abline(v=VAnow+VGnow, lty=2)
  abline(h=VAnow,lty=2)
  abline(h=VAnow-VGnow,lty=2)
  abline(h=VGnow, lty=2)
  legend("topright", legend=c("V_LIS", "V-EIS", "V_AIS", "V_GrIS", "V_all"),
         col=c("red", "blue", "magenta", "green", "black"), lty=1, cex=0.8, text.font=4)