###################################################
### Preparations (load modules, definitions...)
###################################################

import matplotlib

import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import matplotlib._color_data as mcd
from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes, mark_inset


import numpy as np
import math #as ma
import scipy
from scipy import stats, ndimage, interpolate

import pandas as pd
import netCDF4 as nc4
from netCDF4 import Dataset
from netCDF4 import num2date
import pickle
import datetime
import glob
import os
import time

import mpl_toolkits.basemap
from mpl_toolkits.basemap import Basemap, addcyclic, shiftgrid
import pyferret
pyferret.start(quiet=True)

import importlib
#my scripts:
import tools_plot as jplt
importlib.reload(jplt) 
import tools_analysis as jan
importlib.reload(jan) 



def var_pad(var):
    var_p=np.pad(np.pad(var, ((1,1),(0,0)), mode='edge'), ((0,0),(1,1)), mode='wrap')
    return var_p


paths=sorted(glob.glob('../c3a_model_input_output/*Ma*'))


paths_1000ppm=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm'))
paths_1000ppm_S0ini=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_S0ini'))
paths_1000ppm_S0ini_VegFix=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_S0ini_VegFix'))
paths_1000ppm_VegHom=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_VegHom'))
paths_1000ppm_VegTree=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_VegTree'))
paths_1000ppm_VegBare=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_VegBare'))

paths_orb=sorted(glob.glob('../c3a_model_input_output/*Ma*orb*'))
paths_1000ppm_orb=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_orb_*'))
paths_1000ppm_orb_22p0=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_orb_22p0-*'))
paths_1000ppm_orb_24p5=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_orb_24p5-*'))

for ii in paths_1000ppm_orb:
    paths_orb.remove(ii)
paths_SL=sorted(glob.glob('../c3a_model_input_output/*Ma_*ppm_SL*'))
paths_VegHom_SL=sorted(glob.glob('../c3a_model_input_output/*Ma_1000ppm_VegHom_SL*'))
paths_0Ma=sorted(glob.glob('../c3a_model_input_output/000Ma_*'))


paths_225Ma=sorted(glob.glob('../c3a_model_input_output/225Ma_*'))

paths_250ppm=[]
paths_500ppm=[]
paths_1500ppm=[]
paths_2000ppm=[]
paths_4000ppm=[]
paths_1000ppm_25Ma=[]
timeslices_25Ma=np.append([0,65],np.arange(75,250+25,25))
for tt in timeslices_25Ma[1:]:
    age_str='{:03d}'.format(tt)
    paths_250ppm.append(glob.glob('../c3a_model_input_output/'+age_str+'Ma_250ppm')[0])
    paths_500ppm.append(glob.glob('../c3a_model_input_output/'+age_str+'Ma_500ppm')[0])
    paths_1500ppm.append(glob.glob('../c3a_model_input_output/'+age_str+'Ma_1500ppm')[0])
    paths_2000ppm.append(glob.glob('../c3a_model_input_output/'+age_str+'Ma_2000ppm')[0])
    paths_4000ppm.append(glob.glob('../c3a_model_input_output/'+age_str+'Ma_4000ppm')[0])
    
    paths_1000ppm_25Ma.append(glob.glob('../c3a_model_input_output/'+age_str+'Ma_1000ppm')[0])

pCO2_values_df=pd.read_pickle('./meso_BoundaryConditions_pCO2.pkl')
pCO2_values_df
paths_pCO2Foster=[]
paths_pCO2COPSE=[]
for age in pCO2_values_df.index[1:]:
    age_str='{:03d}'.format(age)
    [pCO2_Foster,pCO2_Mills]=pCO2_values_df.loc[age,['pCO2_Foster_round','pCO2_Mills_round']]
    paths_pCO2Foster.append('../c3a_model_input_output/'+age_str+'Ma_'+str(int(pCO2_Foster))+'ppm')
    paths_pCO2COPSE.append('../c3a_model_input_output/'+age_str+'Ma_'+str(int(pCO2_Mills))+'ppm')

paths_problematic=[]
timeslices_problematic=[]  
timeslices_failed=[]+timeslices_problematic 
paths_failed=[]

for tt in timeslices_problematic:
    tt_string='{:03d}'.format(tt)
    paths_problematic.append('../MesoClim/c3a_meso_'+tt_string+'Ma')
    
for tt in timeslices_failed:
    tt_string='{:03d}'.format(tt)
    paths_failed.append('../MesoClim/c3a_meso_'+tt_string+'Ma')
    

paths_problematic.append('../c3a_model_input_output/000Ma_280ppm_ice_runoffuni')   
paths_problematic.append('../c3a_model_input_output/070Ma_1000ppm_runoff')
paths_problematic.append('../c3a_model_input_output/235Ma_1000ppm_stabilize_runoff')
paths_problematic.append('../c3a_model_input_output/235Ma_1000ppm_stabilize_topo')
paths_problematic.append('../c3a_model_input_output/000Ma_280ppm_ice_NoRunoffTuning')


paths_paper_example=['../c3a_model_input_output/075Ma_1000ppm',
                    '../c3a_model_input_output/125Ma_1000ppm',
                    '../c3a_model_input_output/175Ma_1000ppm',
                    '../c3a_model_input_output/225Ma_1000ppm']

paths_paper_appendix=paths[:]
for pp in paths_SL+paths_0Ma+paths_orb+paths_1000ppm_VegTree+paths_1000ppm_VegBare+paths_problematic:
    if pp in paths_paper_appendix:
        paths_paper_appendix.remove(pp)
        
paths_paper=paths[:]
for pp in paths_0Ma+paths_orb+paths_1000ppm_VegTree+paths_1000ppm_VegBare+paths_problematic:
    if pp in paths_paper:
        paths_paper.remove(pp)

    
    
###Colour list
color_dict=mcd.TABLEAU_COLORS #XKCD_COLORS
color_list=list(color_dict.keys())[:]*30
colors_df=pd.DataFrame([], columns=['color'], index=paths)
colors_df.iloc[:,0]=color_list[0:len(paths)]

color_trias=np.array([129.,43.,146.])/255
color_jura=np.array([52.,178.,201.])/255
color_cret=np.array([127.,198.,78.])/255

color_trias_lower=np.array([168.,88.,158.])/255
color_trias_middle=np.array([188.,134.,186.])/255
color_trias_upper=np.array([199.,166.,207.])/255

color_jura_lower=np.array([0.,176.,233.])/255
color_jura_middle=np.array([132.,207.,237.])/255
color_jura_upper=np.array([188.,228.,250.])/255

color_cret_lower=np.array([160.,201.,109.])/255
color_cret_upper=np.array([186.,210.,95.])/255

# see Stadler: Aridity Indexes in Encyclopedia of World Climatology 2005
#https://en.wikipedia.org/wiki/Aridity_index
def func_aridity_idx(prc,radbal):
    L_wv=2264.705 #Latent heat of water vaporization, kJ/kg, 
    #radbal W/m^2=(J/s)/m^2; prc mm/d~kg/d annual mean
    aridity_budyko=100*radbal/(L_wv*prc) # (J/s)/m^2 * 1/(kJ/kg) * 1/(kg/d/m^2) = d/(1000*s) = 24*60*60/1000 = 86.4    
    # radbal/L_wv is potential evaporation (Arora 2002)
    
    return aridity_budyko

   
latsop,lonsop,lonsop_2D,latsop_2D=jan.grid_gen(region='ocean', cell='p')#[0:1+1]
latsot,lonsot,lonsot_2D,latsot_2D=jan.grid_gen(region='ocean', cell='T')#[0:1+1]

latsap,lonsap,lonsap_2D,latsap_2D=jan.grid_gen(region='atmo', cell='p')#[0:1+1]
latsat,lonsat,lonsat_2D,latsat_2D=jan.grid_gen(region='atmo', cell='T')#[0:1+1]


nc_fid = Dataset('./gridarea_mom.nc', 'r')
#nc_attrs, nc_dims, nc_vars = ncdump(nc_fid)
gridcell_areas_mom=nc_fid.variables['cell_area'][:]/1000**2/1000**2 #in the file its in m**2, convert it to 10^6 km^2



timeslices=np.append(0,np.arange(60,255+5,5))
timeslices_25Ma=np.append(np.array([0,65]),timeslices[4::5]) 
timeslices_SL=[100,150,200]


pi_frac=(np.pi/180)*latsot[:]
timeslices_pcolor=np.append(50-2.5,np.arange(60-2.5,255+2.5+5,5))
timeslices_plot=np.append(50,np.arange(60,255+5,5))
timeslices_25Ma_plot=np.append([50,65],np.arange(75,250+25,25))

timeslices_pcolor2D,lats2D =np.meshgrid(-timeslices_pcolor,latsop)
timeslices_scat2D,latsot2D =np.meshgrid(-timeslices_plot,latsot)





################################################################################
### Load data aggragated by script_aggregateData.py
################################################################################

