# Minimum working example for testing the moving rain code. 
# Created by: Jeremy Caves, Dan Ibarra, and Aviv Bachan
# Created on: 3/16/2016
# Last Modified: 10/03/2017 by Dan Ibarra

#setwd("~/Desktop/PlantWxFeedbacks_MS/Plant Wx V5/Code+Misc")
setwd("~/Documents/PlantWx_Submission/PlantWxFeedbacks_MS/Plant Wx V5/Code+Misc")

library(rootSolve)
library(pracma)
library(wesanderson)
library(ggplot2)
library(magrittr)

# graphics.off()  #close plots
# rm(list = ls()) #clear global env
# cat("\014") #clear console

#### Hydrologic Functions ####
myODE <- function(Pe = 1, u = 7, L = 1e7,
                   n = 500, kp = 4e-6, BC_flag_L = 1, BC_flag_R = 3,
                   Win_LBC = 87.5, Win_RBC = 0, omega  = 2.6, 
                   Eo = 3.17e-5 , budyko_flag = 1){
  
  #domain parameters
  dx <- L/n
  x <- seq(dx,L,by = dx)
  
  # D <- (u*dx)/Pe # m2/sec
  D <- (u*L)/Pe # m2/sec
  # print(D)
  
  C_i_minus_one <- (D/(dx^2))+(u/(2*dx))
  C_center <- -(2*D)/(dx^2)
  C_i_plus_one <- (D/(dx^2))-(u/(2*dx))
  
  #choose boundary conditions from the available options
  #BC_flag_L controls left BCs and and BC_flag_R controls right BCs
  # 1 constant flux
  # 2 constant concentration (water vapor content)
  # 3 constant gradient
  
  # Construct the matrix
  A <- diag(x=C_center,nrow=n,ncol=n)
  A[row(A)==col(A)-1] <- C_i_plus_one
  A[row(A)==col(A)+1] <- C_i_minus_one
  
  # Substitute appropriate boundary conditions into the corners of A
  switch(BC_flag_L,
         {#if BC_flag_L = 1 (Robin)
           A[1,1] = -C_i_minus_one},
         {#if BC_flag_L = 2 (Dirichlet)
           A[1,1] = 1 
           A[1,2] = 0
           },
         {#if BC_flag_L = 3 (Neumann)
           A[1,1] = C_center + C_i_minus_one}
         )
  
  switch(BC_flag_R,
         {#if BC_flag_R = 1 (Robin)
           A[n,n] = -C_i_plus_one},
         {#if BC_flag_R = 2 (Dirichlet)
           A[n,n] = 1
           A[n,n-1] = 0
           },
         {#if BC_flag_R = 3 (Neumann)
           A[n,n] = C_center + C_i_plus_one}
         )
  
  #make a forcing vector
  bb = rep(0,n)
  #make a masking vector for the reaction terms f
  #to use with dirichlet conditions
  mask = rep(1,n)
  
  #modify bb
  switch(BC_flag_L,
         {#if BC_flag_L = 1 (Robin)
           bb[1] = Win_LBC/dx},
         {#if BC_flag_L = 2 (Dirichlet)
           bb[1] = -Win_LBC
           mask[1] = 0
           },
         {#if BC_flag_L = 3 (Neumann)
           bb[1] = Win_LBC*dx/D*C_i_minus_one }
  )
  
  switch(BC_flag_R,
         {bb[n] = Win_RBC/dx},
         {bb[n] = -Win_RBC
          mask[n] = 0
         },
         {bb[n] = -Win_RBC*dx/D*C_i_plus_one}
  )
  
  
  #make a function for the Budyko relationship
  ETfun = function(w) 
  {
    switch(budyko_flag,
           { kp*w + Eo - ( (kp*w)^omega + Eo^omega)^(1/omega)  },
           #{ Eo - kp*w*( 1 + ( Eo/(kp*w) )^omega )^(1/omega)  },
           {0}
    )
  }
  
  # Define function to solve
  fun <- function(w) {
    A %*% w + mask*(- kp*w + ETfun(w) ) + bb
  }
  
  #Call multiroot to solve it
  W <- multiroot(f=fun,start=rep(0.1,n), 
                 maxiter = 20000, rtol = 1e-16, atol = 1e-10, positive = T)  
  
  w.out <- W$root
  
  
  # Recalculate W, P, ET, and Runoff
  P <- kp*w.out # Precipitation (kg/m2/sec)
  ET = ETfun(w.out)
  runoff <- (P-ET) # Runoff (kg/m2/sec) 
  #recycling rate
  m <- P/runoff

  #integrate runoff from right to left (assuming advection to the right) 
  cumulative_discharge = rev(cumsum(rev(runoff*dx)))
  # can also do sum(runoff)
  print(cumulative_discharge[1])
  #should be what comes in (Win_LBC).  
   

  # Return Values
  result <- as.data.frame(cbind(x,w.out,P,ET,runoff,m,cumulative_discharge, omega ) )
  colnames(result) <- c("x","W","P","ET","Runoff","m","cumulative_discharge","Omega")
  return(result)

}

#### Weathering Functions ####
# Calculates concentrations
maher <- function(Cmax,q,m=1,Dw=0.03){
  tao <- exp(1)^2
  0.85*Cmax*(tao*m*Dw/q)/(1+tao*m*Dw/q)
}

