#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Oct 18 13:58:07 2021

@author: jene
"""

import glob
# import os
# import datetime as dt
import xarray as xr
# import netCDF4 as nc
# from netCDF4 import MFDataset
import numpy as np
import pickle
from matplotlib import pyplot as plt
import matplotlib as mpl
import matplotlib.path as mpath
import cartopy.crs as ccrs
import cartopy.feature as cfeature
# import warnings
# import scipy.io

def cm_to_inch(value):
    return value/2.54

def calc_spatial_mean(xr_da, weights, lon_name="lon", lat_name="lat"):
    """
    Calculate spatial mean of xarray.DataArray with grid cell weighting.
    
    Parameters
    ----------
    xr_da: xarray.DataArray
        Data to average
    lon_name: str, optional
        Name of x-coordinate
    lat_name: str, optional
        Name of y-coordinate
    radius: float
        Radius of the planet [metres], currently assumed spherical (not important anyway)
    Returns
    -------
    Spatially averaged xarray.DataArray.
    """

    return (xr_da*weights).sum(dim=[lon_name, lat_name],skipna=True)/weights.sum(skipna=True)


def calc_spatial_integral(xr_da, weights, lon_name="lon", lat_name="lat"):
    """
    Calculate spatial integral of xarray.DataArray with grid cell weighting.
    
    Parameters
    ----------
    xr_da: xarray.DataArray
        Data to average
    lon_name: str, optional
        Name of x-coordinate
    lat_name: str, optional
        Name of y-coordinate
    radius: float
        Radius of the planet [metres], currently assumed spherical (not important anyway)
    Returns
    -------
    Spatially averaged xarray.DataArray.
    """
    return (xr_da * weights).sum(dim=[lon_name, lat_name], skipna=True)

def load_data(var, s):
    fname = glob.glob(dataPath + '' + var + '*' + s + '*.dat')[0]
    data = pickle.load(open(fname, 'rb'))
    grid = pickle.load(open(dataPath + 'grid.dat', 'rb'))
    mask = pickle.load(open(dataPath + 'mask.dat', 'rb'))
    grid = grid.where(((grid.lat>60) & (mask>0)))
    # grid.values[np.isnan(grid)] = 0
    
    # attrs = data.attrs
    # data = xr.open_mfdataset('netCDF/' + var + '_' + s + '_*.nc')
    # if var == 'ph':
    #     data = data['phlvl']
    # else:
    #     data = data[var]
    # data.attrs = attrs
    
    
    data = data.where(((data.lat>60) & (mask>0)))
    data = data.groupby('time.year').mean('time')
    # data = data.where(data.year<=2100, drop=True)
    if var == 'fgco2':
        ts = calc_spatial_integral(data, grid, 'x', 'y')
        ts = ts/1e9*(60*60*24*365.25)
    else:
        ts = calc_spatial_mean(data, grid, 'x', 'y')
    return data, ts

def plot_map(ax, data, T, diff = False):
    cmap = mpl.cm.Spectral_r
    # cmap = mpl.cm.bwr
    if data.attrs['var'] == 'sst':
        bounds = np.arange(-5,5.1,.5)
        ticks = np.arange(-5,5.1,2.5)
    if data.attrs['var'] == 'sss':
        bounds = np.arange(-3,3.1,.3)
        ticks = np.arange(-3,3.1,1.5)
    if data.attrs['var'] == 'ph':
        bounds = np.arange(-.5,.51,.05)
        ticks = np.arange(-.5,.51,.25)
    if data.attrs['var'] == 'fgco2':
        # bounds = np.arange(-.5,.51,.1)
        bounds = np.arange(-4,4.1,.5)
        ticks = np.arange(-4,4.1,2)
        data = data/(0.012)*(60*60*24*365.25)
    norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend='neither')
    presData = data.where((data['year']>pres[0]) & (data['year']<pres[1])).mean('year')
    futData= data.where((data['year']>fut[0]) & (data['year']<fut[1])).mean('year')
    timeDiff = futData - presData
    im = ax.pcolormesh(data.lon, data.lat,timeDiff,
                       transform=ccrs.PlateCarree(),
                        cmap=cmap, norm=norm
                       )
    ax.set_extent([-180, 180, 60, 90], ccrs.PlateCarree())
    ax.add_feature(cfeature.LAND, color='lightgray', edgecolor='k')
    ax.coastlines()
    ax.gridlines(draw_labels=False,
              linewidth=.5, color='k', alpha=.5, linestyle=':')
    theta = np.linspace(0, 2*np.pi, 100)
    center, radius = [0.5, 0.5], 0.5
    verts = np.vstack([np.sin(theta), np.cos(theta)]).T
    circle = mpath.Path(verts * radius + center)

    ax.set_boundary(circle, transform=ax.transAxes)
    ax.set_title(T, fontsize=fontsize['tit'])
    return im,ticks

def plot_diff2ref_map(ax, dataRef, data, T):
    cmap = mpl.cm.bwr
    if data.attrs['var'] == 'sst':
        bounds = np.arange(-1.5,1.51,.15)
        ticks = np.arange(-1.5,1.51,.5)
    if data.attrs['var'] == 'sss':
        bounds = np.arange(-1,1.1,.1)
        ticks = np.arange(-.9,.91,.3)
    if data.attrs['var'] == 'ph':
        bounds = np.arange(-.05,.051,.005)
        ticks = np.arange(-.05,.051,.025)
    if data.attrs['var'] == 'fgco2':
        # bounds = np.arange(-.5,.51,.1)
        bounds = np.arange(-1.5,1.51,.15)
        ticks = np.arange(-1.5,1.51,.5)
        data = data/(0.012)*(60*60*24*365.25)
        dataRef = dataRef/(0.012)*(60*60*24*365.25)
    norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend='neither')
    presData = data.where((data['year']>pres[0]) & (data['year']<pres[1])).mean('year')
    futData= data.where((data['year']>fut[0]) & (data['year']<fut[1])).mean('year')
    # timeDiff = futData - presData
    presRefData = dataRef.where((dataRef['year']>pres[0]) & (dataRef['year']<pres[1])).mean('year')
    futRefData= dataRef.where((dataRef['year']>fut[0]) & (dataRef['year']<fut[1])).mean('year')
    # refDiff = futRefData - presRefData
    timeDiff = (futData - presData) - (futRefData - presRefData)
    im = ax.pcolormesh(data.lon, data.lat, timeDiff,
                       transform=ccrs.PlateCarree(),
                        cmap=cmap, norm=norm
                       )
    ax.set_extent([-180, 180, 60, 90], ccrs.PlateCarree())
    ax.add_feature(cfeature.LAND, color='lightgray', edgecolor='k')
    ax.coastlines()
    ax.gridlines(draw_labels=False,
              linewidth=.5, color='k', alpha=.5, linestyle=':')
    theta = np.linspace(0, 2*np.pi, 100)
    center, radius = [0.5, 0.5], 0.5
    verts = np.vstack([np.sin(theta), np.cos(theta)]).T
    circle = mpath.Path(verts * radius + center)

    ax.set_boundary(circle, transform=ax.transAxes)
    ax.set_title(T, fontsize=fontsize['tit'])
    return im,ticks

def plot_20yMap(ax, data, T):
        cmap = mpl.cm.Spectral_r
        if data.attrs['var'] == 'sst':
            bounds = np.arange(-2,8.1,1)
        if data.attrs['var'] == 'sss':
            bounds = np.arange(25,35.1,1)
        if data.attrs['var'] == 'ph':
            bounds = np.arange(7.5,8.01,.05)
        if data.attrs['var'] == 'fgco2':
            bounds = np.arange(-6,0.01,.6)
            data = data/(0.012)*(60*60*24*365.25)
        norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend='both')
        im = ax.pcolormesh(data.lon, data.lat,
                           data.where((data['year']>=mapsYr[0]) & (data['year']<=mapsYr[1])).mean('year'),
                           transform=ccrs.PlateCarree(),
                            cmap=cmap, norm=norm
                           )
        ax.set_extent([-180, 180, 60, 90], ccrs.PlateCarree())
        ax.add_feature(cfeature.LAND, color='lightgray', edgecolor='k')
        ax.coastlines()
        ax.gridlines(draw_labels=False,
                  linewidth=.5, color='k', alpha=.5, linestyle=':')
        theta = np.linspace(0, 2*np.pi, 100)
        center, radius = [0.5, 0.5], 0.5
        verts = np.vstack([np.sin(theta), np.cos(theta)]).T
        circle = mpath.Path(verts * radius + center)
    
        ax.set_boundary(circle, transform=ax.transAxes)
        ax.set_title(T, fontsize=fontsize['tit'])
        return im

fontsize = {'supt':18, 'tit':16, 'axtit':13, 'axtick':11}
figNum =['a) ','b) ','c) ','d) ','e) ','f) ','g) ','h) ','i) ','j) ','k) ',
         'l) ','m) ','n) ','o) ','p) ','q) ','r) ','s) ','t) ','u) ','v) ',
         'w) ','x) ','y) ','z) ', r'$\alpha$)', r'$\beta)$', r'$\gamma)$', r'$\delta)$']
figExt = '.png'

basePath = '/media/ipcc/'
dataPath = basePath + 'NS9252K/KeyCLIM_CMOR/'
dataPath = basePath + 'NS9252K/cmorout/'
dataPath = 'Data/'
refDataPath = basePath + 'NS9034K/CMIP6/'
variable = [
    # 'sst',
    'fgco2',
    # 'sss',
    # 'ph'
    ]
scenario = [
    'cloud2',
    'eddy',
    'iceSheet3',
    'snow'
    ]
grid = pickle.load(open(dataPath + 'grid.dat', 'rb'))

fut = [2071,2100]
pres = [1985,2014]
# mapsYr = [2080,2100]
mapsYr = pres

# ###############################################################################
# ###     Timeseries
# ###############################################################################
# idx = 1
# fig = plt.figure(figsize=(cm_to_inch(29.7),cm_to_inch(21)), dpi=300)
# for var in variable:
#     dataRef, tsRef = load_data(var, 'ref')
#     ax = fig.add_subplot(2,2,idx)
#     ax.plot(tsRef.year, tsRef, label='Reference')
#     for s in scenario:
#         print('Processing ' + var + ', ' + s)
#         data, ts = load_data(var, s)
#         if s == 'ref':
#             ls = 'Reference'
#         if s == 'cloud2':
#             ls = 'Cloud'
#         if s == 'iceSheet3':
#             ls = 'IceSheet'
#         if s == 'snow':
#             ls = 'Snow'
#         if s == 'eddy':
#             ls = 'Eddy'
#         # if data.attrs['var'] == 'fgco2':
#         #     ts = ts/1e9*(60*60*24*365.25)
#         #     tsRef = tsRef/1e9*(60*60*24*365.25)
#         ax.plot(ts.year, ts, label=ls)
#     if data.attrs['var'] == 'fgco2':
#         ax.set_ylabel('CO$_2$ flux [Tg yr${-1}$]', fontsize=fontsize['axtit'])
#         T = 'CO$_2$ flux'
#     if data.attrs['var'] == 'sst':
#         ax.set_ylabel('SST [degC]', fontsize=fontsize['axtit'])
#         T = 'SST'
#     if data.attrs['var'] == 'sss':
#         ax.set_ylabel('SSS [psu]', fontsize=fontsize['axtit'])
#         T = 'SSS'
#     if data.attrs['var'] == 'ph':
#         ax.set_ylabel('pH', fontsize=fontsize['axtit'])
#         T = 'pH'
#     ax.legend()
#     ax.set_xlim([1850,2100])
#     ax.grid(linestyle=':')
#     ax.set_title(figNum[idx-1] + ' ' + T, fontsize=fontsize['tit'])
#     if idx> 2:
#         ax.set_xlabel('Year',fontsize=fontsize['axtit'])
#     idx+=1
# plt.subplots_adjust(wspace=.2,hspace=.3)
# plt.savefig('Plots/timeseries' + figExt)
# plt.show()
# # plt.close()

# # # ###############################################################################
# # # ###     Maps
# # # ###############################################################################
# for var in variable:
#     fig = plt.figure(figsize=(cm_to_inch(29.7),cm_to_inch(21)), dpi=300)
#     gs = fig.add_gridspec(2,6)
#     dataRef, tsRef = load_data(var, 'ref')
#     # ax1 = fig.add_subplot(3,2,1)
#     im = plot_map(fig.add_subplot(gs[0, 0:2],projection=ccrs.NorthPolarStereo()),
#                   dataRef, figNum[0] + ' Reference')
#     idx = 1
#     for s in scenario:
#         print('Processing ' + var + ', ' + s)
#         data, ts = load_data(var, s)
#         if 'ref' in s :
#             ls = 'Reference'
#         if s == 'cloud2':
#             ls = 'Cloud'
#         if s == 'iceSheet3':
#             ls = 'IceSheet'
#         if s == 'snow':
#             ls = 'Snow'
#         if s == 'eddy':
#             ls = 'Eddy'
#         if idx>2:
#             plot_map(fig.add_subplot(gs[1,(idx-3)*3:(idx-2)*3], projection=ccrs.NorthPolarStereo()),
#                   data, figNum[idx] + ' ' + ls)
#         else:
#             plot_map(fig.add_subplot(gs[0,(idx)*2:(idx+1)*2], projection=ccrs.NorthPolarStereo()),
#                   data, figNum[idx] + ' ' + ls)
#         # cbar.set_label('mean(' + variable + ') [' + unit + ']', fontsize=fontsize['axtit'])#, rotation=270, labelpad=+10)
#         # cbar.ax.tick_params(labelsize=fontsize['axtick'])
#         # plt.savefig('Plot/mean_' + Z[3:] + '.png')
#         idx+=1
#     if data.attrs['var'] == 'fgco2':
#         # fig.suptitle(data.attrs['var'],
#         #              fontsize=fontsize['supt'],fontweight="bold")
#         fig.suptitle('$\Delta$CO$_2$ flux', fontsize=fontsize['supt'],fontweight="bold", y=0.95)
#         cbarT = 'CO$_2$ flux [mol C m${-2}$ yr${-1}$]'
#     if data.attrs['var'] == 'sst':
#         fig.suptitle('$\Delta$SST', fontsize=fontsize['supt'],fontweight="bold", y=0.95)
#         cbarT = 'SST [degC]'
#     if data.attrs['var'] == 'sss':
#         fig.suptitle('$\Delta$SSS', fontsize=fontsize['supt'],fontweight="bold", y=0.95)
#         cbarT = 'SSS [psu]'
#     if data.attrs['var'] == 'ph':
#         fig.suptitle('$\Delta$pH', fontsize=fontsize['supt'],fontweight="bold", y=0.95)
#         cbarT = 'pH'
#     cbar_ax = fig.add_axes([0.25, 0.05, 0.5, 0.025])
#     cbar = fig.colorbar(im, orientation="horizontal", cax=cbar_ax)
#     cbar.set_label('$\Delta$' + cbarT, fontsize=fontsize['axtit'])#, rotation=270, labelpad=+10)
#     cbar.ax.tick_params(labelsize=fontsize['axtick'])
#     plt.subplots_adjust(wspace=.1, hspace=.1)
#     plt.savefig('Plots/Maps_' + var + figExt, bbox_inches='tight')
#     plt.show()
#     plt.close()
    

# for var in variable:
#     fig = plt.figure(figsize=(cm_to_inch(29.7),cm_to_inch(21)), dpi=300)
#     gs = fig.add_gridspec(2,6)
#     dataRef, tsRef = load_data(var, 'ref')
#     # ax1 = fig.add_subplot(3,2,1)
#     im = plot_20yMap(fig.add_subplot(gs[0, 0:2],projection=ccrs.NorthPolarStereo()),
#                   dataRef, figNum[0] + ' Reference')
#     idx = 1
#     for s in scenario:
#         print('Processing ' + var + ', ' + s)
#         data, ts = load_data(var, s)
#         # attrs = data.attrs
#         # data = data+dataRef
#         # data.attrs = attrs
#         if 'ref' in s :
#             ls = 'Reference'
#         if s == 'cloud2':
#             ls = 'Cloud'
#         if s == 'iceSheet3':
#             ls = 'IceSheet'
#         if s == 'snow':
#             ls = 'Snow'
#         if s == 'eddy':
#             ls = 'Eddy'
#         if idx>2:
#             plot_20yMap(fig.add_subplot(gs[1,(idx-3)*3:(idx-2)*3], projection=ccrs.NorthPolarStereo()),
#                   data, figNum[idx] + ' ' + ls)
#         else:
#             plot_20yMap(fig.add_subplot(gs[0,(idx)*2:(idx+1)*2], projection=ccrs.NorthPolarStereo()),
#                   data, figNum[idx] + ' ' + ls)
#         # cbar.set_label('mean(' + variable + ') [' + unit + ']', fontsize=fontsize['axtit'])#, rotation=270, labelpad=+10)
#         # cbar.ax.tick_params(labelsize=fontsize['axtick'])
#         # plt.savefig('Plot/mean_' + Z[3:] + '.png')
#         idx+=1
#     if data.attrs['var'] == 'fgco2':
#         # fig.suptitle(data.attrs['var'],
#         #              fontsize=fontsize['supt'],fontweight="bold")
#         fig.suptitle('CO$_2$ flux (' + str(mapsYr[0]) + '-' + str(mapsYr[1]) + ')',
#                       fontsize=fontsize['supt'],fontweight="bold", y=0.95)
#         cbarT = 'CO$_2$ flux [mol C m-2 yr-1]'
#     if data.attrs['var'] == 'sst':
#         fig.suptitle('SST (' + str(mapsYr[0]) + '-' + str(mapsYr[1]) + ')',
#                       fontsize=fontsize['supt'],fontweight="bold", y=0.95)
#         cbarT = 'SST [degC]'
#     if data.attrs['var'] == 'sss':
#         fig.suptitle('SSS (' + str(mapsYr[0]) + '-' + str(mapsYr[1]) + ')',
#                       fontsize=fontsize['supt'],fontweight="bold", y=0.95)
#         cbarT = 'SSS [psu]'
#     if data.attrs['var'] == 'ph':
#         fig.suptitle('pH (' + str(mapsYr[0]) + '-' + str(mapsYr[1]) + ')',
#                       fontsize=fontsize['supt'],fontweight="bold", y=0.95)
#         cbarT = 'pH'
#     cbar_ax = fig.add_axes([0.25, 0.05, 0.5, 0.025])
#     cbar = fig.colorbar(im, orientation="horizontal", cax=cbar_ax)
#     cbar.set_label('' + cbarT, fontsize=fontsize['axtit'])#, rotation=270, labelpad=+10)
#     cbar.ax.tick_params(labelsize=fontsize['axtick'])
#     plt.subplots_adjust(wspace=.1, hspace=.1)
#     plt.savefig('Plots/Maps_' + str(mapsYr[0]) + '-' + str(mapsYr[1]) + '_' + var + figExt, bbox_inches='tight')
#     plt.show()
#     plt.close()
    
# # ###############################################################################
# # ###     Maps
# # ###############################################################################
idxv = 0
fig0 = plt.figure(figsize=(cm_to_inch(28),cm_to_inch(21)), dpi=300)
for var in variable:
    idx = 1
    fig1 = plt.figure(figsize=(cm_to_inch(28),cm_to_inch(21)), dpi=300)
    gs1 = fig1.add_gridspec(2,3)
    # fig2 = plt.figure(figsize=(cm_to_inch(29.7),cm_to_inch(21)), dpi=300)
    # gs2 = fig2.add_gridspec(2,3)
    dataRef, tsRef = load_data(var, 'ref')
    if dataRef.attrs['var'] == 'fgco2':
        # ax0.set_ylabel('CO$_2$ flux [Tg yr$^{-1}$]', fontsize=fontsize['axtit'])
        T = 'CO$_2$ flux'
        fig1.suptitle('$\Delta$CO$_2$ flux (' + str(fut[0]) + '-' + str(fut[1])
                      + ' minus '+ str(pres[0]) + '-' + str(pres[1])
                      + ')', fontsize=fontsize['supt'],fontweight="bold", y=0.95)
        cbarT1 = 'CO$_2$ flux [mol C m$^{-2}$ yr$^{-1}$]'
        # fig2.suptitle('CO$_2$ flux (' + str(mapsYr[0]) + '-' + str(mapsYr[1]) + ')',
        #               fontsize=fontsize['supt'],fontweight="bold", y=0.95)
        cbarT2 = 'CO$_2$ flux [mol C $^{-2}$ yr$^{-1}$]'
    if dataRef.attrs['var'] == 'sst':
        # ax0.set_ylabel('SST [degC]', fontsize=fontsize['axtit'])
        T = 'SST'
        fig1.suptitle('$\Delta$SST (' + str(fut[0]) + '-' + str(fut[1])
                      + ' minus '+ str(pres[0]) + '-' + str(pres[1])
                      + ')', fontsize=fontsize['supt'],fontweight="bold", y=0.95)
        cbarT1 = 'SST [degC]'
        # fig2.suptitle('SST (' + str(mapsYr[0]) + '-' + str(mapsYr[1]) + ')',
        #               fontsize=fontsize['supt'],fontweight="bold", y=0.95)
        cbarT2 = 'SST [degC]'
    if dataRef.attrs['var'] == 'sss':
        # ax0.set_ylabel('SSS [psu]', fontsize=fontsize['axtit'])
        T = 'SSS'
        fig1.suptitle('$\Delta$SSS (' + str(fut[0]) + '-' + str(fut[1])
                      + ' minus '+ str(pres[0]) + '-' + str(pres[1])
                      + ')', fontsize=fontsize['supt'],fontweight="bold", y=0.95)
        cbarT1 = 'SSS [psu]'
        # fig2.suptitle('SSS (' + str(mapsYr[0]) + '-' + str(mapsYr[1]) + ')',
        #               fontsize=fontsize['supt'],fontweight="bold", y=0.95)
        cbarT2 = 'SSS [psu]'
    if dataRef.attrs['var'] == 'ph':
        # ax0.set_ylabel('pH', fontsize=fontsize['axtit'])
        T = 'pH'
        fig1.suptitle('$\Delta$pH (' + str(fut[0]) + '-' + str(fut[1])
                      + ' minus '+ str(pres[0]) + '-' + str(pres[1])
                      + ')', fontsize=fontsize['supt'],fontweight="bold", y=0.95)
        cbarT1 = 'pH'
        # fig2.suptitle('pH (' + str(mapsYr[0]) + '-' + str(mapsYr[1]) + ')',
        #               fontsize=fontsize['supt'],fontweight="bold", y=0.95)
        cbarT2 = 'pH'
    # ax0 = fig0.add_subplot(2,2,idxv+1)
    # ax0.plot(tsRef.year, tsRef, label='Reference')
    ax = fig1.add_subplot(gs1[0, 0],projection=ccrs.NorthPolarStereo())
    ax.set_anchor('C')
    im, ticks = plot_map(ax, dataRef, figNum[0] + ' Reference')
    # cbar_ax = fig1.add_subplot(gs1[0, 0])
    cbar_ax = fig1.add_axes([1/6-.05, 0.5, .25, 0.025])
    cbar_ax.set_anchor('C')
    cbar = fig1.colorbar(im, orientation="horizontal", cax=cbar_ax, ticks = ticks)
    cbar.set_label('$\Delta$' + cbarT1, fontsize=fontsize['axtit'])#, rotation=270, labelpad=+10)
    cbar.ax.tick_params(labelsize=fontsize['axtick'])
    # im = plot_20yMap(fig2.add_subplot(gs2[0, 0:2],projection=ccrs.NorthPolarStereo()),
    #               dataRef, figNum[0] + ' Reference')
    for s in scenario:
        print('Processing ' + var + ', ' + s)
        data, ts = load_data(var, s)
        if 'ref' in s :
            ls = 'Reference'
        if s == 'cloud2':
            ls = 'Cloud'
        if s == 'iceSheet3':
            ls = 'IceSheet'
        if s == 'snow':
            ls = 'Snow'
        if s == 'eddy':
            ls = 'Eddy'
        if idx>2:
            ax = fig1.add_subplot(gs1[1,idx-2], projection=ccrs.NorthPolarStereo())
            # im2 = plot_20yMap(fig2.add_subplot(gs2[1,(idx-3)*3:(idx-2)*3], projection=ccrs.NorthPolarStereo()),
            #       data, figNum[idx] + ' ' + ls)
        else:
            ax = fig1.add_subplot(gs1[0,idx], projection=ccrs.NorthPolarStereo())
            # plot_20yMap(fig2.add_subplot(gs2[0,(idx)*2:(idx+1)*2], projection=ccrs.NorthPolarStereo()),
            #       data, figNum[idx] + ' ' + ls)
        # ax0.plot(ts.year, ts, label=ls)
        ax.set_anchor('C')
        im1, ticks = plot_diff2ref_map(ax, dataRef, data, figNum[idx] + ' ' + ls + ' - Ref.')
        idx+=1
    if data.attrs['var'] == 'fgco2':
        # ax0.set_ylabel('CO$_2$ flux [Tg yr$^{-1}$]', fontsize=fontsize['axtit'])
        T = 'CO$_2$ flux'
        fig1.suptitle('$\Delta$CO$_2$ flux (' + str(fut[0]) + '-' + str(fut[1])
                      + ' minus '+ str(pres[0]) + '-' + str(pres[1])
                      + ')', fontsize=fontsize['supt'],fontweight="bold", y=0.95)
        cbarT1 = 'CO$_2$ flux [mol C m$^{-2}$ yr$^{-1}$]'
        # fig2.suptitle('CO$_2$ flux (' + str(mapsYr[0]) + '-' + str(mapsYr[1]) + ')',
        #               fontsize=fontsize['supt'],fontweight="bold", y=0.95)
        cbarT2 = 'CO$_2$ flux [mol C $^{-2}$ yr$^{-1}$]'
    if data.attrs['var'] == 'sst':
        # ax0.set_ylabel('SST [degC]', fontsize=fontsize['axtit'])
        T = 'SST'
        fig1.suptitle('$\Delta$SST (' + str(fut[0]) + '-' + str(fut[1])
                      + ' minus '+ str(pres[0]) + '-' + str(pres[1])
                      + ')', fontsize=fontsize['supt'],fontweight="bold", y=0.95)
        cbarT1 = 'SST [degC]'
        # fig2.suptitle('SST (' + str(mapsYr[0]) + '-' + str(mapsYr[1]) + ')',
        #               fontsize=fontsize['supt'],fontweight="bold", y=0.95)
        cbarT2 = 'SST [degC]'
    if data.attrs['var'] == 'sss':
        # ax0.set_ylabel('SSS [psu]', fontsize=fontsize['axtit'])
        T = 'SSS'
        fig1.suptitle('$\Delta$SSS (' + str(fut[0]) + '-' + str(fut[1])
                      + ' minus '+ str(pres[0]) + '-' + str(pres[1])
                      + ')', fontsize=fontsize['supt'],fontweight="bold", y=0.95)
        cbarT1 = 'SSS [psu]'
        # fig2.suptitle('SSS (' + str(mapsYr[0]) + '-' + str(mapsYr[1]) + ')',
        #               fontsize=fontsize['supt'],fontweight="bold", y=0.95)
        cbarT2 = 'SSS [psu]'
    if data.attrs['var'] == 'ph':
        # ax0.set_ylabel('pH', fontsize=fontsize['axtit'])
        T = 'pH'
        fig1.suptitle('$\Delta$pH (' + str(fut[0]) + '-' + str(fut[1])
                      + ' minus '+ str(pres[0]) + '-' + str(pres[1])
                      + ')', fontsize=fontsize['supt'],fontweight="bold", y=0.95)
        cbarT1 = 'pH'
        # fig2.suptitle('pH (' + str(mapsYr[0]) + '-' + str(mapsYr[1]) + ')',
        #               fontsize=fontsize['supt'],fontweight="bold", y=0.95)
        cbarT2 = 'pH'
    # ax0.legend()
    # ax0.set_xlim([1850,2100])
    # ax0.grid(linestyle=':')
    # ax0.set_title(figNum[idxv] + ' ' + T, fontsize=fontsize['tit'])
    # if idxv> 1:
    #     ax0.set_xlabel('Year',fontsize=fontsize['axtit'])
    idxv+=1
        
        
    # cbar_ax = fig1.add_axes([0.25, 0.05, 0.5, 0.025])
    cbar_ax = fig1.add_axes([2/3-.25, 0.05, .45, 0.025])
    cbar_ax.set_anchor('C')
    cbar = fig1.colorbar(im1, orientation="horizontal", cax=cbar_ax, ticks = ticks)
    cbar.set_label('$\Delta$' + cbarT1, fontsize=fontsize['axtit'])#, rotation=270, labelpad=+10)
    cbar.ax.tick_params(labelsize=fontsize['axtick'])
    
    # cbar_ax = fig2.add_axes([0.25, 0.05, 0.5, 0.025])
    # cbar = fig2.colorbar(im2, orientation="horizontal", cax=cbar_ax)
    # cbar.set_label('' + cbarT2, fontsize=fontsize['axtit'])#, rotation=270, labelpad=+10)
    # cbar.ax.tick_params(labelsize=fontsize['axtick'])

    # fig1.subplots_adjust(wspace=.1, hspace=.1)
    fig1.savefig('Plots/Maps_' + var + figExt, bbox_inches='tight')
    fig1.show()
    plt.close(fig1)

    # fig2.subplots_adjust(wspace=.1, hspace=.1)
    # fig2.savefig('Plots/Maps_' + str(mapsYr[0]) + '-' + str(mapsYr[1]) + '_' + var + figExt, bbox_inches='tight')
    # fig2.show()
    # plt.close(fig2)
    
fig0.subplots_adjust(wspace=.3,hspace=.3)
fig0.savefig('Plots/timeseries' + figExt)
fig0.show()
plt.close(fig0)