data_aggregated_history_p2_df=pd.read_pickle('./data_aggregated_history_p2_df.pkl')
data_aggregated_p2_arrays_df=pd.read_pickle('./data_aggregated_p2_arrays_df.pkl')
data_aggregated_geography_df=pd.read_pickle('./data_aggregated_geography_df.pkl')
data_aggregated_SeasonalityZonality_df=pd.read_pickle('./data_aggregated_SeasonalityZonality_df.pkl')  






################################################################################
### Figure 1 - GMST Evolution
################################################################################


fig2=plt.figure(figsize=(14.5,10)) 
plt.clf()


gs = GridSpec(3, 1,height_ratios=[1,0.3,0.3]) 

ax = plt.subplot(gs[0,0]);

ax1 = plt.subplot(gs[2,0],sharex=ax);
ax2 = plt.subplot(gs[1,0],sharex=ax);

ax.tick_params(bottom=False,labelbottom=False)


ax1.grid(True,linestyle=':')
ax1_2 = ax1.twinx()

ax2.grid(True,linestyle=':')
ax2_2 = ax2.twinx()
ax2.tick_params(bottom=False,labelbottom=False)

ax.grid(linestyle=':', which='major',zorder=0) #,lw=1.5
ax.grid(linestyle=':', which='minor',lw=0.5,zorder=0)

#matplotlib.rc('font', size=12)
ax.tick_params(labelsize=14)
ax1.tick_params(labelsize=14)
ax1_2.tick_params(labelsize=14)
ax2.tick_params(labelsize=14)
ax2_2.tick_params(labelsize=14)

ax1_2.tick_params(colors='tab:blue', axis='y')
ax2_2.tick_params(colors='tab:blue', axis='y')



ax.plot(-np.asarray(data_aggregated_history_p2_df.loc[paths_1000ppm[:],['age']].values,dtype=float),np.asarray(data_aggregated_history_p2_df.loc[paths_1000ppm[:],['ts_ann_mean']].values,dtype=float),\
        '-o',c='k',markersize=8, label='pCO2_1000ppm',lw=3,zorder=6) #, label='$\mathsf{P_{pCO2\_1000ppm}}$' 
ax.plot(-np.asarray(data_aggregated_history_p2_df.loc[paths_pCO2Foster[:],['age']].values,dtype=float),np.asarray(data_aggregated_history_p2_df.loc[paths_pCO2Foster[:],['ts_ann_mean']].values,dtype=float),\
        '-o',c='darkgreen',markersize=8, label='pCO2_proxy',lw=3,zorder=5) #, label='$\mathsf{P_{pCO2\_proxy}}$'
ax.plot(-np.asarray(data_aggregated_history_p2_df.loc[paths_pCO2COPSE[:],['age']].values,dtype=float),np.asarray(data_aggregated_history_p2_df.loc[paths_pCO2COPSE[:],['ts_ann_mean']].values,dtype=float),\
        '-o',c='tab:brown',markersize=8, label='pCO2_COPSE',lw=3,zorder=4) #, label='$\mathsf{P_{pCO2\_COPSE}}$'


for pp in [paths_250ppm,paths_500ppm,paths_1500ppm,paths_2000ppm,paths_4000ppm]:
    if pp==paths_250ppm:
        
        label='pCO2_const\n (250/500/1500/...'+'\n'+'...2000/4000ppm)\n'
    else:
        label=''
    ax.plot(-np.asarray(data_aggregated_history_p2_df.loc[pp[:],['age']].values,dtype=float),np.asarray(data_aggregated_history_p2_df.loc[pp[:],['ts_ann_mean']].values,dtype=float),\
        '-',marker='o',c='tab:grey',markersize=4, label=label,lw=1,zorder=3)



landfracs_glob=np.asarray(data_aggregated_geography_df.loc[paths_1000ppm[:],['landfrac_glob']].values[:,0],dtype=float)
landfracs_glob_scaled=(100*(1-landfracs_glob[:])-63)*25
landfracs_color='k'
ax1_2.plot(-timeslices[1:],landfracs_glob_scaled[:],'-',c=landfracs_color,lw=2,label='global ocean\narea fraction\n(no axis)',zorder=2) #,label='global land fraction in C3a simulations'
#ax1_2.scatter(-40,landfracs_glob_scaled[0],c=landfracs_color)
idx_min=np.argmin(landfracs_glob)
idx_max=np.argmax(landfracs_glob)
ax1_2.scatter(-timeslices[idx_max+1],landfracs_glob_scaled[idx_max],c=landfracs_color)
ax1_2.scatter(-timeslices[idx_min+1],landfracs_glob_scaled[idx_min],c=landfracs_color)
ax1_2.annotate(str(np.int(np.round(100*(1.-landfracs_glob[idx_min]))))+'%',xy=(-timeslices[idx_min]+1,landfracs_glob_scaled[idx_min]+0),color=landfracs_color,fontsize=12)
ax1_2.annotate(str(np.int(np.round(100*(1.-landfracs_glob[idx_max]))))+'%',xy=(-timeslices[idx_max]+2,landfracs_glob_scaled[idx_max]-10),color=landfracs_color,fontsize=12)

ax.set_xlim(-258,-55)
ax.set_ylim(9,26) #8,24.5
ax.set_yticks(np.arange(10,26+2,2))

ax.set_xticks([-250,-225,-200,-175,-150,-125,-100,-75,-65])

ax.set_xticklabels([250,225,200,175,150,125,100,75,65])

ax.set_ylabel(r'GMST ($\mathrm{^{\circ}C}$)', fontsize=14)

ax.set_title(r'$\bf{(a)}$'+' Simulated Global Mean Surface Air Temperatures (GMST)',loc='left',fontsize=14) #. Dots represent climate states from this study. Broken lines: Other studies.
 
ax.plot([], [], ' ', label='$\mathbf{GMST\ simulated}$'+'\n'+'$\mathbf{here\ for:}$')

h,l = ax.get_legend_handles_labels()

h_sorted=[h[4],h[0],h[1],h[2],h[3]] #h#[]
l_sorted=[l[4],l[0],l[1],l[2],l[3]] #l#[]
legend_ax=ax.legend(h_sorted,l_sorted,ncol=1,loc='upper left',bbox_to_anchor=(1.01, 0.925),frameon=False,fontsize=11) #,labelspacing=0.1





ax1.set_ylim(0,42)
ax1.set_yticks([10,20,30,40])
ax1.set_ylabel(r'GMST ($\mathsf{^{\circ}C}$)',fontsize=14,color='k')
ax1.set_title(r'$\bf{(c)}$'+' Proxy temperature estimates, sea level reconstructions and model global ocean area fraction',loc='left',fontsize=14)
ax1.set_xlabel('Age (Ma)', fontsize=14)


ax1_2.set_ylim(-110,325)
ax1_2.set_yticks([-0,100,200,300])
ax1_2.set_ylabel('Sea Level (m)',fontsize=14,color='tab:blue')
#ax1_2.set_ylabel('glob. land frac. (%)',fontsize=12,color='tab:brown')

h1,l1 = ax1.get_legend_handles_labels()
h1_2,l1_2 = ax1_2.get_legend_handles_labels()
#ax1.legend(h1+h1_2,l1+l1_2,loc=(0.565,0), frameon=False, fontsize=10,ncol=2,bbox_to_anchor=(1.01, 1))# ,loc='lower right'
legend_ax1=ax1.legend(h1+h1_2,l1+l1_2,loc='upper left',bbox_to_anchor=(1.07, 1.1), frameon=False, fontsize=11,ncol=1)

ax2.set_ylabel('$\mathsf{pCO_2}$ (ppm)',fontsize=14)
ax2.set_ylim(0,2500)
ax2.set_yticks([0,500,1000,1500,2000,2500])
#ax2.legend(loc=(0.8,1), frameon=False, fontsize=10) #loc='upper right'

ax2_2.set_ylabel('Ridge Length ($\mathsf{10^3\,km}$)',fontsize=14,color='tab:blue')
#ax2.set_ylim(100,2800)
h2,l2 = ax2.get_legend_handles_labels()
h2_2,l2_2 = ax2_2.get_legend_handles_labels()

legend_ax2=ax2.legend(h2+h2_2,l2+l2_2,loc='upper left',bbox_to_anchor=(1.07, 1.2), frameon=False, fontsize=11,ncol=1) #,labelspacing=0.0

ax2.set_title(r'$\bf{(b)}$'+' Reconstructions of pCO$\mathsf{_2}$ and global mid-ocean ridge length',loc='left',fontsize=14)



jplt.add_geological_timescale(ax1,height=0.11)


ax.grid(linestyle=':', which='major',zorder=0) #,lw=1.5
ax.grid(linestyle=':', which='minor',lw=0.5,zorder=0)



ax.text(-250,24.2,'4000ppm',color='tab:grey',ha='center')
ax.text(-250,22,'2000ppm',color='tab:grey',ha='center')
ax.text(-250,20.7,'1500ppm',color='tab:grey',ha='center')
#ax.text(-62,20,'1000ppm',color='k') #tab:green
ax.text(-250,15.4,'500ppm',color='tab:grey',ha='center')
ax.text(-250,10.8,'250ppm',color='tab:grey',ha='center')