# Calculates fluxes
maher.F <- function(C,q){
  q*C
}

#### Water Plots #### 
# Figure 5 
# First combination LBC = 1, RBC = 3 (Robin to Neumann)
sim13budy15 <- myODE(BC_flag_L = 1, BC_flag_R = 3, omega = 1.5,L = 1e7) # "Minimally Vegetated"
sim13budy17 <- myODE(BC_flag_L = 1, BC_flag_R = 3, omega = 1.72,L = 1e7) # Not used
sim13budy19 <- myODE(BC_flag_L = 1, BC_flag_R = 3, omega = 2.2,L = 1e7) # Gymnosperm
sim13budy26 <- myODE(BC_flag_L = 1, BC_flag_R = 3, omega = 2.5,L = 1e7) # Angiosperm/Modern


par(mfcol=c(3,2))
plot(sim13budy15$x/1000,sim13budy15$W,type="l",ylim=c(0,10),col="DarkRed",las=1,xlab="Distance (km)",ylab="Vapor Content (kg/m2)")
# lines(sim13budy17$x/1000,sim13budy17$W,col="Red")
lines(sim13budy19$x/1000,sim13budy19$W,col="Darkgreen")
lines(sim13budy26$x/1000,sim13budy26$W,col="Blue")
text(x=200,y=9,labels="A",cex=2)

#unitconversion to m/yr
uc = 60*60*24*365/1000

plot(sim13budy15$x/1000,sim13budy15$P*uc,type="l",ylim=c(0,2),col="DarkRed",las=1,xlab="Distance (km)",ylab="P (m/yr)")
# lines(sim13budy17$x/1000,sim13budy17$P*uc,col="Red")
lines(sim13budy19$x/1000,sim13budy19$P*uc,col="Darkgreen")
lines(sim13budy26$x/1000,sim13budy26$P*uc,col="Blue") 
text(x=200,y=1.75,labels="B",cex=2)

plot(sim13budy15$x/1000,sim13budy15$ET*uc,type="l",ylim=c(0,1),col="DarkRed",las=1,xlab="Distance (km)",ylab="ET (m/yr)")
# lines(sim13budy17$x/1000,sim13budy17$ET*uc,col="Red")
lines(sim13budy19$x/1000,sim13budy19$ET*uc,col="Darkgreen")
lines(sim13budy26$x/1000,sim13budy26$ET*uc,col="Blue") 
text(x=200,y=0.9,labels="C",cex=2)

plot(sim13budy15$x/1000,sim13budy15$Runoff*uc,type="l",ylim=c(0,1),col="DarkRed",las=1,xlab="Distance (km)",ylab="Runoff (m/yr)")
# lines(sim13budy17$x/1000,sim13budy17$Runoff*uc,col="Red")
lines(sim13budy19$x/1000,sim13budy19$Runoff*uc,col="Darkgreen")
lines(sim13budy26$x/1000,sim13budy26$Runoff*uc,col="Blue") 
text(x=200,y=0.9,labels="D",cex=2)

plot(sim13budy15$x/1000,sim13budy15$m,type="l",ylim=c(0,7),lty=2,col="DarkRed",las=1,xlab="Distance (km)",ylab="Recycling Ratio (P/Runoff)")
# lines(sim13budy17$x/1000,sim13budy17$m,col="Red")
lines(sim13budy19$x/1000,sim13budy19$m,col="Darkgreen")
lines(sim13budy26$x/1000,sim13budy26$m,col="Blue") 
abline(h=1,lty=1,col="DarkRed") 
text(x=200,y=6.5,labels="E",cex=2)


Eo = 3.17e-5
plot(Eo/sim13budy15$P,sim13budy15$ET/sim13budy15$P,xlim=c(0,5),ylim=c(0,1.05),col="DarkRed",type="l",las=1,xlab="Ep/P",ylab="ET/P",lwd=2)
# lines(Eo/sim13budy17$P,sim13budy17$ET/sim13budy17$P,col="Red",lwd=2)
lines(Eo/sim13budy19$P,sim13budy19$ET/sim13budy19$P,col="DarkGreen",lwd=2)
lines(Eo/sim13budy26$P,sim13budy26$ET/sim13budy26$P,col="Blue",lwd=2) 
lines(c(1,10),c(1,1),lwd=2)
lines(c(-0.2,1),c(-.2,1),lwd=2)

petP = seq(0,10,0.1)
lines.15 = 1 + petP - (1+petP^1.5)^(1/1.5)
# lines.17 = 1 + petP - (1+petP^1.719)^(1/1.719)
lines.19 = 1 + petP - (1+petP^2.2)^(1/2.2)
lines.26 = 1 + petP - (1+petP^2.5)^(1/2.5)
lines(petP,lines.15,col="DarkRed",lty=2)
# lines(petP,lines.17,col="Red",lty=2)
lines(petP,lines.19,col="DarkGreen",lty=2)
lines(petP,lines.26,col="Blue",lty=2)
text(x=0.1,y=1,labels="F",cex=2)

#### Calculate weathering parameters along transect ####
# 1) Upload soil CO2-precipitation relationship
Cotton12=read.csv("Cotton 2012.csv",header=T,skip=1) # Cotton and Sheldon 2012 (GSAB)
head(Cotton12)
dim(Cotton12)

