cmip6_warming_levels

Extract year of passing warming levels and co2 concentration at the time in the cmip6 archive.

# Load packages
library(ncdf4)
library(fields)

## Loading required package: spam

## Loading required package: dotCall64

## Loading required package: grid

## Spam version 2.7-0 (2021-06-25) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction 
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.

## 
## Attaching package: 'spam'

## The following objects are masked from 'package:base':
## 
##     backsolve, forwardsolve

## Loading required package: viridis

## Loading required package: viridisLite

## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'tibble'

## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'pillar'

## See https://github.com/NCAR/Fields for
##  an extensive vignette, other supplements and source code

# r = realization index
# i = initialization index
# p = physics index
# f = forcing index

Helper functions

extractModelnames = function(filenames, startstring, stopstring){
  # Function to extract model names from a vector of cmip filenames
  # that have similar startstring and stopstring

  fs = lapply(filenames,regexpr, pattern = startstring)
  fe = lapply(filenames,regexpr, pattern = stopstring)
  out = rep(NA, length(filenames))

  for(i in 1:length(filenames)){
    out[i] = substr(filenames[i], attr(fs[i][[1]], 'match.length')+1, fe[[i]][1]-1)
  }
  out
}

# rollmean to find the year of passing various thresholds, then index back in to the ssp co2 concentrations
library(zoo)

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

findThres = function(time, x, f = 1/5, thres = 2, y.ix = 1:30){
  # smooth a timeseries and find the point where
  # it crosses a threshold.
  ysmooth = lowess(x = x, f = f)
  ysmooth.start = mean(ysmooth$y[y.ix])
  ysmooth.change = ysmooth$y - ysmooth.start
  out = time[ min(which(ysmooth.change > thres))]
  out
}

findThresRolling = function(time, y, conc, k = 20, thres = 2, bp = 1861:1880){
  # smooth a timeseries with rollmean and find the point where
  # it crosses a threshold.
  # Baseline is calculated using only the years that are in the baseline period
  # - if not all those years are there, the mean is calculated using a shorter
  # baseline period.
  
  require(zoo)
  
  # find the index with the start of the baseline period
  bp.ix = which(round(time) %in% bp)
  bp.mean = mean(y[bp.ix], na.rm = TRUE)
  
  ysmooth = rollmean(y, k = k, na.pad = TRUE)
  #NonNAindex = which(!is.na(ysmooth))
  #firstNonNA = min(NonNAindex)
  #ysmooth.change = ysmooth - ysmooth[firstNonNA]
  ysmooth.change = ysmooth - bp.mean
  y.change = y - bp.mean
  #out = time[ min(which(ysmooth.change > thres))]
  thresyear = tryCatch(time[ min(which(ysmooth.change > thres))],
                     error=function(err) NA)
  
  thresconc = tryCatch(conc[ min(which(ysmooth.change > thres))],
                       error=function(err) NA)
  
  return(list(thresyear=thresyear, thresconc = thresconc, ysmooth = ysmooth,
              ysmooth.change = ysmooth.change,
              y.change = y.change))
}

# data directories
fn.ssp585 = dir('/data/users/hadaw/cmip6/areaavg/tas/', pattern = 'ssp585_r1i1p1')
fn.historical = dir('/data/users/hadaw/cmip6/areaavg/tas/', pattern = 'historical_r1i1p1')

# model name lists
mods.ssp585 = extractModelnames(fn.ssp585, startstring = 'tas_Amon_', stopstring = "_ssp585")
mods.historical = extractModelnames(fn.historical, startstring = 'tas_Amon_',
  stopstring = "_historical")

# Find models common to historical and future projection
common.mods = intersect(mods.ssp585, mods.historical)

nmods = length(common.mods)

fn.ssp585.kept = fn.ssp585[match(common.mods, mods.ssp585)]
fn.historical.kept = fn.historical[ match(common.mods, mods.historical)]

# useful
# https://stackoverflow.com/questions/17215789/extract-a-substring-in-r-according-to-a-pattern
# start with strsplit, pmatch, charmatch

Load and process historical model data

Uses a time series of the last 165 years of the historical runs.

# Historical years
hyrs = 1850:2014
nyr = length(hyrs)

#variable matrix (turn this into a function?
hmat = matrix(nrow = length(fn.historical.kept), ncol = nyr)

