Plot_Probability_Value.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8eeb6d16",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import scipy.io as scio\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import cartopy.feature as cfeature\n",
    "import cartopy.crs as ccrs\n",
    "import netCDF4 as nc\n",
    "import numpy as np\n",
    "from Load_data import Data_from_nc\n",
    "import xarray as xr\n",
    "import math\n",
    "import gc\n",
    "from matplotlib.colors import LinearSegmentedColormap\n",
    "\n",
    "data_land_mask = Data_from_nc('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/surfdata_0.9x1.25_hist_78pfts_CMIP6_simyr1850_c201126.nc')\n",
    "data_mask = data_land_mask.load_variable('PFTDATA_MASK')\n",
    "data_mask = data_mask[29:, :]\n",
    "\n",
    "#create func to compute PR\n",
    "def comp_PR(ctl,scen,q):\n",
    "    #first get quantile from ctl\n",
    "    qctl=ctl.quantile(q,dim='time')\n",
    "    \n",
    "    del ctl\n",
    "    gc.collect()\n",
    "    \n",
    "\n",
    "    #get subsample scenarios\n",
    "    qscen=xr.where(scen>qctl,1,0)\n",
    "    pscen=qscen.sum('time')/len(qscen.time)\n",
    "    \n",
    "    del qscen\n",
    "    gc.collect()\n",
    "\n",
    "    #compute PR\n",
    "    PR_scen=pscen/(1-q)\n",
    "    PR_scen = np.where(data_mask==0, np.nan, PR_scen)\n",
    "    PR_scen = np.where(PR_scen==0, np.nan, PR_scen)\n",
    "    \n",
    "    return PR_scen\n",
    "\n",
    "def get_data_from_mat(file, variable):\n",
    "    var_dict = scio.loadmat(file)\n",
    "    var = var_dict[variable]\n",
    "    var = var[:, :]\n",
    "    return var\n",
    "\n",
    "data_surface = Data_from_nc(\\\n",
    "                            '/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/' + \\\n",
    "                            'surfdata_irrigation_method.nc'\\\n",
    "                           )   #load the data\n",
    "\n",
    "data_irrigation_method = data_surface.load_variable('irrigation_method')\n",
    "\n",
    "data_lon = data_surface.load_variable('LONGXY')\n",
    "data_lon = data_lon[0,:]\n",
    "data_lat = data_surface.load_variable('LATIXY')\n",
    "data_lat = data_lat[29:,0]\n",
    "data_lat = data_lat[:]\n",
    "data_time = np.array(range(262800))\n",
    "\n",
    "def get_data_for_this_study(path, metric):\n",
    "    qctl_IRR_1901_1930_99 = get_data_from_mat(path + \"PercentileValue_IRR_1901_1930_\" + \\\n",
    "                     metric + \"_99.0%.mat\", 'qctl_IRR_1901_1930_99')\n",
    "    qctl_IRR_1901_1930_99= np.array(qctl_IRR_1901_1930_99)\n",
    "    qctl_IRR_1901_1930_995 = get_data_from_mat(path + \"PercentileValue_IRR_1901_1930_\" + \\\n",
    "                     metric + \"_99.5%.mat\", 'qctl_IRR_1901_1930_995')\n",
    "    qctl_IRR_1901_1930_995= np.array(qctl_IRR_1901_1930_995)\n",
    "    qctl_IRR_1901_1930_999 = get_data_from_mat(path + \"PercentileValue_IRR_1901_1930_\" + \\\n",
    "                     metric + \"_99.9%.mat\", 'qctl_IRR_1901_1930_999')\n",
    "    qctl_IRR_1901_1930_999= np.array(qctl_IRR_1901_1930_999)\n",
    "    \n",
    "    qctl_IRR_1985_2014_99 = get_data_from_mat(path + \"PercentileValue_IRR_1985_2014_\" + \\\n",
    "                     metric + \"_99.0%.mat\", 'qctl_IRR_1985_2014_99')\n",
    "    qctl_IRR_1985_2014_99= np.array(qctl_IRR_1985_2014_99)\n",
    "    qctl_IRR_1985_2014_995 = get_data_from_mat(path + \"PercentileValue_IRR_1985_2014_\" + \\\n",
    "                     metric + \"_99.5%.mat\", 'qctl_IRR_1985_2014_995')\n",
    "    qctl_IRR_1985_2014_995= np.array(qctl_IRR_1985_2014_995)\n",
    "    qctl_IRR_1985_2014_999 = get_data_from_mat(path + \"PercentileValue_IRR_1985_2014_\" + \\\n",
    "                     metric + \"_99.9%.mat\", 'qctl_IRR_1985_2014_999')\n",
    "    qctl_IRR_1985_2014_999= np.array(qctl_IRR_1985_2014_999)\n",
    "    \n",
    "    \n",
    "    \n",
    "    qctl_NOI_1901_1930_99 = get_data_from_mat(path + \"PercentileValue_NOI_1901_1930_\" + \\\n",
    "                     metric + \"_99.0%.mat\", 'qctl_NOI_1901_1930_99')\n",
    "    qctl_NOI_1901_1930_99= np.array(qctl_NOI_1901_1930_99)\n",
    "    qctl_NOI_1901_1930_995 = get_data_from_mat(path + \"PercentileValue_NOI_1901_1930_\" + \\\n",
    "                     metric + \"_99.5%.mat\", 'qctl_NOI_1901_1930_995')\n",
    "    qctl_NOI_1901_1930_995= np.array(qctl_NOI_1901_1930_995)\n",
    "    qctl_NOI_1901_1930_999 = get_data_from_mat(path + \"PercentileValue_NOI_1901_1930_\" + \\\n",
    "                     metric + \"_99.9%.mat\", 'qctl_NOI_1901_1930_999')\n",
    "    qctl_NOI_1901_1930_999= np.array(qctl_NOI_1901_1930_999)\n",
    "    \n",
    "    qctl_NOI_1985_2014_99 = get_data_from_mat(path + \"PercentileValue_NOI_1985_2014_\" + \\\n",
    "                     metric + \"_99.0%.mat\", 'qctl_NOI_1985_2014_99')\n",
    "    qctl_NOI_1985_2014_99= np.array(qctl_NOI_1985_2014_99)\n",
    "    qctl_NOI_1985_2014_995 = get_data_from_mat(path + \"PercentileValue_NOI_1985_2014_\" + \\\n",
    "                     metric + \"_99.5%.mat\", 'qctl_NOI_1985_2014_995')\n",
    "    qctl_NOI_1985_2014_995= np.array(qctl_NOI_1985_2014_995)\n",
    "    qctl_NOI_1985_2014_999 = get_data_from_mat(path + \"PercentileValue_NOI_1985_2014_\" + \\\n",
    "                     metric + \"_99.9%.mat\", 'qctl_NOI_1985_2014_999')\n",
    "    qctl_NOI_1985_2014_999= np.array(qctl_NOI_1985_2014_999)\n",
    "    \n",
    "    return qctl_NOI_1901_1930_99, qctl_NOI_1901_1930_995, qctl_NOI_1901_1930_999, qctl_IRR_1901_1930_99, qctl_IRR_1901_1930_995, qctl_IRR_1901_1930_999, qctl_NOI_1985_2014_99, qctl_NOI_1985_2014_995, qctl_NOI_1985_2014_999, qctl_IRR_1985_2014_99, qctl_IRR_1985_2014_995, qctl_IRR_1985_2014_999\n",
    "\n",
    "def plot_PR_absolu(ax, data_xarray, title, metric, perc, levels, ticklabels, unit):\n",
    "    cmap = 'coolwarm'\n",
    "    Blues = mpl.cm.get_cmap('hot_r')\n",
    "    \n",
    "    colors_Blues = [Blues(0.1), Blues(0.3), Blues(0.4), Blues(0.5), Blues(0.6), Blues(0.7), Blues(0.8), Blues(0.9), Blues(1.0)]\n",
    "    cmap = LinearSegmentedColormap.from_list('cmap_name', colors_Blues, N=9)\n",
    "    vmin=1/4\n",
    "    vmax=4\n",
    "    alpha=0.05\n",
    "    levels = levels\n",
    "    im = data_xarray.plot(ax=ax, cmap=cmap,vmin=vmin,vmax=vmax,levels=levels, extend='both', add_colorbar=False, add_labels=False)\n",
    "    cb = plt.colorbar(im,fraction= 0.3, pad= 0.04, extend='both', orientation='horizontal')\n",
    "    cb.ax.tick_params(labelsize=12)\n",
    "    cb.set_ticks(levels)\n",
    "    cb.set_label(label=metric + ' ' + unit, fontsize=14)\n",
    "    cb.set_ticklabels(ticklabels)\n",
    "    ax.coastlines(color='dimgray', linewidth=0.5)\n",
    "    ax.add_feature(cfeature.OCEAN, color='lightgrey')\n",
    "    ax.set_title(title, loc='right', fontsize=14)\n",
    "    ax.set_title(metric + ' ' + perc, loc='left', fontsize=14)\n",
    "    ax.add_feature(cfeature.BORDERS, linestyle='-', linewidth=0.1)\n",
    "    \n",
    "def plot_PR_single(ax, data_xarray, data_signal, title, metric, perc, unit):\n",
    "    cmap = 'coolwarm'\n",
    "    Blues = mpl.cm.get_cmap('coolwarm')\n",
    "    \n",
    "    colors_Blues = [Blues(0.0), Blues(0.1), Blues(0.2), Blues(0.3), 'lightgray', Blues(0.7), Blues(0.8), Blues(0.9), Blues(1.0)]\n",
    "    cmap = LinearSegmentedColormap.from_list('cmap_name', colors_Blues, N=9)\n",
    "    vmin=1/4\n",
    "    vmax=4\n",
    "    alpha=0.05\n",
    "    levels = [-2, -1, -0.5, -0.2, 0.2, 0.5, 1, 2]\n",
    "    im = data_xarray.plot(ax=ax, cmap=cmap,vmin=vmin,vmax=vmax,levels=levels, extend='both', add_colorbar=False, add_labels=False)\n",
    "    cb = plt.colorbar(im,fraction= 0.3, pad= 0.04, extend='both', orientation='horizontal')\n",
    "    cb.ax.tick_params(labelsize=12)\n",
    "    cb.set_label(label = '$\\\\Delta$' + metric + ' ' + unit, fontsize=14)\n",
    "    cb.set_ticks([-2, -1, -0.5, -0.2, 0.2, 0.5, 1, 2])\n",
    "    cb.set_ticklabels([\"-2\", \"-1\", \"-0.5\", \"-0.2\", \"+0.2\", \"+0.5\", \"+1\", \"+2\"])\n",
    "    ax.coastlines(color='dimgray', linewidth=0.5)\n",
    "    ax.add_feature(cfeature.OCEAN, color='lightgrey')\n",
    "    ax.set_title(title, loc='right', fontsize=14)\n",
    "    ax.set_title(metric + ' ' + perc, loc='left', fontsize=14)\n",
    "    ax.add_feature(cfeature.BORDERS, linestyle='-', linewidth=0.1)\n",
    "    ax.contourf(data_lon, data_lat, data_signal, levels=[0, 3.5], colors= 'lightgray')\n",
    "    \n",
    "def plot_PR_single_no_cb(ax, data_xarray, data_signal, title, metric, perc):\n",
    "    cmap = 'coolwarm'\n",
    "    Blues = mpl.cm.get_cmap('coolwarm')\n",
    "    \n",
    "    colors_Blues = [Blues(0.0), Blues(0.1), Blues(0.2), Blues(0.3), 'lightgray', Blues(0.7), Blues(0.8), Blues(0.9), Blues(1.0)]\n",
    "    cmap = LinearSegmentedColormap.from_list('cmap_name', colors_Blues, N=9)\n",
    "    vmin=1/4\n",
    "    vmax=4\n",
    "    alpha=0.05\n",
    "    levels = [-2, -1, -0.5, -0.2, 0.2, 0.5, 1, 2]\n",
    "    im = data_xarray.plot(ax=ax, cmap=cmap,vmin=vmin,vmax=vmax,levels=levels, extend='both', add_colorbar=False, add_labels=False)\n",
    "#     cb = plt.colorbar(im,fraction= 0.3, pad= 0.04, extend='both', orientation='horizontal')\n",
    "#     cb.ax.tick_params(labelsize=12)\n",
    "#     cb.set_ticks([-2, -1, -0.5, -0.2, 0.2, 0.5, 1, 2])\n",
    "#     cb.set_ticklabels([\"-2\", \"-1\", \"-0.5\", \"-0.2\", \"+0.2\", \"+0.5\", \"+1\", \"+2\"])\n",
    "    ax.coastlines(color='dimgray', linewidth=0.5)\n",
    "    ax.add_feature(cfeature.OCEAN, color='lightgrey')\n",
    "    ax.set_title(title, loc='right', fontsize=14)\n",
    "    ax.set_title(metric + ' ' + perc, loc='left', fontsize=14)\n",
    "    ax.add_feature(cfeature.BORDERS, linestyle='-', linewidth=0.1)\n",
    "    ax.contourf(data_lon, data_lat, data_signal, levels=[0, 3.5], colors= 'lightgray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "270938df",
   "metadata": {},
   "outputs": [],
   "source": [
    "def all_for_paper():\n",
    "    \n",
    "    metric = 'TSA'\n",
    "    varname = '$\\mathregular{T_{2m}}$'\n",
    "    unit = '$(^\\circ C)$'\n",
    "    levels = [28+273.15, 30+273.15, 32+273.15, 34+273.15, 36+273.15, 38+273.15, 40+273.15, 42+273.15]\n",
    "    ticklabels = [\"28\", \"30\", \"32\", \"34\", \"36\", \"38\", \"40\", \"42\"]\n",
    "    \n",
    "    qctl_NOI_1901_1930_99_CESM2, qctl_NOI_1901_1930_995_CESM2, qctl_NOI_1901_1930_999_CESM2, qctl_IRR_1901_1930_99_CESM2, qctl_IRR_1901_1930_995_CESM2, qctl_IRR_1901_1930_999_CESM2, qctl_NOI_1985_2014_99_CESM2, qctl_NOI_1985_2014_995_CESM2, qctl_NOI_1985_2014_999_CESM2, qctl_IRR_1985_2014_99_CESM2, qctl_IRR_1985_2014_995_CESM2, qctl_IRR_1985_2014_999_CESM2=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/PR/', metric)\n",
    "    qctl_NOI_1901_1930_99_CESM2_gw, qctl_NOI_1901_1930_995_CESM2_gw, qctl_NOI_1901_1930_999_CESM2_gw, qctl_IRR_1901_1930_99_CESM2_gw, qctl_IRR_1901_1930_995_CESM2_gw, qctl_IRR_1901_1930_999_CESM2_gw, qctl_NOI_1985_2014_99_CESM2_gw, qctl_NOI_1985_2014_995_CESM2_gw, qctl_NOI_1985_2014_999_CESM2_gw, qctl_IRR_1985_2014_99_CESM2_gw, qctl_IRR_1985_2014_995_CESM2_gw, qctl_IRR_1985_2014_999_CESM2_gw=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2_gw/PR/', metric)\n",
    "    qctl_NOI_1901_1930_99_NorESM, qctl_NOI_1901_1930_995_NorESM, qctl_NOI_1901_1930_999_NorESM, qctl_IRR_1901_1930_99_NorESM, qctl_IRR_1901_1930_995_NorESM, qctl_IRR_1901_1930_999_NorESM, qctl_NOI_1985_2014_99_NorESM, qctl_NOI_1985_2014_995_NorESM, qctl_NOI_1985_2014_999_NorESM, qctl_IRR_1985_2014_99_NorESM, qctl_IRR_1985_2014_995_NorESM, qctl_IRR_1985_2014_999_NorESM=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/NorESM/PR/', metric)\n",
    "    qctl_NOI_1901_1930_99_E3SM, qctl_NOI_1901_1930_995_E3SM, qctl_NOI_1901_1930_999_E3SM, qctl_IRR_1901_1930_99_E3SM, qctl_IRR_1901_1930_995_E3SM, qctl_IRR_1901_1930_999_E3SM, qctl_NOI_1985_2014_99_E3SM, qctl_NOI_1985_2014_995_E3SM, qctl_NOI_1985_2014_999_E3SM, qctl_IRR_1985_2014_99_E3SM, qctl_IRR_1985_2014_995_E3SM, qctl_IRR_1985_2014_999_E3SM=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/', metric)\n",
    "    qctl_NOI_1901_1930_99_IPSL, qctl_NOI_1901_1930_995_IPSL, qctl_NOI_1901_1930_999_IPSL, qctl_IRR_1901_1930_99_IPSL, qctl_IRR_1901_1930_995_IPSL, qctl_IRR_1901_1930_999_IPSL, qctl_NOI_1985_2014_99_IPSL, qctl_NOI_1985_2014_995_IPSL, qctl_NOI_1985_2014_999_IPSL, qctl_IRR_1985_2014_99_IPSL, qctl_IRR_1985_2014_995_IPSL, qctl_IRR_1985_2014_999_IPSL=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/IPSL-CM6/PR/', metric)\n",
    "    qctl_NOI_1901_1930_99_CNRM, qctl_NOI_1901_1930_995_CNRM, qctl_NOI_1901_1930_999_CNRM, qctl_IRR_1901_1930_99_CNRM, qctl_IRR_1901_1930_995_CNRM, qctl_IRR_1901_1930_999_CNRM, qctl_NOI_1985_2014_99_CNRM, qctl_NOI_1985_2014_995_CNRM, qctl_NOI_1985_2014_999_CNRM, qctl_IRR_1985_2014_99_CNRM, qctl_IRR_1985_2014_995_CNRM, qctl_IRR_1985_2014_999_CNRM=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/CNRM/', metric)\n",
    "    \n",
    "    qctl_NOI_1901_1930_999 = (qctl_NOI_1901_1930_999_CESM2 + qctl_NOI_1901_1930_999_CESM2_gw + qctl_NOI_1901_1930_999_NorESM + qctl_NOI_1901_1930_999_E3SM + qctl_NOI_1901_1930_999_IPSL + qctl_NOI_1901_1930_999_CNRM)/6\n",
    "    \n",
    "    qctl_NOI_1985_2014_999 = (qctl_NOI_1985_2014_999_CESM2 + qctl_NOI_1985_2014_999_CESM2_gw + qctl_NOI_1985_2014_999_NorESM + qctl_NOI_1985_2014_999_E3SM + qctl_NOI_1985_2014_999_IPSL + qctl_NOI_1985_2014_999_CNRM)/6\n",
    "    \n",
    "    qctl_IRR_1901_1930_999 = (qctl_IRR_1901_1930_999_CESM2 + qctl_IRR_1901_1930_999_CESM2_gw + qctl_IRR_1901_1930_999_NorESM + qctl_IRR_1901_1930_999_E3SM + qctl_IRR_1901_1930_999_IPSL + qctl_IRR_1901_1930_999_CNRM)/6\n",
    "    \n",
    "    qctl_IRR_1985_2014_999 = (qctl_IRR_1985_2014_999_CESM2 + qctl_IRR_1985_2014_999_CESM2_gw + qctl_IRR_1985_2014_999_NorESM + qctl_IRR_1985_2014_999_E3SM + qctl_IRR_1985_2014_999_IPSL + qctl_IRR_1985_2014_999_CNRM)/6\n",
    "    \n",
    "    fig = plt.figure(figsize=(18,16),dpi=300)\n",
    "    fig.subplots_adjust(hspace=0.05, wspace=0.05, left = 0.05, right = 0.95, top = 0.95, bottom = 0.05)\n",
    "    proj=ccrs.PlateCarree()\n",
    "\n",
    "\n",
    "    data_xr_3 = xr.DataArray(qctl_NOI_1901_1930_999, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(431, projection=proj, frameon=False)\n",
    "    plot_PR_absolu(ax, data_xr_3, '1901irr (1901-1930)', varname, '99.9%', levels, ticklabels, unit)\n",
    "    ax.text(0.01, 0.92, 'a', color='dimgrey', fontsize=14, transform=ax.transAxes, weight='bold')\n",
    "    \n",
    "    \n",
    "    \n",
    "    Oth_999_CNRM = qctl_NOI_1985_2014_999_CNRM - qctl_NOI_1901_1930_999_CNRM\n",
    "    Oth_999_NorESM = qctl_NOI_1985_2014_999_NorESM - qctl_NOI_1901_1930_999_NorESM\n",
    "    Oth_999_CESM2 = qctl_NOI_1985_2014_999_CESM2 - qctl_NOI_1901_1930_999_CESM2\n",
    "    Oth_999_CESM2_gw = qctl_NOI_1985_2014_999_CESM2_gw - qctl_NOI_1901_1930_999_CESM2_gw\n",
    "    Oth_999_IPSL = qctl_NOI_1985_2014_999_IPSL - qctl_NOI_1901_1930_999_IPSL\n",
    "    Oth_999_E3SM = qctl_NOI_1985_2014_999_E3SM - qctl_NOI_1901_1930_999_E3SM\n",
    "    \n",
    "    Oth_999_consis_pos = np.zeros([163, 288])\n",
    "    Oth_999_consis_neg = np.zeros([163, 288])\n",
    "    Oth_999_consis_pos[Oth_999_CNRM>=0.2] = Oth_999_consis_pos[Oth_999_CNRM>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_CNRM<=-0.2] = Oth_999_consis_neg[Oth_999_CNRM<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis_pos[Oth_999_NorESM>=0.2] = Oth_999_consis_pos[Oth_999_NorESM>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_NorESM<=-0.2] = Oth_999_consis_neg[Oth_999_NorESM<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis_pos[Oth_999_CESM2>=0.2] = Oth_999_consis_pos[Oth_999_CESM2>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_CESM2<=-0.2] = Oth_999_consis_neg[Oth_999_CESM2<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis_pos[Oth_999_CESM2_gw>=0.2] = Oth_999_consis_pos[Oth_999_CESM2_gw>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_CESM2_gw<=-0.2] = Oth_999_consis_neg[Oth_999_CESM2_gw<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis_pos[Oth_999_IPSL>=0.2] = Oth_999_consis_pos[Oth_999_IPSL>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_IPSL<=-0.2] = Oth_999_consis_neg[Oth_999_IPSL<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis_pos[Oth_999_E3SM>=0.2] = Oth_999_consis_pos[Oth_999_E3SM>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_E3SM<=-0.2] = Oth_999_consis_neg[Oth_999_E3SM<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis = np.maximum(Oth_999_consis_pos, Oth_999_consis_neg)\n",
    "    \n",
    "\n",
    "    data_xr_6 = xr.DataArray(qctl_NOI_1985_2014_999-qctl_NOI_1901_1930_999, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(4,3,7, projection=proj, frameon=False)\n",
    "    plot_PR_single_no_cb(ax, data_xr_6, Oth_999_consis, 'ALL-IE', varname, '99.9%')\n",
    "    ax.text(0.01, 0.92, 'g', color='dimgrey', fontsize=14, transform=ax.transAxes, weight='bold')\n",
    "\n",
    "    \n",
    "    \n",
    "    Irr_999_CNRM = qctl_IRR_1985_2014_999_CNRM - qctl_NOI_1985_2014_999_CNRM\n",
    "    Irr_999_NorESM = qctl_IRR_1985_2014_999_NorESM - qctl_NOI_1985_2014_999_NorESM\n",
    "    Irr_999_CESM2 = qctl_IRR_1985_2014_999_CESM2 - qctl_NOI_1985_2014_999_CESM2\n",
    "    Irr_999_CESM2_gw = qctl_IRR_1985_2014_999_CESM2_gw - qctl_NOI_1985_2014_999_CESM2_gw\n",
    "    Irr_999_IPSL = qctl_IRR_1985_2014_999_IPSL - qctl_NOI_1985_2014_999_IPSL\n",
    "    Irr_999_E3SM = qctl_IRR_1985_2014_999_E3SM - qctl_NOI_1985_2014_999_E3SM\n",
    "    \n",
    "    Irr_999_consis_pos = np.zeros([163, 288])\n",
    "    Irr_999_consis_neg = np.zeros([163, 288])\n",
    "    Irr_999_consis_pos[Irr_999_CNRM>=0.2] = Irr_999_consis_pos[Irr_999_CNRM>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_CNRM<=-0.2] = Irr_999_consis_neg[Irr_999_CNRM<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis_pos[Irr_999_NorESM>=0.2] = Irr_999_consis_pos[Irr_999_NorESM>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_NorESM<=-0.2] = Irr_999_consis_neg[Irr_999_NorESM<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis_pos[Irr_999_CESM2>=0.2] = Irr_999_consis_pos[Irr_999_CESM2>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_CESM2<=-0.2] = Irr_999_consis_neg[Irr_999_CESM2<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis_pos[Irr_999_CESM2_gw>=0.2] = Irr_999_consis_pos[Irr_999_CESM2_gw>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_CESM2_gw<=-0.2] = Irr_999_consis_neg[Irr_999_CESM2_gw<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis_pos[Irr_999_IPSL>=0.2] = Irr_999_consis_pos[Irr_999_IPSL>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_IPSL<=-0.2] = Irr_999_consis_neg[Irr_999_IPSL<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis_pos[Irr_999_E3SM>=0.2] = Irr_999_consis_pos[Irr_999_E3SM>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_E3SM<=-0.2] = Irr_999_consis_neg[Irr_999_E3SM<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis = np.maximum(Irr_999_consis_pos, Irr_999_consis_neg)\n",
    "    \n",
    "\n",
    "    data_xr_9 = xr.DataArray(qctl_IRR_1985_2014_999-qctl_NOI_1985_2014_999, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(4,3,10, projection=proj, frameon=False)\n",
    "    plot_PR_single(ax, data_xr_9, Irr_999_consis, 'IE', varname, '99.9%', unit)\n",
    "    ax.text(0.01, 0.92, 'j', color='dimgrey', fontsize=14, transform=ax.transAxes, weight='bold')\n",
    "    \n",
    "    \n",
    "    \n",
    "    All_999_CNRM = qctl_IRR_1985_2014_999_CNRM - qctl_IRR_1901_1930_999_CNRM\n",
    "    All_999_NorESM = qctl_IRR_1985_2014_999_NorESM - qctl_IRR_1901_1930_999_NorESM\n",
    "    All_999_CESM2 = qctl_IRR_1985_2014_999_CESM2 - qctl_IRR_1901_1930_999_CESM2\n",
    "    All_999_CESM2_gw = qctl_IRR_1985_2014_999_CESM2_gw - qctl_IRR_1901_1930_999_CESM2_gw\n",
    "    All_999_IPSL = qctl_IRR_1985_2014_999_IPSL - qctl_IRR_1901_1930_999_IPSL\n",
    "    All_999_E3SM = qctl_IRR_1985_2014_999_E3SM - qctl_IRR_1901_1930_999_E3SM\n",
    "    \n",
    "    All_999_consis_pos = np.zeros([163, 288])\n",
    "    All_999_consis_neg = np.zeros([163, 288])\n",
    "    All_999_consis_pos[All_999_CNRM>=0.2] = All_999_consis_pos[All_999_CNRM>=0.2]+1\n",
    "    All_999_consis_neg[All_999_CNRM<=-0.2] = All_999_consis_neg[All_999_CNRM<=-0.2]+1\n",
    "    \n",
    "    All_999_consis_pos[All_999_NorESM>=0.2] = All_999_consis_pos[All_999_NorESM>=0.2]+1\n",
    "    All_999_consis_neg[All_999_NorESM<=-0.2] = All_999_consis_neg[All_999_NorESM<=-0.2]+1\n",
    "    \n",
    "    All_999_consis_pos[All_999_CESM2>=0.2] = All_999_consis_pos[All_999_CESM2>=0.2]+1\n",
    "    All_999_consis_neg[All_999_CESM2<=-0.2] = All_999_consis_neg[All_999_CESM2<=-0.2]+1\n",
    "    \n",
    "    All_999_consis_pos[All_999_CESM2_gw>=0.2] = All_999_consis_pos[All_999_CESM2_gw>=0.2]+1\n",
    "    All_999_consis_neg[All_999_CESM2_gw<=-0.2] = All_999_consis_neg[All_999_CESM2_gw<=-0.2]+1\n",
    "    \n",
    "    All_999_consis_pos[All_999_IPSL>=0.2] = All_999_consis_pos[All_999_IPSL>=0.2]+1\n",
    "    All_999_consis_neg[All_999_IPSL<=-0.2] = All_999_consis_neg[All_999_IPSL<=-0.2]+1\n",
    "    \n",
    "    All_999_consis_pos[All_999_E3SM>=0.2] = All_999_consis_pos[All_999_E3SM>=0.2]+1\n",
    "    All_999_consis_neg[All_999_E3SM<=-0.2] = All_999_consis_neg[All_999_E3SM<=-0.2]+1\n",
    "    \n",
    "    All_999_consis = np.maximum(All_999_consis_pos, All_999_consis_neg)\n",
    "    \n",
    "\n",
    "    data_xr_12 = xr.DataArray(qctl_IRR_1985_2014_999-qctl_NOI_1901_1930_999, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(4,3,4, projection=proj, frameon=False)\n",
    "    plot_PR_single_no_cb(ax, data_xr_12, All_999_consis, 'ALL', varname, '99.9%')\n",
    "    ax.text(0.01, 0.92, 'd', color='dimgrey', fontsize=14, transform=ax.transAxes, weight='bold')\n",
    "    \n",
    "    \n",
    "    metric = 'HUMIDEX'\n",
    "    varname = '$\\mathregular{HU}$'\n",
    "    unit = ''\n",
    "    levels = [28, 31, 34, 37, 40, 43, 46, 49]\n",
    "    ticklabels = [\"28\", \"31\", \"34\", \"37\", \"40\", \"43\", \"46\", \"49\"]\n",
    "    \n",
    "    qctl_NOI_1901_1930_99_CESM2, qctl_NOI_1901_1930_995_CESM2, qctl_NOI_1901_1930_999_CESM2, qctl_IRR_1901_1930_99_CESM2, qctl_IRR_1901_1930_995_CESM2, qctl_IRR_1901_1930_999_CESM2, qctl_NOI_1985_2014_99_CESM2, qctl_NOI_1985_2014_995_CESM2, qctl_NOI_1985_2014_999_CESM2, qctl_IRR_1985_2014_99_CESM2, qctl_IRR_1985_2014_995_CESM2, qctl_IRR_1985_2014_999_CESM2=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/PR/', metric)\n",
    "    qctl_NOI_1901_1930_99_CESM2_gw, qctl_NOI_1901_1930_995_CESM2_gw, qctl_NOI_1901_1930_999_CESM2_gw, qctl_IRR_1901_1930_99_CESM2_gw, qctl_IRR_1901_1930_995_CESM2_gw, qctl_IRR_1901_1930_999_CESM2_gw, qctl_NOI_1985_2014_99_CESM2_gw, qctl_NOI_1985_2014_995_CESM2_gw, qctl_NOI_1985_2014_999_CESM2_gw, qctl_IRR_1985_2014_99_CESM2_gw, qctl_IRR_1985_2014_995_CESM2_gw, qctl_IRR_1985_2014_999_CESM2_gw=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2_gw/PR/', metric)\n",
    "    qctl_NOI_1901_1930_99_NorESM, qctl_NOI_1901_1930_995_NorESM, qctl_NOI_1901_1930_999_NorESM, qctl_IRR_1901_1930_99_NorESM, qctl_IRR_1901_1930_995_NorESM, qctl_IRR_1901_1930_999_NorESM, qctl_NOI_1985_2014_99_NorESM, qctl_NOI_1985_2014_995_NorESM, qctl_NOI_1985_2014_999_NorESM, qctl_IRR_1985_2014_99_NorESM, qctl_IRR_1985_2014_995_NorESM, qctl_IRR_1985_2014_999_NorESM=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/NorESM/PR/', metric)\n",
    "    qctl_NOI_1901_1930_99_E3SM, qctl_NOI_1901_1930_995_E3SM, qctl_NOI_1901_1930_999_E3SM, qctl_IRR_1901_1930_99_E3SM, qctl_IRR_1901_1930_995_E3SM, qctl_IRR_1901_1930_999_E3SM, qctl_NOI_1985_2014_99_E3SM, qctl_NOI_1985_2014_995_E3SM, qctl_NOI_1985_2014_999_E3SM, qctl_IRR_1985_2014_99_E3SM, qctl_IRR_1985_2014_995_E3SM, qctl_IRR_1985_2014_999_E3SM=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/', metric)\n",
    "    qctl_NOI_1901_1930_99_IPSL, qctl_NOI_1901_1930_995_IPSL, qctl_NOI_1901_1930_999_IPSL, qctl_IRR_1901_1930_99_IPSL, qctl_IRR_1901_1930_995_IPSL, qctl_IRR_1901_1930_999_IPSL, qctl_NOI_1985_2014_99_IPSL, qctl_NOI_1985_2014_995_IPSL, qctl_NOI_1985_2014_999_IPSL, qctl_IRR_1985_2014_99_IPSL, qctl_IRR_1985_2014_995_IPSL, qctl_IRR_1985_2014_999_IPSL=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/IPSL-CM6/PR/', metric)\n",
    "    qctl_NOI_1901_1930_99_CNRM, qctl_NOI_1901_1930_995_CNRM, qctl_NOI_1901_1930_999_CNRM, qctl_IRR_1901_1930_99_CNRM, qctl_IRR_1901_1930_995_CNRM, qctl_IRR_1901_1930_999_CNRM, qctl_NOI_1985_2014_99_CNRM, qctl_NOI_1985_2014_995_CNRM, qctl_NOI_1985_2014_999_CNRM, qctl_IRR_1985_2014_99_CNRM, qctl_IRR_1985_2014_995_CNRM, qctl_IRR_1985_2014_999_CNRM=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/CNRM/', metric)\n",
    "    \n",
    "    qctl_NOI_1901_1930_999 = (qctl_NOI_1901_1930_999_CESM2 + qctl_NOI_1901_1930_999_CESM2_gw + qctl_NOI_1901_1930_999_NorESM + qctl_NOI_1901_1930_999_E3SM + qctl_NOI_1901_1930_999_IPSL + qctl_NOI_1901_1930_999_CNRM)/6\n",
    "    \n",
    "    qctl_NOI_1985_2014_999 = (qctl_NOI_1985_2014_999_CESM2 + qctl_NOI_1985_2014_999_CESM2_gw + qctl_NOI_1985_2014_999_NorESM + qctl_NOI_1985_2014_999_E3SM + qctl_NOI_1985_2014_999_IPSL + qctl_NOI_1985_2014_999_CNRM)/6\n",
    "    \n",
    "    qctl_IRR_1901_1930_999 = (qctl_IRR_1901_1930_999_CESM2 + qctl_IRR_1901_1930_999_CESM2_gw + qctl_IRR_1901_1930_999_NorESM + qctl_IRR_1901_1930_999_E3SM + qctl_IRR_1901_1930_999_IPSL + qctl_IRR_1901_1930_999_CNRM)/6\n",
    "    \n",
    "    qctl_IRR_1985_2014_999 = (qctl_IRR_1985_2014_999_CESM2 + qctl_IRR_1985_2014_999_CESM2_gw + qctl_IRR_1985_2014_999_NorESM + qctl_IRR_1985_2014_999_E3SM + qctl_IRR_1985_2014_999_IPSL + qctl_IRR_1985_2014_999_CNRM)/6\n",
    "    \n",
    "\n",
    "\n",
    "\n",
    "    data_xr_3 = xr.DataArray(qctl_NOI_1901_1930_999, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(432, projection=proj, frameon=False)\n",
    "    plot_PR_absolu(ax, data_xr_3, '1901irr (1901-1930)', varname, '99.9%', levels, ticklabels, unit)\n",
    "    ax.text(0.01, 0.92, 'b', color='dimgrey', fontsize=14, transform=ax.transAxes, weight='bold')\n",
    "    \n",
    "    \n",
    "    \n",
    "    Oth_999_CNRM = qctl_NOI_1985_2014_999_CNRM - qctl_NOI_1901_1930_999_CNRM\n",
    "    Oth_999_NorESM = qctl_NOI_1985_2014_999_NorESM - qctl_NOI_1901_1930_999_NorESM\n",
    "    Oth_999_CESM2 = qctl_NOI_1985_2014_999_CESM2 - qctl_NOI_1901_1930_999_CESM2\n",
    "    Oth_999_CESM2_gw = qctl_NOI_1985_2014_999_CESM2_gw - qctl_NOI_1901_1930_999_CESM2_gw\n",
    "    Oth_999_IPSL = qctl_NOI_1985_2014_999_IPSL - qctl_NOI_1901_1930_999_IPSL\n",
    "    Oth_999_E3SM = qctl_NOI_1985_2014_999_E3SM - qctl_NOI_1901_1930_999_E3SM\n",
    "    \n",
    "    Oth_999_consis_pos = np.zeros([163, 288])\n",
    "    Oth_999_consis_neg = np.zeros([163, 288])\n",
    "    Oth_999_consis_pos[Oth_999_CNRM>=0.2] = Oth_999_consis_pos[Oth_999_CNRM>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_CNRM<=-0.2] = Oth_999_consis_neg[Oth_999_CNRM<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis_pos[Oth_999_NorESM>=0.2] = Oth_999_consis_pos[Oth_999_NorESM>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_NorESM<=-0.2] = Oth_999_consis_neg[Oth_999_NorESM<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis_pos[Oth_999_CESM2>=0.2] = Oth_999_consis_pos[Oth_999_CESM2>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_CESM2<=-0.2] = Oth_999_consis_neg[Oth_999_CESM2<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis_pos[Oth_999_CESM2_gw>=0.2] = Oth_999_consis_pos[Oth_999_CESM2_gw>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_CESM2_gw<=-0.2] = Oth_999_consis_neg[Oth_999_CESM2_gw<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis_pos[Oth_999_IPSL>=0.2] = Oth_999_consis_pos[Oth_999_IPSL>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_IPSL<=-0.2] = Oth_999_consis_neg[Oth_999_IPSL<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis_pos[Oth_999_E3SM>=0.2] = Oth_999_consis_pos[Oth_999_E3SM>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_E3SM<=-0.2] = Oth_999_consis_neg[Oth_999_E3SM<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis = np.maximum(Oth_999_consis_pos, Oth_999_consis_neg)\n",
    "    \n",
    "\n",
    "    data_xr_6 = xr.DataArray(qctl_NOI_1985_2014_999-qctl_NOI_1901_1930_999, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(4, 3, 8, projection=proj, frameon=False)\n",
    "    plot_PR_single_no_cb(ax, data_xr_6, Oth_999_consis, 'ALL-IE', varname, '99.9%')\n",
    "    ax.text(0.01, 0.92, 'h', color='dimgrey', fontsize=14, transform=ax.transAxes, weight='bold')\n",
    "\n",
    "    \n",
    "    \n",
    "    Irr_999_CNRM = qctl_IRR_1985_2014_999_CNRM - qctl_NOI_1985_2014_999_CNRM\n",
    "    Irr_999_NorESM = qctl_IRR_1985_2014_999_NorESM - qctl_NOI_1985_2014_999_NorESM\n",
    "    Irr_999_CESM2 = qctl_IRR_1985_2014_999_CESM2 - qctl_NOI_1985_2014_999_CESM2\n",
    "    Irr_999_CESM2_gw = qctl_IRR_1985_2014_999_CESM2_gw - qctl_NOI_1985_2014_999_CESM2_gw\n",
    "    Irr_999_IPSL = qctl_IRR_1985_2014_999_IPSL - qctl_NOI_1985_2014_999_IPSL\n",
    "    Irr_999_E3SM = qctl_IRR_1985_2014_999_E3SM - qctl_NOI_1985_2014_999_E3SM\n",
    "    \n",
    "    Irr_999_consis_pos = np.zeros([163, 288])\n",
    "    Irr_999_consis_neg = np.zeros([163, 288])\n",
    "    Irr_999_consis_pos[Irr_999_CNRM>=0.2] = Irr_999_consis_pos[Irr_999_CNRM>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_CNRM<=-0.2] = Irr_999_consis_neg[Irr_999_CNRM<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis_pos[Irr_999_NorESM>=0.2] = Irr_999_consis_pos[Irr_999_NorESM>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_NorESM<=-0.2] = Irr_999_consis_neg[Irr_999_NorESM<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis_pos[Irr_999_CESM2>=0.2] = Irr_999_consis_pos[Irr_999_CESM2>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_CESM2<=-0.2] = Irr_999_consis_neg[Irr_999_CESM2<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis_pos[Irr_999_CESM2_gw>=0.2] = Irr_999_consis_pos[Irr_999_CESM2_gw>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_CESM2_gw<=-0.2] = Irr_999_consis_neg[Irr_999_CESM2_gw<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis_pos[Irr_999_IPSL>=0.2] = Irr_999_consis_pos[Irr_999_IPSL>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_IPSL<=-0.2] = Irr_999_consis_neg[Irr_999_IPSL<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis_pos[Irr_999_E3SM>=0.2] = Irr_999_consis_pos[Irr_999_E3SM>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_E3SM<=-0.2] = Irr_999_consis_neg[Irr_999_E3SM<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis = np.maximum(Irr_999_consis_pos, Irr_999_consis_neg)\n",
    "    \n",
    "\n",
    "    data_xr_9 = xr.DataArray(qctl_IRR_1985_2014_999-qctl_NOI_1985_2014_999, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(4,3,11, projection=proj, frameon=False)\n",
    "    plot_PR_single(ax, data_xr_9, Irr_999_consis, 'IE', varname, '99.9%', unit)\n",
    "    ax.text(0.01, 0.92, 'k', color='dimgrey', fontsize=14, transform=ax.transAxes, weight='bold')\n",
    "    \n",
    "    \n",
    "    \n",
    "    All_999_CNRM = qctl_IRR_1985_2014_999_CNRM - qctl_IRR_1901_1930_999_CNRM\n",
    "    All_999_NorESM = qctl_IRR_1985_2014_999_NorESM - qctl_IRR_1901_1930_999_NorESM\n",
    "    All_999_CESM2 = qctl_IRR_1985_2014_999_CESM2 - qctl_IRR_1901_1930_999_CESM2\n",
    "    All_999_CESM2_gw = qctl_IRR_1985_2014_999_CESM2_gw - qctl_IRR_1901_1930_999_CESM2_gw\n",
    "    All_999_IPSL = qctl_IRR_1985_2014_999_IPSL - qctl_IRR_1901_1930_999_IPSL\n",
    "    All_999_E3SM = qctl_IRR_1985_2014_999_E3SM - qctl_IRR_1901_1930_999_E3SM\n",
    "    \n",
    "    All_999_consis_pos = np.zeros([163, 288])\n",
    "    All_999_consis_neg = np.zeros([163, 288])\n",
    "    All_999_consis_pos[All_999_CNRM>=0.2] = All_999_consis_pos[All_999_CNRM>=0.2]+1\n",
    "    All_999_consis_neg[All_999_CNRM<=-0.2] = All_999_consis_neg[All_999_CNRM<=-0.2]+1\n",
    "    \n",
    "    All_999_consis_pos[All_999_NorESM>=0.2] = All_999_consis_pos[All_999_NorESM>=0.2]+1\n",
    "    All_999_consis_neg[All_999_NorESM<=-0.2] = All_999_consis_neg[All_999_NorESM<=-0.2]+1\n",
    "    \n",
    "    All_999_consis_pos[All_999_CESM2>=0.2] = All_999_consis_pos[All_999_CESM2>=0.2]+1\n",
    "    All_999_consis_neg[All_999_CESM2<=-0.2] = All_999_consis_neg[All_999_CESM2<=-0.2]+1\n",
    "    \n",
    "    All_999_consis_pos[All_999_CESM2_gw>=0.2] = All_999_consis_pos[All_999_CESM2_gw>=0.2]+1\n",
    "    All_999_consis_neg[All_999_CESM2_gw<=-0.2] = All_999_consis_neg[All_999_CESM2_gw<=-0.2]+1\n",
    "    \n",
    "    All_999_consis_pos[All_999_IPSL>=0.2] = All_999_consis_pos[All_999_IPSL>=0.2]+1\n",
    "    All_999_consis_neg[All_999_IPSL<=-0.2] = All_999_consis_neg[All_999_IPSL<=-0.2]+1\n",
    "    \n",
    "    All_999_consis_pos[All_999_E3SM>=0.2] = All_999_consis_pos[All_999_E3SM>=0.2]+1\n",
    "    All_999_consis_neg[All_999_E3SM<=-0.2] = All_999_consis_neg[All_999_E3SM<=-0.2]+1\n",
    "    \n",
    "    All_999_consis = np.maximum(All_999_consis_pos, All_999_consis_neg)\n",
    "    \n",
    "\n",
    "    data_xr_12 = xr.DataArray(qctl_IRR_1985_2014_999-qctl_NOI_1901_1930_999, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(4,3,5, projection=proj, frameon=False)\n",
    "    plot_PR_single_no_cb(ax, data_xr_12, All_999_consis, 'ALL', varname, '99.9%')\n",
    "    ax.text(0.01, 0.92, 'e', color='dimgrey', fontsize=14, transform=ax.transAxes, weight='bold')\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    metric = 'T_w'\n",
    "    varname = '$\\mathregular{T_{w}}$'\n",
    "    unit = '$(^\\circ C)$'\n",
    "    levels = [17, 19, 21, 23, 25, 27, 29, 31]\n",
    "    ticklabels = [\"17\", \"19\", \"24\", \"26\", \"25\", \"27\", \"29\", \"31\"]\n",
    "    \n",
    "    qctl_NOI_1901_1930_99_CESM2, qctl_NOI_1901_1930_995_CESM2, qctl_NOI_1901_1930_999_CESM2, qctl_IRR_1901_1930_99_CESM2, qctl_IRR_1901_1930_995_CESM2, qctl_IRR_1901_1930_999_CESM2, qctl_NOI_1985_2014_99_CESM2, qctl_NOI_1985_2014_995_CESM2, qctl_NOI_1985_2014_999_CESM2, qctl_IRR_1985_2014_99_CESM2, qctl_IRR_1985_2014_995_CESM2, qctl_IRR_1985_2014_999_CESM2=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2/PR/', metric)\n",
    "    qctl_NOI_1901_1930_99_CESM2_gw, qctl_NOI_1901_1930_995_CESM2_gw, qctl_NOI_1901_1930_999_CESM2_gw, qctl_IRR_1901_1930_99_CESM2_gw, qctl_IRR_1901_1930_995_CESM2_gw, qctl_IRR_1901_1930_999_CESM2_gw, qctl_NOI_1985_2014_99_CESM2_gw, qctl_NOI_1985_2014_995_CESM2_gw, qctl_NOI_1985_2014_999_CESM2_gw, qctl_IRR_1985_2014_99_CESM2_gw, qctl_IRR_1985_2014_995_CESM2_gw, qctl_IRR_1985_2014_999_CESM2_gw=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/CESM2_gw/PR/', metric)\n",
    "    qctl_NOI_1901_1930_99_NorESM, qctl_NOI_1901_1930_995_NorESM, qctl_NOI_1901_1930_999_NorESM, qctl_IRR_1901_1930_99_NorESM, qctl_IRR_1901_1930_995_NorESM, qctl_IRR_1901_1930_999_NorESM, qctl_NOI_1985_2014_99_NorESM, qctl_NOI_1985_2014_995_NorESM, qctl_NOI_1985_2014_999_NorESM, qctl_IRR_1985_2014_99_NorESM, qctl_IRR_1985_2014_995_NorESM, qctl_IRR_1985_2014_999_NorESM=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/after_process/NorESM/PR/', metric)\n",
    "    qctl_NOI_1901_1930_99_E3SM, qctl_NOI_1901_1930_995_E3SM, qctl_NOI_1901_1930_999_E3SM, qctl_IRR_1901_1930_99_E3SM, qctl_IRR_1901_1930_995_E3SM, qctl_IRR_1901_1930_999_E3SM, qctl_NOI_1985_2014_99_E3SM, qctl_NOI_1985_2014_995_E3SM, qctl_NOI_1985_2014_999_E3SM, qctl_IRR_1985_2014_99_E3SM, qctl_IRR_1985_2014_995_E3SM, qctl_IRR_1985_2014_999_E3SM=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/E3SM_cdo/', metric)\n",
    "    qctl_NOI_1901_1930_99_IPSL, qctl_NOI_1901_1930_995_IPSL, qctl_NOI_1901_1930_999_IPSL, qctl_IRR_1901_1930_99_IPSL, qctl_IRR_1901_1930_995_IPSL, qctl_IRR_1901_1930_999_IPSL, qctl_NOI_1985_2014_99_IPSL, qctl_NOI_1985_2014_995_IPSL, qctl_NOI_1985_2014_999_IPSL, qctl_IRR_1985_2014_99_IPSL, qctl_IRR_1985_2014_995_IPSL, qctl_IRR_1985_2014_999_IPSL=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/IPSL-CM6/PR/', metric)\n",
    "    qctl_NOI_1901_1930_99_CNRM, qctl_NOI_1901_1930_995_CNRM, qctl_NOI_1901_1930_999_CNRM, qctl_IRR_1901_1930_99_CNRM, qctl_IRR_1901_1930_995_CNRM, qctl_IRR_1901_1930_999_CNRM, qctl_NOI_1985_2014_99_CNRM, qctl_NOI_1985_2014_995_CNRM, qctl_NOI_1985_2014_999_CNRM, qctl_IRR_1985_2014_99_CNRM, qctl_IRR_1985_2014_995_CNRM, qctl_IRR_1985_2014_999_CNRM=get_data_for_this_study('/dodrio/scratch/projects/2022_200/project_output/cesm/yi_yao_IRRMIP/CNRM/', metric)\n",
    "    \n",
    "    qctl_NOI_1901_1930_999 = (qctl_NOI_1901_1930_999_CESM2 + qctl_NOI_1901_1930_999_CESM2_gw + qctl_NOI_1901_1930_999_NorESM + qctl_NOI_1901_1930_999_E3SM + qctl_NOI_1901_1930_999_IPSL + qctl_NOI_1901_1930_999_CNRM)/6\n",
    "    \n",
    "    qctl_NOI_1985_2014_999 = (qctl_NOI_1985_2014_999_CESM2 + qctl_NOI_1985_2014_999_CESM2_gw + qctl_NOI_1985_2014_999_NorESM + qctl_NOI_1985_2014_999_E3SM + qctl_NOI_1985_2014_999_IPSL + qctl_NOI_1985_2014_999_CNRM)/6\n",
    "    \n",
    "    qctl_IRR_1901_1930_999 = (qctl_IRR_1901_1930_999_CESM2 + qctl_IRR_1901_1930_999_CESM2_gw + qctl_IRR_1901_1930_999_NorESM + qctl_IRR_1901_1930_999_E3SM + qctl_IRR_1901_1930_999_IPSL + qctl_IRR_1901_1930_999_CNRM)/6\n",
    "    \n",
    "    qctl_IRR_1985_2014_999 = (qctl_IRR_1985_2014_999_CESM2 + qctl_IRR_1985_2014_999_CESM2_gw + qctl_IRR_1985_2014_999_NorESM + qctl_IRR_1985_2014_999_E3SM + qctl_IRR_1985_2014_999_IPSL + qctl_IRR_1985_2014_999_CNRM)/6\n",
    "    \n",
    "\n",
    "\n",
    "\n",
    "    data_xr_3 = xr.DataArray(qctl_NOI_1901_1930_999, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(433, projection=proj, frameon=False)\n",
    "    plot_PR_absolu(ax, data_xr_3, '1901irr (1901-1930)', varname, '99.9%', levels, ticklabels, unit)\n",
    "    ax.text(0.01, 0.92, 'c', color='dimgrey', fontsize=14, transform=ax.transAxes, weight='bold')\n",
    "    \n",
    "    \n",
    "    \n",
    "    Oth_999_CNRM = qctl_NOI_1985_2014_999_CNRM - qctl_NOI_1901_1930_999_CNRM\n",
    "    Oth_999_NorESM = qctl_NOI_1985_2014_999_NorESM - qctl_NOI_1901_1930_999_NorESM\n",
    "    Oth_999_CESM2 = qctl_NOI_1985_2014_999_CESM2 - qctl_NOI_1901_1930_999_CESM2\n",
    "    Oth_999_CESM2_gw = qctl_NOI_1985_2014_999_CESM2_gw - qctl_NOI_1901_1930_999_CESM2_gw\n",
    "    Oth_999_IPSL = qctl_NOI_1985_2014_999_IPSL - qctl_NOI_1901_1930_999_IPSL\n",
    "    Oth_999_E3SM = qctl_NOI_1985_2014_999_E3SM - qctl_NOI_1901_1930_999_E3SM\n",
    "    \n",
    "    Oth_999_consis_pos = np.zeros([163, 288])\n",
    "    Oth_999_consis_neg = np.zeros([163, 288])\n",
    "    Oth_999_consis_pos[Oth_999_CNRM>=0.2] = Oth_999_consis_pos[Oth_999_CNRM>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_CNRM<=-0.2] = Oth_999_consis_neg[Oth_999_CNRM<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis_pos[Oth_999_NorESM>=0.2] = Oth_999_consis_pos[Oth_999_NorESM>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_NorESM<=-0.2] = Oth_999_consis_neg[Oth_999_NorESM<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis_pos[Oth_999_CESM2>=0.2] = Oth_999_consis_pos[Oth_999_CESM2>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_CESM2<=-0.2] = Oth_999_consis_neg[Oth_999_CESM2<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis_pos[Oth_999_CESM2_gw>=0.2] = Oth_999_consis_pos[Oth_999_CESM2_gw>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_CESM2_gw<=-0.2] = Oth_999_consis_neg[Oth_999_CESM2_gw<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis_pos[Oth_999_IPSL>=0.2] = Oth_999_consis_pos[Oth_999_IPSL>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_IPSL<=-0.2] = Oth_999_consis_neg[Oth_999_IPSL<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis_pos[Oth_999_E3SM>=0.2] = Oth_999_consis_pos[Oth_999_E3SM>=0.2]+1\n",
    "    Oth_999_consis_neg[Oth_999_E3SM<=-0.2] = Oth_999_consis_neg[Oth_999_E3SM<=-0.2]+1\n",
    "    \n",
    "    Oth_999_consis = np.maximum(Oth_999_consis_pos, Oth_999_consis_neg)\n",
    "    \n",
    "\n",
    "    data_xr_6 = xr.DataArray(qctl_NOI_1985_2014_999-qctl_NOI_1901_1930_999, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(4, 3, 9, projection=proj, frameon=False)\n",
    "    plot_PR_single_no_cb(ax, data_xr_6, Oth_999_consis, 'ALL-IE', varname, '99.9%')\n",
    "    ax.text(0.01, 0.92, 'i', color='dimgrey', fontsize=14, transform=ax.transAxes, weight='bold')\n",
    "\n",
    "    \n",
    "    \n",
    "    Irr_999_CNRM = qctl_IRR_1985_2014_999_CNRM - qctl_NOI_1985_2014_999_CNRM\n",
    "    Irr_999_NorESM = qctl_IRR_1985_2014_999_NorESM - qctl_NOI_1985_2014_999_NorESM\n",
    "    Irr_999_CESM2 = qctl_IRR_1985_2014_999_CESM2 - qctl_NOI_1985_2014_999_CESM2\n",
    "    Irr_999_CESM2_gw = qctl_IRR_1985_2014_999_CESM2_gw - qctl_NOI_1985_2014_999_CESM2_gw\n",
    "    Irr_999_IPSL = qctl_IRR_1985_2014_999_IPSL - qctl_NOI_1985_2014_999_IPSL\n",
    "    Irr_999_E3SM = qctl_IRR_1985_2014_999_E3SM - qctl_NOI_1985_2014_999_E3SM\n",
    "    \n",
    "    Irr_999_consis_pos = np.zeros([163, 288])\n",
    "    Irr_999_consis_neg = np.zeros([163, 288])\n",
    "    Irr_999_consis_pos[Irr_999_CNRM>=0.2] = Irr_999_consis_pos[Irr_999_CNRM>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_CNRM<=-0.2] = Irr_999_consis_neg[Irr_999_CNRM<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis_pos[Irr_999_NorESM>=0.2] = Irr_999_consis_pos[Irr_999_NorESM>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_NorESM<=-0.2] = Irr_999_consis_neg[Irr_999_NorESM<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis_pos[Irr_999_CESM2>=0.2] = Irr_999_consis_pos[Irr_999_CESM2>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_CESM2<=-0.2] = Irr_999_consis_neg[Irr_999_CESM2<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis_pos[Irr_999_CESM2_gw>=0.2] = Irr_999_consis_pos[Irr_999_CESM2_gw>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_CESM2_gw<=-0.2] = Irr_999_consis_neg[Irr_999_CESM2_gw<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis_pos[Irr_999_IPSL>=0.2] = Irr_999_consis_pos[Irr_999_IPSL>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_IPSL<=-0.2] = Irr_999_consis_neg[Irr_999_IPSL<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis_pos[Irr_999_E3SM>=0.2] = Irr_999_consis_pos[Irr_999_E3SM>=0.2]+1\n",
    "    Irr_999_consis_neg[Irr_999_E3SM<=-0.2] = Irr_999_consis_neg[Irr_999_E3SM<=-0.2]+1\n",
    "    \n",
    "    Irr_999_consis = np.maximum(Irr_999_consis_pos, Irr_999_consis_neg)\n",
    "    \n",
    "\n",
    "    data_xr_9 = xr.DataArray(qctl_IRR_1985_2014_999-qctl_NOI_1985_2014_999, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(4,3,12, projection=proj, frameon=False)\n",
    "    plot_PR_single(ax, data_xr_9, Irr_999_consis, 'IE', varname, '99.9%', unit)\n",
    "    ax.text(0.01, 0.92, 'l', color='dimgrey', fontsize=14, transform=ax.transAxes, weight='bold')\n",
    "    \n",
    "    \n",
    "    \n",
    "    All_999_CNRM = qctl_IRR_1985_2014_999_CNRM - qctl_IRR_1901_1930_999_CNRM\n",
    "    All_999_NorESM = qctl_IRR_1985_2014_999_NorESM - qctl_IRR_1901_1930_999_NorESM\n",
    "    All_999_CESM2 = qctl_IRR_1985_2014_999_CESM2 - qctl_IRR_1901_1930_999_CESM2\n",
    "    All_999_CESM2_gw = qctl_IRR_1985_2014_999_CESM2_gw - qctl_IRR_1901_1930_999_CESM2_gw\n",
    "    All_999_IPSL = qctl_IRR_1985_2014_999_IPSL - qctl_IRR_1901_1930_999_IPSL\n",
    "    All_999_E3SM = qctl_IRR_1985_2014_999_E3SM - qctl_IRR_1901_1930_999_E3SM\n",
    "    \n",
    "    All_999_consis_pos = np.zeros([163, 288])\n",
    "    All_999_consis_neg = np.zeros([163, 288])\n",
    "    All_999_consis_pos[All_999_CNRM>=0.2] = All_999_consis_pos[All_999_CNRM>=0.2]+1\n",
    "    All_999_consis_neg[All_999_CNRM<=-0.2] = All_999_consis_neg[All_999_CNRM<=-0.2]+1\n",
    "    \n",
    "    All_999_consis_pos[All_999_NorESM>=0.2] = All_999_consis_pos[All_999_NorESM>=0.2]+1\n",
    "    All_999_consis_neg[All_999_NorESM<=-0.2] = All_999_consis_neg[All_999_NorESM<=-0.2]+1\n",
    "    \n",
    "    All_999_consis_pos[All_999_CESM2>=0.2] = All_999_consis_pos[All_999_CESM2>=0.2]+1\n",
    "    All_999_consis_neg[All_999_CESM2<=-0.2] = All_999_consis_neg[All_999_CESM2<=-0.2]+1\n",
    "    \n",
    "    All_999_consis_pos[All_999_CESM2_gw>=0.2] = All_999_consis_pos[All_999_CESM2_gw>=0.2]+1\n",
    "    All_999_consis_neg[All_999_CESM2_gw<=-0.2] = All_999_consis_neg[All_999_CESM2_gw<=-0.2]+1\n",
    "    \n",
    "    All_999_consis_pos[All_999_IPSL>=0.2] = All_999_consis_pos[All_999_IPSL>=0.2]+1\n",
    "    All_999_consis_neg[All_999_IPSL<=-0.2] = All_999_consis_neg[All_999_IPSL<=-0.2]+1\n",
    "    \n",
    "    All_999_consis_pos[All_999_E3SM>=0.2] = All_999_consis_pos[All_999_E3SM>=0.2]+1\n",
    "    All_999_consis_neg[All_999_E3SM<=-0.2] = All_999_consis_neg[All_999_E3SM<=-0.2]+1\n",
    "    \n",
    "    All_999_consis = np.maximum(All_999_consis_pos, All_999_consis_neg)\n",
    "    \n",
    "\n",
    "    data_xr_12 = xr.DataArray(qctl_IRR_1985_2014_999-qctl_NOI_1901_1930_999, coords={'y': data_lat, 'x': data_lon}, dims=[\"y\", \"x\"])\n",
    "    ax = plt.subplot(4,3,6, projection=proj, frameon=False)\n",
    "    plot_PR_single_no_cb(ax, data_xr_12, All_999_consis, 'ALL', varname, '99.9%')\n",
    "    ax.text(0.01, 0.92, 'f', color='dimgrey', fontsize=14, transform=ax.transAxes, weight='bold')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b093b568",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_for_paper()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}