Cotton13=read.csv("Cotton 2013.csv",header=T,skip=1) # Cotton et al., 2013 (Chemical Geology)
head(Cotton13)
dim(Cotton13)

# 2) Calculate MAP-Soil pCO2 Relationships
C13lm <- lm(Cotton13$CO2ppmv~Cotton13$MAP)
C12lm <- lm(Cotton12$CO2ppmv~Cotton12$MAP)

# Combined Relationship
TC <- rbind(Cotton13,Cotton12)
TClm <- lm(TC$CO2ppmv~TC$MAP-1)
summary(TClm)

A.P.CO2scaling <- TClm$coefficients
G.P.CO2scaling <- TClm$coefficients

# 3) Plot these relationships
par(mar=c(5,6,1,1),mfrow=c(1,1))
plot(x=Cotton13$MAP,y=Cotton13$CO2ppmv,bty="n",pch=21,bg="gray",las=1,
     ylab=NA,
     xlab="Mean Annual Precipitation (mm/year)",
     xlim=c(0,1000),ylim=c(0,12000))
mtext(2,line=3.5,text=expression(paste("Soil pCO"[2]," (ppm)")))
points(x=Cotton12$MAP,y=Cotton12$CO2ppmv,pch=21,bg="dodgerblue3")
abline(TClm,lwd=2)

# 4) Calculate Cmax scaling
# Cmax.H <- 250 # Minimum Cmax as in a minimally-vegetated world
CO2.pC <- 1700 # ppm (estimated from Foster et al. Nat Comms 2017)
# Cmax.scaling <- (460-Cmax.H)/(5000)

# Reads in data from GWB simulations--NOT USED
CO2.Cmax <- read.csv("CO2-Cmax2.csv",header=T)
CO2.x <- CO2.Cmax$CO2
Cmax.y <- CO2.Cmax$Cmax
CO2.Cmax.f <- lm(log(Cmax.y)~log(CO2.x))

# plot(CO2.Cmax$CO2,CO2.Cmax$Cmax,log="yx")
# lines(x=CO2.Cmax$CO2,y=exp((predict(CO2.Cmax.f))))

exp(predict(CO2.Cmax.f,data.frame(CO2.x=1000)))

# 5) Create plots for Soil CO2, Cmax, and [HCO3-] along profile

C1.5.m <- maher(Cmax=rep(exp(predict(CO2.Cmax.f,data.frame(CO2.x=CO2.pC))),length(sim13budy15$Runoff)),
                q=sim13budy15$Runoff*uc,
                m=1)

C1.5.m.hydro <- maher(Cmax=exp(predict(CO2.Cmax.f,data.frame(CO2.x=CO2.pC))),
                    q=sim13budy15$Runoff*uc,m=sim13budy15$m,Dw=0.03)
C1.9.m.hydro <- maher(Cmax=exp(predict(CO2.Cmax.f,data.frame(CO2.x=CO2.pC))),
                q=sim13budy19$Runoff*uc,m=sim13budy19$m,Dw=0.03)
C2.6.m.hydro <- maher(Cmax=exp(predict(CO2.Cmax.f,data.frame(CO2.x=CO2.pC))),
                q=sim13budy26$Runoff*uc,m=sim13budy26$m,Dw=0.03)


C1.5.m.rec <- maher(Cmax=exp(predict(CO2.Cmax.f,data.frame(CO2.x=(sim13budy15$P*uc*1000*G.P.CO2scaling)))),
                q=sim13budy15$Runoff*uc,m=sim13budy15$m,Dw=0.03)
# C1.7.m <- maher(Cmax=exp(predict(CO2.Cmax.f,data.frame(CO2.x=(sim13budy17$P*uc*1000*G.P.CO2scaling)))),
#                 q=sim13budy17$Runoff*uc,m=sim13budy17$m,Dw=0.01)
C1.9.m <- maher(Cmax=exp(predict(CO2.Cmax.f,data.frame(CO2.x=(sim13budy19$P*uc*1000*G.P.CO2scaling)))),
                q=sim13budy19$Runoff*uc,m=sim13budy19$m,Dw=0.02652988)
C2.6.m <- maher(Cmax=exp(predict(CO2.Cmax.f,data.frame(CO2.x=(sim13budy26$P*uc*1000*G.P.CO2scaling)))),
                q=sim13budy26$Runoff*uc,m=sim13budy26$m,Dw=0.02576665)

# make last 3 plots of Figure 5
par(mfcol=c(3,2))
plot(x=sim13budy15$x/1000,y=rep(CO2.pC,length(sim13budy15$x)),
     type="l",ylim=c(0,8000),col="Black",las=1,xlab="Distance (km)",
     ylab=expression(paste("Soil CO"[2]," (ppm)")))
text(x=200,y=1100,labels="A",cex=2)

plot(x=sim13budy15$x/1000,y=rep(exp(predict(CO2.Cmax.f,data.frame(CO2.x=CO2.pC))),length(sim13budy15$Runoff)),
     type="l",ylim=c(0,1500),col="Black",las=1,xlab="Distance (km)",
     ylab=expression(paste("[HCO"[3],"-] C"[max]," (",mu,"mol/L)")))