for(i in 1:length(fn.historical.kept)){
# Open and extract data from the netcdf file
  fn = fn.historical.kept[i]
  fnpath = paste0("/data/users/hadaw/cmip6/areaavg/tas/", fn)

  nc = nc_open(fnpath)

  v = ncvar_get(nc, 'tas')
  nc_close(nc)
  
  vtail = tail(v, nyr)
  hmat[i, ] = vtail
}


par(las = 1)
matplot(hyrs, t(hmat), 
        type = 'l', lty = 'solid', main = 'raw model tas', ylab = 'tas (K)')

# Amomalize model runs to the mean of each run
hmat.means = apply(hmat, 1, mean)

# varmat.anom is the variable timeseries matrix anomaly
hmat.anom = sweep(hmat, 1, STATS = hmat.means)

Extract future projection data

Start with a single realisation from each model (r1i1p1), to keep everything the same.
Hadley models use f2 - they don’t have f1

fyrs = 2015:2099
fnyr = length(fyrs)
futuremat = matrix(nrow = length(fn.ssp585.kept), ncol = fnyr)

for(i in 1:length(fn.ssp585.kept)){
# Open and extract data from the netcdf file
  fn = fn.ssp585.kept[i]
  fnpath = paste0("/data/users/hadaw/cmip6/areaavg/tas/", fn)

  nc = nc_open(fnpath)

  v = ncvar_get(nc, 'tas')
  #print(length(v))
  nc_close(nc)
  
  vhead = head(v, fnyr)
 
  futuremat[i, ] = vhead
}


# Amomalize to the mean
# varmat.anom is the variable timeseries matrix anomaly
futuremat.changemat = sweep(futuremat, 1, STATS = futuremat[,1])
future.diff = futuremat.changemat [,ncol(futuremat.changemat )]

full_trajectory <- cbind(hmat, futuremat)
full_years <- c(hyrs,fyrs)

matplot(full_years, t(full_trajectory), type = 'l', lty = 'solid')

Anomalize to early years

anom_years <- 1881:1900
anom_index <- which(full_years%in%anom_years)
historical_mean <- apply(full_trajectory[, anom_index], 1, FUN = mean)

full_trajectory_anomaly = sweep(full_trajectory, 1, STATS = historical_mean)

par(las = 1)
matplot(full_years, t(full_trajectory_anomaly),
        type = 'l', lty = 'solid', col = 'darkgrey',
        bty = 'n',
        main = 'tas anomaly',
        ylab = 'Deg C')

## Load historical and SSP CO2

# Looks like we can get the ssp concentration data from here
# http://greenhousegases.science.unimelb.edu.au/#!/ghg?scenarioid=9&mode=downloads

co2_hist <- read.csv('../data/conc/mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_CMIP_UoM-CMIP-1-2-0_gr1-GMNHSH_0000-2014.csv')

co2_ssp119 <- read.csv('../data/conc/mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-IMAGE-ssp119-1-2-1_gr1-GMNHSH_2015-2500.csv')

co2_ssp126 <- read.csv('../data/conc/mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-IMAGE-ssp126-1-2-1_gr1-GMNHSH_2015-2500.csv')

co2_ssp245 <- read.csv('../data/conc/mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-MESSAGE-GLOBIOM-ssp245-1-2-1_gr1-GMNHSH_2015-2500.csv')

co2_ssp370 <- read.csv('../data/conc/mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-AIM-ssp370-1-2-1_gr1-GMNHSH_2015-2500.csv')

co2_ssp434 <- read.csv('../data/conc/mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-GCAM4-ssp434-1-2-1_gr1-GMNHSH_2015-2500.csv')

co2_ssp460 <- read.csv('../data/conc/mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-GCAM4-ssp460-1-2-1_gr1-GMNHSH_2015-2500.csv')

co2_ssp534 <- read.csv('../data/conc/mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-REMIND-MAGPIE-ssp534-over-1-2-1_gr1-GMNHSH_2015-2500.csv')

co2_ssp585 <- read.csv('../data/conc/mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-REMIND-MAGPIE-ssp585-1-2-1_gr1-GMNHSH_2015-2500.csv')

sort the co2 data into the correct format.

process_co2_data <- function(hist_co2, ssp_co2, years){

  # could get rid of the row names here - they just count rows from 1
  
  co2_hist_ssp <- as.data.frame(rbind(cbind(hist_co2$year, hist_co2$data_mean_global),
                         cbind(ssp_co2$year, ssp_co2$data_mean_global)))
  
  colnames(co2_hist_ssp) <- c('year','data_mean_global')
  
  out <- co2_hist_ssp[which(co2_hist_ssp$year%in%years), ]
  out
}