for llegend in [legend_ax,legend_ax1,legend_ax2]:
    h = llegend.legendHandles
    t = llegend.texts
    renderer = fig2.canvas.get_renderer()

    for i in range(len(h)): #range(6,len(h),1):
        t_str=t[i].get_text()
        if t_str.find('\n')!=-1 and t_str.find('$')==-1:
            if t_str.count('\n')==1:
                fac=1.
            if t_str.count('\n')==2:
                fac=0.85
            hbbox = h[i].get_window_extent(renderer) # bounding box of handle
            tbbox = t[i].get_window_extent(renderer) # bounding box of text
            x = tbbox.x0 # keep default horizontal position    
            y = fac=(hbbox.height - tbbox.height) / 2. + hbbox.y0 # vertically center the  bbox of the text to the bbox of the handle.
            t[i].set_position((x, y)) # set new position of the text

plt.tight_layout(h_pad=-3.25) #h_pad=-1.3
fig2.savefig('./Fig_1_GMST.pdf',bbox_inches='tight', dpi=900)


################################################################################
### Figure 2 - Global Mean Temperature Sensitivity Experiments
################################################################################


fig2=plt.figure(figsize=(8,6)) 
plt.clf()

gs = GridSpec(2, 2,height_ratios=[1,0.4], width_ratios=[1,0.05]) 

ax00 = plt.subplot(gs[1,0]);
#ax00_twin=ax00.twinx();

ax = plt.subplot(gs[0,0]);

for aa in [ax]: 
    aa.tick_params(labelbottom=False)

ax.grid(linestyle=':', which='major',zorder=0) #,lw=1.5
ax.grid(linestyle=':', which='minor',lw=0.5,zorder=0)

#matplotlib.rc('font', size=12)
for aa in [ax,ax00]:
    aa.tick_params(labelsize=12)
    aa.grid(linestyle=':', which='major',zorder=0) #,lw=1.5
    aa.grid(linestyle=':', which='minor',lw=0.5,zorder=0)


ax.plot(-np.asarray(data_aggregated_history_p2_df.loc[paths_1000ppm[1:],['age']].values,dtype=float),np.asarray(data_aggregated_history_p2_df.loc[paths_1000ppm[1:],['ts_ann_mean']].values,dtype=float),\
    '-',marker='.',c='k',markersize=12, label='$\mathtt{{pCO2\_1000ppm}}$',lw=3,zorder=3) #,markeredgecolor='k'
ax.plot(-np.asarray(data_aggregated_history_p2_df.loc[paths_1000ppm_S0ini,['age']].values,dtype=float),np.asarray(data_aggregated_history_p2_df.loc[paths_1000ppm_S0ini,['ts_ann_mean']].values,dtype=float),\
    '-o',c='tab:orange',markersize=8, label='$\mathtt{{S0ini}}$',lw=1,zorder=3) #,markeredgecolor='k'
ax.plot(-np.asarray(data_aggregated_history_p2_df.loc[paths_1000ppm_S0ini_VegFix,['age']].values,dtype=float),np.asarray(data_aggregated_history_p2_df.loc[paths_1000ppm_S0ini_VegFix,['ts_ann_mean']].values,dtype=float),\
    '-o',c='tab:green',markersize=8, label='$\mathtt{{VegFix\_S0ini}}$',lw=1,zorder=3) #,markeredgecolor='k'

ax.plot(-np.asarray(data_aggregated_history_p2_df.loc[paths_1000ppm_VegHom[1:],['age']].values,dtype=float),np.asarray(data_aggregated_history_p2_df.loc[paths_1000ppm_VegHom[1:],['ts_ann_mean']].values,dtype=float),\
    '-o',c='tab:blue',markersize=8, label='$\mathtt{{VegHom}}$',lw=1,zorder=3) #,markeredgecolor='k'


for pp in paths_VegHom_SL:
    #if pp not in paths_problematic:
    if pp not in ['../c3a_meso_timeslice/200Ma_1000ppm_VegHom_SL+80']:
        if pp==paths_VegHom_SL[0]:
            label='$\mathtt{{SL}}$'
        else:
            label=''
        age=-data_aggregated_history_p2_df.loc[pp,['age']][0]
        gmst=data_aggregated_history_p2_df.loc[pp,['ts_ann_mean']][0]
        ax.scatter(age,gmst,c='None',marker='o',edgecolors='tab:blue',linewidths=2.5,zorder=4,label=label)
        #ax.annotate(pp[pp.index('SL'):]+'m',xy=(age,gmst),color='tab:blue',fontsize=12)
        ax.annotate(pp[pp.index('SL')+2:]+'m',xy=(age+2,gmst),color='tab:blue',fontsize=12)
        
for pp in paths_1000ppm_orb:
    if pp==paths_1000ppm_orb_22p0[0]:
        label='$\mathtt{orb}$, obl.=22.0$^{\circ}$'            
    elif pp==paths_1000ppm_orb_24p5[0]:
        label='orb, obl.=24.5$^{\circ}$'
    else:
        label=''

    if pp in paths_1000ppm_orb_22p0:
        marker='s'

    elif pp in paths_1000ppm_orb_24p5:
        marker='D'
    else:
        marker='+'
#    marker='D'
    age=-data_aggregated_history_p2_df.loc[pp,['age']][0]
    gmst=data_aggregated_history_p2_df.loc[pp,['ts_ann_mean']][0]
    ax.scatter(age,gmst,c='None',marker=marker,edgecolors='k',linewidths=2.5,zorder=2,label=label)

        
landfracs_glob=np.asarray(data_aggregated_geography_df.loc[paths_1000ppm[:],['landfrac_glob']].values[:,0],dtype=float)
ax00.plot(-np.asarray(data_aggregated_geography_df.loc[paths_1000ppm[1:],['age']].values),100.*landfracs_glob[1:],'-',c='k',lw=2,zorder=2) #,label='global land fraction',label='global land fraction in C3a simulations'
#ax21.text(0.70,0.8,'global land fraction', fontsize=8,color='m', va='bottom', ha='left',transform=ax01.transAxes,fontweight='normal')
for pp in range(len(paths_VegHom_SL)):
    path=paths_VegHom_SL[pp]
    if path not in ['../c3a_meso_timeslice/200Ma_1000ppm_VegHom_SL+80']:
        age=-data_aggregated_geography_df.loc[path,['age']][0]
        landfrac=data_aggregated_geography_df.loc[path,['landfrac_glob']]
        ax00.scatter(age,100.*data_aggregated_geography_df.loc[path,['landfrac_glob']],c='None',marker='o',edgecolors='tab:blue',linewidths=2.5,zorder=4)
        ax00.annotate(path[path.index('SL')+2:]+'m',xy=(age+3,100.*landfrac-0.5),color='tab:blue',fontsize=12)

    
#ax.set_xlim(-258,-33)
ax.set_xlim(-258,-55)
ax.set_ylim(17,21.6) #8,24.5


ax.legend(ncol=1,loc='upper left',labelspacing=0,frameon=False)

ax.set_ylabel(r'GMST ($\mathrm{^{\circ}C}$)', fontsize=14)

#ax.set_title('Simulated Global Mean Surface Air Temperatures. Dots represent climate states from this study. Broken lines: Other studies. ',loc='left')

ax00.set_ylim((23,40))
ax00.set_yticks([25,30,35,40])
ax00.set_ylabel('Land\nFraction (%)',fontsize=14)
ax00.set_xlabel('Age (Ma)', fontsize=14)
ax00.set_xticks([-250,-225,-200,-175,-150,-125,-100,-75,-65])
ax00.set_xticklabels([250,225,200,175,150,125,100,75,65])


ax00.set_title(r'$\bf{(b)}$ '+'Model global land area fraction',loc='left',fontsize=14)
ax.set_title(r'$\bf{(a)}$ '+'Simulated global mean temperatures',loc='left',fontsize=14)


jplt.add_geological_timescale(ax00,stagenames=False,height=0.11)



plt.tight_layout(h_pad=-0.25) #h_pad=-1.5,
fig2.savefig('./Fig_2_GMST_Sensitivity.pdf',bbox_inches='tight', dpi=900)




################################################################################
### Figure 3 - Seasonality
################################################################################



fig=plt.figure(figsize=(11.5,6)) 
plt.clf()

gs = GridSpec(5, 3,height_ratios=[0.25,0.25,0.15,0.35,0.05], width_ratios=[1,1,0.05]) #, width_ratios=[0.5,1,0.15] , height_ratios=[0.4,1]
#gs = GridSpec(2, 2,height_ratios=[0.5,1], width_ratios=[1,0.05]) 

ax00 = plt.subplot(gs[3:,1]); #[0,0]

ax10 = plt.subplot(gs[0:2+1,1]); #gs[[1:4+1,0]
ax10_cbar = plt.subplot(gs[0:2+1,2]);
ax10_cbar.axis('off')


ax01 = plt.subplot(gs[0:1+1,0]);
ax11 = plt.subplot(gs[2:3+1,0]);
for aa in [ax01,ax11]:
    aa.axis('off')