text(x=200,y=1100,labels="B",cex=2)

plot(x=sim13budy15$x/1000,y=C1.5.m,#sim13budy15$P*uc*1000*P.CO2scaling,
     type="l",ylim=c(0,1200),col="DarkRed",las=1,xlab="Distance (km)",
     ylab=expression(paste("[HCO"[3],"-] (",mu,"mol/L)")))
# lines(sim13budy17$x/1000,C1.7.m,col="Red")
lines(sim13budy15$x/1000,C1.5.m.hydro,col="Darkred",lty=2)
lines(sim13budy19$x/1000,C1.9.m.hydro,col="DarkGreen")
lines(sim13budy26$x/1000,C2.6.m.hydro,col="Blue") 
text(x=200,y=1100,labels="C",cex=2)

plot(x=sim13budy15$x/1000,y=rep(CO2.pC,length(sim13budy15$x)),
     type="l",ylim=c(0,8000),col="DarkRed",las=1,xlab="Distance (km)",
     ylab=expression(paste("Soil CO"[2]," (ppm)")))
# lines(sim13budy17$x/1000,sim13budy17$P*uc*1000*G.P.CO2scaling,col="Red")
lines(sim13budy15$x/1000,sim13budy15$P*uc*1000*A.P.CO2scaling,col="darkred",lty=2)
lines(sim13budy19$x/1000,sim13budy19$P*uc*1000*A.P.CO2scaling,col="DarkGreen")
lines(sim13budy26$x/1000,sim13budy26$P*uc*1000*A.P.CO2scaling,col="Blue")
text(x=200,y=7500,labels="D",cex=2)

plot(x=sim13budy15$x/1000,y=rep(exp(predict(CO2.Cmax.f,data.frame(CO2.x=CO2.pC))),length(sim13budy15$Runoff)),
     type="l",ylim=c(0,1500),col="DarkRed",las=1,xlab="Distance (km)",
     ylab=expression(paste("[HCO"[3],"-] C"[max]," (",mu,"mol/L)")))
# lines(sim13budy17$x/1000,exp(predict(CO2.Cmax.f,data.frame(CO2.x=(sim13budy17$P*uc*1000*G.P.CO2scaling)))),
#       col="Red")
lines(sim13budy15$x/1000,exp(predict(CO2.Cmax.f,data.frame(CO2.x=(sim13budy15$P*uc*1000*G.P.CO2scaling)))),
      col="darkred",lty=2)
lines(sim13budy19$x/1000,exp(predict(CO2.Cmax.f,data.frame(CO2.x=(sim13budy19$P*uc*1000*G.P.CO2scaling)))),
      col="DarkGreen")
lines(sim13budy26$x/1000,exp(predict(CO2.Cmax.f,data.frame(CO2.x=(sim13budy26$P*uc*1000*G.P.CO2scaling)))),
      col="Blue") 
text(x=200,y=1450,labels="E",cex=2)

plot(x=sim13budy15$x/1000,y=C1.5.m,#sim13budy15$P*uc*1000*P.CO2scaling,
     type="l",ylim=c(0,1200),col="DarkRed",las=1,xlab="Distance (km)",
     ylab=expression(paste("[HCO"[3],"-] (",mu,"mol/L)")))
# lines(sim13budy17$x/1000,C1.7.m,col="Red")
lines(sim13budy15$x/1000,C1.5.m.rec,col="Darkred",lty=2)
lines(sim13budy19$x/1000,C1.9.m,col="DarkGreen")
lines(sim13budy26$x/1000,C2.6.m,col="Blue") 
text(x=200,y=1100,labels="F",cex=2)

# Calculate some things related to text
mean(sim13budy15$Runoff)
mean(sim13budy19$Runoff)/mean(sim13budy15$Runoff)
mean(sim13budy26$Runoff)/mean(sim13budy15$Runoff)


# weighted mean recycling factors for normalized flux figure (fig 7)
weighted.mean(sim13budy15$m,sim13budy15$Runoff)
weighted.mean(sim13budy19$m,sim13budy19$Runoff)
weighted.mean(sim13budy26$m,sim13budy26$Runoff)

# Conc for noramlized contour figure (fig 6)
min.conc = weighted.mean(C1.5.m,sim13budy15$Runoff)
min.conc # point A
min.rec.conc = weighted.mean(C1.5.m.rec,sim13budy15$Runoff)
min.rec.conc
min.rec.conc/min.conc # point A'';
gymn.conc = weighted.mean(C1.9.m,sim13budy19$Runoff)
gymn.conc
gymn.conc/min.conc # point B''
angio.conc = weighted.mean(C2.6.m,sim13budy26$Runoff)
angio.conc 
angio.conc/min.conc # point C''
min.rec.conc = weighted.mean(C1.5.m.hydro,sim13budy15$Runoff)
min.rec.conc
min.rec.conc/min.conc # point A';
gymn.conc = weighted.mean(C1.9.m.hydro,sim13budy19$Runoff)
gymn.conc
gymn.conc/min.conc # point B'
angio.conc = weighted.mean(C2.6.m.hydro,sim13budy26$Runoff)
angio.conc 
angio.conc/min.conc # point C'