# could use assign() here


co2_ssp119_const <- process_co2_data(hist_co2 = co2_hist,
                           ssp_co2 = co2_ssp119,
                           years = full_years)

co2_ssp126_const <- process_co2_data(hist_co2 = co2_hist,
                           ssp_co2 = co2_ssp126,
                           years = full_years)
co2_ssp245_const <- process_co2_data(hist_co2 = co2_hist,
                           ssp_co2 = co2_ssp245,
                           years = full_years)
co2_ssp370_const <- process_co2_data(hist_co2 = co2_hist,
                           ssp_co2 = co2_ssp370,
                           years = full_years)

co2_ssp434_const <- process_co2_data(hist_co2 = co2_hist,
                           ssp_co2 = co2_ssp434,
                           years = full_years)

co2_ssp460_const <- process_co2_data(hist_co2 = co2_hist,
                           ssp_co2 = co2_ssp460,
                           years = full_years)
co2_ssp534_const <- process_co2_data(hist_co2 = co2_hist,
                           ssp_co2 = co2_ssp534,
                           years = full_years)

co2_ssp585_const <- process_co2_data(hist_co2 = co2_hist,
                           ssp_co2 = co2_ssp585,
                           years = full_years)

Check CO2 for the SSPs

#doesn't look like we have anything for ssp534

