gchp_vs_gchp_comparison_without_gcpy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare two GCHP datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview of this Notebook\n",
    "\n",
    "* Import dependencies\n",
    "* Define functions: utility, regridding, and plotting\n",
    "* Define data directories and global variables\n",
    "* Compare restart files\n",
    "* Store areas for normalization\n",
    "* Compare SpeciesConc diagnostic collection\n",
    "* Create PDFs for subset of species\n",
    "* Other diagnostic collections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import cubedsphere as cs\n",
    "import xesmf as xe\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "from matplotlib import ticker\n",
    "from matplotlib.backends.backend_pdf import PdfPages\n",
    "from matplotlib.colors import ListedColormap\n",
    "from cartopy import crs\n",
    "from cartopy.mpl.geoaxes import GeoAxes\n",
    "\n",
    "%matplotlib inline\n",
    "import warnings; warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "# Also define colormap. Colormap file source: https://bitbucket.org/gcst/gcpy\n",
    "rgb_WhGrYlRd = np.genfromtxt('/n/home08/elundgren/GC/python/WhGrYlRd.txt',delimiter=' ')\n",
    "WhGrYlRd = ListedColormap(rgb_WhGrYlRd/255.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define general utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_gchp_filepath(outputdir, collection, day, time):\n",
    "    filepath = os.path.join(outputdir, 'GCHP.{}.{}_{}z.nc4'.format(collection,day,time))\n",
    "    return filepath"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_paths( refpath, devpath):\n",
    "    if not os.path.exists(refpath):\n",
    "        print('ERROR! Path 1 does not exist: {}'.format(refpath))\n",
    "    else:\n",
    "        print('Path 1 exists: {}'.format(refpath))\n",
    "    if not os.path.exists(devpath):\n",
    "        print('ERROR! Path 2 does not exist: {}'.format(devpath))\n",
    "    else:\n",
    "        print('Path 2 exists: {}'.format(devpath))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def compare_varnames(refdata, refstr, devdata, devstr):\n",
    "    \n",
    "    # Find common variables in collection by generating lists and list overlap\n",
    "    refvars = [k for k in refdata.data_vars.keys()]\n",
    "    devvars= [k for k in devdata.data_vars.keys()]\n",
    "    commonvars = sorted(list(set(refvars).intersection(set(devvars))))\n",
    "    refonly = [v for v in refvars if v not in devvars]\n",
    "    devonly = [v for v in devvars if v not in refvars]\n",
    "    dimmismatch = [v for v in commonvars if refdata[v].ndim != devdata[v].ndim]\n",
    "    commonvars2D = [v for v in commonvars if refdata[v].ndim == 3]\n",
    "    commonvars3D = [v for v in commonvars if devdata[v].ndim == 4]\n",
    "    \n",
    "    # Print information on common and mismatching variables, as well as dimensions\n",
    "    print('{} common variables ({} are 2-dim and {} are 3-dim)'.format(len(commonvars), len(commonvars2D), len(commonvars3D)))\n",
    "    if len(refonly) > 0:\n",
    "        print('{} variables in {} only (skip)'.format(len(refnly)),refstr)\n",
    "        print('   Variable names: {}'.format(refonly))\n",
    "    if len(devonly) > 0:\n",
    "        print('{} variables in {} only (skip)'.format(len(devonly)),devstr)\n",
    "        print('   Variable names: {}'.format(devonly))\n",
    "    if len(dimmismatch) > 0:\n",
    "        print('{} common variables have different dimensions'.format(len(dimmismatch)))\n",
    "        print('   Variable names: {}'.format(dimmismatch))\n",
    "        \n",
    "    return [commonvars, commonvars2D, commonvars3D]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_stats(refdata, refstr, devdata, devstr, varname):\n",
    "    refvar = refdata[varname]\n",
    "    devvar = devdata[varname]\n",
    "    units = refdata[varname].units\n",
    "    print('Data units: {}'.format(units))\n",
    "    print('Array sizes:')\n",
    "    print('    {}:  {}'.format(refstr,refvar.shape))\n",
    "    print('    {}:  {}'.format(devstr,devvar.shape))\n",
    "    print('Global stats:')\n",
    "    print('  Mean:')\n",
    "    print('    {}:  {}'.format(refstr,np.round(refvar.values.mean(),20)))\n",
    "    print('    {}:  {}'.format(devstr,np.round(devvar.values.mean(),20)))\n",
    "    print('  Min:')\n",
    "    print('    {}:  {}'.format(refstr,np.round(refvar.values.min(),20)))\n",
    "    print('    {}:  {}'.format(devstr,np.round(devvar.values.min(),20)))\n",
    "    print('  Max:')\n",
    "    print('    {}:  {}'.format(refstr,np.round(refvar.values.max(),20)))\n",
    "    print('    {}:  {}'.format(devstr,np.round(devvar.values.max(),20)))\n",
    "    print('  Sum:')\n",
    "    print('    {}:  {}'.format(refstr,np.round(refvar.values.sum(),20)))\n",
    "    print('    {}:  {}'.format(devstr,np.round(devvar.values.sum(),20)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_collection_data(refdir, devdir, collection, day, time):\n",
    "    reffile = get_gchp_filepath(refdir, collection, day, time)\n",
    "    devfile = get_gchp_filepath(devdir, collection, day, time)\n",
    "    check_paths(reffile, devfile)\n",
    "    refdata = xr.open_dataset(reffile)\n",
    "    devdata = xr.open_dataset(devfile)\n",
    "    return [refdata, devdata]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define functions for regridding\n",
    "\n",
    "These functions utilize the xESMF package and the cubedsphere package developed by Jaiwei Zhuang, a graduate student at Harvard University. More specifically I use the cubedsphere package CSGrid class and the xESMF package Regridder class. The latter serves as a wrapper for the ESMF/ESMPy packages. These tools enable simple conservative regridding between lat/lon and cubed sphere horizontal grids as specifically defined by the NASA Global Modeling and Assimilation Office (GMAO). See http://xesmf.readthedocs.io/en/latest/ and https://github.com/JiaweiZhuang/cubedsphere for more information. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_grid_LL(llres):\n",
    "    [dlat,dlon] = list(map(float, llres.split('x')))\n",
    "    lon_b = np.linspace(-180 - dlon/2, 180 - dlon/2, int(360/dlon) + 1, endpoint=True)\n",
    "    lat_b = np.linspace(-90 - dlat/2, 90 + dlat/2, int(180/dlat) + 2, endpoint=True).clip(-90,90)\n",
    "    lat = (lat_b[1:] + lat_b[:-1]) / 2\n",
    "    lon = (lon_b[1:] + lon_b[:-1]) / 2\n",
    "    llgrid = {'lat': lat, \n",
    "              'lon': lon, \n",
    "              'lat_b': lat_b, \n",
    "              'lon_b': lon_b}\n",
    "    return llgrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_grid_CS(csres):\n",
    "    csgrid = cs.csgrid_GMAO(csres)\n",
    "    csgrid_list = [None]*6\n",
    "    for i in range(6):\n",
    "        csgrid_list[i] = {'lat': csgrid['lat'][i], \n",
    "                          'lon': csgrid['lon'][i],\n",
    "                          'lat_b': csgrid['lat_b'][i], \n",
    "                          'lon_b': csgrid['lon_b'][i]}\n",
    "    return [csgrid, csgrid_list]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_regridder_C2L( csres_in, llres_out, weightsdir='.', reuse_weights=False ):\n",
    "    csgrid, csgrid_list = make_grid_CS(csres_in)\n",
    "    llgrid = make_grid_LL(llres_out)\n",
    "    regridder_list = []\n",
    "    for i in range(6):\n",
    "        weightsfile = os.path.join(weightsdir, 'conservative_c{}_{}_{}.nc'.format(str(csres_in), llres_out, str(i)))\n",
    "        regridder = xe.Regridder(csgrid_list[i], llgrid, method='conservative', filename=weightsfile, reuse_weights=reuse_weights)\n",
    "        regridder_list.append(regridder)\n",
    "    return regridder_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_regridder_L2L( llres_in, llres_out, weightsdir='.', reuse_weights=False ):\n",
    "    llgrid_in = make_grid_LL(llres_in)\n",
    "    llgrid_out = make_grid_LL(llres_out)\n",
    "    weightsfile = os.path.join(weightsdir,'conservative_{}_{}.nc'.format(llres_in, llres_out))\n",
    "    regridder = xe.Regridder(llgrid_in, llgrid_out, method='conservative', filename=weightsfile, reuse_weights=reuse_weights)\n",
    "    return regridder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define functions for plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function to plot and compare single vertical level of GCHP\n",
    "\n",
    "Notes: This function currently does not regrid GCHP; therefore the two GCHP files must be at the same cubed sphere resolution. You can use this function to plot interactively or to generate a multi-page pdf of plots. It takes a list of variable names to plot for a single collection only. You can plot for any level and any time slice in the file. By default the colorbars for the plots will have the same range, but you can turn this feature off. Also by default the colorbar of the fractional difference between the model outputs will be limited to +/-2, but you can change this as well via the passed parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_single_level(refdata, refstr, devdata, devstr, varlist=None, weightsdir='.', ilev=0, \n",
    "                         itime=0, savepdf=False, pdfname='gchp_vs_gchp_map.pdf', match_cbar=True, \n",
    "                         full_ratio_range=False, normalize_by_area=False, area_ref=None, \n",
    "                         area_dev=None, check_units=True, flip_vert=False):\n",
    "    \n",
    "    # If no varlist is passed, plot all (surface only for 3D)\n",
    "    if varlist == None:\n",
    "        [varlist, commonvars2D, commonvars3D] = compare_varnames(refdata, devdata)\n",
    "        print('Plotting all common variables (surface only if 3D)')\n",
    "    n_var = len(varlist)\n",
    "\n",
    "    # Get cubed sphere grids and regridder \n",
    "    # for now, do not regrid for gchp vs gchp. Assume same grid.\n",
    "    csres_ref = refdata['lon'].size\n",
    "    [csgrid_ref, csgrid_list_ref] = make_grid_CS(csres_ref)\n",
    "    csres_dev = devdata['lon'].size\n",
    "    [csgrid_dev, csgrid_list_dev] = make_grid_CS(csres_dev)\n",
    "\n",
    "    # Create pdf (if saving)\n",
    "    if savepdf:\n",
    "        print('\\nCreating {} for {} variables'.format(pdfname,n_var))\n",
    "        pdf = PdfPages(pdfname)\n",
    "\n",
    "    # Loop over variables\n",
    "    for ivar in range(n_var):\n",
    "        if savepdf: print('{} '.format(ivar), end='')\n",
    "        varname = varlist[ivar]\n",
    "        \n",
    "        # Do some checks: dimensions and units\n",
    "        varndim_ref = refdata[varname].ndim\n",
    "        varndim_dev = devdata[varname].ndim\n",
    "        if check_units: \n",
    "            assert varndim_ref == varndim_dev, 'Dimensions do not agree for {}!'.format(varname)\n",
    "        units_ref = refdata[varname].units\n",
    "        units_dev = devdata[varname].units\n",
    "        if check_units: \n",
    "            assert units_ref == units_dev, 'Units do not match for {}!'.format(varname)\n",
    "            \n",
    "        # if normalizing by area, adjust units to be per m2, and adjust title string\n",
    "        units = units_ref\n",
    "        varndim = varndim_ref\n",
    "        subtitle_extra = ''\n",
    "                    \n",
    "        # Slice the data\n",
    "        if varndim == 4: \n",
    "            if flip_vert: \n",
    "                ds_ref = refdata[varname].isel(time=itime,lev=71-ilev)\n",
    "                ds_dev = devdata[varname].isel(time=itime,lev=71-ilev)\n",
    "            else: \n",
    "                ds_ref = refdata[varname].isel(time=itime,lev=ilev)\n",
    "                ds_dev = devdata[varname].isel(time=itime,lev=ilev)\n",
    "        elif varndim == 3: \n",
    "            refds = refdata[varname].isel(time=itime)\n",
    "            devds = devdata[varname].isel(time=itime)\n",
    "            \n",
    "        # if normalizing by area, transform on the native grid and adjust units and subtitle string\n",
    "        exclude_list = ['WetLossConvFrac','Prod_','Loss_']\n",
    "        if normalize_by_area and not any(s in varname for s in exclude_list):\n",
    "            ds_ref.values = ds_ref.values / area_ref\n",
    "            ds_dev.values = ds_dev.values / area_dev\n",
    "            units = '{} m-2'.format(units)\n",
    "            subtitle_extra = ', Normalized by Area'\n",
    "            \n",
    "        # Regrid the slices (skip for gchp vs gchp for now)\n",
    "        csdata_ref = ds_ref.data.reshape(6,csres_ref,csres_ref)\n",
    "        csdata_dev = ds_dev.data.reshape(6,csres_dev,csres_dev)\n",
    "        vmin_ref = csdata_ref.min()\n",
    "        vmin_dev = csdata_dev.min()\n",
    "        vmin_cmn = np.min([vmin_ref, vmin_dev])\n",
    "        vmax_ref = csdata_ref.max()\n",
    "        vmax_dev = csdata_dev.max()\n",
    "        vmax_cmn = np.max([vmax_ref, vmax_dev])\n",
    "        if match_cbar: [vmin, vmax] = [vmin_cmn, vmax_cmn]\n",
    "        \n",
    "        # Create 2x2 figure\n",
    "        figs, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, figsize=[12,9], \n",
    "                                                      subplot_kw={'projection': crs.PlateCarree()})\n",
    "        # Give the figure a title\n",
    "        offset = 0.96\n",
    "        fontsize=25\n",
    "        if varndim == 4:\n",
    "            if ilev == 0: levstr = 'Surface'\n",
    "            elif ilev == 22: levstr = '500 hPa'\n",
    "            else: levstr = 'Level ' +  str(ilev-1)\n",
    "            figs.suptitle('{}, {}'.format(varname,levstr), fontsize=fontsize, y=offset)\n",
    "        elif varndim == 3: \n",
    "            figs.suptitle('{}'.format(varname), fontsize=fontsize, y=offset)\n",
    "        else:\n",
    "            print('varndim is 2 for {}! Must be 3 or 4.'.format(varname))\n",
    "            \n",
    "        # Subplot (0,0): Ref\n",
    "        ax0.coastlines()\n",
    "        if not match_cbar: [vmin, vmax] = [vmin_ref, vmax_ref]        \n",
    "        masked_csdata = np.ma.masked_where(np.abs(csgrid_ref['lon'] - 180) < 2, csdata_ref)\n",
    "        for i in range(6):\n",
    "            plot0 = ax0.pcolormesh(csgrid_ref['lon_b'][i,:,:], csgrid_ref['lat_b'][i,:,:], masked_csdata[i,:,:], \n",
    "                                   cmap=WhGrYlRd,vmin=vmin, vmax=vmax)\n",
    "        ax0.set_title('{} (Ref){}\\nC{}'.format(refstr,subtitle_extra,str(csres_ref)))\n",
    "        cb = plt.colorbar(plot0, ax=ax0, orientation='horizontal', pad=0.10)\n",
    "        if (vmax-vmin) < 0.1 or (vmax-vmin) > 100:\n",
    "            cb.locator = ticker.MaxNLocator(nbins=4)\n",
    "            cb.update_ticks()\n",
    "        cb.set_label(units)     \n",
    "        \n",
    "        # Subplot (0,1): Dev\n",
    "        ax1.coastlines()\n",
    "        if not match_cbar: [vmin, vmax] = [vmin_dev, vmax_dev]        \n",
    "        masked_csdata = np.ma.masked_where(np.abs(csgrid_dev['lon'] - 180) < 2, csdata_dev)\n",
    "        for i in range(6):\n",
    "            plot1 = ax1.pcolormesh(csgrid_dev['lon_b'][i,:,:], csgrid_dev['lat_b'][i,:,:], \n",
    "                                   masked_csdata[i,:,:], cmap=WhGrYlRd,vmin=vmin, vmax=vmax)\n",
    "        ax1.set_title('{} (Dev){}\\nC{}'.format(devstr,subtitle_extra,str(csres_dev)))\n",
    "        cb = plt.colorbar(plot1, ax=ax1, orientation='horizontal', pad=0.10)\n",
    "        if (vmax-vmin) < 0.1 or (vmax-vmin) > 100:\n",
    "            cb.locator = ticker.MaxNLocator(nbins=4)\n",
    "            cb.update_ticks()\n",
    "        cb.set_label(units)  \n",
    "            \n",
    "        # Subplot (1,0): Difference\n",
    "        gc_absdiff = csdata_dev - csdata_ref\n",
    "        diffabsmax = max([np.abs(gc_absdiff.min()), np.abs(gc_absdiff.max())])\n",
    "        [vmin, vmax] = [-diffabsmax, diffabsmax]\n",
    "        ax2.coastlines()\n",
    "        # assume the same grid for now in gchp vs gchp\n",
    "        masked_csdata = np.ma.masked_where(np.abs(csgrid_dev['lon'] - 180) < 2, gc_absdiff)\n",
    "        for i in range(6):\n",
    "            plot2 = ax2.pcolormesh(csgrid_dev['lon_b'][i,:,:], csgrid_dev['lat_b'][i,:,:], \n",
    "                                   masked_csdata[i,:,:], cmap=WhGrYlRd,vmin=vmin, vmax=vmax)\n",
    "        ax2.set_title('Difference\\n(Dev - Ref)')   \n",
    "        cb = plt.colorbar(plot2, ax=ax2, orientation='horizontal', pad=0.10)\n",
    "        if (vmax-vmin) < 0.1 or (vmax-vmin) > 100:\n",
    "            cb.locator = ticker.MaxNLocator(nbins=4)\n",
    "            cb.update_ticks()\n",
    "        cb.set_label(units)  \n",
    "    \n",
    "        # Subplot (1,1): Fractional Difference (restrict to +/-2)\n",
    "        gc_fracdiff = (csdata_dev - csdata_ref) / csdata_ref\n",
    "        if full_ratio_range: [vmin, vmax] = [None, None]\n",
    "        else: [vmin, vmax] = [-2, 2]\n",
    "        ax3.coastlines()\n",
    "        # assume the same grid for now in gchp vs gchp\n",
    "        masked_csdata = np.ma.masked_where(np.abs(csgrid_dev['lon'] - 180) < 2, gc_fracdiff)\n",
    "        for i in range(6):\n",
    "            plot3 = ax3.pcolormesh(csgrid_dev['lon_b'][i,:,:], csgrid_dev['lat_b'][i,:,:], \n",
    "                                   masked_csdata[i,:,:], cmap='RdBu_r',vmin=vmin, vmax=vmax)\n",
    "        ax3.set_title('Fractional Difference\\n(Dev - Ref)/Ref')   \n",
    "        cb = plt.colorbar(plot3, ax=ax3, orientation='horizontal', pad=0.10)\n",
    "        cb.set_clim(vmin=vmin, vmax=vmax)\n",
    "        cb.set_label('unitless')     \n",
    "            \n",
    "        if savepdf:    \n",
    "            pdf.savefig(figs)\n",
    "            plt.close(figs)\n",
    "            \n",
    "    if savepdf: pdf.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function to plot and compare zonal means of GCHP\n",
    "\n",
    "Note: Zonal mean is defined as the mean value across a constant latitudinal band. Calculating zonal mean from cubed sphere is not straight-forward. This function therefore regrids model outputs to a common lat-lon grid resolution prior to plotting. Many of the same features available for plotting a single level above are also available for this function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def compare_zonal_mean(refdata, refstr, devdata, devstr, varlist=None, \n",
    "                       weightsdir='.', itime=0, llres_cmp='1x1.25', savepdf=False, \n",
    "                       pdfname='gchp_vs_gchp_map.pdf', match_cbar=True, full_ratio_range=False, \n",
    "                       normalize_by_area=False, area_ref=None, area_dev=None, flip_vert=False):\n",
    "\n",
    "    # If no varlist is passed, plot all 3D variables in the dataset\n",
    "    if varlist == None:\n",
    "        [commonvars, commonvars2D, varlist] = compare_varnames(refdata, devdata)\n",
    "        print('Plotting all 3D variables')\n",
    "    n_var = len(varlist)\n",
    "    \n",
    "    # Get lat-lon grid\n",
    "    llgrid_cmp = make_grid_LL(llres_cmp)\n",
    "\n",
    "    # Get cubed sphere grid and regridder for first data set\n",
    "    csres_ref = refdata['lon'].size\n",
    "    [csgrid_ref, csgrid_list_ref] = make_grid_CS(csres_ref)\n",
    "    cs_regridder_list_ref = make_regridder_C2L(csres_ref, llres_cmp, weightsdir=weightsdir, \n",
    "                                               reuse_weights=True)\n",
    "\n",
    "    # Get cubed sphere grid and regridder for first data set\n",
    "    csres_dev = devdata['lon'].size\n",
    "    [csgrid_dev, csgrid_list_dev] = make_grid_CS(csres_dev)\n",
    "    cs_regridder_list_dev = make_regridder_C2L(csres_dev, llres_cmp, weightsdir=weightsdir, \n",
    "                                               reuse_weights=True)\n",
    "    \n",
    "    # Universal plot setup\n",
    "    xtick_positions = np.arange(-90,91,30)\n",
    "    xticklabels = ['{}$\\degree$'.format(x) for x in xtick_positions]\n",
    "    ytick_positions = np.arange(0,61,20)\n",
    "    yticklabels = [str(y) for y in ytick_positions]\n",
    "    \n",
    "    # Create pdf (if saving)\n",
    "    if savepdf:\n",
    "        print('\\nCreating {} for {} variables'.format(pdfname, n_var))\n",
    "        pdf = PdfPages(pdfname)\n",
    "\n",
    "    # Loop over variables\n",
    "    for ivar in range(n_var):\n",
    "        if savepdf: print('{} '.format(ivar), end='')\n",
    "        varname = varlist[ivar]\n",
    "        \n",
    "        # Do some checks: dimensions and units\n",
    "        varndim_ref = refdata[varname].ndim\n",
    "        varndim_dev = devdata[varname].ndim\n",
    "        nlev = 72\n",
    "        assert varndim_ref == varndim_dev, 'GCHP dimensions do not agree for {}!'.format(varname)\n",
    "        units_ref = refdata[varname].units\n",
    "        units_dev = devdata[varname].units\n",
    "        assert units_ref == units_dev, 'GCHP units do not match for {}!'.format(varname)\n",
    "        \n",
    "        # Set plot extent\n",
    "        extent=(-90,90,0,nlev)\n",
    "        \n",
    "        # if normalizing by area, adjust units to be per m2, and adjust title string\n",
    "        units = units_ref\n",
    "        varndim = varndim_ref\n",
    "        subtitle_extra = ''\n",
    "        \n",
    "        # Slice the data\n",
    "        ds_ref = refdata[varname].isel(time=itime)\n",
    "        ds_dev = devdata[varname].isel(time=itime)\n",
    "\n",
    "        # if normalizing by area, transform on the native grid and adjust units and subtitle string\n",
    "        exclude_list = ['WetLossConvFrac','Prod_','Loss_']\n",
    "        if normalize_by_area and not any(s in varname for s in exclude_list):\n",
    "            ds_ref.values = ds_ref.values / area_ref.values[np.newaxis,:,:]\n",
    "            ds_dev.values = ds_dev.values / area_dev.values[np.newaxis,:,:]\n",
    "            units = '{} m-2'.format(units)\n",
    "            subtitle_extra = ', Normalized by Area'\n",
    "            \n",
    "        # Regrid the slices\n",
    "        if flip_vert: \n",
    "            ds_ref.data = ds_ref.data[::-1,:,:]\n",
    "            ds_dev.data = ds_dev.data[::-1,:,:]\n",
    "        csdata_ref = ds_ref.data.reshape(nlev,6,csres_ref,csres_ref).swapaxes(0,1)\n",
    "        csdata_dev = ds_dev.data.reshape(nlev,6,csres_dev,csres_dev).swapaxes(0,1)\n",
    "        lldata_ref = np.zeros([nlev, llgrid_cmp['lat'].size, llgrid_cmp['lon'].size])\n",
    "        lldata_dev = np.zeros([nlev, llgrid_cmp['lat'].size, llgrid_cmp['lon'].size])\n",
    "        for i in range(6):\n",
    "            regridder_ref = cs_regridder_list_ref[i]\n",
    "            lldata_ref += regridder_ref(csdata_ref[i])\n",
    "            regridder_dev = cs_regridder_list_dev[i]\n",
    "            lldata_dev += regridder_dev(csdata_dev[i])\n",
    "        \n",
    "        # Calculate zonal mean of the regridded data\n",
    "        zm_ref = lldata_ref.mean(axis=2)\n",
    "        zm_dev = lldata_dev.mean(axis=2)\n",
    "            \n",
    "        # Get min and max for colorbar limits\n",
    "        [vmin_ref, vmax_ref] = [zm_ref.min(), zm_ref.max()]\n",
    "        [vmin_dev, vmax_dev] = [zm_dev.min(), zm_dev.max()]\n",
    "        vmin_cmn = np.min([vmin_ref, vmin_dev])\n",
    "        vmax_cmn = np.max([vmax_ref, vmax_dev])\n",
    "        if match_cbar: [vmin, vmax] = [vmin_cmn, vmax_cmn]\n",
    "        \n",
    "        # Create 2x2 figure\n",
    "        figs, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, figsize=[12,12], \n",
    "                                                      subplot_kw={'projection': crs.PlateCarree()})\n",
    "        # Give the page a title\n",
    "        offset = 0.96\n",
    "        fontsize=25\n",
    "        figs.suptitle('{}, Zonal Mean'.format(varname), fontsize=fontsize, y=offset)\n",
    "\n",
    "        # Subplot 0: Ref\n",
    "        if not match_cbar: [vmin, vmax] = [vmin_ref, vmax_ref]\n",
    "        plot0 = ax0.imshow(zm_ref, cmap=WhGrYlRd, extent=extent, vmin=vmin, vmax=vmax)\n",
    "        ax0.set_title('{} (Ref){}\\n{} regridded from {}'.format(refstr, subtitle_extra, \n",
    "                                                                llres_cmp, csres_ref))\n",
    "        ax0.set_aspect('auto')\n",
    "        ax0.set_xticks(xtick_positions)\n",
    "        ax0.set_xticklabels(xticklabels)\n",
    "        ax0.set_yticks(ytick_positions)\n",
    "        ax0.set_yticklabels(yticklabels)\n",
    "        cb = plt.colorbar(plot0, ax=ax0, orientation='horizontal', pad=0.10)\n",
    "        if (vmax-vmin) < 0.001 or (vmax-vmin) > 1000:\n",
    "            cb.locator = ticker.MaxNLocator(nbins=4)\n",
    "            cb.update_ticks()\n",
    "        cb.set_label(units)\n",
    "        \n",
    "        # Subplot 1: Dev\n",
    "        if not match_cbar: [vmin, vmax] = [vmin_dev, vmax_dev]\n",
    "        plot1 = ax1.imshow(zm_dev, cmap=WhGrYlRd, extent=extent, vmin=vmin, vmax=vmax)\n",
    "        ax1.set_title('{} (Dev){}\\n{} regridded from {}'.format(devstr, subtitle_extra, \n",
    "                                                                llres_cmp, csres_dev))\n",
    "        ax1.set_aspect('auto')\n",
    "        ax1.set_xticks(xtick_positions)\n",
    "        ax1.set_xticklabels(xticklabels)\n",
    "        ax1.set_yticks(ytick_positions)\n",
    "        ax1.set_yticklabels(yticklabels)\n",
    "        cb = plt.colorbar(plot1, ax=ax1, orientation='horizontal', pad=0.10)\n",
    "        if (vmax-vmin) < 0.001 or (vmax-vmin) > 1000:\n",
    "            cb.locator = ticker.MaxNLocator(nbins=4)\n",
    "            cb.update_ticks()\n",
    "        cb.set_label(units)\n",
    "            \n",
    "        # Subplot 2: Difference\n",
    "        zm_absdiff = zm_dev - zm_ref\n",
    "        diffabsmax = max([np.abs(zm_absdiff.min()), np.abs(zm_absdiff.max())])\n",
    "        [vmin, vmax] = [-diffabsmax, diffabsmax]\n",
    "        plot2 = ax2.imshow(zm_absdiff, cmap='RdBu_r', extent=extent, vmin=vmin, vmax=vmax)\n",
    "        ax2.set_title('Difference\\n(Dev - Ref)')\n",
    "        ax2.set_aspect('auto')\n",
    "        ax2.set_xticks(xtick_positions)\n",
    "        ax2.set_xticklabels(xticklabels)\n",
    "        ax2.set_yticks(ytick_positions)\n",
    "        ax2.set_yticklabels(yticklabels)\n",
    "        cb = plt.colorbar(plot2, ax=ax2, orientation='horizontal', pad=0.10)\n",
    "        if (vmax-vmin) < 0.001 or (vmax-vmin) > 1000:\n",
    "            cb.locator = ticker.MaxNLocator(nbins=4)\n",
    "            cb.update_ticks()\n",
    "        cb.set_label(units)\n",
    "        \n",
    "        # Subplot 3: Fractional Difference (restrict to +/-2)\n",
    "        zm_fracdiff = (zm_dev - zm_ref) / zm_ref\n",
    "        if full_ratio_range: [vmin, vmax] = [None, None]\n",
    "        else: [vmin, vmax] = [-2, 2]\n",
    "        plot3 = ax3.imshow(zm_fracdiff, vmin=vmin, vmax=vmax, cmap='RdBu_r', extent=extent)\n",
    "        ax3.set_title('Fractional Difference\\n(Dev - Ref)/Ref')\n",
    "        ax3.set_aspect('auto')\n",
    "        ax3.set_xticks(xtick_positions)\n",
    "        ax3.set_xticklabels(xticklabels)\n",
    "        ax3.set_yticks(ytick_positions)\n",
    "        ax3.set_yticklabels(yticklabels)\n",
    "        cb = plt.colorbar(plot3, ax=ax3, orientation='horizontal', pad=0.10)\n",
    "        cb.set_clim(vmin=vmin, vmax=vmax)\n",
    "        cb.set_label('unitless')      \n",
    "            \n",
    "        if savepdf:    \n",
    "            pdf.savefig(figs)\n",
    "            plt.close(figs)\n",
    "            \n",
    "    if savepdf: pdf.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define directories and global information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Path 1 exists: /n/home08/elundgren/GC/testruns/12.0.3/gchp_standard_ref/OutputDir\n",
      "Path 2 exists: /n/home08/elundgren/GC/testruns/12.0.3/gchp_standard_dev/OutputDir\n"
     ]
    }
   ],
   "source": [
    "# Shared high-level directory\n",
    "testdir = '/n/home08/elundgren/GC/testruns/12.0.3'\n",
    "\n",
    "# Ref and dev run output directories\n",
    "refdir = os.path.join(testdir,'gchp_standard_ref/OutputDir')\n",
    "devdir = os.path.join(testdir,'gchp_standard_dev/OutputDir')\n",
    "\n",
    "# Ref and dev strings (for use in plots and messages)\n",
    "refstr='12.0.2'\n",
    "devstr='12.0.3'\n",
    "\n",
    "# Set directory to store plots\n",
    "plotsdir = os.path.join(testdir,'plots')\n",
    "\n",
    "# Set directory to store regridding weights\n",
    "weightsdir = '/n/home08/elundgren/GC/python/regrid_weights'\n",
    "\n",
    "# Define cubed sphere resolutions (used for restart file comparison)\n",
    "csres_ref = 48\n",
    "csres_dev = 48\n",
    "\n",
    "# Check that paths exist\n",
    "check_paths(refdir, devdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Print Ref netcdf filenames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GCHP.Emissions.20160716_1200z.nc4\n",
      "GCHP.LevelEdgeDiags.20160716_1200z.nc4\n",
      "GCHP.Aerosols.20160716_1200z.nc4\n",
      "GCHP.JValues.20160716_1200z.nc4\n",
      "GCHP.ConcAfterChem.20160716_1200z.nc4\n",
      "GCHP.WetLossConv.20160716_1200z.nc4\n",
      "GCHP.StateMet_inst.20160801_0000z.nc4\n",
      "GCHP.DryDep.20160716_1200z.nc4\n",
      "GCHP.ProdLoss.20160716_1200z.nc4\n",
      "GCHP.SpeciesConc_inst.20160801_0000z.nc4\n",
      "GCHP.WetLossLS.20160716_1200z.nc4\n",
      "GCHP.SpeciesConc_avg.20160716_1200z.nc4\n",
      "GCHP.CloudConvFlux.20160716_1200z.nc4\n",
      "GCHP.StateMet_avg.20160716_1200z.nc4\n",
      "GCHP.StateChm.20160716_1200z.nc4\n",
      "GCHP.AerosolMass.20160716_1200z.nc4\n"
     ]
    }
   ],
   "source": [
    "reffiles = [k for k in os.listdir(refdir) if '.nc' in k]\n",
    "for k in reffiles:\n",
    "    print(k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Print Dev netcdf filenames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GCHP.Emissions.20160716_1200z.nc4\n",
      "GCHP.LevelEdgeDiags.20160716_1200z.nc4\n",
      "GCHP.Aerosols.20160716_1200z.nc4\n",
      "GCHP.JValues.20160716_1200z.nc4\n",
      "GCHP.ConcAfterChem.20160716_1200z.nc4\n",
      "GCHP.WetLossConv.20160716_1200z.nc4\n",
      "GCHP.StateMet_inst.20160801_0000z.nc4\n",
      "GCHP.DryDep.20160716_1200z.nc4\n",
      "GCHP.ProdLoss.20160716_1200z.nc4\n",
      "GCHP.SpeciesConc_inst.20160801_0000z.nc4\n",
      "GCHP.WetLossLS.20160716_1200z.nc4\n",
      "GCHP.SpeciesConc_avg.20160716_1200z.nc4\n",
      "GCHP.CloudConvFlux.20160716_1200z.nc4\n",
      "GCHP.StateMet_avg.20160716_1200z.nc4\n",
      "GCHP.StateChm.20160716_1200z.nc4\n",
      "GCHP.AerosolMass.20160716_1200z.nc4\n"
     ]
    }
   ],
   "source": [
    "devfiles = [k for k in os.listdir(devdir) if '.nc' in k]\n",
    "for k in devfiles:\n",
    "    print(k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create Regridding Weights (executing this section is not required)\n",
    "\n",
    "The xESMF Regridder object can create weights for regridding from one grid resolution to another, or it can be provided weights that are pre-generated to save time. Generating weights therefore only needs to be done once for each combination of input and output grid type and resolution. This can be done by calling the make_regridder functions and passing reuse_weights=False. Once the weights exist, all subsequent calls to make a regridder can be passed reuse_weights=True along with the location of the pre-existing weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is here only as an example. Plotting will create the weights on-the-fly if they do not yet exist.\n",
    "#regridder_C2L_list = make_regridder_C2L(48, '1x1.25', weightsdir=weightsdir, reuse_weights=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare Restart Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Path 1 exists: /n/home08/elundgren/GC/testruns/12.0.3/gchp_standard_ref/OutputDir/../initial_GEOSChem_rst.c48_standard.nc\n",
      "Path 2 exists: /n/home08/elundgren/GC/testruns/12.0.3/gchp_standard_dev/OutputDir/../initial_GEOSChem_rst.c48_standard.nc\n"
     ]
    }
   ],
   "source": [
    "reffile_rst = os.path.join(refdir,'..','initial_GEOSChem_rst.c{}_standard.nc'.format(csres_ref))\n",
    "devfile_rst = os.path.join(devdir,'..','initial_GEOSChem_rst.c{}_standard.nc'.format(csres_dev))\n",
    "check_paths(reffile_rst, devfile_rst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "277 common variables (0 are 2-dim and 277 are 3-dim)\n"
     ]
    }
   ],
   "source": [
    "refdata_rst = xr.open_dataset(reffile_rst)\n",
    "devdata_rst = xr.open_dataset(devfile_rst)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(refdata_rst, refstr, devdata_rst, devstr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data units: mol mol-1\n",
      "Array sizes:\n",
      "    12.0.2:  (1, 72, 288, 48)\n",
      "    12.0.3:  (1, 72, 288, 48)\n",
      "Global stats:\n",
      "  Mean:\n",
      "    12.0.2:  1.1734857530143472e-10\n",
      "    12.0.3:  1.1734857530143472e-10\n",
      "  Min:\n",
      "    12.0.2:  0.0\n",
      "    12.0.3:  0.0\n",
      "  Max:\n",
      "    12.0.2:  2.2552185119373291e-10\n",
      "    12.0.3:  2.2552185119373291e-10\n",
      "  Sum:\n",
      "    12.0.2:  0.00011680032912408933\n",
      "    12.0.3:  0.00011680032912408933\n"
     ]
    }
   ],
   "source": [
    "get_stats(refdata_rst, refstr, devdata_rst, devstr, 'SPC_CFC11')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing restarts requires flipping GCHP in the vertical. Do this by passing flip_vert=True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#compare_zonal_mean(refdata_rst, refstr, devdata_rst, devstr,\n",
    "#                   varlist=['SPC_O3'], weightsdir=weightsdir, flip_vert=True )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#compare_single_level(refdata_rst, refstr, devdata_rst, devstr,\n",
    "#                     varlist=['SPC_O3'], ilev=0, weightsdir=weightsdir, flip_vert=True )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Store Grid Areas for Normalizing by Area\n",
    "\n",
    "Eventually this will be replaced by retrieving areas from xESMF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Path 1 exists: /n/home08/elundgren/GC/testruns/12.0.3/gchp_standard_ref/OutputDir/GCHP.StateMet_avg.20160716_1200z.nc4\n",
      "Path 2 exists: /n/home08/elundgren/GC/testruns/12.0.3/gchp_standard_dev/OutputDir/GCHP.StateMet_avg.20160716_1200z.nc4\n"
     ]
    }
   ],
   "source": [
    "day = '20160716'\n",
    "time = '1200'\n",
    "[refdata, devdata] = get_collection_data(refdir, devdir, 'StateMet_avg', day, time)\n",
    "varname = 'Met_AREAM2'\n",
    "ref_area = refdata[varname].isel(lev=0,time=0)\n",
    "dev_area = devdata[varname].isel(lev=0,time=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspect time-averaged species concentration diagnostic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Path 1 exists: /n/home08/elundgren/GC/testruns/12.0.3/gchp_standard_ref/OutputDir/GCHP.SpeciesConc_avg.20160716_1200z.nc4\n",
      "Path 2 exists: /n/home08/elundgren/GC/testruns/12.0.3/gchp_standard_dev/OutputDir/GCHP.SpeciesConc_avg.20160716_1200z.nc4\n"
     ]
    }
   ],
   "source": [
    "day = '20160716'\n",
    "time = '1200'\n",
    "collection = 'SpeciesConc_avg'\n",
    "[refdata, devdata] = get_collection_data(refdir, devdir, collection, day, time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "161 common variables (0 are 2-dim and 161 are 3-dim)\n"
     ]
    }
   ],
   "source": [
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(refdata, refstr, devdata, devstr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['SpeciesConc_CFC11', 'SpeciesConc_CFC113', 'SpeciesConc_CFC114', 'SpeciesConc_CFC115', 'SpeciesConc_CFC12', 'SpeciesConc_HCFC123', 'SpeciesConc_HCFC141b', 'SpeciesConc_HCFC142b', 'SpeciesConc_HCFC22']\n"
     ]
    }
   ],
   "source": [
    "CFCs=[k for k in commonvars if 'CFC' in k]\n",
    "print(CFCs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data units: mol mol-1 dry\n",
      "Array sizes:\n",
      "    12.0.2:  (1, 72, 288, 48)\n",
      "    12.0.3:  (1, 72, 288, 48)\n",
      "Global stats:\n",
      "  Mean:\n",
      "    12.0.2:  5.978070816148318e-11\n",
      "    12.0.3:  1.149659673016501e-10\n",
      "  Min:\n",
      "    12.0.2:  0.0\n",
      "    12.0.3:  0.0\n",
      "  Max:\n",
      "    12.0.2:  2.0787468968386236e-10\n",
      "    12.0.3:  2.1939146332972115e-10\n",
      "  Sum:\n",
      "    12.0.2:  5.9501409850781783e-05\n",
      "    12.0.3:  0.00011442884715506807\n"
     ]
    }
   ],
   "source": [
    "varname='SpeciesConc_CFC11'\n",
    "get_stats(refdata, refstr, devdata, devstr, varname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#compare_single_level(refdata, refstr, devdata, devstr, varlist=[varname], weightsdir=weightsdir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compare_zonal_mean(refdata, refstr, devdata, devstr, varlist=[varname], weightsdir=weightsdir )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create PDF of plots for list of species"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define subset of variables in files\n",
    "#varlist = commonvars3D\n",
    "varlist = [varname]\n",
    "#varlist = [k for k in commonvars3D if 'CFC' in k]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Creating /n/home08/elundgren/GC/testruns/12.0.3/plots/SpeciesConc_CFC11_Surface.pdf for 1 variables\n",
      "0 "
     ]
    }
   ],
   "source": [
    "# Surface\n",
    "#pdfname = os.path.join(plotsdir,'{}_Surface.pdf'.format(ref_collection))\n",
    "pdfname = os.path.join(plotsdir,'{}_Surface.pdf'.format(varname))\n",
    "#pdfname = os.path.join(plotsdir,'{}_Surface.pdf'.format('CFCs'))\n",
    "compare_single_level(refdata, refstr, devdata, devstr, varlist=varlist, ilev=0,\n",
    "                     weightsdir=weightsdir, savepdf=True, pdfname=pdfname )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Creating /n/home08/elundgren/GC/testruns/12.0.3/plots/SpeciesConc_CFC11_ZonalMean.pdf for 1 variables\n",
      "0 "
     ]
    }
   ],
   "source": [
    "# Zonal mean\n",
    "#pdfname = os.path.join(plotsdir,'{}_ZonalMean.pdf'.format(ref_collection))\n",
    "pdfname = os.path.join(plotsdir,'{}_ZonalMean.pdf'.format(varname))\n",
    "#pdfname = os.path.join(plotsdir,'{}_ZonalMean.pdf'.format('CFCs'))\n",
    "compare_zonal_mean(refdata, refstr, devdata, devstr, varlist=varlist, \n",
    "                   weightsdir=weightsdir, savepdf=True, pdfname=pdfname )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other diagnostic collections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['GCHP.Emissions.20160716_1200z.nc4',\n",
       " 'GCHP.LevelEdgeDiags.20160716_1200z.nc4',\n",
       " 'GCHP.Aerosols.20160716_1200z.nc4',\n",
       " 'GCHP.JValues.20160716_1200z.nc4',\n",
       " 'GCHP.ConcAfterChem.20160716_1200z.nc4',\n",
       " 'GCHP.WetLossConv.20160716_1200z.nc4',\n",
       " 'GCHP.StateMet_inst.20160801_0000z.nc4',\n",
       " 'GCHP.DryDep.20160716_1200z.nc4',\n",
       " 'GCHP.ProdLoss.20160716_1200z.nc4',\n",
       " 'GCHP.SpeciesConc_inst.20160801_0000z.nc4',\n",
       " 'GCHP.WetLossLS.20160716_1200z.nc4',\n",
       " 'GCHP.SpeciesConc_avg.20160716_1200z.nc4',\n",
       " 'GCHP.CloudConvFlux.20160716_1200z.nc4',\n",
       " 'GCHP.StateMet_avg.20160716_1200z.nc4',\n",
       " 'GCHP.StateChm.20160716_1200z.nc4',\n",
       " 'GCHP.AerosolMass.20160716_1200z.nc4']"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Inspect the available files again, but only show those that match\n",
    "commonfiles = [k for k in reffiles if k in devfiles]\n",
    "commonfiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Template for inspecting any collection:\n",
    "\n",
    "# To specify file:\n",
    "#day = {edit}\n",
    "#time = {edit}\n",
    "#collection = {edit}\n",
    "\n",
    "# For retrieving the data:\n",
    "#[refdata, devdata] = get_collection_data(refdir, devdir, collection, day, time)\n",
    "\n",
    "# For inspecting interactively:\n",
    "#refdata\n",
    "#devdata\n",
    "#[commonvars, commonvars2D, commonvars3D] = compare_varnames(refdata, refstr, devdata, devstr)\n",
    "# varname = {edit}\n",
    "#get_stats(refdata, refstr, devdata, devstr, varname)\n",
    "# varlist = {edit}\n",
    "#compare_single_level(refdata, refstr, devdata, devstr, varlist=[varname], weightsdir=weightsdir)\n",
    "#compare_zonal_mean(refdata, refstr, devdata, devstr, varlist=[varname], weightsdir=weightsdir )\n",
    "\n",
    "# For creating PDFs:\n",
    "#prefix = {edit}\n",
    "#pdfname = os.path.join(plotsdir,'{}_Surface.pdf'.format(prefix))\n",
    "#compare_single_level(refdata, refstr, devdata, devstr, varlist=varlist, ilev=0,\n",
    "#                     weightsdir=weightsdir, savepdf=True, pdfname=pdfname )\n",
    "#pdfname = os.path.join(plotsdir,'{}_ZonalMean.pdf'.format(prefix))\n",
    "#compare_zonal_mean(refdata, refstr, devdata, devstr, varlist=varlist, \n",
    "#                   weightsdir=weightsdir, savepdf=True, pdfname=pdfname )\n",
    "\n",
    "# For other function options use shift-tab while hovering over the function name in a cell that calls it"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}