# weighted mean Cmax concentrations
min.Cmax = weighted.mean(exp(predict(CO2.Cmax.f,data.frame(CO2.x=(sim13budy15$P*uc*1000*G.P.CO2scaling)))),sim13budy15$Runoff)
Gym.Cmax = weighted.mean(exp(predict(CO2.Cmax.f,data.frame(CO2.x=(sim13budy19$P*uc*1000*G.P.CO2scaling)))),sim13budy19$Runoff)
Mod.Cmax = weighted.mean(exp(predict(CO2.Cmax.f,data.frame(CO2.x=(sim13budy26$P*uc*1000*G.P.CO2scaling)))),sim13budy26$Runoff)

min.Cmax/exp(predict(CO2.Cmax.f,data.frame(CO2.x=CO2.pC)))
Gym.Cmax/exp(predict(CO2.Cmax.f,data.frame(CO2.x=CO2.pC)))
Mod.Cmax/exp(predict(CO2.Cmax.f,data.frame(CO2.x=CO2.pC)))


# Calculate effective Dw values from Cmax increase
ModernAvg.Dw = 0.03
Minrec.Dw = min.Cmax/exp(predict(CO2.Cmax.f,data.frame(CO2.x=CO2.pC)))
Gym.Dw = ModernAvg.Dw *(min.Cmax/Gym.Cmax)
Mod.Dw = ModernAvg.Dw *(min.Cmax/Mod.Cmax)

# New Contour Plot ####
# By Dan May 23 2017
Dw_Range = c(seq(0.003,0.01,0.0001),seq(0.011,0.1,0.001),seq(0.11,0.5,0.01))
mm = 61.0168 # mol/g
tau = exp(2)
omega_Range <- seq(1,3.5,0.025)

BareCumFlux = matrix(NA,ncol=length(omega_Range),nrow=length(Dw_Range))
PlantCumFlux = matrix(NA,ncol=length(omega_Range),nrow=length(Dw_Range))
PlantThermoFlux = matrix(NA,ncol=length(omega_Range),nrow=length(Dw_Range))

BareMeanConc= matrix(NA,ncol=length(omega_Range),nrow=length(Dw_Range))
PlantMeanConc = matrix(NA,ncol=length(omega_Range),nrow=length(Dw_Range))
PlantThermoConc = matrix(NA,ncol=length(omega_Range),nrow=length(Dw_Range))

for (i in 1:length(omega_Range)){
  fluxes <- myODE(BC_flag_L = 1, BC_flag_R = 3, omega = omega_Range[i] ,L = 1e7)
  
  for (j in 1:length(Dw_Range)){
    
    
    BareConc <- maher(Cmax=rep(exp(predict(CO2.Cmax.f,data.frame(CO2.x=CO2.pC))),length(sim13budy15$Runoff)),
                      q=fluxes$Runoff*uc,
                      m=1,Dw=Dw_Range[j])
    PlantConc <- maher(Cmax=exp(predict(CO2.Cmax.f,data.frame(CO2.x=(fluxes$P*uc*1000*G.P.CO2scaling)))),
                       q=fluxes$Runoff*uc,
                       m=fluxes$m,Dw=Dw_Range[j])
    
    PlantThermo <- maher(Cmax=exp(predict(CO2.Cmax.f,data.frame(CO2.x=(fluxes$P*uc*1000*G.P.CO2scaling)))),
                         q=fluxes$Runoff*uc,
                         m=1,Dw=Dw_Range[j])

    BareMeanConc[j,i] <- weighted.mean(BareConc,fluxes$Runoff)
    PlantMeanConc[j,i] <- weighted.mean(PlantConc,fluxes$Runoff)
    PlantThermoConc[j,i] <- weighted.mean(PlantThermo,fluxes$Runoff)
    
    BareCumFlux[j,i] = sum(fluxes$Runoff*uc*BareConc)
    PlantCumFlux[j,i] = sum(fluxes$Runoff*uc*PlantConc)
    PlantThermoFlux[j,i] = sum(fluxes$Runoff*uc*PlantThermo)
  }
}

# make plot of Contours
Dw_RangeLog = log10(Dw_Range)
Omega_RangeLog = log10(omega_Range)
colorwes=wes_palette(name="GrandBudapest2",n=15,type="continuous")

NormalizeConc = BareMeanConc[91,21]
# NormalizeFlux = BareCumFlux[11,11]

dev.off()

contourlevels = c(0,3.5)
levels.set = seq(0,3.5,0.25)

filled.contour(y=Omega_RangeLog,x=Dw_RangeLog,z=BareMeanConc/NormalizeConc,levels=levels.set, color=function(x)terrain.colors(x),
               xlab="Dw (m/yr)",ylab="Omega (w)",zlim=contourlevels,xlim=c(log10(c(0.003,0.5))),
               plot.axes={box(which="plot",lwd=3); abline(h=log10(c(1.5,2.2,2.6)),lwd=2,col="black",lty=1); abline(v=log10(c(0.01,0.03,0.3)),lwd=2,col="black",lty=1);axis(1,at=log10(c(0.003,0.005,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5)),lwd=3,labels = c(0.003,0.005,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5));
                 axis(2,at=log10(c(1,1.5,2,2.5,3,3.5,4,5)),labels=c(1,1.5,2,2.5,3,3.5,4,5),lwd=3);
                 contour(y=Omega_RangeLog,x=Dw_RangeLog,z=BareMeanConc/NormalizeConc,method="edge",add=T,axes=T,lty=4,lwd=1.5,labcex=1.25,levels=levels.set,zlim=contourlevels); points(x=log10(0.01),y=250,pch=15,cex=3); title("No Recycling pCO2 = 1700ppm");
                 points(log10(0.03),log10(1.5),pch=15,cex=2.5)})