cbPal <- rev(c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00"))#, "#CC79A7"))

ssplist <- c('119', '126', '245', '370', '434', '460', '585')

co2_list <- paste0('co2_ssp', ssplist, '_const')

# run through each of the SSPs

plot(co2_ssp585_const, type = 'n', bty = 'n')

for(i in 1:length(co2_list)){
  co2 <- get(co2_list[i])
  lines(co2, col = cbPal[i], lwd = 1.5)
}

Load historical tas

getTAShist <- function(fn.historical.kept, hyrs = 1850:2014){
  
  nyr = length(hyrs)

  #variable matrix (turn this into a function?
  hmat = matrix(nrow = length(fn.historical.kept), ncol = nyr)


  for(i in 1:length(fn.historical.kept)){
  # Open and extract data from the netcdf file
    fn = fn.historical.kept[i]
    fnpath = paste0("/data/users/hadaw/cmip6/areaavg/tas/", fn)

    nc = nc_open(fnpath)

    v = ncvar_get(nc, 'tas')
    nc_close(nc)
  
    vtail = tail(v, nyr)
    hmat[i, ] = vtail
  }
  hmat
}

load tas for each SSP

getTASssp <- function(ssp, mods.historical, fyrs = 2015:2099){
  
  fn.ssp <- dir('/data/users/hadaw/cmip6/areaavg/tas/', pattern = paste0('ssp',ssp,'_r1i1p1'))
  mods.ssp = extractModelnames(fn.ssp, startstring = 'tas_Amon_', stopstring = paste0('_ssp', ssp))
  
  common.mods = intersect(mods.ssp, mods.historical)
  nmods = length(common.mods)
  fn.ssp.kept = fn.ssp[match(common.mods, mods.ssp)]
  
  fnyr = length(fyrs)
  futuremat = matrix(nrow = length(fn.ssp.kept), ncol = fnyr)
  
  for(i in 1:length(fn.ssp.kept)){
    # Open and extract data from the netcdf file
    fn = fn.ssp.kept[i]
    fnpath = paste0("/data/users/hadaw/cmip6/areaavg/tas/", fn)
    
    nc = nc_open(fnpath)
    
    v = ncvar_get(nc, 'tas')
    #print(length(v))
    nc_close(nc)
    
    vhead = head(v, fnyr)
    vl <- length(vhead)
    
    futuremat[i,1:vl ] = vhead
  }
  
  return(list(tas = futuremat, mods = common.mods, fn.ssp = fn.ssp, mods = mods.ssp))
}

Generate full trajectories of tas, splicing historical with ssp

for(i in 1:length(ssplist)){
  print(i)
  
  fyrs <- 2015:2099
  hyrs <- 1850:2014
  TASssp <- getTASssp(ssp = ssplist[i], mods.historical = mods.historical, fyrs = 2015:2099)
  
  # find a way to match what's coming out here
  
  # data directories
  #fn.ssp = TASssp$fn.ssp
  #fn.historical = dir('/data/users/hadaw/cmip6/areaavg/tas/', pattern = 'historical_r1i1p1')

# model name lists
  mods.ssp = TASssp$mods
  
  mods.historical = extractModelnames(fn.historical, startstring = 'tas_Amon_',
  stopstring = "_historical")

# Find models common to historical and future projection
  common.mods = intersect(mods.ssp, mods.historical)

  nmods = length(common.mods)

  #fn.ssp.kept = fn.ssp[match(common.mods, mods.ssp)]
  fn.historical.kept = fn.historical[ match(common.mods, mods.historical)]
  
  TAShist <- getTAShist(fn.historical.kept)

# hmat will change each time  
  full_trajectory_ssp <- cbind(TAShist, TASssp$tas)
  
  assign(paste0('full_trajectory_tas_ssp',ssplist[i]), full_trajectory_ssp)
  
  assign(paste0('models_tas_ssp',ssplist[i]), TASssp$mods)
  
}

## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7

Year of passing 1.5 and 2 degrees for SSP585

# re-write this as a function
passyears_ssp585 <- as.data.frame(matrix(NA, nrow = length(common.mods), ncol = 5))
colnames(passyears_ssp585) <- c('modname', "year_1.5","conc_1.5", "year_2", "conc_2")

passyears_ssp585[,1] <- common.mods 

for(i in 1:length(common.mods)){
  
  out1_5 <- findThresRolling(time = co2_ssp585_const$year, y = full_trajectory[i,],
                         conc = co2_ssp585_const$data_mean_global, k = 20, thres = 1.5, bp = 1861:1880)
  
  passyears_ssp585[i, "year_1.5"] <- out1_5$thresyear
  passyears_ssp585[i, "conc_1.5"] <- out1_5$thresconc
  
  
  out2 <- findThresRolling(time = co2_ssp585_const$year, y = full_trajectory[i,],
                         conc = co2_ssp585_const$data_mean_global, k = 20, thres = 2, bp = 1861:1880)
  
  passyears_ssp585[i,"year_2"] <- out2$thresyear
  passyears_ssp585[i,"conc_2"] <- out2$thresconc
}


generatePassyears <- function(mods, co2, tas, k = 20, bp = 1861:1880){
  
  # mods is a vector of model names, length the same as rows in tas
  
  stopifnot(length(mods)==nrow(tas))
  
  passyears <- as.data.frame(matrix(NA, nrow = length(mods), ncol = 9))
  
  colnames(passyears) <- c('modname', "year_1.5","conc_1.5", "year_2", "conc_2", "year_3", "conc_3", "year_4", "conc_4")

  passyears[,1] <- mods 
  
  for(i in 1:length(mods)){
    
    out1_5 <- findThresRolling(time = co2$year, y = tas[i,],
                               conc = co2$data_mean_global, k = k, thres = 1.5, bp = bp)
    
    passyears[i, "year_1.5"] <- out1_5$thresyear
    passyears[i, "conc_1.5"] <- out1_5$thresconc
    
    out2 <- findThresRolling(time = co2$year, y = tas[i,],
                             conc = co2$data_mean_global, k = k, thres = 2, bp = bp)
    
    passyears[i,"year_2"] <- out2$thresyear
    passyears[i,"conc_2"] <- out2$thresconc
    
     out3 <- findThresRolling(time = co2$year, y = tas[i,],
                              conc = co2$data_mean_global, k = k, thres = 3, bp = bp)
    
     passyears[i,"year_3"] <- out3$thresyear
     passyears[i,"conc_3"] <- out3$thresconc
    
     out4 <- findThresRolling(time = co2$year, y = tas[i,],
                              conc = co2$data_mean_global, k = k, thres = 4, bp = bp)
    
     passyears[i,"year_4"] <- out4$thresyear
     passyears[i,"conc_4"] <- out4$thresconc
    
  }
  
  passyears
  
}

test <- findThresRolling(time = co2_ssp585_const$year, y = full_trajectory_tas_ssp119[1, ], conc = co2$data_mean_global,
                         k = 20, thres = 3, bp = 1861:1880)

## Warning in min(which(ysmooth.change > thres)): no non-missing arguments to min;
## returning Inf

## Warning in min(which(ysmooth.change > thres)): no non-missing arguments to min;
## returning Inf

plotRollingThresh <- function(fTRobj){
  
  
  plot(fTRobj$ysmooth.change, type = 'l')
  lines(fTRobj$y.change, type = 'l', col = 'grey')
  
  
}


plotRollingThresh(test)

Calculate the years the SSPs pass warming levels

for(i in 1:length(co2_list)){
  
  passyears <- generatePassyears(mods = get(paste0('models_tas_ssp',ssplist[i])), co2 = get(co2_list[i]), 
                                 tas = get(paste0('full_trajectory_tas_ssp', ssplist[i])))
  
  assign(paste0('passyears_ssp',ssplist[i]), passyears) 
  
}

Plot the years the SSPs pass warming levels

for(i in 1:length(ssplist)){
  
  anom_years <- 1881:1900
  anom_index <- which(full_years%in%anom_years)
  
  full_trajectory <- get(paste0('full_trajectory_tas_ssp',ssplist[i]))
  historical_mean <- apply(full_trajectory[, anom_index], 1, FUN = mean)
  
  full_trajectory_anomaly = sweep(full_trajectory, 1, STATS = historical_mean)
  
  par(las = 1)
  matplot(full_years, t(full_trajectory_anomaly),
        type = 'l', lty = 'solid', col = 'darkgrey',
        bty = 'n',
        main = paste0('SSP',ssplist[i]),
        ylab = 'Deg C')
  
 passyears = get(paste0('passyears_ssp',ssplist[i]))

 points(passyears[,'year_1.5'],rep(1.5, length(get(paste0('models_tas_ssp', ssplist[i])))), pch = 20, col = 'skyblue2')
 points(passyears[,'year_2'], rep(2, length(get(paste0('models_tas_ssp', ssplist[i])))), pch = 20, col = 'tomato2')
 points(passyears[,'year_3'], rep(3, length(get(paste0('models_tas_ssp', ssplist[i])))), pch = 20, col = 'darkorchid')
 points(passyears[,'year_4'], rep(4, length(get(paste0('models_tas_ssp', ssplist[i])))), pch = 20, col = 'darkturquoise')
  
}

Check it works with the smoothing

Save the CMIP6 year of passing thresholds and accompanying CO2 concentration as RData and CSVs.

passyears_list <- as.list(ls(pattern = 'passyears_ssp'))


for(i in passyears_list){
  
  passyears = get(i)
  save(passyears, file = paste0(i,'.RData'))
  
  write.table(passyears, file = paste0(i,'.csv'))
  
}

Load CMIP5 passyears

rcp85_passyears_2 <- read.table('~/at2degreesgit/RCP85_combined_2deg.txt')

rcp60_passyears_2 <- read.table('~/at2degreesgit/RCP60_combined_2deg.txt')

rcp45_passyears_2 <- read.table('~/at2degreesgit/RCP45_combined_2deg.txt')

rcp26_passyears_2 <- read.table('~/at2degreesgit/RCP26_combined_2deg.txt')

#plot co2 concentrations
par(mar = c(5,6,3,3), las = 1)
plot(passyears_ssp119$conc_2, rep(1, length(passyears_ssp119$conc_2)),  xlim = c(370,700), ylim = c(0,8),
     type = 'n', axes = FALSE, ylab = '', xlab = 'CO2 concentration at 2 degrees')
abline(h = 1:7, lty = 'dotted')

ssplist_ordered <- c('ssp119','ssp126', 'ssp434', 'ssp245', 'ssp460', 'ssp370', 'ssp585' )
axis(2, at = 1:7, labels = ssplist_ordered)
axis(1)

pch = 21
cex = 1.2
col = 'tomato4'
bg = 'tomato'
points(passyears_ssp119$conc_2, rep(1, length(passyears_ssp119$conc_2)), col = col, bg = bg, pch = pch, cex = cex)

points(rcp26_passyears_2$CO2_ppmv, rep(1.8, length(rcp26_passyears_2$CO2_ppmv)), col = 'skyblue3', pch = 21, cex = cex, bg = 'skyblue2')
points(passyears_ssp126$conc_2, rep(2, length(passyears_ssp126$conc_2)), col = col, bg = bg, pch = pch, cex = cex)

points(passyears_ssp434$conc_2, rep(3, length(passyears_ssp434$conc_2)), col = col, bg = bg, pch = pch, cex = cex)

points(rcp45_passyears_2$CO2_ppmv, rep(3.8, length(rcp45_passyears_2$CO2_ppmv)), col = 'skyblue3', pch = 21, cex = cex,bg = 'skyblue2')
points(passyears_ssp245$conc_2, rep(4, length(passyears_ssp245$conc_2)), col = col, bg = bg, pch = pch, cex = cex)

points(rcp60_passyears_2$CO2_ppmv, rep(4.8, length(rcp60_passyears_2$CO2_ppmv)), col = 'skyblue3', pch = 21, cex = cex, bg = 'skyblue2')
points(passyears_ssp460$conc_2, rep(5, length(passyears_ssp460$conc_2)), col = col, bg = bg, pch = pch, cex = cex)

points(passyears_ssp370$conc_2, rep(6, length(passyears_ssp370$conc_2)), col = col, bg = bg, pch = pch, cex = cex)

points(rcp85_passyears_2$CO2_ppmv, rep(6.8, length(rcp85_passyears_2$CO2_ppmv)), col = 'skyblue3', pch = 21, cex = cex, bg = 'skyblue2')
points(passyears_ssp585$conc_2, rep(7, length(passyears_ssp585$conc_2)), col = col, bg = bg, pch = pch, cex = cex)

mtext('CO2 concentration passing 2 degrees', side = 3, adj = 0, line = -1, cex = 1.3)
legend('bottomright', legend = c('CMIP6','CMIP5'), col = c(col, 'skyblue2'), pch = 21, pt.bg = c(bg, 'skyblue3'), bty = 'n')

rcp85_passyears_1_5 <- read.table('~/at2degreesgit/RCP85_combined_1_5_deg.txt')

rcp60_passyears_1_5 <- read.table('~/at2degreesgit/RCP60_combined_1_5_deg.txt')

rcp45_passyears_1_5 <- read.table('~/at2degreesgit/RCP45_combined_1_5_deg.txt')

rcp26_passyears_1_5 <- read.table('~/at2degreesgit/RCP26_combined_1_5_deg.txt')

#plot co2 concentrations
par(mar = c(5,6,3,3), las = 1)
plot(passyears_ssp119$conc_1.5, rep(1, length(passyears_ssp119$conc_1.5)),  xlim = c(370,700), ylim = c(0,8),
     type = 'n', axes = FALSE, ylab = '', xlab = 'CO2 concentration at 1.5 degrees')
abline(h = 1:7, lty = 'dotted')

ssplist_ordered <- c('ssp119','ssp126', 'ssp434', 'ssp245', 'ssp460', 'ssp370', 'ssp585' )
axis(2, at = 1:7, labels = ssplist_ordered)
axis(1)

pch = 21
cex = 1.2
col = 'tomato4'
bg = 'tomato'
points(passyears_ssp119$conc_1.5, rep(1, length(passyears_ssp119$conc_1.5)), col = col, bg = bg, pch = pch, cex = cex)

points(rcp26_passyears_1_5$CO2_ppmv, rep(1.8, length(rcp26_passyears_1_5$CO2_ppmv)), col = 'skyblue3', pch = 21, cex = cex, bg = 'skyblue2')
points(passyears_ssp126$conc_1.5, rep(2, length(passyears_ssp126$conc_1.5)), col = col, bg = bg, pch = pch, cex = cex)

points(passyears_ssp434$conc_1.5, rep(3, length(passyears_ssp434$conc_1.5)), col = col, bg = bg, pch = pch, cex = cex)

points(rcp45_passyears_1_5$CO2_ppmv, rep(3.8, length(rcp45_passyears_1_5$CO2_ppmv)), col = 'skyblue3', pch = 21, cex = cex,bg = 'skyblue2')
points(passyears_ssp245$conc_1.5, rep(4, length(passyears_ssp245$conc_1.5)), col = col, bg = bg, pch = pch, cex = cex)

points(rcp60_passyears_1_5$CO2_ppmv, rep(4.8, length(rcp60_passyears_1_5$CO2_ppmv)), col = 'skyblue3', pch = 21, cex = cex, bg = 'skyblue2')
points(passyears_ssp460$conc_1.5, rep(5, length(passyears_ssp460$conc_1.5)), col = col, bg = bg, pch = pch, cex = cex)

points(passyears_ssp370$conc_1.5, rep(6, length(passyears_ssp370$conc_1.5)), col = col, bg = bg, pch = pch, cex = cex)

points(rcp85_passyears_1_5$CO2_ppmv, rep(6.8, length(rcp85_passyears_1_5$CO2_ppmv)), col = 'skyblue3', pch = 21, cex = cex, bg = 'skyblue2')
points(passyears_ssp585$conc_1.5, rep(7, length(passyears_ssp585$conc_1.5)), col = col, bg = bg, pch = pch, cex = cex)

mtext('CO2 concentration passing 1.5 degrees', side = 3, adj = 0, line = -1, cex = 1.3)
legend('bottomright', legend = c('CMIP6','CMIP5'), col = c(col, 'skyblue2'), pch = 21, pt.bg = c(bg, 'skyblue3'), bty = 'n')