cafepp_eg3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# e.g. 3: Examine 3d winds for jet structures...\n",
    "\n",
    "# last edited 25 July 2018"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
       "    return false;\n",
       "}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
    "    return false;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BEGIN\n",
      "sys.version= 3.6.6 |Anaconda, Inc.| (default, Jun 28 2018, 17:14:51) \n",
      "[GCC 7.2.0]\n",
      "hostname= oa-35-cdc\n",
      "END\n"
     ]
    }
   ],
   "source": [
    "print('BEGIN')\n",
    "\n",
    "###\n",
    "\n",
    "import socket\n",
    "import re\n",
    "import os\n",
    "import sys\n",
    "import inspect\n",
    "import pickle\n",
    "import glob\n",
    "\n",
    "print('sys.version=',sys.version)\n",
    "\n",
    "hostname=socket.gethostname()\n",
    "\n",
    "print('hostname=',hostname)\n",
    "\n",
    "__file__='jupyter_notebook' #this can be deleted when written to a python script and loaded as module.\n",
    "\n",
    "import cartopy.crs as ccrs\n",
    "from cartopy.util import add_cyclic_point\n",
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "from matplotlib.mlab import griddata\n",
    "\n",
    "from decadal_diag import \\\n",
    "  fractional_year_from_num2date, \\\n",
    "  generate_daily_month_indices, \\\n",
    "  shade_2d_latlon, \\\n",
    "  get_timestamp_number, \\\n",
    "  shade_2d_curvilinear, \\\n",
    "  shade_2d_simple, \\\n",
    "  nino_indices, \\\n",
    "  plot_xy\n",
    "\n",
    "from n_data_funcs import n_data_funcs\n",
    "  \n",
    "cafe_experiment_titles = {4:'v0_forecast', 1:'v1_forecast', 3:'v2_forecast', 9:'v1', \\\n",
    "  5:'v2', 6:'v3', 7:'aug17p', 2:'v1_assim', 10:'etkf', 11:'nov17n', 12:'jul18b', 13:'v0'}\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "CRED = '\\033[91m'\n",
    "CEND = '\\033[0m'\n",
    "\n",
    "output_name='eg3.pkl'\n",
    "input_name=output_name\n",
    "\n",
    "print('END')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'tuple'>\n",
      "2\n"
     ]
    }
   ],
   "source": [
    "a=()\n",
    "\n",
    "#a=(np.zeros(10),)\n",
    "\n",
    "a=a+(np.zeros(10),)\n",
    "a=a+(np.zeros(10),)\n",
    "\n",
    "print(type(a))\n",
    "\n",
    "print(len(a))\n",
    "\n",
    "b,c=a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BEGIN\n",
      "sys.version= 3.6.6 |Anaconda, Inc.| (default, Jun 28 2018, 17:14:51) \n",
      "[GCC 7.2.0]\n",
      "hostname= oa-35-cdc\n",
      "this is vm32\n",
      "Current Working Directory= /OSM/CBR/OA_DCFP/work/col414/cafepp\n",
      "\u001b[91mProcessing v1 control...\u001b[0m\n",
      "n_data_funcs.init: no ensembles in this example.\n",
      "calculate_filedatetime_info: case 1: no ensembles, one input file.\n",
      "calculate_quantity: Turning on diagnostics.\n",
      "calculate_quantity: Diag= True\n",
      "calculate_quantity: key,value= quantity latlon_region\n",
      "calculate_quantity: calculating a quantity...\n",
      "calculate_quantity: value_split= ['latlon_region']\n",
      "calculate_quantity: latlon_region chosen.\n",
      "n_data_funcs.calculate_quantity: Using whole lat/lon region or that specified in instance.\n",
      "cafe_control_v1_monthly_quantity.shape= (372, 24, 90, 144)\n",
      "Diag= False\n",
      "monthly_clim_anom: All years have 12 months.\n",
      "n_data_funcs.init: no ensembles in this example.\n",
      "calculate_filedatetime_info: case 1: no ensembles, one input file.\n",
      "calculate_quantity: Turning on diagnostics.\n",
      "calculate_quantity: Diag= True\n",
      "calculate_quantity: key,value= quantity latlon_region\n",
      "calculate_quantity: calculating a quantity...\n",
      "calculate_quantity: value_split= ['latlon_region']\n",
      "calculate_quantity: latlon_region chosen.\n",
      "n_data_funcs.calculate_quantity: Using whole lat/lon region or that specified in instance.\n",
      "cafe_control_v1_monthly_quantity.shape= (372, 24, 90, 144)\n",
      "Diag= False\n",
      "monthly_clim_anom: All years have 12 months.\n"
     ]
    },
    {
     "ename": "SystemExit",
     "evalue": "STOP!:jupyter_notebook line number: 266",
     "output_type": "error",
     "traceback": [
      "An exception has occurred, use %tb to see the full traceback.\n",
      "\u001b[0;31mSystemExit\u001b[0m\u001b[0;31m:\u001b[0m STOP!:jupyter_notebook line number: 266\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/OSM/CBR/OA_DCFP/apps/col414/anaconda3/envs/cafepp_36_new/lib/python3.6/site-packages/IPython/core/interactiveshell.py:2971: UserWarning: To exit: use 'exit', 'quit', or Ctrl-D.\n",
      "  warn(\"To exit: use 'exit', 'quit', or Ctrl-D.\", stacklevel=1)\n"
     ]
    }
   ],
   "source": [
    "print('BEGIN')\n",
    "\n",
    "import socket\n",
    "import re\n",
    "import os\n",
    "import sys\n",
    "\n",
    "print('sys.version=',sys.version)\n",
    "\n",
    "hostname=socket.gethostname()\n",
    "\n",
    "print('hostname=',hostname)\n",
    "\n",
    "if(re.match('raijin',hostname)):\n",
    "  print('this is rajin')\n",
    "  rundir='/short/v14/mac599/cafepp/rundir'\n",
    "  topdir=''\n",
    "elif(re.match('oa-3.-cdc',hostname)):\n",
    "  print('this is vm32')\n",
    "  rundir='/OSM/CBR/OA_DCFP/work/col414/cafepp'\n",
    "  topdir='/OSM/CBR/OA_DCFP/data/CAFEPP/CMIP6'\n",
    "  topdir_ncep='/OSM/CBR/OA_DCFP/data/CAFEPP/'\n",
    "else:\n",
    "  raise SystemExit('host not known:'+__file__+' line number: '+str(inspect.stack()[0][2]))\n",
    "\n",
    "os.chdir(rundir)\n",
    "\n",
    "os.chdir('/OSM/CBR/OA_DCFP/work/col414/cafepp')\n",
    "\n",
    "print('Current Working Directory=',os.getcwd())\n",
    "\n",
    "SWITCH_OFF=True #temporary for easy skipping of this section.\n",
    "SWITCH_OFF=False\n",
    "\n",
    "if(not SWITCH_OFF):\n",
    "  \n",
    "  max_ensemble=11\n",
    "  max_ensemble=3 #first 2 ensembles\n",
    "#===============================================================================\n",
    "\n",
    "  \n",
    "  SWITCH_OFF=True #temporary for easy skipping of this section.\n",
    "  #SWITCH_OFF=False\n",
    "\n",
    "  if(not SWITCH_OFF):\n",
    "    print(CRED+'Processing v0 forecast...'+CEND)\n",
    "\n",
    "    pval=4\n",
    "    cafe_vars=['ua','ta']\n",
    "\n",
    "    cafe_forecast_v0_quantity_monthly_vars=()\n",
    "    cafe_forecast_v0_quantity_monthlyclimatology_vars=()\n",
    "\n",
    "    for cafe_cnt,cafe_var in enumerate(cafe_vars):\n",
    "\n",
    "      cafe_forecast_v0_allfiles=[]\n",
    "      for ensemble in range(2,max_ensemble+1):\n",
    "        cafe_forecast_v0_idir='/OSM/CBR/OA_DCFP/data/CAFEPP/CMIP6/CMIP6/CMIP/CSIRO/CAFE-1-0/historical/r'+str(ensemble)+'i1p'+str(pval)+'f1/day/'+cafe_var+'/gn/v20171025'\n",
    "        cafe_forecast_v0_allfiles.append(cafe_forecast_v0_idir+'/'+cafe_var+'_day_historical_CAFE-1-0_r'+str(ensemble)+'i1p'+str(pval)+'f1_gn_20070101-20121231.nc')\n",
    "\n",
    "      print('cafe_forecast_v0_allfiles=',cafe_forecast_v0_allfiles)\n",
    "\n",
    "      cafe_forecast_v0_files=[]\n",
    "      for cnt,cafe_forecast_v0_file in enumerate(cafe_forecast_v0_allfiles):\n",
    "        cafe_forecast_v0_files.append([cafe_forecast_v0_file])\n",
    "\n",
    "      print('cafe_forecast_v0_files=',cafe_forecast_v0_files)\n",
    "\n",
    "      #raise SystemExit('STOP!:'+__file__+' line number: '+str(inspect.stack()[0][2]))\n",
    "\n",
    "      cbeg_cafe_forecast_v0=2007\n",
    "      cend_cafe_forecast_v0=2012\n",
    "      abeg_cafe_forecast_v0=2007\n",
    "      aend_cafe_forecast_v0=2012\n",
    "\n",
    "      grid_label='cafe'\n",
    "\n",
    "      cafe_forecast_v0=n_data_funcs(Diag=False,input_files=cafe_forecast_v0_files,input_var_name=cafe_var)\n",
    "\n",
    "      cafe_forecast_v0.calculate_filedatetime_info(Diag=False,calendar='julian')\n",
    "\n",
    "      cafe_forecast_v0.get_latlon_info(Diag=False)\n",
    "\n",
    "      cafe_forecast_v0.get_lev_info(lev='plev')\n",
    "\n",
    "      cafe_forecast_v0_daily_quantity=cafe_forecast_v0.calculate_quantity('dummy',Diag=True,quantity='latlon_region')\n",
    "\n",
    "      print('cafe_forecast_v0_daily_quantity.shape=',cafe_forecast_v0_daily_quantity.shape)\n",
    "\n",
    "      cafe_forecast_v0.daily_monthly_indices_info(timesep=1) #would be 1 for daily data and 24 for hourly data.\n",
    "\n",
    "      cafe_forecast_v0_monthly_from_daily=cafe_forecast_v0.daily_to_monthly(EndOption=1,input=cafe_forecast_v0_daily_quantity)\n",
    "\n",
    "      cafe_forecast_v0_quantity_monthlyclimatology, cafe_forecast_v0_quantity_monthly = cafe_forecast_v0.monthly_clim_anom( \\\n",
    "        Diag=False, input=cafe_forecast_v0_monthly_from_daily, AnnOut=True, ZeroClim=True, \\\n",
    "        cbeg=cbeg_cafe_forecast_v0, cend=cend_cafe_forecast_v0, abeg=abeg_cafe_forecast_v0, aend=aend_cafe_forecast_v0)\n",
    "\n",
    "    #   cafe_forecast_v0_quantity_monthlyclimatology, cafe_forecast_v0_quantity_monthlyanomaly=cafe_forecast_v0.monthly_clim_anom( \\\n",
    "    #     Diag=False,input=cafe_forecast_v0_monthly_from_daily, AnnOut=True, ZeroClim=False, \\\n",
    "    #     cbeg=cbeg_cafe_forecast_v0, cend=cend_cafe_forecast_v0, abeg=abeg_cafe_forecast_v0, aend=aend_cafe_forecast_v0)\n",
    "\n",
    "      cafe_forecast_v0_quantity_monthlyclimatology_vars = cafe_forecast_v0_quantity_monthlyclimatology_vars + (cafe_forecast_v0_quantity_monthlyclimatology,)\n",
    "      cafe_forecast_v0_quantity_monthly_vars = cafe_forecast_v0_quantity_monthly_vars + (cafe_forecast_v0_quantity_monthly,)\n",
    "    \n",
    "    raise SystemExit('STOP!:'+__file__+' line number: '+str(inspect.stack()[0][2]))\n",
    "\n",
    "#===============================================================================\n",
    "\n",
    "  \n",
    "  SWITCH_OFF=True #temporary for easy skipping of this section.\n",
    "  #SWITCH_OFF=False\n",
    "\n",
    "  if(not SWITCH_OFF):\n",
    "    print(CRED+'Processing v1 forecast...'+CEND)\n",
    "\n",
    "    pval=1\n",
    "    cafe_vars=['ua','ta']\n",
    "\n",
    "    cafe_forecast_v1_quantity_monthly_vars=()\n",
    "    cafe_forecast_v1_quantity_monthlyclimatology_vars=()\n",
    "    \n",
    "    for cafe_cnt,cafe_var in enumerate(cafe_vars):\n",
    "\n",
    "      cafe_forecast_v1_allfiles=[]\n",
    "      for ensemble in range(2,max_ensemble+1):\n",
    "        cafe_forecast_v1_idir='/OSM/CBR/OA_DCFP/data/CAFEPP/CMIP6/CMIP6/CMIP/CSIRO/CAFE-1-0/historical/r'+str(ensemble)+'i1p'+str(pval)+'f1/day/'+cafe_var+'/gn/v20171025'\n",
    "        cafe_forecast_v1_allfiles.append(cafe_forecast_v1_idir+'/'+cafe_var+'_day_historical_CAFE-1-0_r'+str(ensemble)+'i1p'+str(pval)+'f1_gn_20070101-20121231.nc')\n",
    "\n",
    "      print('cafe_forecast_v1_allfiles=',cafe_forecast_v1_allfiles)\n",
    "\n",
    "      cafe_forecast_v1_files=[]\n",
    "      for cnt,cafe_forecast_v1_file in enumerate(cafe_forecast_v1_allfiles):\n",
    "        cafe_forecast_v1_files.append([cafe_forecast_v1_file])\n",
    "\n",
    "      print('cafe_forecast_v1_files=',cafe_forecast_v1_files)\n",
    "\n",
    "      cbeg_cafe_forecast_v1=2007\n",
    "      cend_cafe_forecast_v1=2012\n",
    "      abeg_cafe_forecast_v1=2007\n",
    "      aend_cafe_forecast_v1=2012\n",
    "\n",
    "      grid_label='cafe'\n",
    "\n",
    "      cafe_forecast_v1=n_data_funcs(Diag=False,input_files=cafe_forecast_v1_files,input_var_name=cafe_var)\n",
    "\n",
    "      cafe_forecast_v1.calculate_filedatetime_info(Diag=False,calendar='julian')\n",
    "\n",
    "      cafe_forecast_v1.get_latlon_info(Diag=False)\n",
    "\n",
    "      cafe_forecast_v1.get_lev_info(lev='plev')\n",
    "\n",
    "      cafe_forecast_v1_daily_quantity=cafe_forecast_v1.calculate_quantity('dummy',Diag=True,quantity='latlon_region')\n",
    "\n",
    "      print('cafe_forecast_v1_daily_quantity.shape=',cafe_forecast_v1_daily_quantity.shape)\n",
    "\n",
    "      cafe_forecast_v1.daily_monthly_indices_info(timesep=1) #would be 1 for daily data and 24 for hourly data.\n",
    "\n",
    "      cafe_forecast_v1_monthly_from_daily=cafe_forecast_v1.daily_to_monthly(EndOption=1,input=cafe_forecast_v1_daily_quantity)\n",
    "\n",
    "      cafe_forecast_v1_quantity_monthlyclimatology, cafe_forecast_v1_quantity_monthly = cafe_forecast_v1.monthly_clim_anom( \\\n",
    "        Diag=False, input=cafe_forecast_v1_monthly_from_daily, AnnOut=True, ZeroClim=True, \\\n",
    "        cbeg=cbeg_cafe_forecast_v1, cend=cend_cafe_forecast_v1, abeg=abeg_cafe_forecast_v1, aend=aend_cafe_forecast_v1)    \n",
    "\n",
    "    #   cafe_forecast_v1_quantity_monthlyclimatology, cafe_forecast_v1_quantity_monthlyanomaly=cafe_forecast_v1.monthly_clim_anom( \\\n",
    "    #     Diag=False,input=cafe_forecast_v1_monthly_from_daily, AnnOut=True, ZeroClim=False, \\\n",
    "    #     cbeg=cbeg_cafe_forecast_v1, cend=cend_cafe_forecast_v1, abeg=abeg_cafe_forecast_v1, aend=aend_cafe_forecast_v1)\n",
    "\n",
    "      cafe_forecast_v1_quantity_monthlyclimatology_vars = cafe_forecast_v1_quantity_monthlyclimatology_vars + (cafe_forecast_v1_quantity_monthlyclimatology,)\n",
    "      cafe_forecast_v1_quantity_monthly_vars = cafe_forecast_v1_quantity_monthly_vars + (cafe_forecast_v1_quantity_monthly,)\n",
    "      \n",
    "    #raise SystemExit('STOP!:'+__file__+' line number: '+str(inspect.stack()[0][2]))\n",
    "\n",
    "#===============================================================================\n",
    "\n",
    "  \n",
    "  SWITCH_OFF=True #temporary for easy skipping of this section.\n",
    "  #SWITCH_OFF=False\n",
    "\n",
    "  if(not SWITCH_OFF):\n",
    "    print(CRED+'Processing v0 control...'+CEND)\n",
    "\n",
    "    pval=13\n",
    "    cafe_vars=['ua','ta']\n",
    "\n",
    "    cafe_control_v0_quantity_monthly_vars=()\n",
    "    cafe_control_v0_quantity_monthlyclimatology_vars=()\n",
    "    \n",
    "    for cafe_cnt,cafe_var in enumerate(cafe_vars):\n",
    "    \n",
    "      cafe_control_v0_idir='/OSM/CBR/OA_DCFP/data/CAFEPP/CMIP6/CMIP6/CMIP/CSIRO/CAFE-1-0/piControl/r1i1p'+str(pval)+'f1/Amon/'+cafe_var+'/gn/v20171025'\n",
    "      cafe_control_v0_files=[[cafe_control_v0_idir+'/'+cafe_var+'_Amon_piControl_CAFE-1-0_r1i1p'+str(pval)+'f1_gn_047001-050012.nc']]\n",
    "\n",
    "      cbeg_cafe_control_v0=471\n",
    "      cend_cafe_control_v0=500\n",
    "      abeg_cafe_control_v0=471\n",
    "      aend_cafe_control_v0=500\n",
    "\n",
    "      grid_label='cafe'\n",
    "\n",
    "      cafe_control_v0=n_data_funcs(Diag=False,input_files=cafe_control_v0_files,input_var_name=cafe_var)\n",
    "\n",
    "      cafe_control_v0.calculate_filedatetime_info(Diag=False,calendar='noleap')\n",
    "\n",
    "      cafe_control_v0.get_latlon_info(Diag=False)\n",
    "\n",
    "      cafe_control_v0.get_lev_info(lev='plev')\n",
    "\n",
    "      cafe_control_v0_monthly_quantity=cafe_control_v0.calculate_quantity('dummy',Diag=True,quantity='latlon_region')\n",
    "\n",
    "      print('cafe_control_v0_monthly_quantity.shape=',cafe_control_v0_monthly_quantity.shape)\n",
    "\n",
    "      cafe_control_v0_quantity_monthlyclimatology, cafe_control_v0_quantity_monthly = cafe_control_v0.monthly_clim_anom( \\\n",
    "        Diag=False, input=cafe_control_v0_monthly_quantity, AnnOut=True, ZeroClim=True, \\\n",
    "        cbeg=cbeg_cafe_control_v0, cend=cend_cafe_control_v0, abeg=abeg_cafe_control_v0, aend=aend_cafe_control_v0)    \n",
    "\n",
    "    #   cafe_control_v0_quantity_monthlyclimatology, cafe_control_v0_quantity_monthlyanomaly=cafe_control_v0.monthly_clim_anom( \\\n",
    "    #     Diag=False,input=cafe_control_v0_monthly_quantity, AnnOut=True, ZeroClim=False, \\\n",
    "    #     cbeg=cbeg_cafe_control_v0, cend=cend_cafe_control_v0, abeg=abeg_cafe_control_v0, aend=aend_cafe_control_v0)\n",
    "\n",
    "      cafe_control_v0_quantity_monthlyclimatology_vars = cafe_control_v0_quantity_monthlyclimatology_vars + (cafe_control_v0_quantity_monthlyclimatology,)\n",
    "      cafe_fcontrol_v0_quantity_monthly_vars = cafe_control_v0_quantity_monthly_vars + (cafe_control_v0_quantity_monthly,)\n",
    "\n",
    "    raise SystemExit('STOP!:'+__file__+' line number: '+str(inspect.stack()[0][2]))\n",
    "\n",
    "#===============================================================================\n",
    "  \n",
    "  SWITCH_OFF=True #temporary for easy skipping of this section.\n",
    "  SWITCH_OFF=False\n",
    "\n",
    "  if(not SWITCH_OFF):\n",
    "    print(CRED+'Processing v1 control...'+CEND)\n",
    "\n",
    "    pval=9\n",
    "    cafe_vars=['ua','ta']\n",
    "\n",
    "    cafe_control_v1_quantity_monthly_vars=()\n",
    "    cafe_control_v1_quantity_monthlyclimatology_vars=()\n",
    "    \n",
    "    for cafe_cnt,cafe_var in enumerate(cafe_vars):\n",
    "\n",
    "      cafe_control_v1_idir='/OSM/CBR/OA_DCFP/data/CAFEPP/CMIP6/CMIP6/CMIP/CSIRO/CAFE-1-0/piControl/r1i1p'+str(pval)+'f1/Amon/'+cafe_var+'/gn/v20171025'\n",
    "      cafe_control_v1_files=[[cafe_control_v1_idir+'/'+cafe_var+'_Amon_piControl_CAFE-1-0_r1i1p'+str(pval)+'f1_gn_047001-050012.nc']]\n",
    "\n",
    "      cbeg_cafe_control_v1=471\n",
    "      cend_cafe_control_v1=500\n",
    "      abeg_cafe_control_v1=471\n",
    "      aend_cafe_control_v1=500\n",
    "\n",
    "      grid_label='cafe'\n",
    "\n",
    "      cafe_control_v1=n_data_funcs(Diag=False,input_files=cafe_control_v1_files,input_var_name=cafe_var)\n",
    "\n",
    "      cafe_control_v1.calculate_filedatetime_info(Diag=False,calendar='noleap')\n",
    "\n",
    "      cafe_control_v1.get_latlon_info(Diag=False)\n",
    "\n",
    "      cafe_control_v1.get_lev_info(lev='plev')\n",
    "\n",
    "      cafe_control_v1_monthly_quantity=cafe_control_v1.calculate_quantity('dummy',Diag=True,quantity='latlon_region')\n",
    "\n",
    "      print('cafe_control_v1_monthly_quantity.shape=',cafe_control_v1_monthly_quantity.shape)\n",
    "\n",
    "      cafe_control_v1_quantity_monthlyclimatology, cafe_control_v1_quantity_monthly = cafe_control_v1.monthly_clim_anom( \\\n",
    "        Diag=False, input=cafe_control_v1_monthly_quantity, AnnOut=True, ZeroClim=True, \\\n",
    "        cbeg=cbeg_cafe_control_v1, cend=cend_cafe_control_v1, abeg=abeg_cafe_control_v1, aend=aend_cafe_control_v1)    \n",
    "\n",
    "    #   cafe_control_v1_quantity_monthlyclimatology, cafe_control_v1_quantity_monthlyanomaly=cafe_control_v1.monthly_clim_anom( \\\n",
    "    #     Diag=False,input=cafe_control_v1_monthly_quantity, AnnOut=True, ZeroClim=False, \\\n",
    "    #     cbeg=cbeg_cafe_control_v1, cend=cend_cafe_control_v1, abeg=abeg_cafe_control_v1, aend=aend_cafe_control_v1)\n",
    "\n",
    "      cafe_control_v1_quantity_monthlyclimatology_vars = cafe_control_v1_quantity_monthlyclimatology_vars + (cafe_control_v1_quantity_monthlyclimatology,)\n",
    "      cafe_fcontrol_v1_quantity_monthly_vars = cafe_control_v1_quantity_monthly_vars + (cafe_control_v1_quantity_monthly,)\n",
    "\n",
    "    raise SystemExit('STOP!:'+__file__+' line number: '+str(inspect.stack()[0][2]))\n",
    "\n",
    "#===============================================================================\n",
    "\n",
    "  \n",
    "  SWITCH_OFF=True #temporary for easy skipping of this section.\n",
    "  #SWITCH_OFF=False\n",
    "  print(CRED+'Processing v2 control...'+CEND)\n",
    "\n",
    "  if(not SWITCH_OFF):\n",
    "    pval=5\n",
    "    cafe_vars=['ua','ta']\n",
    "\n",
    "    cafe_control_v2_quantity_monthly_vars=()\n",
    "    cafe_control_v2_quantity_monthlyclimatology_vars=()\n",
    "\n",
    "    for cafe_cnt,cafe_var in enumerate(cafe_vars):\n",
    "\n",
    "      cafe_control_v2_idir='/OSM/CBR/OA_DCFP/data/CAFEPP/CMIP6/CMIP6/CMIP/CSIRO/CAFE-1-0/piControl/r1i1p'+str(pval)+'f1/Amon/'+cafe_var+'/gn/v20171025'\n",
    "      cafe_control_v2_files=[[cafe_control_v2_idir+'/'+cafe_var+'_Amon_piControl_CAFE-1-0_r1i1p'+str(pval)+'f1_gn_047001-050012.nc']]\n",
    "\n",
    "      cbeg_cafe_control_v2=471\n",
    "      cend_cafe_control_v2=500\n",
    "      abeg_cafe_control_v2=471\n",
    "      aend_cafe_control_v2=500\n",
    "\n",
    "      grid_label='cafe'\n",
    "\n",
    "      cafe_control_v2=n_data_funcs(Diag=False,input_files=cafe_control_v2_files,input_var_name=cafe_var)\n",
    "\n",
    "      cafe_control_v2.calculate_filedatetime_info(Diag=False,calendar='noleap')\n",
    "\n",
    "      cafe_control_v2.get_latlon_info(Diag=False)\n",
    "\n",
    "      cafe_control_v2.get_lev_info(lev='plev')\n",
    "\n",
    "      cafe_control_v2_monthly_quantity=cafe_control_v2.calculate_quantity('dummy',Diag=True,quantity='latlon_region')\n",
    "\n",
    "      print('cafe_control_v2_monthly_quantity.shape=',cafe_control_v2_monthly_quantity.shape)\n",
    "\n",
    "      cafe_control_v2_quantity_monthlyclimatology, cafe_control_v2_quantity_monthly = cafe_control_v2.monthly_clim_anom( \\\n",
    "        Diag=False, input=cafe_control_v2_monthly_quantity, AnnOut=True, ZeroClim=True, \\\n",
    "        cbeg=cbeg_cafe_control_v2, cend=cend_cafe_control_v2, abeg=abeg_cafe_control_v2, aend=aend_cafe_control_v2)    \n",
    "\n",
    "    #   cafe_control_v2_quantity_monthlyclimatology, cafe_control_v2_quantity_monthlyanomaly=cafe_control_v2.monthly_clim_anom( \\\n",
    "    #     Diag=False,input=cafe_control_v2_monthly_quantity, AnnOut=True, ZeroClim=False, \\\n",
    "    #     cbeg=cbeg_cafe_control_v2, cend=cend_cafe_control_v2, abeg=abeg_cafe_control_v2, aend=aend_cafe_control_v2)\n",
    "\n",
    "      cafe_control_v2_quantity_monthlyclimatology_vars = cafe_control_v2_quantity_monthlyclimatology_vars + (cafe_control_v2_quantity_monthlyclimatology,)\n",
    "      cafe_fcontrol_v2_quantity_monthly_vars = cafe_control_v2_quantity_monthly_vars + (cafe_control_v2_quantity_monthly,)\n",
    "      \n",
    "    raise SystemExit('STOP!:'+__file__+' line number: '+str(inspect.stack()[0][2]))\n",
    "\n",
    "#===============================================================================\n",
    "\n",
    "#Clobber=False\n",
    "Clobber=True\n",
    "\n",
    "#output_type='pdf'\n",
    "#output_type='png'\n",
    "\n",
    "#title=True\n",
    "#title=False #turn off title across whole set.\n",
    "\n",
    "if(Clobber):\n",
    "\n",
    "  print(CRED+'overwriting output PKL file and regenerating...'+CEND)\n",
    "\n",
    "  cafe_control_v0_lev=cafe_control_v0.lev\n",
    "  cafe_control_v0_lat=cafe_control_v0.lat\n",
    "  \n",
    "  cafe_control_v1_lev=cafe_control_v1.lev\n",
    "  cafe_control_v1_lat=cafe_control_v1.lat\n",
    "  \n",
    "  cafe_control_v2_lev=cafe_control_v2.lev\n",
    "  cafe_control_v2_lat=cafe_control_v2.lat\n",
    "  \n",
    "  cafe_forecast_v0_lev=cafe_forecast_v0.lev\n",
    "  cafe_forecast_v0_lat=cafe_forecast_v0.lat\n",
    "  \n",
    "  cafe_forecast_v1_lev=cafe_forecast_v1.lev\n",
    "  cafe_forecast_v1_lat=cafe_forecast_v1.lat\n",
    "  \n",
    "  pkl_objects=( \\\n",
    "    cafe_vars, \\\n",
    "    cafe_control_v0_quantity_monthlyclimatology_vars, cafe_control_v0_lev, cafe_control_v0_lat, \\\n",
    "    cafe_control_v1_quantity_monthlyclimatology_vars, cafe_control_v1_lev, cafe_control_v1_lat, \\\n",
    "    cafe_control_v2_quantity_monthlyclimatology_vars, cafe_control_v2_lev, cafe_control_v2_lat, \\\n",
    "    cafe_forecast_v0_quantity_monthlyclimatology_vars, cafe_forecast_v0_lev, cafe_forecast_v0_lat, \\\n",
    "    cafe_forecast_v1_quantity_monthlyclimatology_vars, cafe_forecast_v1_lev, cafe_forecast_v1_lat, \\\n",
    "    )\n",
    "  \n",
    "if((os.path.exists(output_name) and Clobber) or (not os.path.exists(output_name))):\n",
    "  print(CRED+'Pkl file exists and deleting...'+CEND)\n",
    "  \n",
    "  if(os.path.exists(output_name)): os.remove(output_name)\n",
    "\n",
    "  pickling_out = open(output_name, \"wb\")\n",
    "  pickle.dump( pkl_objects, pickling_out)\n",
    "  pickling_out.close()\n",
    "\n",
    "if(os.path.exists(output_name) and not Clobber):\n",
    "  print(CRED+'Pkl file exists and keeping...'+CEND)\n",
    "    \n",
    "print('END')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BEGIN\n",
      "\u001b[91moverwriting output PKL file and regenerating...\u001b[0m\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'cafe_control_v0' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32mjupyter_notebook\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     14\u001b[0m   \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mCRED\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m'overwriting output PKL file and regenerating...'\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mCEND\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     15\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 16\u001b[0;31m   \u001b[0mcafe_control_v0_lev\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcafe_control_v0\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlev\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     17\u001b[0m   \u001b[0mcafe_control_v0_lat\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcafe_control_v0\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlat\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     18\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'cafe_control_v0' is not defined"
     ]
    }
   ],
   "source": [
    "print('BEGIN')\n",
    "\n",
    "Clobber=False\n",
    "Clobber=True\n",
    "\n",
    "output_type='pdf'\n",
    "output_type='png'\n",
    "\n",
    "#title=True\n",
    "#title=False #turn off title across whole set.\n",
    "\n",
    "if(Clobber):\n",
    "\n",
    "  print(CRED+'overwriting output PKL file and regenerating...'+CEND)\n",
    "\n",
    "  cafe_control_v0_lev=cafe_control_v0.lev\n",
    "  cafe_control_v0_lat=cafe_control_v0.lat\n",
    "  \n",
    "  cafe_control_v1_lev=cafe_control_v1.lev\n",
    "  cafe_control_v1_lat=cafe_control_v1.lat\n",
    "  \n",
    "  cafe_control_v2_lev=cafe_control_v2.lev\n",
    "  cafe_control_v2_lat=cafe_control_v2.lat\n",
    "  \n",
    "  cafe_forecast_v0_lev=cafe_forecast_v0.lev\n",
    "  cafe_forecast_v0_lat=cafe_forecast_v0.lat\n",
    "  \n",
    "  cafe_forecast_v1_lev=cafe_forecast_v1.lev\n",
    "  cafe_forecast_v1_lat=cafe_forecast_v1.lat\n",
    "  \n",
    "  pkl_objects=( \\\n",
    "    cafe_vars, \\\n",
    "    cafe_control_v0_quantity_monthlyclimatology_vars, cafe_control_v0_lev, cafe_control_v0_lat, \\\n",
    "    cafe_control_v1_quantity_monthlyclimatology_vars, cafe_control_v1_lev, cafe_control_v1_lat, \\\n",
    "    cafe_control_v2_quantity_monthlyclimatology_vars, cafe_control_v2_lev, cafe_control_v2_lat, \\\n",
    "    cafe_forecast_v0_quantity_monthlyclimatology_vars, cafe_forecast_v0_lev, cafe_forecast_v0_lat, \\\n",
    "    cafe_forecast_v1_quantity_monthlyclimatology_vars, cafe_forecast_v1_lev, cafe_forecast_v1_lat, \\\n",
    "    )\n",
    "  \n",
    "if((os.path.exists(output_name) and Clobber) or (not os.path.exists(output_name))):\n",
    "  print(CRED+'Pkl file exists and deleting...'+CEND)\n",
    "  \n",
    "  if(os.path.exists(output_name)): os.remove(output_name)\n",
    "\n",
    "  pickling_out = open(output_name, \"wb\")\n",
    "  pickle.dump( pkl_objects, pickling_out)\n",
    "  pickling_out.close()\n",
    "\n",
    "if(os.path.exists(output_name) and not Clobber):\n",
    "  print(CRED+'Pkl file exists and keeping...'+CEND)\n",
    "    \n",
    "print('END')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[91mPkl read from new or existing...create different output names compared to input...\u001b[0m\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'cafe_contv0_data' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32mjupyter_notebook\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      9\u001b[0m   \u001b[0mpickling_in\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m \u001b[0mcafe_contv0_data_var1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcafe_contv0_data_var2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcafe_contv0_data\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     12\u001b[0m \u001b[0mcafe_contv1_data_var1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcafe_contv0_data_var2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcafe_contv1_data\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     13\u001b[0m \u001b[0mcafe_contv2_data_var1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcafe_contv0_data_var2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcafe_contv2_data\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'cafe_contv0_data' is not defined"
     ]
    }
   ],
   "source": [
    "print(CRED+'Pkl read from new or existing...create different output names compared to input...'+CEND)\n",
    "\n",
    "# output_type='pdf'\n",
    "# output_type='png'\n",
    "\n",
    "if(os.path.exists(output_name)):\n",
    "  pickling_in = open(input_name, \"rb\")\n",
    "  cafe_vars, \\\n",
    "  cafe_contv0_data, cafe_contv0_lev, cafe_contv0_lat, \\\n",
    "  cafe_contv1_data, cafe_contv1_lev, cafe_contv1_lat, \\\n",
    "  cafe_contv2_data, cafe_contv2_lev, cafe_contv2_lat, \\\n",
    "  cafe_forcv0_data, cafe_forcv0_lev, cafe_forcv0_lat, \\\n",
    "  cafe_forcv1_data, cafe_forcv1_lev, cafe_forcv1_lat \\\n",
    "  = pickle.load(pickling_in)\n",
    "  pickling_in.close()  \n",
    "\n",
    "cafe_contv0_data_var1,cafe_contv0_data_var2 = cafe_contv0_data\n",
    "cafe_contv1_data_var1,cafe_contv0_data_var2 = cafe_contv1_data\n",
    "cafe_contv2_data_var1,cafe_contv0_data_var2 = cafe_contv2_data\n",
    "cafe_forcv0_data_var1,cafe_forcv0_data_var2 = cafe_forcv0_data\n",
    "cafe_forcv1_data_var1,cafe_forcv1_data_var2 = cafe_forcv1_data\n",
    "\n",
    "#print(cafe_contv1_data.shape)\n",
    "\n",
    "#https://matplotlib.org/api/_as_gen/matplotlib.pyplot.yscale.html\n",
    "  \n",
    "xysize=(12,8)\n",
    "xysize=(6,4)\n",
    "\n",
    "clevs_full=[-48,-36,-24,-12,0,12,24,36,48]\n",
    "clevs_diff=[-4,-3,-2,-1,0,1,2,3,4]\n",
    "\n",
    "shade_2d_simple(np.average(cafe_contv0_data_var1,axis=-1), \\\n",
    "  title='cont v0 zonal winds', units='m/s', \\\n",
    "  reverse_yaxis=True, yscale='log', \\\n",
    "  xyvals=(cafe_contv0_lat, cafe_contv0_lev/100.), \\\n",
    "  xysize=xysize, xlab='Latitude', ylab='Pressure (hPa)', ytik=[10,100,500,1000], grid=True, clevs=clevs_full, n_contour=[0])\n",
    "#   output=('xsec',output_type))\n",
    "\n",
    "shade_2d_simple(np.average(cafe_contv1_data_var1,axis=-1), \\\n",
    "  title='cont v1 zonal winds', units='m/s', \\\n",
    "  reverse_yaxis=True, yscale='log', \\\n",
    "  xyvals=(cafe_contv1_lat, cafe_contv1_lev/100.), \\\n",
    "  xysize=xysize, xlab='Latitude', ylab='Pressure (hPa)', ytik=[10,100,500,1000], grid=True, clevs=clevs_full, n_contour=[0])\n",
    "\n",
    "shade_2d_simple(np.average(cafe_contv2_data_var1,axis=-1), \\\n",
    "  title='cont v2 zonal winds', units='m/s', \\\n",
    "  reverse_yaxis=True, yscale='log', \\\n",
    "  xyvals=(cafe_contv2_lat, cafe_contv2_lev/100.), \\\n",
    "  xysize=xysize, xlab='Latitude', ylab='Pressure (hPa)', ytik=[10,100,500,1000], grid=True, clevs=clevs_full, n_contour=[0])\n",
    "\n",
    "#===============================================================================\n",
    "\n",
    "# shade_2d_simple(np.average(cafe_contv2_data,axis=-1)-np.average(cafe_contv1_data,axis=-1), \\\n",
    "#   title='cont v2-v1 zonal winds', units='m/s', \\\n",
    "#   reverse_yaxis=True, yscale='log', \\\n",
    "#   xyvals=(cafe_contv2_lat, cafe_contv2_lev/100.), \\\n",
    "#   xysize=xysize, xlab='Latitude', ylab='Pressure (hPa)', ytik=[10,100,500,1000], grid=True, clevs=clevs_diff, n_contour=[0])\n",
    "\n",
    "#===============================================================================\n",
    "\n",
    "shade_2d_simple(np.average(np.average(cafe_forcv0_data_var1,axis=-1),axis=0), \\\n",
    "  title='forc v0 zonal winds', units='m/s', \\\n",
    "  reverse_yaxis=True, yscale='log', \\\n",
    "  xyvals=(cafe_forcv0_lat, cafe_forcv0_lev/100.), \\\n",
    "  xysize=xysize, xlab='Latitude', ylab='Pressure (hPa)', ytik=[10,100,500,1000], grid=True, clevs=clevs_full, n_contour=[0])\n",
    "\n",
    "shade_2d_simple(np.average(np.average(cafe_forcv1_data_var1,axis=-1),axis=0), \\\n",
    "  title='forc v1 zonal winds', units='m/s', \\\n",
    "  reverse_yaxis=True, yscale='log', \\\n",
    "  xyvals=(cafe_forcv1_lat, cafe_forcv1_lev/100.), \\\n",
    "  xysize=xysize, xlab='Latitude', ylab='Pressure (hPa)', ytik=[10,100,500,1000], grid=True, clevs=clevs_full, n_contour=[0])\n",
    "\n",
    "shade_2d_simple(np.average(np.average(cafe_forcv1_data_var1,axis=-1),axis=0)-np.average(cafe_contv1_data_var1,axis=-1), \\\n",
    "  title='forc v1-cont v1 zonal winds', units='m/s', \\\n",
    "  reverse_yaxis=True, yscale='log', \\\n",
    "  xyvals=(cafe_forcv1_lat, cafe_forcv1_lev/100.), \\\n",
    "  xysize=xysize, xlab='Latitude', ylab='Pressure (hPa)', ytik=[10,100,500,1000], grid=True, clevs=clevs_diff, \\\n",
    "  contour_array=(cafe_forcv1_lat, cafe_forcv1_lev/100., np.average(cafe_contv0_data_var1,axis=-1)))\n",
    "\n",
    "shade_2d_simple(np.average(cafe_contv0_data_var1,axis=-1), \\\n",
    "  title='cont v1 zonal winds', units='m/s', \\\n",
    "  reverse_yaxis=True, yscale='log', \\\n",
    "  xyvals=(cafe_forcv0_lat, cafe_forcv0_lev/100.), \\\n",
    "  xysize=xysize, xlab='Latitude', ylab='Pressure (hPa)', ytik=[10,100,500,1000], grid=True, clevs=clevs_full, \\\n",
    "  contour_array=(cafe_forcv1_lat, cafe_forcv1_lev/100., np.average(cafe_contv0_data_var2,axis=-1)-273.15))\n",
    "\n",
    "# shade_2d_simple(np.average(cafe_contv1_data,axis=2), \\\n",
    "#   title='zonal winds', units='m/s', \\\n",
    "#   reverse_yaxis=True, yscale='linear', \\\n",
    "#   xyvals=(cafe_control_file_var_lat, cafe_control_file_var_lev/100.), \\\n",
    "#   xysize=xysize, xlab='Latitude', ylab='Pressure (hPa)') #, xlim=(-90,90),ylim=(0,200))\n",
    "\n",
    "# shade_2d_simple(np.average(cafe_contv1_data,axis=2), \\\n",
    "#   title='zonal winds', units='m/s', \\\n",
    "#   reverse_yaxis=True, yscale='symlog', \\\n",
    "#   xyvals=(cafe_control_file_var_lat, cafe_control_file_var_lev/100.), \\\n",
    "#   xysize=xysize, xlab='Latitude', ylab='Pressure (hPa)')\n",
    "\n",
    "# shade_2d_simple(np.average(cafe_contv1_data,axis=2), \\\n",
    "#   title='zonal winds', units='m/s', \\\n",
    "#   reverse_yaxis=True, yscale='log', \\\n",
    "#   xyvals=(cafe_control_file_var_lat, cafe_control_file_var_lev/100.), \\\n",
    "#   xysize=xysize, xlab='Latitude', ylab='Pressure (hPa)', ylim=(10,1000), \\\n",
    "#   ytik=[10,100,500,1000])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}