filled.contour(y=Omega_RangeLog,x=Dw_RangeLog,z=PlantMeanConc/NormalizeConc,levels=levels.set,color=function(x)terrain.colors(x),
               xlab="Dw (m/yr)",ylab="Omega (w)",zlim=contourlevels,#xlim=c(log10(c(0.003,0.5))),
               plot.axes={box(which="plot",lwd=3); abline(h=log10(c(1.5,2.2,2.6)),lwd=2,col="black",lty=1); abline(v=log10(c(0.01,0.03,0.3)),lwd=2,col="black",lty=1);axis(1,at=log10(c(0.003,0.005,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5)),lwd=3,labels = c(0.003,0.005,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5)); 
                 axis(2,at=log10(c(1,1.5,2,2.5,3,3.5,4,5)),labels=c(1,1.5,2,2.5,3,3.5,4,5),lwd=3);
                 contour(y=Omega_RangeLog,x=Dw_RangeLog,z=PlantMeanConc/NormalizeConc,method="edge",add=T,axes=T,lty=4,lwd=1.5,labcex=1.25,levels=levels.set,zlim=contourlevels); points(x=log10(0.01),y=250,pch=15,cex=3); title("With Plant Recycling and Soil CO2-CMax Feedbacks");
                 points(log10(c(0.03,Gym.Dw,Mod.Dw)),log10(c(1.5,2.2,2.6)),pch=15,cex=2.5);
                 lines(log10(c(0.03,Gym.Dw,Mod.Dw)),log10(c(1.5,2.2,2.6)),lwd=5,lty=3)})
                 # points(log10(c(0.03,Gym.Dw,Mod.Dw)*10),log10(c(1.5,1.9,2.6)),pch=15,cex=2.5,col="grey");
                 # lines(log10(c(0.03,Gym.Dw,Mod.Dw)*10),log10(c(1.5,1.9,2.6)),lwd=5,lty=3,col="grey")
                 # points(log10(c(0.03,Gym.Dw,Mod.Dw)/3),log10(c(1.5,1.9,2.6)),pch=15,cex=2.5,col="grey");
                 # lines(log10(c(0.03,Gym.Dw,Mod.Dw)/3),log10(c(1.5,1.9,2.6)),lwd=5,lty=3,col="grey")


filled.contour(y=Omega_RangeLog,x=Dw_RangeLog,z=PlantThermoConc/NormalizeConc,levels=levels.set,color=function(x)terrain.colors(x),
               xlab="Dw (m/yr)",ylab="Omega (w)",zlim=contourlevels,#xlim=c(log10(c(0.003,0.5))),
               plot.axes={box(which="plot",lwd=3); abline(h=log10(c(1.5,2.2,2.6)),lwd=2,col="black",lty=1); abline(v=log10(c(0.01,0.03,0.3)),lwd=2,col="black",lty=1);axis(1,at=log10(c(0.003,0.005,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5)),lwd=3,labels = c(0.003,0.005,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5)); 
                 axis(2,at=log10(c(1,1.5,2,2.5,3,3.5,4,5)),labels=c(1,1.5,2,2.5,3,3.5,4,5),lwd=3);
                 contour(y=Omega_RangeLog,x=Dw_RangeLog,z=PlantThermoConc/NormalizeConc,method="edge",add=T,axes=T,lty=4,lwd=1.5,labcex=1.25,levels=levels.set,zlim=contourlevels); points(x=log10(0.01),y=250,pch=15,cex=3); title("Without Plant Recycling but with Soil CO2-CMax Feedbacks");
                 points(log10(c(0.03,0.03,0.03)),log10(c(1.5,2.2,2.6)),pch=15,cex=2.5);
                 lines(log10(c(0.03,0.03,0.03)),log10(c(1.5,2.2,2.6)),lwd=5,lty=3)})



# points(log10(c(0.03,Gym.Dw,Mod.Dw)*10),log10(c(1.5,1.9,2.6)),pch=15,cex=2.5,col="grey");
# lines(log10(c(0.03,Gym.Dw,Mod.Dw)*10),log10(c(1.5,1.9,2.6)),lwd=5,lty=3,col="grey")
# points(log10(c(0.03,Gym.Dw,Mod.Dw)/3),log10(c(1.5,1.9,2.6)),pch=15,cex=2.5,col="grey");
# lines(log10(c(0.03,Gym.Dw,Mod.Dw)/3),log10(c(1.5,1.9,2.6)),lwd=5,lty=3,col="grey")