ax11_cbar = plt.subplot(gs[4,0]);
ax11_cbar.axis('off')   
    
  
    
for aa in [ax10]: #,ax01,ax20,ax21,ax20 ,ax21 ,ax10 ,ax11
    aa.tick_params(labelbottom=False)

ax00.plot(-timeslices_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm',['ts_variations_globalmean']][0][:],'-',c='k',lw=2,zorder=4,label='pCO2_1000ppm')
ax00.plot(-timeslices_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_pCO2Foster',['ts_variations_globalmean']][0][:],'--',c='k',lw=2,zorder=4,label='pCO2_proxy')
ax00.plot(-timeslices_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_pCO2COPSE',['ts_variations_globalmean']][0][:],':',c='k',lw=2,zorder=4,label='pCO2_COPSE')

landareas_zonal=data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm',['landareas_zonal']][0]


vmin=0; vmax=7;
a10=ax10.pcolor(timeslices_pcolor2D[:,1:],lats2D[:,1:],landareas_zonal[:,:],cmap='Greys',zorder=1,vmin=vmin,vmax=vmax) #,vmax=25
cbar10=fig.colorbar(a10, ax=ax10_cbar, orientation='vertical', boundaries=np.arange(vmin,vmax+0.5,0.5),ticks=np.arange(vmin,vmax+1,1),fraction=1, pad=0.04) #,extend='max' extend=extend, shrink=0.621, pad=0.02

vmin=0; vmax=40;

for pp in range(2):
    if pp==0:
        path='../c3a_model_input_output/225Ma_1000ppm'
        ax_path=ax01
        age=225
        
    if pp==1:
        path='../c3a_model_input_output/075Ma_1000ppm'
        ax_path=ax11
        age=75  
    continental_mask=data_aggregated_geography_df.loc[path,['continental_mask']][0]
    ts_season_arr=data_aggregated_p2_arrays_df.loc[path, ['ts_seasonal_arr']].values[0]
    ts_season_variation=np.max(ts_season_arr,axis=2)-np.min(ts_season_arr,axis=2) 
    ts_season_variation_intp=jan.interpol_atm_to_ocn_grid(ts_season_variation)
    ts_season_variation_intp_masked=np.ma.MaskedArray(ts_season_variation_intp,mask=continental_mask)
    a11=jplt.plot_world(np.ma.MaskedArray(ts_season_variation_intp_masked), varname='', lim_l=vmin, lim_u=vmax, cbar_delta=1, cbar_tick_freq=2, units='mm/d',interpol_var=False,\
                         projection='robin', title_opt='off', axes='off', anno='off',latlabels=[True,False,False,False], continents='on',continents_timeslice=age,plates='off',\
                         colourmap='RdBu_r',extend='max', var_digits=0,\
                         cont='off', vec='off',scat='off',hatches='off',fig_handle=[fig,ax_path])[0] #[1]                         text_corner=str(age)+'Ma',text_corner_posi=[0.1,0.55],\

cbar11=fig.colorbar(a11, ax=ax11_cbar, orientation='horizontal',fraction=1.1, boundaries=np.arange(vmin,vmax+1,1),ticks=np.arange(vmin,vmax+5,5)) #, ticks=cbar_ticks, cax=ax1,pad=-2 , shrink=0.9
#cbar11.ax.set_title('$\mathsf{\Delta SAT_{zonal}}$ ($\mathsf{^{\circ}C}$)',fontsize=10) #,y=-0.2
cbar11.ax.set_xlabel('$\mathsf{\Delta SAT_{season}}$ ($\mathsf{^{\circ}C}$)',fontsize=14)
cbar11.ax.tick_params(labelsize=12) #,labeltop=True,labelbottom=False

cont_color_list=['None','tab:red','tab:orange','tab:blue','tab:cyan','magenta'] #,'w'
cont_levels_list=[0,2.5,10,20,30,40]
cont=ax10.contour(timeslices_scat2D[:,1:],latsot2D[:,1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm',['ts_season_variation_zonalmean']][0][:,0:],
                  levels=cont_levels_list,linewidths=2,colors=cont_color_list,zorder=3)

fmt='{}{}'.format('%g', '$\mathsf{^{\circ}C}$')
ax10.clabel(cont, fontsize=12, inline=True,fmt=fmt)
for ccll in range(len(cont_levels_list)):
    cbar11.ax.plot([(cont_levels_list[ccll]-vmin)/(vmax-vmin),(cont_levels_list[ccll]-vmin)/(vmax-vmin)],[0,100],c=cont_color_list[ccll],lw=3,zorder=10)
#cbar11.ax.plot([0,1],[0,1],c='k',lw=3,zorder=10)



ax00.plot(-timeslices_25Ma_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm_VegHom',['ts_variations_globalmean']][0][:],'-',c='tab:blue',lw=1,zorder=4,label='VegHom')
ax00.plot(-timeslices_25Ma_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm_S0ini_VegFix',['ts_variations_globalmean']][0],'-',c='limegreen',lw=1,zorder=4,label='VegFix_S0ini')
ax00.plot(-timeslices_25Ma_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm_S0ini',['ts_variations_globalmean']][0],'-',c='tab:orange',lw=1,zorder=4,label='S0ini')
ax00.plot(-timeslices_25Ma_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_500ppm',['ts_variations_globalmean']][0][:],'v',c='k',markerfacecolor='tab:cyan',markeredgewidth=1.5,lw=2,zorder=4,label='500ppm')
ax00.plot(-timeslices_25Ma_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_2000ppm',['ts_variations_globalmean']][0][:],'^',c='k',markerfacecolor='tab:red',markeredgewidth=1.5,lw=2,zorder=4, label='2000ppm')

for pp in range(len(paths_VegHom_SL)):
    #if paths_VegHom_SL[pp] not in paths_problematic: 
    if paths_VegHom_SL[pp] not in ['../c3a_model_input_output/200Ma_1000ppm_VegHom_SL+80']:
        if pp==0:
            label='SL'
        else:
            label=''
        ax00.plot(-data_aggregated_history_p2_df.loc[paths_VegHom_SL[pp],['age']][0],data_aggregated_SeasonalityZonality_df.loc['paths_VegHom_SL',['ts_variations_globalmean']][0][pp],'o',\
                  c='tab:blue',markerfacecolor='None',markeredgewidth=2,lw=2,zorder=6,label=label)

for pp in range(len(paths_1000ppm_orb)):
        path=paths_1000ppm_orb[pp]

       
        if path==paths_1000ppm_orb_22p0[0]:
            label='$\mathtt{orb}$, 22.0$^{\circ}$'            
        elif path==paths_1000ppm_orb_24p5[0]:
            label='$\mathtt{orb}$, 24.5$^{\circ}$'
        else:
            label=''
        if path in paths_1000ppm_orb_22p0:
            marker='s'
            
        elif path in paths_1000ppm_orb_24p5:
            marker='D'
        else:
            marker='+'
        ax00.plot(-data_aggregated_history_p2_df.loc[path,['age']][0],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm_orb',['ts_variations_globalmean']][0][pp],linestyle='None',marker=marker,\
                  c='k',markerfacecolor='None',markeredgewidth=2,lw=2,zorder=1,label=label)


    
for aa in [ax00,ax10]: #,ax01,ax11,ax21,ax31,ax20,ax30,ax30,ax31
    aa.grid(True,linestyle=':')
    aa.set_xticks(np.append(np.arange(-250,-75+25,25),[-65,-50]))
    #aa.set_xlim((-257.5,-47.5))
    aa.set_xlim((-257.5,-57.5))
    
for aa in [ax10]: #,ax11,ax30,ax31,ax30,ax31
    aa.set_ylim((-90,90))
    aa.set_yticks(np.arange(-90,90+30,30))
    
for aa in [ax10]:      #,ax30,ax30
    aa.set_yticklabels(['90$^{\circ}$S','60$^{\circ}$S','30$^{\circ}$S','0$^{\circ}$','30$^{\circ}$N','60$^{\circ}$N','90$^{\circ}$N'])


for aa in [ax00]:     #,ax11 ax30,ax31,ax30,ax31

    aa.set_xticklabels(np.append(np.arange(250,75-25,-25),[65,0]))
# for aa in [ax30,ax31]:
    aa.set_xlabel('Age (Ma)',fontsize=14)

for aa in [ax00,ax10]:    #,ax01,ax11,ax20,ax30,ax21,ax31,ax31 ,ax30
    aa.tick_params(labelsize=12)

for cbar in [cbar10]: #,cbar30
    cbar.ax.set_ylabel('Land Area $\mathsf{(10^6\,km^2)}$',fontsize=14)
    cbar.ax.tick_params(labelsize=12)
    


ax00.set_ylim((11.5,25.5))
ax00.set_yticks([15,20,25])

ax00.set_ylabel('$\mathsf{\Delta SAT_{season}\ (^{\circ}C)}$',fontsize=14)


ax00.legend(loc=(1.0,0.00),frameon=False,labelspacing=0.1,ncol=1,fontsize=8) #loc='upper right'loc=(0.7,0.73)



ax00.set_title(r'$\bf{(d)}$ '+'Seasonal SAT contrasts (global avg.)',loc='left',fontsize=14)
ax10.set_title(r'$\bf{(c)}$ '+'Seasonal SAT contrasts (zonal avg., contours)',loc='left',fontsize=14)
ax01.set_title(r'$\bf{(a)}$'+' 225Ma', loc='left',fontsize=14)
ax11.set_title(r'$\bf{(b)}$'+' 75Ma', loc='left',fontsize=14)




for aa in [ax00]: #ax30,ax31,,ax11
    jplt.add_geological_timescale(aa,height=0.085,fontsize=8,stagenames=False,textoffset=0.015) #,fontweight='normal'

plt.tight_layout(w_pad=-7,h_pad=-0.75) #w_pad=-4,h_pad=-0.3

fig.savefig('./Fig_3_SAT_seasonality.pdf',bbox_inches='tight', dpi=900)










################################################################################
### Figure 4 - Meridional Temperature Gradient
################################################################################



fig=plt.figure(figsize=(6.4,8)) 
plt.clf()

gs = GridSpec(3, 2,height_ratios=[0.9,0.5,0.5], width_ratios=[1,0.05]) 

ax00 = plt.subplot(gs[1,0]);
ax00_SP = plt.subplot(gs[2,0],sharey=ax00);
#ax00_twin=ax00.twinx();

ax10 = plt.subplot(gs[0,0]);
ax10_cbar = plt.subplot(gs[0,1]);
ax10_cbar.axis('off')



for aa in [ax10,ax00]: #,ax01,ax20,ax21,ax20 ,ax21 ,ax10 ,ax11
    aa.tick_params(labelbottom=False)


landareas_zonal=data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm',['landareas_zonal']][0]

vmin=0; vmax=7
a10=ax10.pcolor(timeslices_pcolor2D[:,1:],lats2D[:,1:],landareas_zonal[:,:],cmap='Greys',zorder=1,vmin=vmin,vmax=vmax) #,vmax=25
cbar10=fig.colorbar(a10, ax=ax10_cbar, orientation='vertical', boundaries=np.arange(vmin,vmax+0.5,0.5),ticks=np.arange(vmin,vmax+1,1),fraction=1, pad=0.04) #,extend='max' extend=extend, shrink=0.621, pad=0.02



cont_color_list=['magenta','tab:red','tab:orange','tab:blue','tab:cyan']#['tab:cyan','tab:blue','gold','tab:orange','tab:red'] # ,'tab:orange''tab:green'] #] #,'tab:cyan','w']
#['magenta','tab:red','tab:orange','tab:blue','tab:cyan']
cont_levels_list=[-5,0,5,15,25]
#cont_levels_list=[0,5,10,15]
cont=ax10.contour(timeslices_scat2D[:,1:],latsot2D[:,1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm',['ts_anns_air_zonalmean']][0][:,0:],levels=cont_levels_list,linewidths=2,colors=cont_color_list,zorder=3)
#ax10.contour(timeslices_scat2D,latsot2D[:,0:],monsoon_seasonality_zonality_df.loc['paths_pCO2Foster',['ts_anns_zonality']][0][:,0:],levels=cont_levels_list,linewidths=2,linestyles='dashed',colors=cont_color_list,zorder=3)
#ax10.contour(timeslices_scat2D,latsot2D[:,0:],monsoon_seasonality_zonality_df.loc['paths_pCO2COPSE',['ts_anns_zonality']][0][:,0:],levels=cont_levels_list,linewidths=2,linestyles='dotted',colors=cont_color_list,zorder=3)
fmt='{}{}'.format('%d', '$\mathsf{^{\circ}C}$')
ax10.clabel(cont, fontsize=14, inline=True,fmt=fmt) #, fmt=cont_label_fmt

for ii in range(2):
    if ii==0:
        aa=ax00
    elif ii==1:
        aa=ax00_SP
    aa.plot(-timeslices_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm',['ts_grads']][0][ii,:],'-',c='k',lw=2,zorder=4,label='pCO2_1000ppm')
    aa.plot(-timeslices_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_pCO2Foster',['ts_grads']][0][ii,:],'--',c='k',lw=2,zorder=4,label='pCO2_proxy')
    aa.plot(-timeslices_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_pCO2COPSE',['ts_grads']][0][ii,:],':',c='k',lw=2,zorder=4,label='pCO2_COPSE')


    aa.plot(-timeslices_25Ma_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm_VegHom',['ts_grads']][0][ii,:],'-',c='tab:blue',lw=1,zorder=4,label='VegHom')

    aa.plot(-timeslices_25Ma_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm_S0ini_VegFix',['ts_grads']][0][ii,:],'-',c='limegreen',lw=1,zorder=4,label='VegFix_S0ini')

    aa.plot(-timeslices_25Ma_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm_S0ini',['ts_grads']][0][ii,:],'-',c='tab:orange',lw=1,zorder=4,label='S0ini')

    aa.plot(-timeslices_25Ma_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_500ppm',['ts_grads']][0][ii,:],'v',c='k',markerfacecolor='tab:cyan',markeredgewidth=1.5,lw=2,zorder=4,label='500ppm')
    aa.plot(-timeslices_25Ma_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_2000ppm',['ts_grads']][0][ii,:],'^',c='k',markerfacecolor='tab:red',markeredgewidth=1.5,lw=2,zorder=4,label='2000ppm')

    for pp in range(len(paths_VegHom_SL)):
        #if paths_VegHom_SL[pp] not in paths_problematic: 
        if paths_VegHom_SL[pp] not in ['../c3a_model_input_output/200Ma_1000ppm_VegHom_SL+80']:
            if pp==0:
                label='SL'
            else:
                label=''
            aa.plot(-data_aggregated_history_p2_df.loc[paths_VegHom_SL[pp],['age']][0],data_aggregated_SeasonalityZonality_df.loc['paths_VegHom_SL',['ts_grads']][0][ii,pp],'o',\
                      c='tab:blue',markerfacecolor='None',markeredgewidth=2,lw=2,zorder=6,label=label)
    for pp in range(len(paths_1000ppm_orb)):
        path=paths_1000ppm_orb[pp]

        if path==paths_1000ppm_orb_22p0[0]:
            label='orb, 22.0'            
        elif path==paths_1000ppm_orb_24p5[0]:
            label='orb, 24.5'
        else:
            label=''
        if path in paths_1000ppm_orb_22p0:
            marker='s'            
        elif path in paths_1000ppm_orb_24p5:
            marker='D'
        else:
            marker='+'
        aa.plot(-data_aggregated_history_p2_df.loc[path,['age']][0],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm_orb',['ts_grads']][0][ii,pp],marker=marker,\
                  c='k',markerfacecolor='None',markeredgewidth=2,linestyle='None',zorder=1,label=label)


    
    
for aa in [ax00,ax00_SP,ax10]: #,ax01,ax11,ax21,ax31,ax20,ax30,ax30,ax31
    aa.grid(True,linestyle=':')
    aa.set_xticks(np.append(np.arange(-250,-75+25,25),[-65,-50]))
    #aa.set_xlim((-257.5,-47.5))
    aa.set_xlim((-257.5,-57.5))
    
for aa in [ax10]: #,ax11,ax30,ax31,ax30,ax31
    aa.set_ylim((-90,90))
    aa.set_yticks(np.arange(-90,90+30,30))
    
for aa in [ax10]:      #,ax30,ax30
    aa.set_yticklabels(['90$^{\circ}$S','60$^{\circ}$S','30$^{\circ}$S','0$^{\circ}$','30$^{\circ}$N','60$^{\circ}$N','90$^{\circ}$N'])
    #aa.set_ylabel('latitude $\mathsf{\phi\ (^{\circ})}$',fontsize=14)

for aa in [ax00_SP]:     #,ax11ax30,ax31,ax30,ax31
    #aa.set_xticklabels(np.append(np.arange(-250,-75+25,25),[-65,0]))
    aa.set_xticklabels(np.append(np.arange(250,75-25,-25),[65,0]))


for aa in [ax00,ax00_SP,ax10]:    #,ax01,ax11,ax20,ax30,ax21,ax31,ax31 ,ax30
    aa.tick_params(labelsize=12)

for cbar in [cbar10]: #,cbar30
    cbar.ax.set_ylabel('Land Area $\mathsf{(10^6\,km^2)}$',fontsize=14)
    cbar.ax.tick_params(labelsize=12)
    

ax00.set_ylim((21,36))



ax00.set_ylabel('$\mathsf{\Delta SAT_{merid,NH}\ (^{\circ}C)}$',fontsize=14)
ax00_SP.set_ylabel('$\mathsf{\Delta SAT_{merid,SH}\ (^{\circ}C)}$',fontsize=14)

ax00_SP.set_xlabel('Age (Ma)',fontsize=14)
    

ax10.set_title(r'$\bf{(a)}$ '+'SAT (zonal average, contours)',loc='left',fontsize=14)
ax00.set_title(r'$\bf{(b)}$ '+'Northern meridional SAT gradient',loc='left',fontsize=14)
ax00_SP.set_title(r'$\bf{(c)}$ '+'Southern meridional SAT gradient',loc='left',fontsize=14)


ax00_SP.legend(loc=(1.00,0.00),frameon=False,labelspacing=0.2,ncol=1,fontsize=8,borderaxespad=0) #loc='upper right'loc=(0.7,0.73)
#ax00.legend(loc=(1,0.1),frameon=False,labelspacing=0.1,ncol=1,fontsize=8)

for aa in [ax00_SP]: #,ax11ax30,ax31,
    jplt.add_geological_timescale(aa,height=0.085,fontsize=8,stagenames=False) #,fontweight='normal'

plt.tight_layout(w_pad=-8,h_pad=-0.4) #w_pad=-4,h_pad=-0.3 w_pad=-4,h_pad=-0.5

fig.savefig('./Fig_4_SAT_meridionalgradient.pdf',bbox_inches='tight', dpi=900)








################################################################################
### Figure 5 - Zonality
################################################################################




fig=plt.figure(figsize=(11.5,6)) 
plt.clf()

gs = GridSpec(5, 3,height_ratios=[0.25,0.25,0.15,0.35,0.05], width_ratios=[1,1,0.05]) #, width_ratios=[0.5,1,0.15] , height_ratios=[0.4,1]


ax00 = plt.subplot(gs[3:,1]); #[0,0]

ax10 = plt.subplot(gs[0:2+1,1]); #gs[[1:4+1,0]
ax10_cbar = plt.subplot(gs[0:2+1,2]);
ax10_cbar.axis('off')


ax01 = plt.subplot(gs[0:1+1,0]);
ax11 = plt.subplot(gs[2:3+1,0]);
for aa in [ax01,ax11]:
    aa.axis('off')
ax11_cbar = plt.subplot(gs[4,0]);
ax11_cbar.axis('off') 
    


for aa in [ax10]: #,ax01,ax20,ax21,ax20 ,ax21 ,ax10 ,ax11
    aa.tick_params(labelbottom=False)



ax00.plot(-timeslices_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm',['ts_zonality_globalmean']][0][:],'-',c='k',lw=2,zorder=4,label='pCO2_1000ppm')
ax00.plot(-timeslices_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_pCO2Foster',['ts_zonality_globalmean']][0][:],'--',c='k',lw=2,zorder=4,label='pCO2_proxy')
ax00.plot(-timeslices_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_pCO2COPSE',['ts_zonality_globalmean']][0][:],':',c='k',lw=2,zorder=4,label='pCO2_COPSE')

data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm',['prc_anns_zonality']][0]


vmin=0; vmax=4
a10=ax10.pcolor(timeslices_pcolor2D[:,1:],lats2D[:,1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm',['ts_anns_zonality']][0][:,:],cmap='RdBu_r',zorder=1,vmin=vmin,vmax=vmax) #,vmax=25
#ax10.pcolor(timeslices_pcolor2D[:,0:1+1]+5,lats2D[:,0:1+1],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm',['ts_anns_zonality']][0][:,0:1+1],cmap='RdBu_r',zorder=1,vmin=vmin,vmax=vmax) #,vmax=25
cbar10=fig.colorbar(a10, ax=ax10_cbar, orientation='vertical',fraction=1.1, pad=0.04, boundaries=np.arange(vmin,vmax+0.25,0.25),ticks=np.arange(vmin,vmax+1,1)) #

cont_color_list=['None','tab:red','tab:blue'] # ,'tab:orange''tab:green'] #] #,'tab:cyan','w']
cont_levels_list=[0,1,2]
#cont_levels_list=[0,5,10,15]
ax10.contour(timeslices_scat2D[:,1:],latsot2D[:,1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm',['ts_anns_zonality']][0][:,0:],levels=cont_levels_list,linewidths=2,colors=cont_color_list,zorder=3)

for ccll in range(len(cont_levels_list)):
    cbar10.ax.plot([0,1],[(cont_levels_list[ccll]-vmin)/(vmax-vmin),(cont_levels_list[ccll]-vmin)/(vmax-vmin)],c=cont_color_list[ccll],lw=3)




ax00.plot(-timeslices_25Ma_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm_VegHom',['ts_zonality_globalmean']][0][:],'-',c='tab:blue',lw=1,zorder=4,label='VegHom')

ax00.plot(-timeslices_25Ma_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm_S0ini_VegFix',['ts_zonality_globalmean']][0],'-',c='limegreen',lw=1,zorder=4,label='VegFix_S0ini')

ax00.plot(-timeslices_25Ma_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm_S0ini',['ts_zonality_globalmean']][0],'-',c='tab:orange',lw=1,zorder=4,label='S0ini')

ax00.plot(-timeslices_25Ma_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_500ppm',['ts_zonality_globalmean']][0][:],'v',c='k',markerfacecolor='tab:cyan',markeredgewidth=1.5,lw=2,zorder=4,label='500ppm')
ax00.plot(-timeslices_25Ma_plot[1:],data_aggregated_SeasonalityZonality_df.loc['paths_2000ppm',['ts_zonality_globalmean']][0][:],'^',c='k',markerfacecolor='tab:red',markeredgewidth=1.5,lw=2,zorder=4,label='2000ppm')



for pp in range(len(paths_VegHom_SL)):
    #if paths_VegHom_SL[pp] not in paths_problematic: 
    if paths_VegHom_SL[pp] not in ['../c3a_model_input_output/200Ma_1000ppm_VegHom_SL+80']:
        if pp==0:
            label='SL'
        else:
            label=''
        ax00.plot(-data_aggregated_history_p2_df.loc[paths_VegHom_SL[pp],['age']][0],data_aggregated_SeasonalityZonality_df.loc['paths_VegHom_SL',['ts_zonality_globalmean']][0][pp],'o',\
                  c='tab:blue',markerfacecolor='None',markeredgewidth=2,lw=2,zorder=6,label=label)


for pp in range(len(paths_1000ppm_orb)):
        path=paths_1000ppm_orb[pp]
        if pp==0:
            label='orb'
        else:
            label=''

        marker='D'
        ax00.plot(-data_aggregated_history_p2_df.loc[path,['age']][0],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm_orb',['ts_zonality_globalmean']][0][pp],'D',\
                  c='k',markerfacecolor='None',markeredgewidth=2,lw=2,zorder=1,label=label)

for pp in range(2):
    if pp==0:
        path='../c3a_model_input_output/225Ma_1000ppm'
        ax_path=ax01
        age=225
    if pp==1:
        path='../c3a_model_input_output/075Ma_1000ppm'
        ax_path=ax11
        age=75  
    ts_ann_arr=data_aggregated_p2_arrays_df.loc[path, ['ts_ann_arr']].values[0]
    continental_mask=data_aggregated_geography_df.loc[path,['continental_mask']][0]
    ts_ann_intp=jan.interpol_atm_to_ocn_grid(ts_ann_arr)
    ts_ann_intp_masked=np.ma.MaskedArray(ts_ann_intp,mask=continental_mask)
    ts_anns_zonalmean=np.ma.mean(ts_ann_intp_masked,axis=1)
    ts_anns_zonality=ts_ann_intp_masked-np.transpose(np.tile(ts_anns_zonalmean,(96,1)))
    a11=jplt.plot_world(np.ma.MaskedArray(ts_anns_zonality), varname='prc_jja', lim_l=-4, lim_u=4, cbar_delta=0.5, cbar_tick_freq=2, units='mm/d',interpol_var=False,\
                         projection='robin', title_opt='off', axes='off', anno='off',latlabels=[True,False,False,False], continents='on',continents_timeslice=age,plates='off',\
                         colourmap='RdBu_r',extend='max', var_digits=0,\
                         cont='off', vec='off',scat='off',hatches='off',fig_handle=[fig,ax_path])[0] #[1] ,text_corner=str(age)+'Ma'
#, text_corner=str(age)+'Ma'
cbar11=fig.colorbar(a11, ax=ax11_cbar, orientation='horizontal',fraction=1, boundaries=np.arange(-4,4+0.5,0.5),ticks=np.arange(-4,4+1,1)) #, ticks=cbar_ticks, cax=ax1,pad=-2 , shrink=0.9
#cbar11.ax.set_title('$\mathsf{\Delta SAT_{zonal}}$ ($\mathsf{^{\circ}C}$)',fontsize=10) #,y=-0.2
cbar11.ax.set_xlabel('$\mathsf{\Delta SAT_{zonal}}$ ($\mathsf{^{\circ}C}$)',fontsize=14)
cbar11.ax.tick_params(labelsize=12) #,labeltop=True,labelbottom=False


    
for aa in [ax00,ax10]: #,ax01,ax11,ax21,ax31,ax20,ax30,ax30,ax31
    aa.grid(True,linestyle=':')
    aa.set_xticks(np.append(np.arange(-250,-75+25,25),[-65,-50]))

    aa.set_xlim((-257.5,-57.5))
    
for aa in [ax10]: #,ax11,ax30,ax31,ax30,ax31
    aa.set_ylim((-90,90))
    aa.set_yticks(np.arange(-90,90+30,30))
    
for aa in [ax10]:      #,ax30,ax30
    aa.set_yticklabels(['90$^{\circ}$S','60$^{\circ}$S','30$^{\circ}$S','0$^{\circ}$','30$^{\circ}$N','60$^{\circ}$N','90$^{\circ}$N'])
    #aa.set_ylabel('latitude $\mathsf{\phi\ (^{\circ})}$',fontsize=14)

for aa in [ax00]:     #,ax11ax30,ax31,ax30,ax31

    aa.set_xticklabels(np.append(np.arange(250,75-25,-25),[65,0]))


for aa in [ax00,ax10]:    #,ax01,ax11,ax20,ax30,ax21,ax31,ax31 ,ax30
    aa.tick_params(labelsize=12)

for cbar in [cbar10]: #,cbar30
    #cbar.ax.set_ylabel('$\mathsf{AZ_{SAT,zonal}\ (^{\circ}C)}$',fontsize=14)
    cbar.ax.set_ylabel('$\mathsf{|\Delta SAT_{zonal}|\ (^{\circ}C)}$',fontsize=14)
    cbar.ax.tick_params(labelsize=12)
    


ax00.set_ylim((0.7,2.5))

ax00.set_ylabel('$\mathsf{|\Delta SAT_{zonal}|\ (^{\circ}C)}$',fontsize=14)

ax00.set_xlabel('Age (Ma)',fontsize=14)
    


ax00.set_title(r'$\bf{(d)}$ '+'Zonal SAT contrasts (global average)',loc='left',fontsize=14)
ax10.set_title(r'$\bf{(c)}$ '+'Zonal SAT contrasts (zonal average)',loc='left',fontsize=14)
ax01.set_title(r'$\bf{(a)}$'+' 225Ma', loc='left',fontsize=14)
ax11.set_title(r'$\bf{(b)}$'+' 75Ma', loc='left',fontsize=14)



ax00.legend(loc=(1.0,0.00),frameon=False,labelspacing=0.2,ncol=1,fontsize=8,borderaxespad=0) #loc='upper right'loc=(0.7,0.73)

for aa in [ax00]: #,ax11ax30,ax31,
    jplt.add_geological_timescale(aa,height=0.085,fontsize=8,stagenames=False,textoffset=0.015) #,fontweight='normal'
plt.tight_layout(w_pad=-7,h_pad=-0.75) #w_pad=-4,h_pad=-0.3 w_pad=-4,h_pad=-0.5



fig.savefig('./Fig_5_SAT_zonality.pdf',bbox_inches='tight', dpi=900)














################################################################################
### Figure 6 - Maps with SAT and aridity
################################################################################




interpol_opt=True
plot_opt_legend=False

#fig=plt.figure(figsize=(6,12.2)) 
fig=plt.figure(figsize=(12,7.5)) 
plt.clf()
#gs = GridSpec(5, 1,height_ratios=[1,1,1,1,0.24]) #, width_ratios=[0.5,1,0.15] , height_ratios=[0.4,1]
gs = GridSpec(3, 4,height_ratios=[1,1,0.24], width_ratios=[0.025,1,1,0.025])


ax0 = plt.subplot(gs[0,0:1+1]);
ax0.axis('off')
ax1 = plt.subplot(gs[0,2:3+1]);
ax1.axis('off')
ax2 = plt.subplot(gs[1,0:1+1]);
ax2.axis('off')
ax3 = plt.subplot(gs[1,2:3+1]);
ax3.axis('off')
ax4 = plt.subplot(gs[2,2]);
ax4.axis('off')
    


axs=[ax0,ax1,ax2,ax3] #,ax4
panel_nrs=['(a)','(b)','(c)','(d)']


for ii in range(len(paths_paper_example[:])): #['../c3a_model_input_output/000Ma_280ppm_ice'] : #
    path=paths_paper_example[ii]
    ax=axs[::-1][ii]

    
    age=data_aggregated_history_p2_df.loc[path,['age']][0]
    age_str='{:03d}'.format(age)
    
    
    for pp in [path,'../c3a_model_input_output/'+age_str+'Ma_500ppm']:
        prc_ann_arr=data_aggregated_p2_arrays_df.loc[pp, ['prc_ann_arr']].values[0] 
        rbsr_ann_arr=data_aggregated_p2_arrays_df.loc[pp, ['rbsr_ann']][0]
        prc_ann_intp=jan.interpol_atm_to_ocn_grid(prc_ann_arr) 
        rbsr_ann_intp=jan.interpol_atm_to_ocn_grid(rbsr_ann_arr) 
        aridity_continents=func_aridity_idx(prc_ann_intp,rbsr_ann_intp)
        if pp==path:
            aridity=np.copy(aridity_continents)
        elif pp=='../c3a_model_input_output/'+age_str+'Ma_500ppm':
            aridity_500ppm=np.copy(aridity_continents)

    continental_mask=data_aggregated_geography_df.loc[path,['continental_mask']][0]

    if ii in [0,2]:
        latlabels=[False,True,False,False]
    else:
        latlabels=[False,False,False,False]
    m = Basemap(projection='robin', llcrnrlat=-90, urcrnrlat=90,llcrnrlon=-180, urcrnrlon=180, resolution='c', lon_0=0, ax=ax)
    m.drawmapboundary()  
    m.drawparallels(np.arange(-90.,120.,30.), linewidth=0.5,labels=latlabels)
    m.drawmeridians(np.arange(0.,420.,60.), linewidth=0.5)

    xpam, ypam = m(lonsap_2D, latsap_2D)
    xom, yom = m(lonsot_2D, latsot_2D)
    xpom, ypom = m(lonsop_2D, latsop_2D)
    xtam, ytam = m(lonsat_2D, latsat_2D)


    continental_mask=data_aggregated_geography_df.loc[path,['continental_mask']][0]
    #m.readshapefile('../c3a_meso_reconstructions/plate_polygons/'+age_str+'Ma_reconstruction_polygons','polygons', color='tab:grey', linewidth=1)
    ax.contour(xom,yom,continental_mask,levels=[0.5], linewidths=1.75, colors='k')
    ts_ann_arr=data_aggregated_p2_arrays_df.loc[path, ['ts_ann_arr']].values[0]
#     ts_season_variation=np.max(ts_season_arr,axis=2)-np.min(ts_season_arr,axis=2) 
    ts_ann_arr_intp=jan.interpol_atm_to_ocn_grid(ts_ann_arr)
#     ts_season_variation_intp_masked=np.ma.MaskedArray(ts_season_variation_intp,mask=continental_mask)
    vmin=-10; vmax=35;
    cont_levels_list=[-5,0,5,15,25]
    cont_color_list=['magenta','tab:red','tab:orange','tab:blue','tab:cyan']
    a1=ax.pcolormesh(xpom, ypom, ts_ann_arr_intp, cmap='RdBu_r',vmin=vmin, vmax=vmax, snap=True) 
    cont=ax.contour(xom,yom,ts_ann_arr_intp,levels=cont_levels_list, linewidths=1.75, colors=cont_color_list)

    
    matplotlib.rcParams['hatch.linewidth'] = 1.25
    matplotlib.rcParams['hatch.color']='magenta' #whitesmoke
    #if interpol_opt==True:
    ax.contourf(xom, yom, np.ma.masked_array(aridity,mask=continental_mask), hatches=['///',None], levels=[1.5,100], colors='none', zorder=5) #alpha=0. 
    ax.contourf(xom, yom, np.ma.masked_array(aridity_500ppm,mask=continental_mask), hatches=['\\\\\\',None], levels=[1.5,100], colors='none', zorder=5) #alpha=0.


    
cbar=fig.colorbar(mappable=a1, ax=ax4, boundaries=np.arange(vmin,vmax+1,1), ticks=np.arange(vmin,vmax+5,5), orientation='horizontal',fraction=0.65,pad=-0.0) #, cax=ax1,pad=-2 , shrink=0.9

for ccll in range(len(cont_levels_list)):
    cbar.ax.plot([(cont_levels_list[ccll]-vmin)/(vmax-vmin),(cont_levels_list[ccll]-vmin)/(vmax-vmin)],[0,1],c=cont_color_list[ccll],lw=3)

cbar.ax.set_title('$\mathsf{SAT_{ann} \ (^{\circ}C)}$',fontsize=14)
cbar.ax.tick_params(labelsize=12)


plt.tight_layout(h_pad=-0,w_pad=-1.3)  #w_pad=4 w_pad=-14,h_pad=-2

fig.savefig('./Fig_6_Maps_SAT_Aridity.pdf',bbox_inches='tight', pad_inches=0.0, dpi=900)  














################################################################################
### Figure 7 - Aridity
################################################################################


fig=plt.figure(figsize=(12,6), constrained_layout=True) 
plt.clf()

gs = GridSpec(3, 2,height_ratios=[0.6,0.25,0.6], width_ratios=[1,1])  #[0.7,0.25,0.55]

ax00 = plt.subplot(gs[2,0]);


ax10 = plt.subplot(gs[0:1+1,0]);
ax10_cbar=fig.add_axes([0.49,0.46,0.048,0.45])

ax10_cbar.axis('off')


ax11 = plt.subplot(gs[0:1+1,1]);
ax01 = plt.subplot(gs[2,1],sharex=ax11);


for aa in [ax01,ax11]: #,ax21,ax20,ax21,ax20 ,ax21 ,ax10 ,ax11
    aa.tick_params(right=True,labelright=True,left=False,labelleft=False)
    aa.yaxis.set_label_position("right")

vmin=0; vmax=7.0;
a10=ax10.pcolor(timeslices_pcolor2D[:,1:],lats2D[:,1:],data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm',['aridity_1p5_zonalarea']][0][:,:],cmap='RdBu_r',zorder=1,vmin=vmin,vmax=vmax) #,vmax=25
cbar10=plt.colorbar(a10, ax=ax10_cbar, orientation='vertical', boundaries=np.arange(vmin,vmax+0.5,0.5),ticks=np.arange(vmin,vmax+1,1),fraction=1) #,extend='max' extend=extend, shrink=0.621, pad=0.02


cont_color_list=['None','tab:red','tab:blue'] #,'w'
cont_levels_list=[0,2,4]


for aa in [ax00,ax01]:
    if aa==ax00:
        var_plot='aridity_1p5_globarea'
        var_fac=100.
    if aa==ax01:
        var_plot='aridity_continents_globalmean'
        var_fac=1.
    aa.plot(-np.asarray(data_aggregated_history_p2_df.loc[paths_pCO2COPSE[:],['age']].values,dtype=float),var_fac*data_aggregated_SeasonalityZonality_df.loc['paths_pCO2COPSE',[var_plot]][0][:],\
        ':',c='k',markersize=5, label='pCO2_COPSE',lw=2,zorder=3) #,markeredgecolor='k'
    aa.plot(-np.asarray(data_aggregated_history_p2_df.loc[paths_pCO2Foster[:],['age']].values,dtype=float),var_fac*data_aggregated_SeasonalityZonality_df.loc['paths_pCO2Foster',[var_plot]][0][:],\
        '--',c='k',markersize=5, label='pCO2_proxy',lw=2,zorder=3) #,markeredgecolor='k'
    aa.plot(-np.asarray(data_aggregated_history_p2_df.loc[paths_1000ppm[:],['age']].values,dtype=float),var_fac*data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm',[var_plot]][0][:],\
        '-',c='k',markersize=5, label='pCO2_1000ppm',lw=2,zorder=3) #,markeredgecolor='k'
    aa.plot(-np.asarray(data_aggregated_history_p2_df.loc[paths_2000ppm[:],['age']].values,dtype=float),var_fac*data_aggregated_SeasonalityZonality_df.loc['paths_2000ppm',[var_plot]][0][:],\
        c='None',marker='v',markersize=5,markerfacecolor='r',markeredgecolor='k', label='2000ppm',linestyle='None',zorder=3)
    aa.plot(-np.asarray(data_aggregated_history_p2_df.loc[paths_500ppm[:],['age']].values,dtype=float),var_fac*data_aggregated_SeasonalityZonality_df.loc['paths_500ppm',[var_plot]][0][:],\
        c='None',marker='^',markersize=5,markerfacecolor='cyan',markeredgecolor='k', label='500ppm',linestyle='None',zorder=3)
    aa.plot(-np.asarray(data_aggregated_history_p2_df.loc[paths_1000ppm_S0ini[:],['age']].values,dtype=float),var_fac*data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm_S0ini',[var_plot]][0][:],\
        c='tab:orange',markersize=5, label='S0ini',lw=1,zorder=3) #,markeredgecolor='k'
    aa.plot(-np.asarray(data_aggregated_history_p2_df.loc[paths_1000ppm_S0ini_VegFix[:],['age']].values,dtype=float),var_fac*data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm_S0ini_VegFix',[var_plot]][0][:],\
        c='limegreen',markersize=5, label='VegFix_S0ini',lw=1,zorder=3) #,markeredgecolor='k'
    aa.plot(-np.asarray(data_aggregated_history_p2_df.loc[paths_1000ppm_VegHom[:],['age']].values,dtype=float),var_fac*data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm_VegHom',[var_plot]][0][:],\
        '-',c='tab:blue',marker='',markersize=5, label='VegHom',lw=1,zorder=3)

    for pp in range(len(paths_VegHom_SL)):
        path=paths_VegHom_SL[pp]
        if path not in ['../c3a_model_input_output/200Ma_1000ppm_VegHom_SL+80']:
            if pp==0:
                label='SL'
            else:
                label=''
            age=-data_aggregated_history_p2_df.loc[path,['age']][0]    
            aridity_continental_globalmean=var_fac*data_aggregated_SeasonalityZonality_df.loc['paths_VegHom_SL',[var_plot]][0][pp]
            #aridity_continental_globmean=data_aggregated_history_p2_df.loc[pp,['aridity_continents_globmean']][0]
            aa.scatter(age,aridity_continental_globalmean,c='None',marker='o',edgecolors='tab:blue',linewidths=2.5,zorder=4,label=label)
            #ax00.annotate(path[path.index('SL'):]+'m',xy=(age+3,aridity_continental_globalmean),color='k',fontsize=8)
            
    for pp in range(len(paths_1000ppm_orb)):
            path=paths_1000ppm_orb[pp]
            if path==paths_1000ppm_orb_22p0[0]:
                label='orb, 22.0'            
            elif path==paths_1000ppm_orb_24p5[0]:
                label='orb, 24.5'
            else:
                label=''
            if path in paths_1000ppm_orb_22p0:
                marker='s'

            elif path in paths_1000ppm_orb_24p5:
                marker='D'
            else:
                marker='+'
            age=-data_aggregated_history_p2_df.loc[path,['age']][0]   
            aridity_continental_globalmean=var_fac*data_aggregated_SeasonalityZonality_df.loc['paths_1000ppm_orb',[var_plot]][0][pp]

paths_pp=paths_1000ppm
evap_NrCell_perlats=np.ma.empty((48,len(paths_pp)))

    
for aa in [ax00,ax10,ax11,ax01]: #,ax21,ax30,ax31
    aa.grid(True,linestyle=':')
    aa.set_xticks(np.append(np.arange(-250,-75+25,25),[-65,-50]))
    #aa.set_xlim((-257.5,-47.5))
    aa.set_xlim((-257.5,-57.5))
    
for aa in [ax10]: #,ax30,ax31
    aa.set_ylim((-100,90))
    aa.set_yticks(np.arange(-90,90+30,30))
    
for aa in [ax10]:      #,ax30
    aa.set_yticklabels(['90$^{\circ}$S','60$^{\circ}$S','30$^{\circ}$S','0$^{\circ}$','30$^{\circ}$N','60$^{\circ}$N','90$^{\circ}$N'])

for aa in [ax00,ax01]:     #ax30,ax31
     aa.set_xticklabels(np.append(np.arange(250,75-25,-25),[65,0]))


for aa in [ax00,ax10,ax01,ax11]:    #,ax21,ax31 ,ax30
    aa.tick_params(labelsize=12)

for cbar in [cbar10]:
    cbar.ax.set_ylabel('$\mathsf{A_{land,dry}}$ ($\mathsf{10^6\,km^2}$)',fontsize=14)
    #cbar.ax.set_ylabel('AI (a.u.)',fontsize=14)
    cbar.ax.tick_params(labelsize=12)

#ax01.set_ylim(7,40)
ax11.set_ylim((-1.5,23))
ax11.set_ylabel('Evaporite Volume ($\mathsf{10^3\,km^3}$)',fontsize=14)

ax01.set_ylim((0.91,1.75))

ax00.set_ylabel('$\mathsf{f_{land,dry}}$ (%)',fontsize=14)

ax00.set_xlabel('Age (Ma)',fontsize=14)
ax01.set_xlabel('Age (Ma)',fontsize=14)
ax01.set_ylabel('AI (a.u.)',fontsize=14)

for aa in [ax10,ax11]:
    aa.tick_params(labelbottom=False)


ax00.set_ylim((16,60))
ax00.set_yticks([20,30,40,50,60])



    


ax10.set_title(r'$\bf{(a)}$ '+'Dry land area (zonal sum)', loc='left',fontsize=14)
ax00.set_title(r'$\bf{(b)}$ '+'Global dry land area fraction', loc='left',fontsize=14)
ax01.set_title(r'$\bf{(d)}$ '+'Global aridity index',loc='left',fontsize=14)
ax11.set_title(r'$\bf{(c)}$ '+'Coal/evaporite accumulation',loc='left',fontsize=14)


for aa in [ax00]:
    jplt.add_geological_timescale(aa,height=0.085,fontsize=8,stagenames=False,textoffset=0.015) #,fontweight='normal'
jplt.add_geological_timescale(ax01,height=0.085,fontsize=8,stagenames=False,textoffset=0.015)


fig.savefig('./Fig_7_aridity.pdf',bbox_inches='tight', dpi=900)