# filled.contour(y=Omega_RangeLog,x=Dw_RangeLog,z=PlantCumFlux/NormalizeFlux,nlevels=contourlevels[2]*2, color=function(x)terrain.colors(x),
#                xlab="Dw (m/yr)",ylab="Omega (w)",zlim=contourlevels,
#                plot.axes={box(which="plot",lwd=3); abline(h=log10(c(1.529,1.719,1.913,2.6)),lwd=3,col="black",lty=1); abline(v=log10(c(0.01,0.03,0.3)),lwd=3,col="black",lty=1);axis(1,at=log10(c(0.005,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5)),lwd=3,labels = c(0.005,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5)); 
#                  axis(2,at=log10(c(1,2,3,4,5)),labels=c(1,2,3,4,5),lwd=3);
#                  contour(y=Omega_RangeLog,x=Dw_RangeLog,z=PlantCumFlux/NormalizeFlux,method="edge",add=T,axes=T,lty=4,lwd=1.5,labcex=1.5,nlevels=contourlevels[2]*2,zlim=contourlevels); points(x=log10(0.01),y=250,pch=15,cex=3); title("Vegetated World (Cumulative Flux)")})



sim13budy15 <- myODE(BC_flag_L = 1, BC_flag_R = 1, omega = 1.5,L = 1e7, Win_RBC = -34.97106) # "Minimally Vegetated"
sim13budy17 <- myODE(BC_flag_L = 1, BC_flag_R = 1, omega = 1.72,L = 1e7, Win_RBC = -34.97106) # Not used
sim13budy19 <- myODE(BC_flag_L = 1, BC_flag_R = 1, omega = 2.2,L = 1e7, Win_RBC = -34.97106) # Gymnosperm
sim13budy26 <- myODE(BC_flag_L = 1, BC_flag_R = 1, omega = 2.5,L = 1e7, Win_RBC = -34.97106) # Angiosperm/Modern

par(mfcol=c(3,2))
plot(sim13budy15$x/1000,sim13budy15$W,type="l",ylim=c(0,10),col="DarkRed",las=1,xlab="Distance (km)",ylab="Vapor Content (kg/m2)")
# lines(sim13budy17$x/1000,sim13budy17$W,col="Red")
lines(sim13budy19$x/1000,sim13budy19$W,col="Darkgreen")
lines(sim13budy26$x/1000,sim13budy26$W,col="Blue")
text(x=200,y=9,labels="A",cex=2)



par(mfcol=c(3,2))
plot(sim13budy15$x/1000,sim13budy15$W,type="l",ylim=c(0,10),col="DarkRed",las=1,xlab="Distance (km)",ylab="Vapor Content (kg/m2)")
# lines(sim13budy17$x/1000,sim13budy17$W,col="Red")
lines(sim13budy19$x/1000,sim13budy19$W,col="Darkgreen")
lines(sim13budy26$x/1000,sim13budy26$W,col="Blue")
text(x=200,y=9,labels="A",cex=2)

#unitconversion to m/yr
uc = 60*60*24*365/1000

plot(sim13budy15$x/1000,sim13budy15$P*uc,type="l",ylim=c(0,1.5),col="DarkRed",las=1,xlab="Distance (km)",ylab="P (m/yr)")
# lines(sim13budy17$x/1000,sim13budy17$P*uc,col="Red")
lines(sim13budy19$x/1000,sim13budy19$P*uc,col="Darkgreen")
lines(sim13budy26$x/1000,sim13budy26$P*uc,col="Blue") 
text(x=200,y=1.75,labels="B",cex=2)

plot(sim13budy15$x/1000,sim13budy15$ET*uc,type="l",ylim=c(0,1),col="DarkRed",las=1,xlab="Distance (km)",ylab="ET (m/yr)")
# lines(sim13budy17$x/1000,sim13budy17$ET*uc,col="Red")
lines(sim13budy19$x/1000,sim13budy19$ET*uc,col="Darkgreen")
lines(sim13budy26$x/1000,sim13budy26$ET*uc,col="Blue") 
text(x=200,y=0.9,labels="C",cex=2)

plot(sim13budy15$x/1000,sim13budy15$Runoff*uc,type="l",ylim=c(0,0.5),col="DarkRed",las=1,xlab="Distance (km)",ylab="Runoff (m/yr)")
# lines(sim13budy17$x/1000,sim13budy17$Runoff*uc,col="Red")
lines(sim13budy19$x/1000,sim13budy19$Runoff*uc,col="Darkgreen")
lines(sim13budy26$x/1000,sim13budy26$Runoff*uc,col="Blue") 
text(x=200,y=0.9,labels="D",cex=2)


#### Calculate weathering parameters along transect ####
# 1) Upload soil CO2-precipitation relationship
Cotton12=read.csv("Cotton 2012.csv",header=T,skip=1) # Cotton and Sheldon 2012 (GSAB)
head(Cotton12)
dim(Cotton12)

Cotton13=read.csv("Cotton 2013.csv",header=T,skip=1) # Cotton et al., 2013 (Chemical Geology)
head(Cotton13)
dim(Cotton13)

# 2) Calculate MAP-Soil pCO2 Relationships
C13lm <- lm(Cotton13$CO2ppmv~Cotton13$MAP)
C12lm <- lm(Cotton12$CO2ppmv~Cotton12$MAP)

# Combined Relationship
TC <- rbind(Cotton13,Cotton12)
TClm <- lm(TC$CO2ppmv~TC$MAP-1)
summary(TClm)

A.P.CO2scaling <- TClm$coefficients
G.P.CO2scaling <- TClm$coefficients

# 3) Plot these relationships
# par(mar=c(5,6,1,1),mfrow=c(1,1))
# plot(x=Cotton13$MAP,y=Cotton13$CO2ppmv,bty="n",pch=21,bg="gray",las=1,
#      ylab=NA,
#      xlab="Mean Annual Precipitation (mm/year)",
#      xlim=c(0,1000),ylim=c(0,12000))
# mtext(2,line=3.5,text=expression(paste("Soil pCO"[2]," (ppm)")))
# points(x=Cotton12$MAP,y=Cotton12$CO2ppmv,pch=21,bg="dodgerblue3")
# abline(TClm,lwd=2)

# 4) Calculate Cmax scaling
# Cmax.H <- 250 # Minimum Cmax as in a minimally-vegetated world
CO2.pC <- 1700 # ppm (estimated from Foster et al. Nat Comms 2017)
# Cmax.scaling <- (460-Cmax.H)/(5000)

# Reads in data from GWB simulations--NOT USED
CO2.Cmax <- read.csv("CO2-Cmax2.csv",header=T)
CO2.x <- CO2.Cmax$CO2
Cmax.y <- CO2.Cmax$Cmax
CO2.Cmax.f <- lm(log(Cmax.y)~log(CO2.x))

# plot(CO2.Cmax$CO2,CO2.Cmax$Cmax,log="yx")
# lines(x=CO2.Cmax$CO2,y=exp((predict(CO2.Cmax.f))))

exp(predict(CO2.Cmax.f,data.frame(CO2.x=1000)))

# 5) Create plots for Soil CO2, Cmax, and [HCO3-] along profile

C1.5.m <- maher(Cmax=rep(exp(predict(CO2.Cmax.f,data.frame(CO2.x=CO2.pC))),length(sim13budy15$Runoff)),
                q=sim13budy15$Runoff*uc,
                m=1)

C1.5.m.hydro <- maher(Cmax=exp(predict(CO2.Cmax.f,data.frame(CO2.x=CO2.pC))),
                      q=sim13budy15$Runoff*uc,m=sim13budy15$m,Dw=0.3)
C1.9.m.hydro <- maher(Cmax=exp(predict(CO2.Cmax.f,data.frame(CO2.x=CO2.pC))),
                      q=sim13budy19$Runoff*uc,m=sim13budy19$m,Dw=0.3)
C2.6.m.hydro <- maher(Cmax=exp(predict(CO2.Cmax.f,data.frame(CO2.x=CO2.pC))),
                      q=sim13budy26$Runoff*uc,m=sim13budy26$m,Dw=0.3)


C1.5.m.rec <- maher(Cmax=exp(predict(CO2.Cmax.f,data.frame(CO2.x=(sim13budy15$P*uc*1000*G.P.CO2scaling)))),
                    q=sim13budy15$Runoff*uc,m=sim13budy15$m,Dw=0.3)
# C1.7.m <- maher(Cmax=exp(predict(CO2.Cmax.f,data.frame(CO2.x=(sim13budy17$P*uc*1000*G.P.CO2scaling)))),
#                 q=sim13budy17$Runoff*uc,m=sim13budy17$m,Dw=0.01)
C1.9.m <- maher(Cmax=exp(predict(CO2.Cmax.f,data.frame(CO2.x=(sim13budy19$P*uc*1000*G.P.CO2scaling)))),
                q=sim13budy19$Runoff*uc,m=sim13budy19$m,Dw=0.2652988)
C2.6.m <- maher(Cmax=exp(predict(CO2.Cmax.f,data.frame(CO2.x=(sim13budy26$P*uc*1000*G.P.CO2scaling)))),
                q=sim13budy26$Runoff*uc,m=sim13budy26$m,Dw=0.2576665)

plot(x=sim13budy15$x/1000,y=rep(CO2.pC,length(sim13budy15$x)),
     type="l",ylim=c(0,8000),col="DarkRed",las=1,xlab="Distance (km)",
     ylab=expression(paste("Soil CO"[2]," (ppm)")))
# lines(sim13budy17$x/1000,sim13budy17$P*uc*1000*G.P.CO2scaling,col="Red")
lines(sim13budy15$x/1000,sim13budy15$P*uc*1000*A.P.CO2scaling,col="darkred",lty=2)
lines(sim13budy19$x/1000,sim13budy19$P*uc*1000*A.P.CO2scaling,col="DarkGreen")
lines(sim13budy26$x/1000,sim13budy26$P*uc*1000*A.P.CO2scaling,col="Blue")
text(x=200,y=7500,labels="D",cex=2)


plot(x=sim13budy15$x/1000,y=C1.5.m,#sim13budy15$P*uc*1000*P.CO2scaling,
     type="l",ylim=c(0,1500),col="DarkRed",las=1,xlab="Distance (km)",
     ylab=expression(paste("[HCO"[3],"-] (",mu,"mol/L)")))
# lines(sim13budy17$x/1000,C1.7.m,col="Red")
lines(sim13budy15$x/1000,C1.5.m.rec,col="Darkred",lty=2)
lines(sim13budy19$x/1000,C1.9.m,col="DarkGreen")
lines(sim13budy26$x/1000,C2.6.m,col="Blue") 
text(x=200,y=1100,labels="F",cex=2)