cafepp_monthly_nino34_from_sst.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 12,
        "hidden": false,
        "row": 0,
        "width": 12
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "source": [
    "# CAFEPP (utilising cafepp.py)\n",
    "\n",
    "***\n",
    "\n",
    "## This example will produce a set of monthly SST files across the V1 assimilation run (interpolated onto a 1x1 latxlon grid), and upon success, generate an average over the nino3.4 region and plot the time-series.\n",
    "\n",
    "***\n",
    "\n",
    "# Various settings required:\n",
    "\n",
    "## BATCH determines whether it will be sent to the queue via qsub command or run interactively.\n",
    "\n",
    "## BACKGROUND is alternative to BATCH on systems where there is no queue system so that multiple jobs can be run in quick succession.\n",
    "\n",
    "## CLEAN determines whether the run directory will be emptied prior to processing.\n",
    "\n",
    "Choose to run either 1. assimilation 2. forecast, then continue with plotting options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 4,
        "hidden": false,
        "row": 12,
        "width": 4
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BEGIN\n",
      "sys.version= 2.7.14 |Anaconda, Inc.| (default, Dec  7 2017, 17:05:42) \n",
      "[GCC 7.2.0]\n",
      "hostname= oa-32-cdc\n",
      "Removing run directory and reestablish it.\n",
      "Running in directory /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_20\n",
      "Using script cafepp_monthly_forecast.py\n",
      "END\n"
     ]
    }
   ],
   "source": [
    "print('BEGIN')\n",
    "\n",
    "###\n",
    "\n",
    "from __future__ import print_function #this is to allow print(,file=xxx) feature\n",
    "\n",
    "import datetime\n",
    "import shutil\n",
    "import os\n",
    "import sys\n",
    "import socket\n",
    "import re\n",
    "import inspect\n",
    "\n",
    "__file__='jupyter_notebook' #this can be deleted when written to a python script and loaded as module.\n",
    "\n",
    "print('sys.version=',sys.version)\n",
    "hostname=socket.gethostname()\n",
    "print('hostname=',hostname)\n",
    "\n",
    "#dvar='tos' #variable to generate as well as further proces and plot.\n",
    "#dvar='pr'\n",
    "#dvar='psl'\n",
    "#dvar='siconc'\n",
    "#dvar='msftyz' #was called msftyyz\n",
    "#dvar='ts'\n",
    "#dvar='thetao'\n",
    "#dvar='tos'\n",
    "#dvar='siconc'\n",
    "\n",
    "BATCH=True #submit to queue\n",
    "BATCH=False #run interactively but in a batch temporary area.\n",
    "\n",
    "BACKGROUND=True #run in background\n",
    "BACKGROUND=False #run in foreground (if not BATCH=True)\n",
    "\n",
    "CLEAN=False #don't remove rundir, just use it.\n",
    "CLEAN=True #remove rundir and recreate it.\n",
    "\n",
    "if(BATCH and BACKGROUND):\n",
    "  raise SystemExit('Cannot have both BATCH and BACKGROUND:'+__file__+' line number: '+str(inspect.stack()[0][2]))\n",
    "  \n",
    "if(BATCH):\n",
    "  print('Submitting to queue.')\n",
    "  \n",
    "if(BACKGROUND):\n",
    "  print('Putting in background.')\n",
    "\n",
    "if(CLEAN):\n",
    "  print('Removing run directory and reestablish it.')\n",
    "\n",
    "if(re.match('raijin',hostname)):\n",
    "  cmipdir='/short/v14/mac599' #this might be different to predir for other users.\n",
    "  predir='/short/v14/mac599' #this is directory area for temporary cafepp files.\n",
    "elif(re.match('oa-3.-cdc',hostname)):\n",
    "  cmipdir='/OSM/CBR/OA_DCFP/data/CAFEPP/CMIP6' #this might be different to predir for other users.\n",
    "  predir='/OSM/CBR/OA_DCFP/data/col414' #this is directory area for temporary cafepp files.\n",
    "elif(re.match('tube-hba',hostname)):\n",
    "  cmipdir='/OSM/HBA/OA_DECADAL_CLIMATE/work/col414' #this might be different to predir for other users.\n",
    "  predir='/OSM/HBA/OA_DECADAL_CLIMATE/work/col414' #this is directory area for temporary cafepp files.\n",
    "else:\n",
    "  raise SystemExit('hostname not known:'+__file__+' line number: '+str(inspect.stack()[0][2])) \n",
    "\n",
    "topdir=predir+'/'+'cafepp'\n",
    "#script='cafepp_monthly_assimilation.py'\n",
    "script='cafepp_monthly_forecast.py'\n",
    "#script='cafepp_monthly_forecast_v2.py'\n",
    "#script='cafepp_monthly_control.py'\n",
    "\n",
    "# rundir=topdir+'/'+'rundir20171124143249'\n",
    "# rundir=topdir+'/'+'rundir20171128165302'\n",
    "# rundir=topdir+'/'+'rundir'+datetime.datetime.today().strftime('%Y%m%d%H%M%S')\n",
    "#rundir=topdir+'/'+'rundir' #temporary running directory\n",
    "rundir=topdir+'/'+'rundir_20' #temporary running directory\n",
    "print('Running in directory '+rundir)\n",
    "print('Using script '+script)\n",
    "\n",
    "if(re.match('raijin',hostname)):\n",
    "  srcdir='/home/599/mac599/decadal' #location of main cafepp code\n",
    "  prodir=srcdir+'/paper_analysis' #project directory\n",
    "elif(re.match('oa-3.-cdc',hostname)):\n",
    "  srcdir='/OSM/CBR/OA_DCFP/work/col414/cafepp' #'/home/599/mac599/decadal' #location of main cafepp code\n",
    "  prodir='/OSM/CBR/OA_DCFP/work/col414/cafepp' #srcdir+'/paper_analysis' #project directory\n",
    "elif(re.match('tube-hba',hostname)):\n",
    "  srcdir='/OSM/HBA/OA_DECADAL_CLIMATE/work/col414/cafepp' #'/home/599/mac599/decadal' #location of main cafepp code\n",
    "  prodir='/OSM/HBA/OA_DECADAL_CLIMATE/work/col414/cafepp' #srcdir+'/paper_analysis' #project directory  \n",
    "else:\n",
    "  raise SystemExit('host not known:'+__file__+' line number: '+str(inspect.stack()[0][2]))\n",
    "\n",
    "print('END')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 4,
        "height": 4,
        "hidden": false,
        "row": 12,
        "width": 4
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "source": [
    "# Here we make directories, copy across necessary JSON and python code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 8,
        "height": 15,
        "hidden": false,
        "row": 12,
        "width": 4
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BEGIN\n",
      "copying json files, generating symlinks, cmor tables, queue script, if necessary.\n",
      "Copying /OSM/CBR/OA_DCFP/work/col414/cafepp/decadal_diag.py to  /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_20\n",
      "Copying /OSM/CBR/OA_DCFP/work/col414/cafepp/decadal_diag_py2.py to  /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_20\n",
      "Copying /OSM/CBR/OA_DCFP/work/col414/cafepp/decadal_diag_py3.py to  /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_20\n",
      "Copying /OSM/CBR/OA_DCFP/work/col414/cafepp/cafepp.py to  /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_20\n",
      "Copying /OSM/CBR/OA_DCFP/work/col414/cafepp/app_funcs.py to  /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_20\n",
      "Copying /OSM/CBR/OA_DCFP/work/col414/cafepp/ProcTime.py to  /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_20\n",
      "Copying /OSM/CBR/OA_DCFP/work/col414/cafepp/cafepp_monthly_forecast.py to  /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_20\n",
      "Copying /OSM/CBR/OA_DCFP/work/col414/cafepp/JsonTemplates/cafepp_experiments.json to  /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_20\n",
      "Copying /OSM/CBR/OA_DCFP/work/col414/cafepp/JsonTemplates/cafepp_csiro-gfdl.json to  /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_20\n",
      "Copying /OSM/CBR/OA_DCFP/work/col414/cafepp/JsonTemplates/cafepp_vars.json to  /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_20\n",
      "Copying /OSM/CBR/OA_DCFP/work/col414/cafepp/cafepp_monthly_assimilation.json /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_20/JsonTemplates/cafepp_monthly_assimilation.json\n",
      "prodir= /OSM/CBR/OA_DCFP/work/col414/cafepp\n",
      "rundir= /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_20\n",
      "END\n"
     ]
    }
   ],
   "source": [
    "print('BEGIN')\n",
    "\n",
    "print('copying json files, generating symlinks, cmor tables, queue script, if necessary.')\n",
    "if(os.path.exists(rundir) and CLEAN):\n",
    "  shutil.rmtree(rundir)\n",
    "  os.mkdir(rundir)\n",
    "  os.mkdir(rundir+'/'+'JsonTemplates') # exist_ok=True only python3\n",
    "elif(not os.path.exists(rundir)):\n",
    "  os.mkdir(rundir)\n",
    "  os.mkdir(rundir+'/'+'JsonTemplates') # exist_ok=True only python3\n",
    "\n",
    "if(os.path.isdir(cmipdir+'/'+'CMIP6')):\n",
    "  #print('hello')\n",
    "  os.symlink(cmipdir+'/'+'CMIP6',rundir+'/'+'CMIP6')\n",
    "else:\n",
    "  print('there')\n",
    "  os.mkdir(cmipdir)\n",
    "os.symlink(prodir+'/'+'TablesTemplates',rundir+'/'+'TablesTemplates')\n",
    "os.symlink(prodir+'/'+'cmip6-cmor-tables',rundir+'/'+'cmip6-cmor-tables')\n",
    "\n",
    "if(BATCH or BACKGROUND):\n",
    "  print('Copying and editing '+prodir+'/'+'qjob.csh to '+rundir+'/'+'qjob.csh')\n",
    "  ifh=open(prodir+'/'+'qjob.csh')\n",
    "  ofh=open(rundir+'/'+'qjob.csh','w')\n",
    "  for i,line in enumerate(ifh):\n",
    "    line=line.replace('CAFEPP_SCRIPT','./'+script+' RUNDIR')\n",
    "    line=line.replace('RUNDIR',rundir)\n",
    "    line=line.replace('CONDA_SOURCE','. /short/v14/mac599/anaconda3/etc/profile.d/conda.sh')\n",
    "    line=line.replace('CONDA_ACTIVATE','conda activate cafepp_27_scipy')\n",
    "    print(line,file=ofh,end='')\n",
    "  ifh.close()\n",
    "  ofh.close()\n",
    "\n",
    "vector_string=['decadal_diag.py','decadal_diag_py2.py','decadal_diag_py3.py','cafepp.py','app_funcs.py','ProcTime.py']\n",
    "for i,file_now in enumerate(vector_string):\n",
    "  print('Copying '+srcdir+'/'+file_now+' to ',rundir)\n",
    "  shutil.copyfile(srcdir+'/'+file_now,rundir+'/'+file_now)\n",
    "  \n",
    "vector_string=[script]\n",
    "for i,file_now in enumerate(vector_string):\n",
    "  print('Copying '+prodir+'/'+file_now+' to ',rundir)\n",
    "  shutil.copyfile(prodir+'/'+file_now,rundir+'/'+file_now)\n",
    "  os.chmod(rundir+'/'+script,500) #make executable incase want to run it interacively from terminal.\n",
    "  \n",
    "vector_string=[]\n",
    "#vector_string.append(script) #may need to do edits?\n",
    "#vector_string.append('cafepp_monthly_assimilation.json')\n",
    "vector_string.append('JsonTemplates'+'/'+'cafepp_experiments.json')\n",
    "vector_string.append('JsonTemplates'+'/'+'cafepp_csiro-gfdl.json')\n",
    "vector_string.append('JsonTemplates'+'/'+'cafepp_vars.json')\n",
    "#vector_string.append('cafepp_monthly_forecast.py')\n",
    "for i,file_now in enumerate(vector_string):\n",
    "  print('Copying '+prodir+'/'+file_now+' to ',rundir)\n",
    "  shutil.copyfile(prodir+'/'+file_now,rundir+'/'+file_now)\n",
    "  \n",
    "print('Copying '+prodir+'/'+'cafepp_monthly_assimilation.json',rundir+'/'+'JsonTemplates/cafepp_monthly_assimilation.json')\n",
    "shutil.copyfile(prodir+'/'+'cafepp_monthly_assimilation.json',rundir+'/'+'JsonTemplates/cafepp_monthly_assimilation.json')\n",
    "  \n",
    "if(not os.path.exists(rundir+'/'+'CMIP5')):\n",
    "  if(re.match('raijin',hostname)):\n",
    "    #print('No need to set CMIP5 symbolic link on raijin.')\n",
    "    os.symlink('/g/data/p66/mac599/CMIP5',rundir+'/'+'CMIP5') \n",
    "  elif(re.match('oa-3.-cdc',hostname)):\n",
    "    os.symlink('/OSM/CBR/OA_DCFP/data/CAFEPP/g/data/p66/mac599/CMIP5',rundir+'/'+'CMIP5')\n",
    "  elif(re.match('tube-hba',hostname)):\n",
    "    os.symlink('/OSM/HBA/OA_DECADAL_CLIMATE/work/col414/CMIP5',rundir+'/'+'CMIP5')\n",
    "  else:\n",
    "    raise SystemExit('Dont know how to set CMIP5 symbolic link on this machine.'+__file__+' line number: '+str(inspect.stack()[0][2]))\n",
    "\n",
    "#   srcdir='/home/599/mac599/decadal' #location of main cafepp code\n",
    "#   prodir=srcdir+'/paper_analysis' #project directory\n",
    " \n",
    "print('prodir=',prodir)\n",
    "print('rundir=',rundir)\n",
    "\n",
    "#break\n",
    "print('END')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 6,
        "hidden": false,
        "row": 27,
        "width": 12
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "source": [
    "# We establish where to execute the job.\n",
    "\n",
    "## We import a function that is relatively simple to loop over all necessary years, months, ensembles as required.\n",
    "## Different applications will require a different module to be written (often small and relatively simple)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 4,
        "hidden": false,
        "row": 16,
        "width": 4
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BEGIN\n",
      "Current Working Directory= /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_20\n",
      ".\n",
      "you have not specified the environment variable: 'CDAT_LOCATION' , trying to import cdms2 anyway\n",
      "END\n"
     ]
    }
   ],
   "source": [
    "print('BEGIN')\n",
    "\n",
    "os.chdir(rundir)\n",
    "print('Current Working Directory=',os.getcwd())\n",
    "\n",
    "os.environ['APP_OUTPATH'] = '.'\n",
    "print(os.getenv('APP_OUTPATH'))\n",
    "\n",
    "import getpass\n",
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "import os\n",
    "from time import strftime\n",
    "import netCDF4\n",
    "from math import radians, cos, sin, asin, sqrt\n",
    "import seawater\n",
    "#import sys\n",
    "import getopt\n",
    "import string\n",
    "from decadal_diag import MustHaveAllLevs,diag_acc_drake,diag_acc_africa,diag_mozmbq,diag_aabw,diag_nadw,\\\n",
    "diag_pp,diag_nflux,diag_ep,diag_ssh,diag_moc,diag_moc_atlantic,diag_moc_pacific,diag_moc_indian,\\\n",
    "diag_shice_cover,diag_nhice_cover,diag_nino34,xtra_nino34,init_data,sum_data,avg_data,filemonth_index,\\\n",
    "data_wavg,time_avg,diag_nhblocking_index,diag_rws5,finish,diag_msftyz,make_mask3D,diag_mfo,transPort,\\\n",
    "diag_rws500,create_odirs,create_ofils,diag_iod,diag_iod,xtra_iod,atmos_vertical_interpolate,diag_isothetaoNc,\\\n",
    "calc_iso_surface,calc_isoN,grab_var_meta,diag_psl,diag_hfls,diag_heat_content,diag_salt_content,\\\n",
    "diag_north_heat_trans,diag_north_salt_trans,ocean_vertical_interpolate,diag_thetao0to80m,diag_varNl,\\\n",
    "uncomment_json,process_json,modify_json,get_daily_indices_for_monthlyave,diag_maxdTbydz,diag_depmaxdTbydz,\\\n",
    "diag_dTbydz,shade_2d_simple,shade_2d_latlon,diag_zmld_boyer,zmld_boyer,sigmatheta,diag_zmld_so,\\\n",
    "zmld_so,diag_spice,spice,diag_bigthetao,diag_soabs,diag_spiciness,diag_potrho,fractional_year_from_num2date,\\\n",
    "new_monthly_array_shape,restrict_input_files,get_timestamp_number\n",
    "from decadal_diag import data_wavg_ProcTime,cmor_file_parts,cmor_directory_parts,cmor_ripf_parts,file_spec_summary,generate_daily_month_indices\n",
    "from decadal_diag import plot_map_box,plot_xy_climatology,plot_xy_ensemble,calculate_monthly_climatology_anomaly_from_monthly\n",
    "\n",
    "if(sys.version_info.major==2):\n",
    "  from decadal_diag_py2 import shade_2d_polar\n",
    "elif(sys.version_info.major==3):\n",
    "  from decadal_diag_py3 import shade_2d_polar\n",
    "else:\n",
    "  raise SystemExit('Not ready for this version of python:'+__file__+' line number: '+str(inspect.stack()[0][2]))\n",
    "    \n",
    "import cmor\n",
    "import cdtime\n",
    "from app_funcs import *\n",
    "import json\n",
    "import pprint\n",
    "from datetime import date\n",
    "import filecmp\n",
    "from shutil import copyfile\n",
    "import cdms2\n",
    "import inspect\n",
    "import socket\n",
    "import glob\n",
    "from matplotlib.mlab import griddata\n",
    "import scipy.sparse as sps\n",
    "import cartopy.crs as ccrs\n",
    "from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter\n",
    "from cartopy.util import add_cyclic_point\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['mathtext.default'] = 'regular'\n",
    "import matplotlib.pyplot as plt\n",
    "from gridfill import fill as poisson_fill\n",
    "\n",
    "#from ProcTime import ProcTime\n",
    "\n",
    "#del(cafepp_monthly_forecast)\n",
    "#if('cafepp_monthly_forecast' not in sys.modules):\n",
    "#  print('no')\n",
    "#else:\n",
    "#  print('yes')\n",
    "\n",
    "if(script=='cafepp_monthly_assimilation.py'):\n",
    "  import cafepp_monthly_assimilation #this script/function needs to be edited to limit variable.\n",
    "elif(script=='cafepp_monthly_forecast.py'):\n",
    "  import cafepp_monthly_forecast #this script/function needs to be edited to limit year,month,ensemble range,variable if desired. Will need to restart kernel for it to take effect.\n",
    "elif(script=='cafepp_monthly_control.py'):\n",
    "  import cafepp_monthly_control\n",
    "\n",
    "#print(sys.modules)\n",
    "\n",
    "print('END')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 4,
        "height": 4,
        "hidden": false,
        "row": 16,
        "width": 4
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "source": [
    "## information common to processing all runs..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 4,
        "hidden": false,
        "row": 20,
        "width": 4
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BEGIN.\n",
      "indices_nino= ['nino34', 'nino3', 'nino4', 'nino1+2']\n",
      "indices_i_gr2= [[190, 239], [210, 269], [150, 209], [270, 280]]\n",
      "END.\n"
     ]
    }
   ],
   "source": [
    "print('BEGIN.')\n",
    "\n",
    "import math\n",
    "\n",
    "#raise Exception('STOP!')\n",
    "\n",
    "#del(get_timestamp_number)\n",
    "#from decadal_diag import get_timestamp_number\n",
    "\n",
    "nmy=12\n",
    "\n",
    "rad = 4.0*math.atan(1.0)/180.0\n",
    "#print('rad=',rad)\n",
    "\n",
    "lons_nino = [[190, 240, 240, 190, 190], [210, 270, 270, 210, 210], [160, 210, 210, 160, 160], [270, 280, 280, 270, 270]] # indices_nino=['nino34','nino3','nino4','nino1+2']\n",
    "lats_nino = [[-5, -5, 5, 5, -5], [-5.1, -5.1, 5.1, 5.1, -5.1], [-5, -5, 5, 5, -5], [-10, -10, 0, 0, -10]] # 170W-120W:190-240, 150W-90W:210-270, 160E-150W:160-210, 90W-80W:270-280 (10S-0 latitude, others all 5S to 5N)\n",
    "indices_nino=['nino34','nino3','nino4','nino1+2'] # 170W-120W:190-240, 150W-90W:210-270, 160E-150W:160-210, 90W-80W:270-280 (10S-0 latitude, others all 5S to 5N)\n",
    "indices_label=['Ni$\\~{n}$o3.4','Ni$\\~{n}$o3','Ni$\\~{n}$o4','Ni$\\~{n}$o1+2']\n",
    "  \n",
    "#if(cafe_grid=='gr2'):\n",
    "#elif(cafe_grid=='gn'):\n",
    "\n",
    "indices_i_gr2,indices_j_gr2=[[190,239],[210,269],[150,209],[270,280]],[[85,94],[85,94],[85,94],[80,89]] #check this, whether I need +1 also...what about fractional cells?\n",
    "indices_i_gn,indices_j_gn=[[110,159],[130,189],[80,129],[190,199]],[[122,151],[122,151],[122,151],[107,136]] #check this, whether I need +1 also...what about fractional cells?\n",
    "indices_i_ncep2,indices_j_ncep2=[[102,128],[112,144],[86,112],[144,149]], [[45,50],[45,50],[45,50],[42,46]]\n",
    "\n",
    "#indices_select=array('i',[0,1,2,3])\n",
    "#indices_select=[np.int(0)]\n",
    "#indices_select=range(3)\n",
    "#print(type(indices_select))\n",
    "#print('indices_select=',indices_select)00\n",
    "#used for speeding things up, either can choose a single indice or all - don't know why can't use : in eval for e.g.\n",
    "sss='-1' #last\n",
    "sss='0' #first\n",
    "sss='ALL'\n",
    "\n",
    "#raise Exception('STOP!')\n",
    "\n",
    "if(sss!='ALL'):\n",
    "  indices_nino=[indices_nino[eval(sss)]]\n",
    "  indices_label=[indices_label[eval(sss)]]\n",
    "  indices_i_gr2=[indices_i_gr2[eval(sss)]]\n",
    "  indices_j_gr2=[indices_j_gr2[eval(sss)]]\n",
    "  indices_i_gn=[indices_i_gn[eval(sss)]]\n",
    "  indices_j_gn=[indices_j_gn[eval(sss)]]\n",
    "  indices_i_ncep2=[indices_i_ncep2[eval(sss)]]\n",
    "  indices_j_ncep2=[indices_j_ncep2[eval(sss)]]\n",
    "print('indices_nino=',indices_nino)\n",
    "\n",
    "print('indices_i_gr2=',indices_i_gr2)\n",
    "\n",
    "nindices_nino=len(indices_nino)\n",
    "\n",
    "#else:\n",
    "#  raise Exception('STOP!')\n",
    "#raise Exception('STOP!')\n",
    "#del(get_timestamp_number)\n",
    "#from decadal_diag import get_timestamp_number\n",
    "\n",
    "rad = 4.0*math.atan(1.0)/180.0\n",
    "#print('rad=',rad)\n",
    "\n",
    "print('END.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 4,
        "hidden": false,
        "row": 33,
        "width": 12
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "source": [
    "# 1. Here we are procesing an assimilation, set script='cafepp_monthly_assimilation.py' above.\n",
    "\n",
    "## Either submit it to the queue or run interactively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 4,
        "height": 4,
        "hidden": false,
        "row": 20,
        "width": 4
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "outputs": [],
   "source": [
    "print('BEGIN')\n",
    "\n",
    "SWITCH_OFF=True #temporary for easy skipping of this section.\n",
    "SWITCH_OFF=False\n",
    "\n",
    "if(not SWITCH_OFF):\n",
    "  print('prodir=',prodir)\n",
    "  \n",
    "  #dvar='t16d'\n",
    "  #dvar='t20d'\n",
    "  dvar='t22d'\n",
    "  \n",
    "  ybeg=2002 #2002\n",
    "  yend=2016 #2016\n",
    "  mbeg=1 #1\n",
    "  mend=12 #6\n",
    "  \n",
    "  NoClobber=True #do not clobber\n",
    "  #NoClobber=False #clobber\n",
    "  \n",
    "  cafe_experiment='coupled_da/OUTPUT-2step-nobreeding-carbon2'\n",
    "  \n",
    "  line_kwargs='rundir='+rundir+ \\\n",
    "    ' dvar='+dvar+ \\\n",
    "    ' ybeg='+str(ybeg)+ \\\n",
    "    ' yend='+str(yend)+ \\\n",
    "    ' mbeg='+str(mbeg)+ \\\n",
    "    ' mend='+str(mend)+ \\\n",
    "    ' NoClobber='+str(NoClobber)+ \\\n",
    "    ' cafe_experiment='+str(cafe_experiment)\n",
    "\n",
    "  if(script=='cafepp_monthly_assimilation.py'):\n",
    "    \n",
    "    if(BATCH or BACKGROUND):\n",
    "      os.chmod(script,500)\n",
    "      \n",
    "      if(BACKGROUND):\n",
    "        os.chmod('qjob.csh',500)\n",
    "      \n",
    "      print('Copying and editing '+prodir+'/'+'qjob.csh to '+rundir+'/'+'qjob.csh')\n",
    "      ifh=open(prodir+'/'+'qjob.csh')\n",
    "      ofh=open(rundir+'/'+'qjob.csh','w')\n",
    "      for i,line in enumerate(ifh):\n",
    "        line=line.replace('CAFEPP_SCRIPT','./'+script+' '+line_kwargs)\n",
    "        line=line.replace('RUNDIR',rundir)\n",
    "        line=line.replace('CONDA_SOURCE','. /OSM/CBR/OA_DCFP/apps/col414/anaconda3/etc/profile.d/conda.sh')\n",
    "        line=line.replace('CONDA_ACTIVATE','conda activate cafepp_27_scipy')\n",
    "        print(line,file=ofh,end='')\n",
    "      ifh.close()\n",
    "      ofh.close()\n",
    "      \n",
    "      #raise Exception('STOP!')\n",
    "      \n",
    "      if(BATCH):\n",
    "        os.system('qsub '+rundir+'/'+'qjob.csh')\n",
    "      else:\n",
    "        print('Current Working Directory=',os.getcwd())\n",
    "        \n",
    "        import commands\n",
    "        ret = commands.getoutput('nohup /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_16/qjob.csh > /dev/null 2>&1 &')\n",
    "        print(ret)\n",
    "\n",
    "        raise Exception('STOP!')\n",
    "        \n",
    "    else:\n",
    "      print('Current Working Directory=',os.getcwd())\n",
    "      os.chdir(rundir)\n",
    "      \n",
    "      kwargs=dict(x.split('=', 1) for x in line_kwargs.split(' '))\n",
    "      cafepp_monthly_assimilation.main(**kwargs)\n",
    "  \n",
    "print('END')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 5,
        "hidden": false,
        "row": 37,
        "width": 12
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "source": [
    "\n",
    "# 2. Here we are processing multi-ensemble forecast, set script='cafepp_monthly_forecast.py' above.\n",
    "## Either submit it to the queue or run interactively.\n",
    "\n",
    "queued version needs more work...\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "code_folding": [],
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 4,
        "hidden": false,
        "row": 42,
        "width": 4
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BEGIN\n",
      "Current Working Directory= /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_20\n",
      "kwargs= {'StopAfterOne': 'False', 'dvar': 't20d', 'ebeg': '1', 'ybeg_first': '2002', 'mbeg_first': '2', 'yend': '2016', 'mbeg_norm': '1', 'eend': '11', 'num_months_truncate': '60', 'mend_norm': '12', 'ybeg': '2015', 'rundir': '/OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_20', 'NoClobber': 'True', 'cafe_experiment': 'v1_forecast', 'yend_last': '2016', 'mend_last': '5'}\n"
     ]
    },
    {
     "ename": "Exception",
     "evalue": "STOP!",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m\u001b[0m",
      "\u001b[0;31mException\u001b[0mTraceback (most recent call last)",
      "\u001b[0;32mjupyter_notebook\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     96\u001b[0m       \u001b[0mkwargs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msplit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'='\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mx\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mline_kwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msplit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m' '\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     97\u001b[0m       \u001b[0;32mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'kwargs='\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 98\u001b[0;31m       \u001b[0;32mraise\u001b[0m \u001b[0mException\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'STOP!'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     99\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    100\u001b[0m       \u001b[0mcafepp_monthly_forecast\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmain\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mException\u001b[0m: STOP!"
     ]
    }
   ],
   "source": [
    "print('BEGIN')\n",
    "\n",
    "#from importlib import reload\n",
    "#import importlib\n",
    "#print(dir(importlib))\n",
    "#cafepp_monthly_forecast=importlib.import_module(cafepp_monthly_forecast)\n",
    "#raise Exception('STOP!')\n",
    "\n",
    "SWITCH_OFF=True #temporary for easy skipping of this section.\n",
    "SWITCH_OFF=False\n",
    "\n",
    "if(not SWITCH_OFF):\n",
    "  #print('prodir=',prodir)\n",
    "  \n",
    "#   dvar='tos'\n",
    "#   dvar='msftyz'\n",
    "#   dvar='siconc'\n",
    "#   dvar='ua'\n",
    "#   dvar='ta'\n",
    "#   dvar='zg'\n",
    "#   dvar='t16d'\n",
    "  dvar='t20d'\n",
    "#   dvar='t22d'\n",
    "#   dvar='thetao'\n",
    "#   dvar='tauu'\n",
    "\n",
    "  StopAfterOne='False'\n",
    "  #StopAfterOne='True'\n",
    "\n",
    "  NoClobber=True #do not clobber\n",
    "  #NoClobber=False #clobber\n",
    "\n",
    "  cafe_experiment='v1_forecast'\n",
    "  #cafe_experiment='v0_forecast'\n",
    "  \n",
    "  ###\n",
    "  \n",
    "  ybeg=2015 #2002\n",
    "  yend=2016 #2016\n",
    "  mbeg_norm=1 #1\n",
    "  mend_norm=12 #12\n",
    "  ybeg_first=2002\n",
    "  yend_last=2016\n",
    "  mbeg_first=2\n",
    "  mend_last=5\n",
    "  ebeg=1 #1\n",
    "  eend=11 #11\n",
    "  num_months_truncate=60\n",
    "  #num_months_truncate=72\n",
    "  #num_months_truncate=600\n",
    "\n",
    "  line_kwargs='StopAfterOne='+StopAfterOne+' num_months_truncate='+str(num_months_truncate)+' NoClobber='+str(NoClobber)+' rundir='+rundir+ \\\n",
    "    ' dvar='+dvar+' cafe_experiment='+cafe_experiment+ \\\n",
    "    ' ybeg='+str(ybeg)+' yend='+str(yend)+ \\\n",
    "    ' mbeg_norm='+str(mbeg_norm)+' mend_norm='+str(mend_norm)+ \\\n",
    "    ' ybeg_first='+str(ybeg_first)+' yend_last='+str(yend_last)+ \\\n",
    "    ' mbeg_first='+str(mbeg_first)+' mend_last='+str(mend_last)+ \\\n",
    "    ' ebeg='+str(ebeg)+' eend='+str(eend)  \n",
    "  \n",
    "  if(script=='cafepp_monthly_forecast.py'):\n",
    "\n",
    "    if(BATCH or BACKGROUND):\n",
    "      os.chmod(script,500)\n",
    "      \n",
    "      if(BACKGROUND):\n",
    "        os.chmod('qjob.csh',500)\n",
    "      \n",
    "      print('Copying and editing '+prodir+'/'+'qjob.csh to '+rundir+'/'+'qjob.csh')\n",
    "      ifh=open(prodir+'/'+'qjob.csh')\n",
    "      ofh=open(rundir+'/'+'qjob.csh','w')\n",
    "      for i,line in enumerate(ifh):\n",
    "        line=line.replace('CAFEPP_SCRIPT','./'+script+' '+line_kwargs)\n",
    "        line=line.replace('RUNDIR',rundir)\n",
    "        line=line.replace('CONDA_SOURCE','. /OSM/CBR/OA_DCFP/apps/col414/anaconda3/etc/profile.d/conda.sh')\n",
    "        line=line.replace('CONDA_ACTIVATE','conda activate cafepp_27_scipy')\n",
    "        print(line,file=ofh,end='')\n",
    "      ifh.close()\n",
    "      ofh.close()\n",
    "      \n",
    "      #raise Exception('STOP!')\n",
    "      \n",
    "      if(BATCH):\n",
    "        os.system('qsub '+rundir+'/'+'qjob.csh')\n",
    "      else:\n",
    "        print('Current Working Directory=',os.getcwd())\n",
    "        \n",
    "#         import commands\n",
    "#         ret = commands.getoutput('nohup /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_16/qjob.csh > /dev/null 2>&1 &')\n",
    "#         print(ret)\n",
    "\n",
    "#         import subprocess\n",
    "#         j=subprocess.call(['/OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_16/qjob.csh'], shell=True)\n",
    "#         print(j)\n",
    "\n",
    "        #os.system(''./qjob.csh')\n",
    "        raise Exception('STOP!')\n",
    "        \n",
    "    else:\n",
    "      print('Current Working Directory=',os.getcwd())\n",
    "      os.chdir(rundir)\n",
    "      \n",
    "      kwargs=dict(x.split('=', 1) for x in line_kwargs.split(' '))\n",
    "      print('kwargs=',kwargs)\n",
    "      raise Exception('STOP!')\n",
    "      \n",
    "      cafepp_monthly_forecast.main(**kwargs)\n",
    "  \n",
    "#      cafepp_monthly_forecast.main(rundir, dvar=dvar, cafe_experiment='v1_forecast', \\\n",
    "#        ybeg=2002, yend=2016, mbeg=1, mend=12, \\\n",
    "#        ybeg_first=2002, yend_last=2016, mbeg_first=2, mend_last=5, \\\n",
    "#        ebeg=1, eend=11)\n",
    "    \n",
    "#      cafepp_monthly_forecast.main(rundir, dvar=dvar, cafe_experiment=cafe_experiment, \\\n",
    "#        ybeg=2015, yend=2015, mbeg=1, mend=1, \\\n",
    "#        ybeg_first=2002, yend_last=2016, mbeg_first=2, mend_last=5, \\\n",
    "#        ebeg=1, eend=11)  \n",
    "  \n",
    "print('END')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print('BEGIN.')\n",
    "\n",
    "SWITCH_OFF=True #temporary for easy skipping of this section.\n",
    "#SWITCH_OFF=False\n",
    "\n",
    "if(not SWITCH_OFF):\n",
    "  \n",
    "  import math\n",
    "  os.chdir(rundir)\n",
    "\n",
    "  colors=['black','red','green','blue','orange','brown','purple','pink']\n",
    "  \n",
    "  #forc_nens=11\n",
    "  grid_label='gr2'\n",
    "  grid_label='gn'\n",
    "  forc_files_stringA='CMIP6/CMIP/CSIRO/CAFE-1-0/historical/r[1-9]i1p1f1/Omon/'+dvar+'/'+grid_label+'/v20171025/'+dvar+'_Omon_historical_CAFE-1-0_r[1-9]i1p1f1_'+grid_label+'_2011??-??????.nc'\n",
    "  forc_files_stringB='CMIP6/CMIP/CSIRO/CAFE-1-0/historical/r[1-9]i1p1f1/Omon/'+dvar+'/'+grid_label+'/v20171025/'+dvar+'_Omon_historical_CAFE-1-0_r[1-9]i1p1f1_'+grid_label+'_2012??-??????.nc'\n",
    "  forc_files_stringC='CMIP6/CMIP/CSIRO/CAFE-1-0/historical/r[1-9]i1p1f1/Omon/'+dvar+'/'+grid_label+'/v20171025/'+dvar+'_Omon_historical_CAFE-1-0_r[1-9]i1p1f1_'+grid_label+'_2013??-??????.nc'\n",
    "  forc_files_stringD='CMIP6/CMIP/CSIRO/CAFE-1-0/historical/r[1-9]i1p1f1/Omon/'+dvar+'/'+grid_label+'/v20171025/'+dvar+'_Omon_historical_CAFE-1-0_r[1-9]i1p1f1_'+grid_label+'_2014??-??????.nc'\n",
    "\n",
    "  forc_files_stringA=''\n",
    "  forc_files_stringD=''\n",
    "\n",
    "  forc_files=sorted(restrict_input_files(glob.glob(forc_files_stringA)+glob.glob(forc_files_stringB)+glob.glob(forc_files_stringC)+glob.glob(forc_files_stringD),24,72)) #later want to split these into different ensembles. Will assume same input files for all valid ensembles.\n",
    "\n",
    "  nforc_files=len(forc_files)\n",
    "  print('forc_files[0]=',forc_files[0])\n",
    "  #print('forc input_files ('+str(nforc_files)+') =',forc_files)\n",
    "\n",
    "  datetime_all,datetime_uniq,ripf_all,ripf_uniq=file_spec_summary(forc_files,True)\n",
    "\n",
    "  #raise Exception('STOP!')\n",
    "print('END.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 4,
        "hidden": false,
        "row": 46,
        "width": 12
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "source": [
    "# 3. Here we are processing uninitialised control experiment, set script='cafepp_monthly_control.py' above\n",
    "\n",
    "Either submit it to the queue or run interactive."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# turn off scrolling on output cells..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%%javascript\n",
    "#IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
    "#    return false;\n",
    "#}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 4,
        "height": 4,
        "hidden": false,
        "row": 42,
        "width": 4
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "outputs": [],
   "source": [
    "print('BEGIN')\n",
    "\n",
    "# %load_ext autoreload\n",
    "# %autoreload 2\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "SWITCH_OFF=True #temporary for easy skipping of this section.\n",
    "SWITCH_OFF=False\n",
    "\n",
    "#cafe_experiment='v0'\n",
    "#cafe_experiment='v1'\n",
    "cafe_experiment='v2'\n",
    "# cafe_experiment='v3'\n",
    "#cafe_experiment='nov17n' #ybeg,yend=200,239\n",
    "#cafe_experiment='jul18b' #ybeg,yend=1980,2009\n",
    "  \n",
    "#cafe_experiment='v0' #t_surf in atmos_month_*nc files.\n",
    "ybeg=491\n",
    "yend=500\n",
    "mbeg=1\n",
    "mend=12\n",
    "\n",
    "###\n",
    "\n",
    "#dvar='so'\n",
    "dvar='thetao'\n",
    "#dvar='tos'\n",
    "#dvar='tas'\n",
    "#dvar='siconc'\n",
    "#dvar='msftyz' #was called msftyyz\n",
    "#dvar='ua'\n",
    "#dvar='ta'\n",
    "#dvar='zg'\n",
    "#dvar='t16d'\n",
    "#dvar='t20d'\n",
    "#dvar='t22d'\n",
    "\n",
    "NoClobber=True #do not clobber\n",
    "NoClobber=False #clobber\n",
    "\n",
    "line_kwargs= \\\n",
    "  'NoClobber='+str(NoClobber)+ \\\n",
    "  ' rundir='+rundir+ \\\n",
    "  ' dvar='+dvar+ \\\n",
    "  ' cafe_experiment='+cafe_experiment+ \\\n",
    "  ' ybeg='+str(ybeg)+ \\\n",
    "  ' yend='+str(yend)+ \\\n",
    "  ' mbeg='+str(mbeg)+ \\\n",
    "  ' mend='+str(mend)\n",
    "\n",
    "if(not SWITCH_OFF):\n",
    "  \n",
    "  print('prodir=',prodir)\n",
    "  \n",
    "  if(script=='cafepp_monthly_control.py'):\n",
    "    if(BATCH or BACKGROUND):\n",
    "      os.chmod(script,500)\n",
    "      \n",
    "      if(BACKGROUND):\n",
    "        os.chmod('qjob.csh',500)\n",
    "    \n",
    "      print('Copying and editing '+prodir+'/'+'qjob.csh to '+rundir+'/'+'qjob.csh')\n",
    "      ifh=open(prodir+'/'+'qjob.csh')\n",
    "      ofh=open(rundir+'/'+'qjob.csh','w')\n",
    "      for i,line in enumerate(ifh):\n",
    "        line=line.replace('CAFEPP_SCRIPT','./'+script+' '+line_kwargs)\n",
    "        line=line.replace('RUNDIR',rundir)\n",
    "#         line=line.replace('CONDA_SOURCE','. /short/v14/mac599/anaconda3/etc/profile.d/conda.sh') #raijin\n",
    "        line=line.replace('CONDA_SOURCE','. /OSM/CBR/OA_DCFP/apps/col414/anaconda3/etc/profile.d/conda.sh') #vm's\n",
    "        line=line.replace('CONDA_ACTIVATE','conda activate cafepp_27_scipy')\n",
    "        print(line,file=ofh,end='')\n",
    "      ifh.close()\n",
    "      ofh.close()\n",
    "  \n",
    "      if(BATCH):\n",
    "        #os.system('qsub '+rundir+'/'+'qjob.csh')\n",
    "        pass\n",
    "      else:\n",
    "        print('Current Working Directory=',os.getcwd())\n",
    "        \n",
    "#         import commands\n",
    "#         ret = commands.getoutput('nohup /OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_16/qjob.csh > /dev/null 2>&1 &')\n",
    "#         print(ret)\n",
    "\n",
    "#         import subprocess\n",
    "#         j=subprocess.call(['/OSM/CBR/OA_DCFP/data/col414/cafepp/rundir_16/qjob.csh'], shell=True)\n",
    "#         print(j)\n",
    "\n",
    "        #os.system(''./qjob.csh')\n",
    "  \n",
    "        raise Exception('STOP!')        \n",
    "\n",
    "    else:\n",
    "      print('Current Working Directory=',os.getcwd())\n",
    "      os.chdir(rundir)\n",
    "\n",
    "      kwargs=dict(x.split('=', 1) for x in line_kwargs.split(' '))\n",
    "      cafepp_monthly_control.main(**kwargs)\n",
    "  \n",
    "print('END')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "raise Exception('STOP!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 8,
        "height": 4,
        "hidden": false,
        "row": 42,
        "width": 4
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "source": [
    "# Process/Plot control run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 24,
        "hidden": false,
        "row": 50,
        "width": 4
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "outputs": [],
   "source": [
    "print('BEGIN.')\n",
    "\n",
    "pval=5\n",
    "#pval=1\n",
    "\n",
    "run='historial'\n",
    "run='piControl'\n",
    "\n",
    "cont_files=['CMIP6/CMIP/CSIRO/CAFE-1-0/'+run+'/r1i1p'+str(pval)+'f1/Omon/tos/gn/v20171025/tos_Omon_'+run+'_CAFE-1-0_r1i1p'+str(pval)+'f1_gn_009101-010012.nc', \\\n",
    "            'CMIP6/CMIP/CSIRO/CAFE-1-0/'+run+'/r1i1p'+str(pval)+'f1/Omon/tos/gn/v20171025/tos_Omon_'+run+'_CAFE-1-0_r1i1p'+str(pval)+'f1_gn_019101-020012.nc', \\\n",
    "            'CMIP6/CMIP/CSIRO/CAFE-1-0/'+run+'/r1i1p'+str(pval)+'f1/Omon/tos/gn/v20171025/tos_Omon_'+run+'_CAFE-1-0_r1i1p'+str(pval)+'f1_gn_029101-030012.nc', \\\n",
    "            'CMIP6/CMIP/CSIRO/CAFE-1-0/'+run+'/r1i1p'+str(pval)+'f1/Omon/tos/gn/v20171025/tos_Omon_'+run+'_CAFE-1-0_r1i1p'+str(pval)+'f1_gn_039101-040012.nc', \\\n",
    "            'CMIP6/CMIP/CSIRO/CAFE-1-0/'+run+'/r1i1p'+str(pval)+'f1/Omon/tos/gn/v20171025/tos_Omon_'+run+'_CAFE-1-0_r1i1p'+str(pval)+'f1_gn_049101-050012.nc']\n",
    "\n",
    "#cont_files=cont_files[0:2]\n",
    "\n",
    "cont_nexps=len(cont_files)\n",
    "\n",
    "for cont_cnt,cont_file in enumerate(cont_files):\n",
    "  \n",
    "  print('cont input file=',cont_file)\n",
    "\n",
    "  cont_ifh=netCDF4.Dataset(cont_file)\n",
    "\n",
    "  cont_lat=cont_ifh.variables['latitude'][:,0]\n",
    "  cont_lon=cont_ifh.variables['longitude'][0,:]\n",
    "\n",
    "  #forc_clat=np.zeros((nforc_files),dtype=float) #if latitude range changes will need to have custom clat.\n",
    "  cont_clat=np.cos(cont_lat[:]*rad)\n",
    "\n",
    "  cont_time=cont_ifh.variables['time']\n",
    "\n",
    "  cont_date_time_stamp=netCDF4.num2date(cont_time[:],cont_time.units,cont_time.calendar)\n",
    "\n",
    "  cont_ntime=len(cont_time)\n",
    "\n",
    "  '''\n",
    "  print('date_time_stamp=',date_time_stamp)\n",
    "  num_stamp=netCDF4.date2num(date_time_stamp,time.units,time.calendar)\n",
    "  print('num_stamp=',num_stamp)\n",
    "  print('year_fraction=',year_fraction)\n",
    "  '''\n",
    "  \n",
    "  if(cont_cnt==0):\n",
    "    cont_year_fraction_monthly=ma.zeros((cont_ntime,cont_nexps),dtype=float)\n",
    "    cont_nino_monthly=ma.zeros((cont_ntime,nindices_nino,cont_nexps),dtype=float)    \n",
    "  cont_year_fraction_monthly[:,cont_cnt]=fractional_year_from_num2date(cont_date_time_stamp,cont_time.calendar)\n",
    "\n",
    "  #print('type(cont_date_time_stamp)=',type(cont_date_time_stamp))\n",
    "  #print('cont_date_time_stamp=',cont_date_time_stamp)\n",
    "\n",
    "  #raise Exception('STOP!')\n",
    "  for cont_k,indice in enumerate(indices_nino):\n",
    "    #cont_k=indices_nino.index(indice)\n",
    "    cont_imin,cont_imax=indices_i_gn[indices_nino.index(indice)][0],indices_i_gn[indices_nino.index(indice)][1]\n",
    "    cont_jmin,cont_jmax=indices_j_gn[indices_nino.index(indice)][0],indices_j_gn[indices_nino.index(indice)][1]\n",
    "    #print('indices.index(indice),cont_imin,cont_imax,cont_jmin,cont_jmax=',indices.index(indice),cont_imin,cont_imax,cont_jmin,cont_jmax)\n",
    "    #print('cont_lat,cont_lon=',cont_lat[cont_jmin:cont_jmax],cont_lon[cont_imin:cont_imax])\n",
    "\n",
    "    cont_nino_monthly[:,cont_k,cont_cnt]=np.average(np.average(cont_ifh.variables[dvar][:,cont_jmin:cont_jmax,cont_imin:cont_imax],axis=1,weights=cont_clat[cont_jmin:cont_jmax]),axis=1)\n",
    "  cont_ifh.close()\n",
    "  #indices_i_gn,indices_j_gn=[[122,152],[122,152],[122,152]],[[110,160],[110,160],[110,160]] #check this, whether I need +1 also...\n",
    "  #raise Exception('STOP!')\n",
    "  #print('cont_year_fraction=',cont_year_fraction_monthly)\n",
    "\n",
    "  %matplotlib inline\n",
    "\n",
    "  fix,ax=plt.subplots()\n",
    "  for cont_k,indice in enumerate(indices_nino):\n",
    "    ax.plot(cont_year_fraction_monthly[:,cont_cnt],cont_nino_monthly[:,indices_nino.index(indice),cont_cnt],label=indices_label[indices_nino.index(indices_nino[cont_k])])\n",
    "  legend=ax.legend(loc='lower right',shadow=False,fontsize='x-large')\n",
    "  plt.title('Monthly values (CONTROL run)')\n",
    "  plt.xlabel('Year')\n",
    "  plt.ylabel('$^o$C')\n",
    "  #plt.show()\n",
    "\n",
    "#raise Exception('STOP!')\n",
    "\n",
    "  cont_ybeg=cont_date_time_stamp[0].year\n",
    "  cont_yend=cont_date_time_stamp[-1].year\n",
    "\n",
    "  cont_ydiff=cont_yend-cont_ybeg+1\n",
    "  nmy=12\n",
    "\n",
    "  cont_MissingMonths=False\n",
    "\n",
    "  cont_first_month=cont_date_time_stamp[0].month\n",
    "  cont_last_month=cont_date_time_stamp[-1].month\n",
    "\n",
    "  cont_missing_months_beg,cont_missing_months_end=0,0\n",
    "\n",
    "  cont_cyear_beg_skip,cont_cyear_end_skip=0,1\n",
    "\n",
    "  if(cont_first_month!=1):\n",
    "    cont_missing_months_beg=12-cont_first_month\n",
    "    cont_cyear_beg_skip=1\n",
    "    cont_MissingMonths=True\n",
    "\n",
    "  if(cont_last_month!=12):\n",
    "    cont_missing_months_end=12-cont_last_month\n",
    "    cont_cyear_end_skip=2\n",
    "    cont_MissingMonths=True\n",
    "  \n",
    "  if(cont_MissingMonths):\n",
    "    print('There are missing months in the set. '+str(cont_missing_months_beg)+' at beginning and '+str(cont_missing_months_end)+' at end.')\n",
    "    print('Currently years with missing months are not used in generating long term monthly climatology.')\n",
    "    print('And missing months will be set to missing in the final time-series.')\n",
    "    cont_ts_beg,cont_ts_end,cont_ts_avg,cont_dt_beg,cont_dt_end,cont_dt_avg=get_timestamp_number(cont_ybeg,cont_yend,1,12,cont_time.units,cont_time.calendar)\n",
    "    cont_year_fraction_monthly_full=cont_fractional_year_from_num2date(cont_ts_avg,cont_time.calendar)\n",
    "  \n",
    "    cont_nino_monthly_full=ma.masked_all(nindices_nino,cont_ydiff*nmy,dtype=float) #ensure missing months are masked out.\n",
    "    cont_last_month_index=cont_ydiff*nmy-cont_last_month\n",
    "    cont_nino_monthly_full[cont_first_month-1:cont_last_month_index,:]=cont_nino_monthly\n",
    "  else:\n",
    "    print('All years have 12 months.')\n",
    "  \n",
    "    cont_nino_monthly_full=cont_nino_monthly\n",
    "    cont_year_fraction_monthly_full=cont_year_fraction_monthly\n",
    "    \n",
    "#print('cont_nino_monthly_full.shape=',cont_nino_monthly_full.shape)\n",
    "\n",
    "cont_nino_climatology,cont_nino_monthly_anomaly = calculate_monthly_climatology_anomaly_from_monthly(cont_nino_monthly_full)\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "for plot_this_exp,cont_file in enumerate(cont_files):\n",
    "  cont_zero=np.zeros(len(cont_year_fraction_monthly_full[:,plot_this_exp]))\n",
    "  fix,ax=plt.subplots()\n",
    "  for cont_k,indice in enumerate(indices_nino):\n",
    "    ax.plot(cont_year_fraction_monthly_full[:,plot_this_exp],cont_nino_monthly_anomaly[:,cont_k,plot_this_exp],label=indices_label[indices_nino.index(indices_nino[cont_k])])\n",
    "  plt.plot(cont_year_fraction_monthly_full[:,plot_this_exp],cont_zero,color='black')\n",
    "\n",
    "  legend=ax.legend(loc='lower right',shadow=False,fontsize='x-large')\n",
    "  plt.title('Monthly anomalies (CONTROL run)')\n",
    "  plt.xlabel('Year')\n",
    "  plt.ylabel('$^o$C')\n",
    "  plt.show()\n",
    "\n",
    "print('END.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 7,
        "hidden": false,
        "row": 74,
        "width": 12
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "source": [
    "# Now to plot data, this will depend on success of previous steps producing necessary outputs with cafepp. Note that we have data on the native grid, where we can discover the points that are in the nino3.4 region.\n",
    "\n",
    "# Process/Plot forecast run, each forecast is a seperate run.\n",
    "\n",
    "Find minimum and maximum year/month from the set of input files\n",
    "Create a dummy time stamp that spans the entire set\n",
    "Read in data into a multi dimensional array (indice, ensemble, forecast, time)\n",
    "Have a set of indices which describe the beginning and ending of each forecast, for quick retrieval\n",
    "\n",
    "Would like to have anomalies relative to various climatologies, these could include different year periods, experiments to form climatology or even zero climatology."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 4,
        "height": 12,
        "hidden": false,
        "row": 50,
        "width": 4
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print('BEGIN.')\n",
    "\n",
    "import math\n",
    "os.chdir(rundir)\n",
    "\n",
    "colors=['black','red','green','blue','orange','brown','purple','pink']\n",
    "  \n",
    "#forc_nens=11\n",
    "grid_label='gr2'\n",
    "grid_label='gn'\n",
    "forc_files_stringA='CMIP6/CMIP/CSIRO/CAFE-1-0/historical/r[1-3]i1p1f1/Omon/'+dvar+'/'+grid_label+'/v20171025/'+dvar+'_Omon_historical_CAFE-1-0_r[1-3]i1p1f1_'+grid_label+'_2011??-??????.nc'\n",
    "forc_files_stringB='CMIP6/CMIP/CSIRO/CAFE-1-0/historical/r[1-3]i1p1f1/Omon/'+dvar+'/'+grid_label+'/v20171025/'+dvar+'_Omon_historical_CAFE-1-0_r[1-3]i1p1f1_'+grid_label+'_2012??-??????.nc'\n",
    "forc_files_stringC='CMIP6/CMIP/CSIRO/CAFE-1-0/historical/r[1-3]i1p1f1/Omon/'+dvar+'/'+grid_label+'/v20171025/'+dvar+'_Omon_historical_CAFE-1-0_r[1-3]i1p1f1_'+grid_label+'_2013??-??????.nc'\n",
    "forc_files_stringD='CMIP6/CMIP/CSIRO/CAFE-1-0/historical/r[1-3]i1p1f1/Omon/'+dvar+'/'+grid_label+'/v20171025/'+dvar+'_Omon_historical_CAFE-1-0_r[1-3]i1p1f1_'+grid_label+'_2014??-??????.nc'\n",
    "\n",
    "forc_files_stringA=''\n",
    "forc_files_stringD=''\n",
    "\n",
    "forc_files=sorted(restrict_input_files(glob.glob(forc_files_stringA)+glob.glob(forc_files_stringB)+glob.glob(forc_files_stringC)+glob.glob(forc_files_stringD),24,72)) #later want to split these into different ensembles. Will assume same input files for all valid ensembles.\n",
    "\n",
    "nforc_files=len(forc_files)\n",
    "print('forc_files[0]=',forc_files[0])\n",
    "#print('forc input_files ('+str(nforc_files)+') =',forc_files)\n",
    "\n",
    "datetime_all,datetime_uniq,ripf_all,ripf_uniq=file_spec_summary(forc_files,False)\n",
    "\n",
    "#print('len(datetime_uniq)=',len(datetime_uniq))\n",
    "#print('len(datetime_all)=',len(datetime_all))\n",
    "#print('len(ripf_uniq)=',len(ripf_uniq))\n",
    "#print('len(ripf_all)=',len(ripf_all))\n",
    "\n",
    "#print(ripf_all)\n",
    "\n",
    "#if(len(datetime_all)!=(len(datetime_uniq)*len(ripf_uniq))):\n",
    "#  print('Possible issue with datetime_uniq,ripf_uniq.')\n",
    "  \n",
    "#if(len(datetime_all)!=len(ripf_all)):\n",
    "#  print('Possible issue with datetime/ripf.')\n",
    "\n",
    "#raise Exception('STOP!')\n",
    "\n",
    "'''for forc_i,forc_input_file in enumerate(forc_files):\n",
    "  print('forc_i,forc_input_file=',forc_i,forc_input_file)\n",
    "\n",
    "  forc_input_file_tail=os.path.basename(forc_input_file)\n",
    "  forc_var,forc_table,forc_experiment,forc_model,forc_ripf,forc_grid,forc_datetime=cmor_file_parts(forc_input_file_tail)\n",
    "  forc_rval,forc_ival,forc_pval,forc_fval=cmor_ripf_parts(forc_ripf)\n",
    "  forc_input_file_head=string.split(forc_input_file,sep=forc_input_file_tail)[0]\n",
    "  forc_version,forc_grid,forc_var,forc_table,forc_ripf,forc_experiment,forc_model,forc_institution,forc_activity,forc_cmip=cmor_directory_parts(forc_input_file_head)\n",
    "\n",
    "  #print('forc_input_file_tail=',forc_input_file_tail)\n",
    "  #print('forc_input_file_head=',forc_input_file_head)\n",
    "\n",
    "  #print('forc_var,forc_table,forc_experiment,forc_model,forc_ripf,forc_grid,forc_datetime=',forc_var,forc_table,forc_experiment,forc_model,forc_ripf,forc_grid,forc_datetime)\n",
    "  \n",
    "  #raise Exception('STOP!')\n",
    "  #print('forc_version,forc_grid,forc_var,forc_table,forc_ripf,forc_experiment,forc_model,forc_institution,forc_activity,forc_cmip=',forc_version,forc_grid,forc_var,forc_table,forc_ripf,forc_experiment,forc_model,forc_institution,forc_activity,forc_cmip)\n",
    "  \n",
    "  #print('forc_ripf=',forc_ripf)\n",
    "  #print('forc_rval,forc_ival,forc_pval,forc_fval',forc_rval,forc_ival,forc_pval,forc_fval)\n",
    "  \n",
    "  #raise Exception('STOP!')'''\n",
    "  \n",
    "#raise Exception('STOP!')\n",
    "\n",
    "forc_ens=[]\n",
    "forc_ifhs=[]\n",
    "forc_date_time_stamps=[]\n",
    "forc_times=[]\n",
    "forc_mins=[]\n",
    "forc_maxs=[]\n",
    "for forc_i,forc_input_file in enumerate(forc_files):\n",
    "  #print('forc_i,forc_input_file=',forc_i,forc_input_file)\n",
    "  forc_input_file_tail=os.path.basename(forc_input_file)\n",
    "  forc_input_file_head=string.split(forc_input_file,sep=forc_input_file_tail)[0]\n",
    "  forc_var,forc_table,forc_experiment,forc_model,forc_ripf,forc_grid,forc_datetime=cmor_file_parts(forc_input_file_tail)\n",
    "  forc_rval,forc_ival,forc_pval,forc_fval=cmor_ripf_parts(forc_ripf)\n",
    "  forc_version,forc_grid,forc_var,forc_table,forc_ripf,forc_experiment,forc_model,forc_institution,forc_activity,forc_cmip=cmor_directory_parts(forc_input_file_head)\n",
    "  forc_ens.append(forc_rval)\n",
    "  #print('forc_i=',forc_i)\n",
    "  forc_ifhs.append(netCDF4.Dataset(forc_input_file))\n",
    "  forc_times.append(forc_ifhs[forc_i].variables['time'])\n",
    "  \n",
    "  #print(forc_ifhs[forc_i].variables['time'][:])\n",
    "  \n",
    "  #raise Exception('STOP!')                 \n",
    "\n",
    "  if(forc_i>0):  #check that same calendar/units are in use, otherwise would need to convert to common one.\n",
    "    if(forc_times[forc_i].calendar!=forc_times[forc_i-1].calendar): raise SystemExit('calendars not matching:'+__file__+' line number: '+str(inspect.stack()[0][2]))\n",
    "    if(forc_times[forc_i].units!=forc_times[forc_i-1].units): raise SystemExit('units not matching:'+__file__+' line number: '+str(inspect.stack()[0][2]))\n",
    "  forc_mins.append(min(forc_times[forc_i]))\n",
    "  forc_maxs.append(max(forc_times[forc_i]))\n",
    "  forc_date_time_stamps.append(netCDF4.num2date(forc_times[forc_i][:],forc_times[forc_i].units,forc_times[forc_i].calendar))\n",
    "  if(forc_i==0):\n",
    "    npforc_date_time_stamps=netCDF4.num2date(forc_times[forc_i][:],forc_times[forc_i].units,forc_times[forc_i].calendar)\n",
    "  else:\n",
    "    npforc_date_time_stamps=np.append(npforc_date_time_stamps,netCDF4.num2date(forc_times[forc_i][:],forc_times[forc_i].units,forc_times[forc_i].calendar))\n",
    "\n",
    "np.set_printoptions(threshold='nan') #will print out whole array\n",
    "np.set_printoptions(threshold=1000)\n",
    "\n",
    "#print('npforc_date_time_stamps=',npforc_date_time_stamps)\n",
    "#print('npforc_date_time_stamps.shape=',npforc_date_time_stamps.shape)\n",
    "\n",
    "#raise Exception('STOP!')                 \n",
    "#print('type(forc_date_time_stamps)=',type(forc_date_time_stamps))\n",
    "\n",
    "#cont_date_time_stamp=netCDF4.num2date(cont_time[:],cont_time.units,cont_time.calendar)\n",
    "\n",
    "#forc_year_fraction_monthly=fractional_year_from_num2date(npforc_date_time_stamps,forc_times[0].calendar) #get working...\n",
    "\n",
    "#print('forc_times=',forc_times[:])\n",
    "#print('len(forc_times)=',len(forc_times))\n",
    "\n",
    "forc_nptimes=np.array(forc_times)\n",
    "#print('forc_nptimes=',forc_nptimes[:])\n",
    "#print('type(forc_nptimes)=',type(forc_nptimes))\n",
    "\n",
    "forc_npftimes=forc_nptimes.flatten()\n",
    "#print('forc_npftimes=',forc_npftimes[:])\n",
    "\n",
    "#print(forc_ens)\n",
    "#print('forc_year_fraction_monthly=',forc_year_fraction_monthly[:])\n",
    "\n",
    "#raise Exception('STOP!')\n",
    "\n",
    "year_min=netCDF4.num2date(min(forc_mins),forc_times[0].units,forc_times[0].calendar)\n",
    "year_max=netCDF4.num2date(max(forc_maxs),forc_times[0].units,forc_times[0].calendar) #not year maximum as need to consider end of each experiment.\n",
    "\n",
    "print('Forecast start times go from ',min(forc_mins),'to',max(forc_maxs),'or',year_min,'to',year_max)\n",
    "print('Ensembles go from ',min(forc_ens),'to',max(forc_ens)) #later on restrict to only valid ensembles to keep arrays compact eg 1,2,3,11 but not 1-11.\n",
    "\n",
    "npforc_ens=np.array(list(set(forc_ens)))\n",
    "#print(npforc_ens)\n",
    "\n",
    "forc_nens=len(npforc_ens)\n",
    "#print(forc_nens)\n",
    "#raise Exception('STOP!')\n",
    "\n",
    "forc_ybeg,forc_yend=netCDF4.num2date(min(forc_mins),forc_times[0].units,forc_times[0].calendar).year,netCDF4.num2date(max(forc_maxs),forc_times[0].units,forc_times[0].calendar).year\n",
    "forc_mbeg,forc_mend=netCDF4.num2date(min(forc_mins),forc_times[0].units,forc_times[0].calendar).month,netCDF4.num2date(max(forc_maxs),forc_times[0].units,forc_times[0].calendar).month\n",
    "\n",
    "forc_j,forc_k,forc_date,forc_m,forc_n,forc_time=get_timestamp_number(forc_ybeg,forc_yend,forc_mbeg,forc_mend,forc_times[0].units,forc_times[0].calendar)\n",
    "\n",
    "forc_ntime=len(forc_time)\n",
    "\n",
    "forc_year_fraction_monthly=fractional_year_from_num2date(forc_date,forc_times[0].calendar)\n",
    "\n",
    "#print('forc_ntime=',forc_ntime)\n",
    "#print('forc_date=',forc_date)\n",
    "#print('forc_time',forc_time)\n",
    "#print('forc_year_fraction_monthly',forc_year_fraction_monthly)\n",
    "\n",
    "forc_beg,forc_end=np.zeros(nforc_files,dtype=int),np.zeros(nforc_files,dtype=int) #these indices define beg/end time indices.\n",
    "\n",
    "forc_beg_cnt,forc_end_cnt=np.zeros((forc_nens,(nforc_files/forc_nens)),dtype=int),np.zeros((forc_nens,(nforc_files/forc_nens)),dtype=int) #these indices define beg/end time indices.\n",
    "\n",
    "#forc_nino_monthly=ma.zeros((nindices_nino,forc_nens,nforc_files,forc_ntime),dtype=float) #set them all to missing value, then assign the forecasts to the various segments in the array.\\\\\n",
    "#forc_nino_monthly=ma.masked_equal( np.zeros((nindices_nino,forc_nens,(nforc_files/forc_nens),forc_ntime),dtype=float), 0.0) #set them all to missing value, then assign the forecasts to the various segments in the array.\\\\\n",
    "forc_nino_monthly=ma.masked_equal( np.zeros((forc_ntime,nindices_nino,forc_nens,(nforc_files/forc_nens)),dtype=float), 0.0) #set them all to missing value, then assign the forecasts to the various segments in the array.\\\\\n",
    "\n",
    "#print(forc_nino_monthly.view(ma.MaskedArray))\n",
    "\n",
    "#raise Exception('STOP!')\n",
    "\n",
    "forc_check=ma.zeros((nindices_nino,forc_nens,nforc_files),dtype=int) #can use this to see how the arrays are being populated or not.\n",
    "\n",
    "#print('forc_nino_monthly.shape=',forc_nino_monthly.shape)\n",
    "\n",
    "#forc_beg_date=date_time_stamps[0][0] #this would be for total set of times.\n",
    "#forc_end_date=date_time_stamps[-1][-1] #this would be for total set of times.\n",
    "\n",
    "forc_lat=forc_ifhs[0].variables['latitude'][:,0]\n",
    "forc_lon=forc_ifhs[0].variables['longitude'][0,:]\n",
    "\n",
    "#forc_clat=np.zeros((nforc_files),dtype=float) #if latitude range changes will need to have custom clat.\n",
    "forc_clat=np.cos(forc_lat[:]*rad)\n",
    "\n",
    "#print(indices_i[0][0],indices_i[0][1])\n",
    "\n",
    "forc_ens_cnt=np.zeros(forc_nens,dtype=int)\n",
    "\n",
    "for forc_i,forc_input_file in enumerate(forc_files):\n",
    "  #print('forc_i,forc_input_file=',forc_i,forc_input_file)\n",
    "  forc_beg_date=forc_date_time_stamps[forc_i][0]\n",
    "  forc_end_date=forc_date_time_stamps[forc_i][-1]\n",
    "  #print('forc_beg_date,forc_end_date=',forc_beg_date,forc_end_date)\n",
    "\n",
    "  #print(type(forc_ens[forc_i]))\n",
    "  #raise Exception('STOP!')\n",
    "  \n",
    "  forc_begend_time=[netCDF4.date2num(forc_beg_date,forc_times[i].units,forc_times[i].calendar), netCDF4.date2num(forc_end_date,forc_times[i].units,forc_times[i].calendar)]\n",
    "  forc_loc_beg,forc_loc_end=np.where(forc_time[:]==forc_begend_time[0],1,0),np.where(forc_time[:]==forc_begend_time[1],1,0)\n",
    "  \n",
    "  #print('forc_loc_beg,end=',forc_loc_beg,forc_loc_end)\n",
    "  #raise Exception('STOP!')\n",
    "  \n",
    "  forc_beg[forc_i],forc_end[forc_i]=np.argmax(forc_loc_beg),np.argmax(forc_loc_end)\n",
    "  \n",
    "  #print('forc_beg[forc_i],end[forc_i]=',forc_beg[forc_i],forc_end[forc_i])\n",
    "  #raise Exception('STOP!')\n",
    "  forc_ens_cnt[forc_ens[forc_i]-1]+=1 #this is used to put each forecast in a unique ensemble slot in the array, the array should end up with equal values if the experiment is consistent & complete.\n",
    "  #print('forc_ens[forc_i]-1,forc_ens_cnt[forc_ens[forc_i]-1]-1=',forc_ens[forc_i]-1,forc_ens_cnt[forc_ens[forc_i]-1]-1)\n",
    "\n",
    "  forc_beg_cnt[forc_ens[forc_i]-1,forc_ens_cnt[forc_ens[forc_i]-1]-1], forc_end_cnt[forc_ens[forc_i]-1,forc_ens_cnt[forc_ens[forc_i]-1]-1]=np.argmax(forc_loc_beg),np.argmax(forc_loc_end)\n",
    "\n",
    "  for forc_k,indice in enumerate(indices_nino):\n",
    "\n",
    "    forc_imin,forc_imax=indices_i_gn[forc_k][0],indices_i_gn[forc_k][1]\n",
    "    forc_jmin,forc_jmax=indices_j_gn[forc_k][0],indices_j_gn[forc_k][1]\n",
    "    #print('forc_k,forc_imin,forc_imax,forc_jmin,forc_jmax=',forc_k,forc_imin,forc_imax,forc_jmin,forc_jmax)\n",
    "    \n",
    "    forc_check[forc_k,forc_ens[forc_i]-1,forc_i] = forc_check[forc_k,forc_ens[forc_i]-1,forc_i] + 1 #can use this for checking. For example, onece an array is set to 1 it should not be set/reset again (no overlap).\n",
    "    \n",
    "    #print('forc_nino_monthly.shape=',forc_nino_monthly.shape)\n",
    "    \n",
    "    forc_nino_monthly[forc_beg[forc_i]:forc_end[forc_i]+1,forc_k,forc_ens[forc_i]-1,forc_ens_cnt[forc_ens[forc_i]-1]-1]=np.average(np.average(forc_ifhs[forc_i].variables[dvar][:,forc_jmin:forc_jmax+1,forc_imin:forc_imax+1],axis=1,weights=forc_clat[forc_jmin:forc_jmax+1]),axis=1)\n",
    "    #forc_nino_monthly[forc_k,forc_ens[forc_i]-1,forc_ens_cnt[forc_ens[forc_i]-1]-1,forc_beg[forc_i]:forc_end[forc_i]+1]=np.average(np.average(forc_ifhs[forc_i].variables[dvar][:,forc_jmin:forc_jmax+1,forc_imin:forc_imax+1],axis=1,weights=forc_clat[forc_jmin:forc_jmax+1]),axis=1)\n",
    "    #nb = raw_input('Press enter: ')\n",
    "    #raise Exception('STOP!')\n",
    "    \n",
    "#print('forc_check=',forc_check)\n",
    "\n",
    "#print('forc_ens_cnt=',forc_ens_cnt)\n",
    "\n",
    "#print('forc_beg=',forc_beg)\n",
    "#print('forc_end=',forc_end)\n",
    "\n",
    "#print('forc_beg_cnt=',forc_beg_cnt)\n",
    "#print('forc_end_cnt=',forc_end_cnt)\n",
    "\n",
    "#raise Exception('STOP!')\n",
    "\n",
    "#forc_nino_anomaly=forc_nino_monthly.copy() #make a copy in preparation for generating anomalies.\n",
    "\n",
    "forc_nino_anomaly_step1=np.expand_dims(forc_nino_monthly,-1)\n",
    "\n",
    "forc_nino_anomaly=np.tile(forc_nino_anomaly_step1,(cont_nexps))\n",
    "\n",
    "print('forc_nino_anomaly.shape=',forc_nino_anomaly.shape)\n",
    "\n",
    "#raise Exception('STOP!')\n",
    "\n",
    "forc_years,forc_months,forc_days,forc_hours=[],[],[],[]\n",
    "forc_yearsM1,forc_monthsM1,forc_daysM1,forc_hoursM1=[],[],[],[]\n",
    "for forc_i,forc_date_now in enumerate(forc_date):\n",
    "  forc_years.append(forc_date_now.year)\n",
    "  forc_months.append(forc_date_now.month)\n",
    "  forc_days.append(forc_date_now.day)\n",
    "  forc_hours.append(forc_date_now.hour)\n",
    "  forc_yearsM1.append(forc_date_now.year-1)\n",
    "  forc_monthsM1.append(forc_date_now.month-1)\n",
    "  forc_daysM1.append(forc_date_now.day-1)\n",
    "  forc_hoursM1.append(forc_date_now.hour-1)\n",
    "\n",
    "plot_index='nino34'\n",
    "#plot_index='nino3'\n",
    "#plot_index='nino4'\n",
    "#plot_index='nino1+2'\n",
    "\n",
    "forc_ens_cnt=np.zeros(forc_nens,dtype=int)\n",
    "\n",
    "#some of these things could be put into a function.\n",
    "for forc_i,forc_input_file in enumerate(forc_files):\n",
    "  #print('forc_input_file=',forc_input_file)\n",
    "  forc_year_beg,forc_year_end=forc_years[forc_beg[forc_i]],forc_years[forc_end[forc_i]]\n",
    "  forc_month_beg,forc_month_end=forc_months[forc_beg[forc_i]],forc_months[forc_end[forc_i]]\n",
    "  #print('forc_i,forc_year_beg,forc_year_end,forc_month_beg,forc_month_end=',forc_year_beg,forc_year_end,forc_month_beg,forc_month_end)\n",
    "  forc_mbeg,forc_mend=1,12\n",
    "  forc_ybeg,forc_yend=forc_year_beg,forc_year_end\n",
    "  if(forc_month_beg!=1): forc_ybeg,forc_mbeg=forc_year_beg+1,1\n",
    "  if(forc_month_end!=12): forc_yend,forc_mbeg=forc_year_end-1,12\n",
    "  #print('forc_ybeg,forc_yend,forc_mbeg,forc_mend=',forc_ybeg,forc_yend,forc_mbeg,forc_mend)\n",
    "  #print('forc_nino_monthly[0,forc_ens[forc_i],forc_i,forc_beg[forc_i]:forc_end[forc_i]]=',forc_nino34_monthly[0,0,forc_i,forc_beg[forc_i]:forc_end[forc_i]])\n",
    "  #print('forc_years[forc_beg[forc_i]:forc_end[forc_i]+1]=',forc_years[forc_beg[forc_i]:forc_end[forc_i]+1])\n",
    "  #print('forc_months[forc_beg[forc_i]:forc_end[forc_i]+1]=',forc_months[forc_beg[forc_i]:forc_end[forc_i]+1]) #just need to find first january and last december in this list to extract the indices...\n",
    "  #print('cont_nino_climatology.shape=',cont_nino_climatology.shape) #cont_climatology calculated in previous cell.\n",
    "  #print('cont_nino_climatology=',cont_nino_climatology)\n",
    "  #print(cont_nino_climatology[0,forc_months[forc_beg[forc_i]]-1])\n",
    "  #print(cont_nino_climatology[0,forc_monthsM1[forc_beg[forc_i]:forc_end[forc_i]+1]])\n",
    "  forc_ens_cnt[forc_ens[forc_i]-1]+=1\n",
    "  \n",
    "  for forc_k,indice in enumerate(indices_nino):\n",
    "    #None\n",
    "    #pretty cool, each month of climatology (12) will be mapped according to values of forc_monthsM1\n",
    "    #forc_nino_anomaly[forc_k,forc_ens[forc_i]-1,forc_ens_cnt[forc_ens[forc_i]-1]-1,forc_beg[forc_i]:forc_end[forc_i]+1] = forc_nino_anomaly[forc_k,forc_ens[forc_i]-1,forc_ens_cnt[forc_ens[forc_i]-1]-1,forc_beg[forc_i]:forc_end[forc_i]+1] - cont_nino_climatology[forc_k,forc_monthsM1[forc_beg[forc_i]:forc_end[forc_i]+1]]\n",
    "    for clim_cnt in range(cont_nexps):\n",
    "      forc_nino_anomaly[forc_beg[forc_i]:forc_end[forc_i]+1,forc_k,forc_ens[forc_i]-1,forc_ens_cnt[forc_ens[forc_i]-1]-1,clim_cnt] = \\\n",
    "      forc_nino_anomaly[forc_beg[forc_i]:forc_end[forc_i]+1,forc_k,forc_ens[forc_i]-1,forc_ens_cnt[forc_ens[forc_i]-1]-1,clim_cnt] - cont_nino_climatology[forc_monthsM1[forc_beg[forc_i]:forc_end[forc_i]+1],forc_k,clim_cnt]\n",
    "  \n",
    "  #raise Exception('STOP!')\n",
    "  \n",
    "  forc_year_fraction_monthly_now=forc_year_fraction_monthly[forc_beg[forc_i]:forc_end[forc_i]+1]\n",
    "  #forc_zero=np.zeros(len(forc_year_fraction_monthly_now))\n",
    "  #print('forc_ens[forc_i]=',forc_ens[forc_i])\n",
    "  \n",
    "  plt.plot(forc_year_fraction_monthly_now,forc_nino_anomaly[forc_beg[forc_i]:forc_end[forc_i]+1,indices_nino.index(plot_index),forc_ens[forc_i]-1,forc_ens_cnt[forc_ens[forc_i]-1]-1,0],color=colors[forc_ens[forc_i]-1])\n",
    "  #plt.plot(forc_year_fraction_monthly_now,forc_nino_monthly[indices_nino.index(plot_index),forc_ens[forc_i]-1,forc_ens_cnt[forc_ens[forc_i]-1]-1,forc_beg[forc_i]:forc_end[forc_i]+1],color=colors[forc_ens[forc_i]-1])\n",
    "\n",
    "#  print('Times go from ',min(forc_mins),'to',max(forc_maxs),'or',\n",
    "#    netCDF4.num2date(min(forc_mins),forc_times[0].units,forc_times[0].calendar),'to',netCDF4.num2date(max(forc_maxs),forc_times[0].units,forc_times[0].calendar))\n",
    "\n",
    "#year_min=netCDF4.num2date(min(forc_mins),forc_times[0].units,forc_times[0].calendar)\n",
    "#year_max=netCDF4.num2date(max(forc_maxs),forc_times[0].units,forc_times[0].calendar)\n",
    "all_years=np.linspace(year_min.year,year_max.year,20)\n",
    "forc_zero=np.zeros(len(all_years))\n",
    "plt.plot(all_years,forc_zero,color='black',linestyle=':')\n",
    "\n",
    "plt.title('Monthly anomalies')\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel(indices_label[indices_nino.index(plot_index)]+' ($^o$C)')\n",
    "plt.show()\n",
    "\n",
    "print('END')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Process observational data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print('BEGIN')\n",
    "\n",
    "import matplotlib.ticker as ticker\n",
    "\n",
    "nmy=12\n",
    "\n",
    "#del(generate_daily_month_indices)\n",
    "#from decadal_diag import generate_daily_month_indices\n",
    "\n",
    "ivar_str_ncep2='skt'\n",
    "#ybeg_data_ncep2=2000 #2000;assume 12 months in every year.\n",
    "#yend_data_ncep2=2016 #2016;assume 12 months in every year.\n",
    "\n",
    "grid_label='ncep2'\n",
    "ncep2_files_string='/short/v14/mac599/ncep2/skt.sfc.gauss.????.nc'\n",
    "ncep2_files=sorted(glob.glob(ncep2_files_string))\n",
    "nncep2_files=len(ncep2_files)\n",
    "print('ncep2 input_files ('+str(nncep2_files)+') =',ncep2_files)\n",
    "\n",
    "ifhN_ncep2=netCDF4.MFDataset(ncep2_files)\n",
    "ifh0_ncep2=netCDF4.MFDataset(ncep2_files[0])\n",
    "\n",
    "time_daily_ncep2=ifhN_ncep2.variables['time']\n",
    "#lat=ifh0.variables['lat']\n",
    "#lon=ifh0.variables['lon']\n",
    "\n",
    "lat_ncep2=ifh0_ncep2.variables['lat'][:]\n",
    "lon_ncep2=ifh0_ncep2.variables['lon'][:]\n",
    "\n",
    "ntime_daily_ncep2=len(time_daily_ncep2)\n",
    "#print('lat=',lat)\n",
    "\n",
    "ifh0_ncep2.close()\n",
    "\n",
    "#print('rad=',rad)\n",
    "clat_ncep2=np.cos(lat_ncep2[:]*rad)\n",
    "\n",
    "#time_daily_ncep2_calendar='julian'\n",
    "time_daily_ncep2_calendar='proleptic_gregorian'\n",
    "\n",
    "date_time_stamp_daily_ncep2=netCDF4.num2date(time_daily_ncep2[:],time_daily_ncep2.units,time_daily_ncep2_calendar)\n",
    "\n",
    "year_fraction_daily_ncep2=fractional_year_from_num2date(date_time_stamp_daily_ncep2,time_daily_ncep2_calendar)\n",
    "\n",
    "nino_daily_ncep2=ma.zeros((ntime_daily_ncep2,nindices_nino),dtype=float)\n",
    "#raise Exception('STOP!')\n",
    "for ncep2_k,indice in enumerate(indices_nino):\n",
    "  #ncep2_k=indices_nino.index(indice)\n",
    "  ncep2_imin,ncep2_imax=indices_i_ncep2[indices_nino.index(indice)][0],indices_i_ncep2[indices_nino.index(indice)][1]\n",
    "  ncep2_jmin,ncep2_jmax=indices_j_ncep2[indices_nino.index(indice)][0],indices_j_ncep2[indices_nino.index(indice)][1]\n",
    "  #print('indices.index(indice),imin_ncep2,imax_ncep2,jmin_ncep2,jmax_ncep2=',indices.index(indice),imin_ncep2,imax_ncep2,jmin_ncep2,jmax_ncep2)\n",
    "  #print('lat_ncep2,lon_ncep2=',lat_ncep2[jmin_ncep2:jmax_ncep2],lon_ncep2[imin_ncep2:imax_ncep2])\n",
    "\n",
    "  nino_daily_ncep2[:,ncep2_k]=np.average(np.average(ifhN_ncep2.variables['skt'][:,ncep2_jmin:ncep2_jmax,ncep2_imin:ncep2_imax],axis=1,weights=clat_ncep2[ncep2_jmin:ncep2_jmax]),axis=1)\n",
    "\n",
    "#indices_i_gn,indices_j_gn=[[122,152],[122,152],[122,152]],[[110,160],[110,160],[110,160]] #check this, whether I need +1 also...\n",
    "#raise Exception('STOP!')\n",
    "#print('year_fraction_ncep2=',year_fraction_monthly_ncep2)\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "fix,ax=plt.subplots()\n",
    "for ncep2_k,indice in enumerate(indices_nino):\n",
    "  ax.plot(year_fraction_daily_ncep2,nino_daily_ncep2[:,indices_nino.index(indice)]-273.15,label=indices_label[indices_nino.index(indices_nino[ncep2_k])])\n",
    "  ax.xaxis.set_major_locator(ticker.MultipleLocator(2.0))\n",
    "legend=ax.legend(loc='lower right',shadow=False,fontsize='x-large')\n",
    "plt.title('Daily values (NCEP2)')\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel('$^o$C')\n",
    "plt.grid(True,linestyle='-')\n",
    "plt.show()\n",
    "\n",
    "what_to_keep='all_times'\n",
    "#what_to_keep='all_but_first'\n",
    "#what_to_keep='all_but_last'\n",
    "#what_to_keep='all_but_first_and_last'\n",
    "\n",
    "if(what_to_keep=='all_times'):\n",
    "  #all times kept:\n",
    "  date_time_stamp_daily_new_ncep2=date_time_stamp_daily_ncep2\n",
    "  nino_daily_new_ncep2=nino_daily_ncep2\n",
    "  year_fraction_daily_new_ncep2=year_fraction_daily_ncep2\n",
    "  \n",
    "elif(what_to_keep=='all_but_first'):\n",
    "  #remove first day:\n",
    "  date_time_stamp_daily_new_ncep2=date_time_stamp_daily_ncep2[1::]\n",
    "  nino_daily_new_ncep2=nino_daily_ncep2[1::,:]\n",
    "  year_fraction_daily_new_ncep2=year_fraction_daily_ncep2[1::]\n",
    "  \n",
    "elif(what_to_keep=='all_but_last'):\n",
    "  #remove last day:\n",
    "  date_time_stamp_daily_new_ncep2=date_time_stamp_daily_ncep2[0:-1]\n",
    "  nino_daily_new_ncep2=nino_daily_ncep2[0:-1,:]\n",
    "  year_fraction_daily_new_ncep2=year_fraction_daily_ncep2[0:-1]\n",
    "\n",
    "elif(what_to_keep=='all_but_first_and_last'):\n",
    "  #remove first and last day:\n",
    "  date_time_stamp_daily_new_ncep2=date_time_stamp_daily_ncep2[1:-1]\n",
    "  nino_daily_new_ncep2=nino_daily_ncep2[1:-1,:]\n",
    "  year_fraction_daily_new_ncep2=year_fraction_daily_ncep2[1:-1]\n",
    "  \n",
    "else:\n",
    "  raise Exception('That what_to_keep option doesnt exist.')\n",
    "\n",
    "#print('len(date_time_stamp_daily_new_ncep2)=',len(date_time_stamp_daily_new_ncep2))\n",
    "\n",
    "daily_month_indice_beg_ncep2,daily_month_indice_end_ncep2,daily_year_beg_ncep2,daily_year_end_ncep2,daily_month_beg_ncep2,daily_month_end_ncep2,daily_day_beg_ncep2,daily_day_end_ncep2,beg_month_partial_ncep2,end_month_partial_ncep2 = \\\n",
    "  generate_daily_month_indices(date_time_stamp_daily_new_ncep2,time_daily_ncep2.units,time_daily_ncep2_calendar,24)\n",
    "\n",
    "#raise Exception('STOP!')\n",
    "\n",
    "nino_monthly_ncep2=ma.zeros((len(daily_month_indice_beg_ncep2),nindices_nino),dtype=float)\n",
    "\n",
    "#raise Exception('STOP!')\n",
    "\n",
    "for k,indice in enumerate(indices_nino):\n",
    "  for month in range(0,len(daily_month_indice_beg_ncep2)):\n",
    "    nino_monthly_ncep2[month,k]=np.average(nino_daily_new_ncep2[daily_month_indice_beg_ncep2[month]:daily_month_indice_end_ncep2[month]+1,k])\n",
    "\n",
    "#raise Exception('STOP!')\n",
    "\n",
    "year_fraction_monthly_ncep2=ma.zeros(len(daily_month_indice_beg_ncep2))\n",
    "for month in range(0,len(daily_month_indice_beg_ncep2)):\n",
    "  year_fraction_monthly_ncep2[month]=np.average(year_fraction_daily_new_ncep2[daily_month_indice_beg_ncep2[month]:daily_month_indice_end_ncep2[month]+1])\n",
    "\n",
    "num_stamp_new_ncep2=netCDF4.date2num(date_time_stamp_daily_new_ncep2,time_daily_ncep2.units,time_daily_ncep2_calendar)\n",
    "\n",
    "num_stamp_monthly_ncep2=np.zeros(len(daily_month_indice_beg_ncep2))\n",
    "#print('num_stamp_monthly_ncep2.shape=',num_stamp_monthly_ncep2.shape)\n",
    "\n",
    "for month in range(0,len(daily_month_indice_beg_ncep2)):\n",
    "  num_stamp_monthly_ncep2[month]=np.average(num_stamp_new_ncep2[daily_month_indice_beg_ncep2[month]:daily_month_indice_end_ncep2[month]+1])\n",
    "\n",
    "date_time_stamp_monthly_ncep2=netCDF4.num2date(num_stamp_monthly_ncep2,time_daily_ncep2.units,time_daily_ncep2_calendar)\n",
    "\n",
    "#raise Exception('STOP!')\n",
    "  \n",
    "EndOption=1 #use months with no days missing at begin and end months (i.e. discard partial months if they exist).\n",
    "#EndOption=2 #use months with days missing at both begin and end months (i.e. don't discard partial months if they exist @ beg/end).\n",
    "\n",
    "print('beg,end_month_partial_ncep2=',beg_month_partial_ncep2,end_month_partial_ncep2)\n",
    "\n",
    "if(EndOption==1):\n",
    "  print('Discarding beg&/end month if they exist.')\n",
    "  if(beg_month_partial_ncep2 or end_month_partial_ncep2):\n",
    "    if(beg_month_partial_ncep2 and end_month_partial_ncep2):\n",
    "      print('type#1')\n",
    "      nino_monthly_new_ncep2=nino_monthly_ncep2[1:-1,:]\n",
    "      year_fraction_monthly_new_ncep2=year_fraction_monthly_ncep2[1:-1]\n",
    "      date_time_stamp_monthly_new_ncep2=date_time_stamp_monthly_ncep2[1:-1]\n",
    "      \n",
    "    elif(not beg_month_partial_ncep2 and end_month_partial_ncep2):\n",
    "      print('type#2')\n",
    "      nino_monthly_new_ncep2=nino_monthly_ncep2[0:-1,:]\n",
    "      year_fraction_monthly_new_ncep2=year_fraction_monthly_ncep2[0:-1]\n",
    "      date_time_stamp_monthly_new_ncep2=date_time_stamp_monthly_ncep2[0:-1]\n",
    "\n",
    "    elif(beg_month_partial_ncep2 and not end_month_partial_ncep2):\n",
    "      print('type#3')\n",
    "      nino_monthly_newnino_ncep2=nino_monthly_ncep2[1::,:]\n",
    "      year_fraction_monthly_new_ncep2=year_fraction_monthly_ncep2[1::]\n",
    "      date_time_stamp_monthly_new_ncep2=date_time_stamp_monthly_ncep2[1::]\n",
    "\n",
    "    else:\n",
    "      raise SystemExit('Shouldnt get here:'+__file__+' line number: '+str(inspect.stack()[0][2]))\n",
    "  else:\n",
    "    print('type#4')\n",
    "    nino_monthly_new_ncep2=nino_monthly_ncep2\n",
    "    year_fraction_monthly_new_ncep2=year_fraction_monthly_ncep2\n",
    "    date_time_stamp_monthly_new_ncep2=date_time_stamp_monthly_ncep2\n",
    "    \n",
    "elif(EndOption==2):\n",
    "  print('Keeping beg/end month or both.')\n",
    "  nino_monthly_new_ncep2=nino_monthly_ncep2\n",
    "  year_fraction_monthly_new_ncep2=year_fraction_monthly_ncep2\n",
    "  date_time_stamp_monthly_new_ncep2=date_time_stamp_monthly_ncep2\n",
    "  \n",
    "else:\n",
    "  raise SystemExit('EndOption can be only 1 or 2:'+__file__+' line number: '+str(inspect.stack()[0][2]))\n",
    "  \n",
    "#raise Exception('STOP!')\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "#print(nino_monthly_new)\n",
    "\n",
    "#print('year_fraction_monthly_new=',year_fraction_monthly_new)\n",
    "\n",
    "plot_index='nino34'\n",
    "#plot_index='nino1+2'\n",
    "#plot_index='nino4'\n",
    "\n",
    "fix,ax=plt.subplots()\n",
    "ax.plot(year_fraction_monthly_new_ncep2,nino_monthly_new_ncep2[:,indices_nino.index(plot_index)]-273.15,color='blue',label='mon_from_day')\n",
    "legend=ax.legend(loc='lower left',shadow=False,fontsize='large')\n",
    "ax.xaxis.set_major_locator(ticker.MultipleLocator(2.0))\n",
    "plt.title('Monthly values (NCEP2)')\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel(indices_label[indices_nino.index(plot_index)] +'($^o$C)')\n",
    "plt.grid(True,linestyle='-')\n",
    "plt.show()\n",
    "\n",
    "#raise Exception('STOP!')\n",
    "\n",
    "#print('date_time_stamp_monthly_new_ncep2=',date_time_stamp_monthly_new_ncep2)\n",
    "ybeg_ncep2=date_time_stamp_monthly_new_ncep2[0].year\n",
    "yend_ncep2=date_time_stamp_monthly_new_ncep2[-1].year\n",
    "\n",
    "ydiff_monthly_ncep2=yend_ncep2-ybeg_ncep2+1\n",
    "\n",
    "MissingMonths_ncep2=False\n",
    "\n",
    "first_month_ncep2=date_time_stamp_monthly_new_ncep2[0].month\n",
    "last_month_ncep2=date_time_stamp_monthly_new_ncep2[-1].month\n",
    "\n",
    "#print('ybeg,yend_ncep2=',ybeg_ncep2,yend_ncep2)\n",
    "\n",
    "#print('first,last_month_ncep2=',first_month_ncep2,last_month_ncep2)\n",
    "\n",
    "missing_months_beg_ncep2,missing_months_end_ncep2=0,0\n",
    "\n",
    "cyear_beg_skip_ncep2,cyear_end_skip_ncep2=0,1\n",
    "\n",
    "if(first_month_ncep2!=1):\n",
    "  missing_months_beg_ncep2=12-first_month_ncep2\n",
    "  cyear_beg_skip_ncep2=1\n",
    "  MissingMonths_ncep2=True\n",
    "\n",
    "if(last_month_ncep2!=12):\n",
    "  missing_months_end_ncep2=12-last_month_ncep2\n",
    "  cyear_end_skip_ncep2=2\n",
    "  MissingMonths_ncep2=True\n",
    "  \n",
    "if(MissingMonths_ncep2):\n",
    "  print('There are missing months in the set. '+str(missing_months_beg_ncep2)+' at beginning and '+str(missing_months_end_ncep2)+' at end.')\n",
    "  print('Currently years with missing months are not used in generating long term monthly climatology.')\n",
    "  print('And missing months will be set to missing in the final time-series.')\n",
    "  ts_beg_ncep2,ts_end_ncep2,ts_avg_ncep2,dt_beg_ncep2,dt_end_ncep2,dt_avg_ncep2=get_timestamp_number(ybeg_ncep2,yend_ncep2,1,12,time_ncep2.units,time_ncep2_calendar)\n",
    "  year_fraction_monthly_full_ncep2=fractional_year_from_num2date(ts_avg_ncep2,time_ncep2_calendar)\n",
    "  \n",
    "  nino_monthly_full_ncep2=ma.masked_all((ydiff_monthly_ncep2*nmy,nindices_nino),dtype=float) #ensure missing months are masked out.\n",
    "  last_month_index_ncep2=ydiff_monthly_ncep2*nmy-last_month_ncep2\n",
    "  nino_monthly_full_ncep2[first_month_ncep2-1:last_month_index_ncep2,:]=nino_monthly_ncep2\n",
    "else:\n",
    "  print('All years have 12 months.')\n",
    "  \n",
    "  nino_monthly_full_ncep2=nino_monthly_ncep2.copy()\n",
    "  year_fraction_monthly_full_ncep2=year_fraction_monthly_ncep2.copy()\n",
    "\n",
    "  print('nino_monthly_full_ncep2.shape=',nino_monthly_full_ncep2.shape)\n",
    "\n",
    "nino_climatology_ncep2,nino_monthly_anomaly_ncep2 = calculate_monthly_climatology_anomaly_from_monthly(nino_monthly_full_ncep2)\n",
    "\n",
    "#raise Exception('STOP!')\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "zero_ncep2=np.zeros(len(year_fraction_monthly_full_ncep2))\n",
    "fix,ax=plt.subplots()\n",
    "for ncep2_k,indice in enumerate(indices_nino):\n",
    "  ax.plot(year_fraction_monthly_full_ncep2,nino_monthly_anomaly_ncep2[:,ncep2_k],label=indices_label[indices_nino.index(indices_nino[ncep2_k])])\n",
    "plt.plot(year_fraction_monthly_full_ncep2,zero_ncep2,color='black')\n",
    "\n",
    "legend=ax.legend(loc='upper left',shadow=False,fontsize='x-large')\n",
    "ax.xaxis.set_major_locator(ticker.MultipleLocator(2.0))\n",
    "plt.title('Monthly anomalies (NCEP2)')\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel('$^o$C')\n",
    "plt.grid(True,linestyle='-')\n",
    "plt.show()\n",
    "\n",
    "print('END')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 8,
        "height": 18,
        "hidden": false,
        "row": 50,
        "width": 4
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "#https://www.youtube.com/watch?v=6SHnmho7zCs\n",
    "\n",
    "def plot_map_box(ind,indices_label,indices_nino,lats_nino,lons_nino):\n",
    "  \n",
    "  colors=['black','red','green','blue','orange','brown','purple','pink']\n",
    "\n",
    "  ind_list=[]\n",
    "  for f,ff in enumerate(ind):\n",
    "    ind_list.append(indices_nino.index(ind[f])) \n",
    "  \n",
    "  map = Basemap(projection='cyl', \n",
    "  llcrnrlat=-30, llcrnrlon=100, \n",
    "  urcrnrlat=30, urcrnrlon=320,\n",
    "  lat_0=0, lon_0=180)\n",
    "\n",
    "  map.drawmapboundary() #fill_color='aqua'\n",
    "  map.fillcontinents() #color='coral',lake_color='aqua'\n",
    "  map.drawcoastlines()\n",
    "\n",
    "  #print('len(lons_nino)=',len(lons_nino))\n",
    "\n",
    "  for i,ii in enumerate(ind_list):\n",
    "    x,y = map(lons_nino[ii], lats_nino[ii])\n",
    "    map.plot(x, y, marker=None,color=colors[i],linewidth=1)\n",
    "    xmid,ymid = map((lons_nino[ii][0]+lons_nino[ii][1])/2,(lats_nino[ii][0]+lats_nino[ii][2])/2)\n",
    "    plt.text(xmid,ymid,indices_nino[ii],fontsize=6,horizontalalignment='left',verticalalignment='center')  \n",
    "\n",
    "  plt.show()\n",
    "\n",
    "def plot_xy_climatology(ind,obs,mod,clim,indices_label,obs_toggle):\n",
    "\n",
    "  colors=['black','red','green','blue','orange','brown','purple','pink']\n",
    "\n",
    "  ind_list=[]\n",
    "  for f,ff in enumerate(ind):\n",
    "    ind_list.append(indices_nino.index(ind[f])) \n",
    "\n",
    "  #print('indices_label=',indices_label)\n",
    "  #print('ind=',ind)\n",
    "  #print('ind_list=',ind_list)\n",
    "  \n",
    "  fix,ax=plt.subplots()\n",
    "  if(obs_toggle):\n",
    "    plt.title('Monthly climatology (obs:dashed)')\n",
    "  else:\n",
    "    plt.title('Monthly climatology')\n",
    "\n",
    "  plt.xlabel('Month')\n",
    "  plt.ylabel('$^o$C')\n",
    "  plt.xticks(range(12), ['J','F','M','A','M','J','J','A','S','O','N','D'], rotation='horizontal')\n",
    "\n",
    "  for j in ind_list:\n",
    "    ax.plot(range(12),mod[:,j,clim],color=colors[j],label=indices_label[j])\n",
    "    if(obs_toggle):\n",
    "      ax.plot(range(12),obs[:,j]-273.15,color=colors[j],linestyle=':')\n",
    "  plt.grid(True,linestyle='-')\n",
    "  legend=ax.legend(loc='lower right',shadow=False,fontsize='small',title='model indice')\n",
    "  plt.show()\n",
    "  \n",
    "def plot_xy_ensemble(mod_month,mod_anom,mod_time,obs_month,obs_anom,obs_time,indices_nino,ensembles, \\\n",
    "  ind,ens,forc,clim,indices_label,forc_beg,forc_end,forc_beg_cnt,forc_end_cnt,files,anom_toggle, \\\n",
    "  ensemble_average_toggle,obs_toggle,print_file):\n",
    "  #print('data.shape=',data.shape)\n",
    "  #print('ind,indices_nino=',ind,indices_nino)\n",
    "  #print('ens,ensembles=',ens,ensembles)\n",
    "  #x=range(len(data[0,0,0,:]))\n",
    "  #print('data.shape=',data.shape)\n",
    "  #print('forc_beg=',forc_beg)\n",
    "  #print('x=',x)\n",
    "  #plt.plot(time,data[ind,ens,0,0:24])\n",
    "  \n",
    "  colors=['black','red','green','blue','orange','brown','purple','pink']\n",
    "\n",
    "  #print('Current Working Directory=',os.getcwd())\n",
    "  \n",
    "  #print('ensemble_average_toggle=',ensemble_average_toggle)\n",
    "\n",
    "  #if(print_file==None):\n",
    "  #  print('No print file.')\n",
    "  #else:\n",
    "  #  print('print_file=',print_file)\n",
    "  #print('jjj=',jjj)\n",
    "  if(anom_toggle):\n",
    "    mod_data=mod_anom[:,:,:,:,clim]\n",
    "    obs_data=obs_anom\n",
    "  else:\n",
    "    mod_data=mod_month\n",
    "    obs_data=obs_month-273.15\n",
    "  #print('obs_data=',obs_data)\n",
    "  \n",
    "  ens_list=[]\n",
    "  for e,ee in enumerate(ens):\n",
    "    ens_list.append(ensembles.index(ens[e]))\n",
    "    \n",
    "  ind_list=[]\n",
    "  for f,ff in enumerate(ind):\n",
    "    ind_list.append(indices_nino.index(ind[f])) \n",
    "    \n",
    "  fix,ax=plt.subplots()\n",
    "  #for i in range(len(files)/len(ensembles)):\n",
    "  for i in forc:\n",
    "    #print('i=',i)\n",
    "    ival=int(i.split(':')[0])\n",
    "    #print('ival=',ival)\n",
    "    #time_now=time[forc_beg[i]:forc_end[i]+1]\n",
    "    for j in ind_list:\n",
    "      #print('i,j,k=',i,j,k)\n",
    "      for k in ens_list:\n",
    "          mod_time_now_test=mod_time[forc_beg_cnt[k,ival]:forc_end_cnt[k,ival]+1]\n",
    "      #print('time_now=',time_now)\n",
    "          if(i==forc[0] and j==ind_list[0]):\n",
    "            legend_label=str(k+1)\n",
    "          else:\n",
    "            legend_label=None\n",
    "          ax.plot(mod_time_now_test,mod_data[forc_beg_cnt[k,ival]:forc_end_cnt[k,ival]+1,j,k,ival],color=colors[k],label=legend_label) #indice, ensemble, forecast, date/time\n",
    "          if(k==ens_list[-1] and ensemble_average_toggle):\n",
    "            #xxx=np.average(data[j,ens_list,ival,forc_beg_cnt[k,ival]:forc_end_cnt[k,ival]+1],axis=0)\n",
    "            #print(xxx.shape)\n",
    "            ax.plot(mod_time_now_test,np.average(mod_data[forc_beg_cnt[k,ival]:forc_end_cnt[k,ival]+1,j,ens_list,ival],axis=0),color='pink',linewidth=5,label='EA')\n",
    "          if(k==ens_list[-1] and i==forc[-1] and obs_toggle):\n",
    "            #print('hello')\n",
    "            ax.plot(obs_time,obs_data[:,j],color='lightblue',label='ncep2')\n",
    "  legend=ax.legend(loc='lower right',shadow=False,fontsize='small',title='ensemble')\n",
    "  if(anom_toggle):\n",
    "    plt.title('Monthly anomalies')\n",
    "    plt.ylim([-6,6])\n",
    "  else:\n",
    "    plt.title('Monthly values')\n",
    "  plt.xlabel('Year')\n",
    "  #plt.ylabel(indices_label[ind]+' ($^o$C)')\n",
    "  plt.ylabel('$^o$C')\n",
    "  #plt.ylim([-6,6])\n",
    "  #plt.xlim([x[0],x[-1]])\n",
    "  plt.xlim([mod_time[0],mod_time[-1]])\n",
    "  plt.grid(True,linestyle='-')\n",
    "  if(anom_toggle):\n",
    "   cont_zero=np.zeros(len(mod_time))\n",
    "   plt.plot(mod_time,cont_zero,color='black',linestyle=':')\n",
    "  \n",
    "  if(print_file==None or print_file==''):\n",
    "    None\n",
    "  else:\n",
    "    print('Print to file ',os.getcwd()+'/'+print_file+'.png')\n",
    "    dummy=plt.savefig(print_file)\n",
    "    print_file=None\n",
    "  #return()\n",
    "  plt.show()\n",
    "\n",
    "#print('datetime_uniq=',datetime_uniq+range(24))\n",
    "datetime_uniq_xtra=[str(i)+': '+datetime_uniq[i] for i,x in enumerate(datetime_uniq)]\n",
    "#print('datetime_uniq_xtra=',datetime_uniq_xtra)\n",
    "#print('len(datetime_uniq_xtra)=',len(datetime_uniq_xtra))\n",
    "\n",
    "clim_list=[int(i) for i in range(cont_nexps)]\n",
    "\n",
    "ensembles=[str(x) for x in npforc_ens]\n",
    "indice_nino_multi = widgets.SelectMultiple(description='nino indice',options=indices_nino,value=indices_nino)\n",
    "ensemble_multi = widgets.SelectMultiple(description='ensemble',options=ensembles,value=ensembles)\n",
    "forecast_multi = widgets.SelectMultiple(description='forecast',options=datetime_uniq_xtra,value=datetime_uniq_xtra,rows=10)\n",
    "climatology_select = widgets.Select(description='climatology',options=clim_list,value=clim_list[0])\n",
    "anom_toggle=widgets.ToggleButton(description='Show Anoms.',value=False,button_style='danger')\n",
    "\n",
    "obs_toggle=widgets.ToggleButton(description='Show Obs.',value=False,button_style='danger')\n",
    "\n",
    "ensemble_average_toggle=widgets.ToggleButton(description='Show E.A.',value=False,button_style='danger')\n",
    "style = {'description_width': 'initial'}\n",
    "print_file=widgets.Text(value=None,placeholder='Print file name (delete to stop printing it):',continuous_update=False,style=style)\n",
    "\n",
    "#container=widgets.HBox([indice_nino_multi,ensemble_multi,forecast_multi,anom_toggle])\n",
    "#display(container)\n",
    "  \n",
    "#raise Exception('STOP!')\n",
    "#print(print_file.keys)\n",
    "#print(anom_toggle.keys)\n",
    "\n",
    "mod_month_to_plot=forc_nino_monthly.copy()\n",
    "mod_anom_to_plot=forc_nino_anomaly.copy()\n",
    "\n",
    "obs_month_to_plot=nino_monthly_full_ncep2.copy()\n",
    "obs_anom_to_plot=nino_monthly_anomaly_ncep2.copy()\n",
    "\n",
    "dummy1=widgets.interact(plot_xy_ensemble, mod_time=widgets.fixed(forc_year_fraction_monthly), mod_month=widgets.fixed(mod_month_to_plot), mod_anom=widgets.fixed(mod_anom_to_plot), \\\n",
    "  obs_time=widgets.fixed(year_fraction_monthly_full_ncep2), obs_month=widgets.fixed(obs_month_to_plot), obs_anom=widgets.fixed(obs_anom_to_plot), \\\n",
    "  indices_label=widgets.fixed(indices_label), forc_beg=widgets.fixed(forc_beg), forc_end=widgets.fixed(forc_end), \\\n",
    "  forc_beg_cnt=widgets.fixed(forc_beg_cnt), forc_end_cnt=widgets.fixed(forc_end_cnt), files=widgets.fixed(forc_files), \\\n",
    "  ensembles=widgets.fixed(ensembles), indices_nino=widgets.fixed(indices_nino), \\\n",
    "  ind=indice_nino_multi, ens=ensemble_multi, forc=forecast_multi, clim=climatology_select, anom_toggle=anom_toggle, ensemble_average_toggle=ensemble_average_toggle, obs_toggle=obs_toggle, print_file=print_file)\n",
    "\n",
    "dummy2=widgets.interact(plot_xy_climatology, \\\n",
    "  obs=widgets.fixed(nino_climatology_ncep2), mod=widgets.fixed(cont_nino_climatology), \\\n",
    "  indices_label=widgets.fixed(indices_label), \\\n",
    "  ind=indice_nino_multi, clim=climatology_select, obs_toggle=obs_toggle)\n",
    "\n",
    "dummy3=widgets.interact(plot_map_box, indices_label=widgets.fixed(indices_label), indices_nino=widgets.fixed(indices_nino), lats_nino=widgets.fixed(lats_nino), lons_nino=widgets.fixed(lons_nino), \\\n",
    "  ind=indice_nino_multi, obs_toggle=obs_toggle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 4,
        "height": 4,
        "hidden": false,
        "row": 62,
        "width": 4
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "outputs": [],
   "source": [
    "for i in [0,3,5]:\n",
    "  print(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "hidden": true
       },
       "report_default": {
        "hidden": true
       }
      }
     }
    }
   },
   "outputs": [],
   "source": [
    "from __future__ import print_ function\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#%matplotlib nbagg\n",
    "#%matplotlib notebook\n",
    "%matplotlib inline\n",
    "\n",
    "print('forc_nino_anomaly.shape=',forc_nino_anomaly.shape) #indice, ensemble, forecast, date/time\n",
    "\n",
    "#indices_nino=['nino34','nino3','nino4']\n",
    "\n",
    "print('indices_nino=',indices_nino)\n",
    "\n",
    "indices_nino_dict={}\n",
    "for i,d in enumerate(indices_nino):\n",
    "  indices_nino_dict[d]=i\n",
    "\n",
    "print('indices_nino_dict=',indices_nino_dict)\n",
    "\n",
    "ensembles=[str(x) for x in npforc_ens]\n",
    "\n",
    "print('ensembles=',ensembles)\n",
    "\n",
    "print('forc_nens=',forc_nens)\n",
    "\n",
    "#fig, ax = plt.subplots(1,figsize=(10,4))\n",
    "#plt.suptitle('Monthly anomalies')\n",
    "\n",
    "def test_plot(ensembles, indices_nino):\n",
    "  #from IPython.display import display\n",
    "  print('ensembles,indces_nino=',ensembles,indices_nino)\n",
    "  #display(ensemble_range, indice_nino_range)\n",
    "  #ax.clear()\n",
    "  #all_years=np.linspace(2001,2020,20)\n",
    "  #line=all_years-2001+.1\n",
    "  #t = np.arange(0., 5., 0.2)\n",
    "  #ax.plot( t, t**2, 'bs')\n",
    "  #ax.plot(all_years,line,color='black',linestyle=':')\n",
    "  #ax.legend(loc=1)\n",
    "  #plt.show()\n",
    "  plt.figure(2)\n",
    "  x=np.linspace(-10,10,num=100)\n",
    "  plt.plot(x,x**ensembles+ensembles*10)\n",
    "  plt.show()\n",
    "\n",
    "\n",
    "#ensemble_range = widgets.SelectMultiple(description='ensemble',options=ensembles,value=ensembles,disabled=False)\n",
    "\n",
    "#indice_nino_range = widgets.SelectMultiple(description='nino indice',options=indices_nino,value=indices_nino,disabled=False)\n",
    "\n",
    "#display(ensemble_range, indice_nino_range)\n",
    "\n",
    "#interactive_plot=widgets.interactive(interactive_plot, ensembles=(0,forc_nens-1,1), indices_nino={'nino34':0, 'nino3':1, 'nino4':2})\n",
    "testit=widgets.interactive(test_plot, ensembles=(0,forc_nens-1,1), indices_nino=indices_nino_dict)\n",
    "\n",
    "#print(testit.value)\n",
    "\n",
    "#del(interactive_plot)\n",
    "#ensemble_range.observe(update_plot, names='value')\n",
    "\n",
    "#indice_nino_range.observe(update_plot, names='value')\n",
    "\n",
    "testit\n",
    "#interactive_plot\n",
    "\n",
    "#forc_nino_anomaly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "hidden": true
       },
       "report_default": {
        "hidden": true
       }
      }
     }
    }
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "print('npforc_ens=',npforc_ens)\n",
    "\n",
    "ensembles=[str(x) for x in npforc_ens]\n",
    "print('ensembles=',ensembles)\n",
    "\n",
    "import ipywidgets as widgets\n",
    "\n",
    "#print(widgets.Button.on_click.__doc__)\n",
    "\n",
    "ensemble_range = widgets.SelectMultiple(description='ensemble',options=ensembles,value=ensembles)\n",
    "\n",
    "#print(ensemble_range.keys)\n",
    "\n",
    "#indice_nino_range = widgets.SelectMultiple(description='nino indice',options=['nino34','nino3'],value=['nino34'])\n",
    "indice_nino_range = widgets.SelectMultiple(description='nino indice',options=indices_nino,value=indices_nino)\n",
    "\n",
    "container=widgets.HBox([ensemble_range,indice_nino_range])\n",
    "\n",
    "display(container)\n",
    "#display(ensemble_range, indice_nino_range)\n",
    "#display(indice_nino_range)\n",
    "\n",
    "from IPython.display import display\n",
    "\n",
    "#def plotit():\n",
    "fig, ax = plt.subplots(1,figsize=(10,4))\n",
    "\n",
    "def on_value_change(change):\n",
    "  #ax.clear()\n",
    "  plt.figure(2)\n",
    "  x=np.linspace(-10,10,num=100)\n",
    "  plt.plot(x,x**3+1*10)\n",
    "  plt.show()\n",
    "  print('hello')\n",
    "  print('old',change['old'])\n",
    "  print('new',change['new'])\n",
    "  print('name',change['name'])\n",
    "\n",
    "ensemble_range.observe(on_value_change, names='value')\n",
    "\n",
    "indice_nino_range.observe(on_value_change, names='value')\n",
    "\n",
    "print('ensemble_range=',ensemble_range.value)\n",
    "\n",
    "print('indice_nino_range=',indice_nino_range.value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "hidden": true
       },
       "report_default": {
        "hidden": true
       }
      }
     }
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from ipywidgets import interactive\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "def f(m, b):\n",
    "    plt.figure(2)\n",
    "    x = np.linspace(-10, 10, num=1000)\n",
    "    plt.plot(x, m * x + b)\n",
    "    plt.ylim(-5, 5)\n",
    "    plt.show()\n",
    "\n",
    "interactive_plot = interactive(f, m=(-2.0, 2.0), b=(-3, 3, 0.5))\n",
    "output = interactive_plot.children[-1]\n",
    "output.layout.height = '350px'\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "hidden": true
       },
       "report_default": {
        "hidden": true
       }
      }
     }
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "a=[0,1,2,3,4,0,1,2,3]\n",
    "b=np.array([10,20,30,40,50])\n",
    "print(b[a])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 4,
        "height": 4,
        "hidden": false,
        "row": 66,
        "width": 4
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "source": [
    "# Process/Plot assimilation run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "hidden": true
       },
       "report_default": {
        "hidden": true
       }
      }
     }
    }
   },
   "outputs": [],
   "source": [
    "print('BEGIN.')\n",
    "\n",
    "import math\n",
    "\n",
    "rad = 4.0*math.atan(1.0)/180.0\n",
    "#print('rad=',rad)\n",
    "\n",
    "os.chdir(rundir)\n",
    "\n",
    "grid_label='gr2'\n",
    "grid_label='gn'\n",
    "assim_files_string='CMIP6/CMIP/CSIRO/CAFE-1-0/historical/r1i1p2f1/Omon/'+dvar+'/'+grid_label+'/v20171025/'+dvar+'_Omon_historical_CAFE-1-0_r1i1p2f1_'+grid_label+'_??????-??????.nc'\n",
    "assim_files=glob.glob(assim_files_string)\n",
    "\n",
    "print('assim input_files ('+str(len(assim_files))+') =',assim_files)\n",
    "\n",
    "assim_input_file=assim_files[5]\n",
    "\n",
    "#print(datetime.datetime.fromtimestamp(os.stat(input_file).st_mtime))\n",
    "print('assim input file=',assim_input_file,'time=',datetime.datetime.fromtimestamp(os.stat(assim_input_file).st_mtime))\n",
    "\n",
    "assim_ifh0=netCDF4.Dataset(assim_input_file) #choose if more than 1.\n",
    "\n",
    "assim_time=assim_ifh0.variables['time']\n",
    "assim_lat=assim_ifh0.variables['latitude'][:,0]\n",
    "assim_lon=assim_ifh0.variables['longitude'][0,:]\n",
    "assim_clat=np.cos(assim_lat[:]*rad)\n",
    "\n",
    "#print('lat=',lat[122:152])\n",
    "#print('lon=',lon[110:160])\n",
    "#print('clat=',clat)\n",
    "#break\n",
    "\n",
    "assim_nino34_monthly=np.average(np.average(assim_ifh0.variables[dvar][:,122:152,110:160],axis=1,weights=assim_clat[122:152]),axis=1) #need to add in area weighting strictly\n",
    "\n",
    "'''\n",
    "print('time.units=',time.units)\n",
    "print('time.calendar=',time.calendar)\n",
    "print('time=',time)\n",
    "print('time[:]=',time[:])\n",
    "'''\n",
    "\n",
    "assim_date_time_stamp=netCDF4.num2date(assim_time[:],assim_time.units,assim_time.calendar)\n",
    "\n",
    "'''\n",
    "print('date_time_stamp=',date_time_stamp)\n",
    "num_stamp=netCDF4.date2num(date_time_stamp,time.units,time.calendar)\n",
    "print('num_stamp=',num_stamp)\n",
    "print('year_fraction=',year_fraction)\n",
    "'''\n",
    "assim_year_fraction_monthly=fractional_year_from_num2date(assim_date_time_stamp,assim_time.calendar)\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "plt.title('Monthly values')\n",
    "plt.plot(assim_year_fraction_monthly,assim_nino34_monthly[:])\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel('Ni$\\~{n}$o34 ($^o$C)')\n",
    "plt.show()\n",
    "print('END')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 8,
        "height": 4,
        "hidden": false,
        "row": 68,
        "width": 4
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "source": [
    "In the following monthly values will be calculated and then plotted."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 4,
        "hidden": false,
        "row": 81,
        "width": 12
       },
       "report_default": {
        "hidden": false
       }
      }
     }
    }
   },
   "source": [
    "From the monthly values now we can calculate long term monthly climatologies and the month-to-month anomalies, and then plot the time-series. These are generally of more interest than the full monthly values.\n",
    "\n",
    "Want to generalise climatology, anomaly generation so that any data product can be processed (e.g. (time), (time, lat, lon), (time, depth, lat, lon), (time, lev, lat, lon).\n",
    "\n",
    "Also want to pad out years with missing months."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "hidden": true
       },
       "report_default": {
        "hidden": true
       }
      }
     }
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print('BEGIN')\n",
    "\n",
    "ybeg=date_time_stamp[0].year\n",
    "yend=date_time_stamp[-1].year\n",
    "\n",
    "ydiff=yend-ybeg+1\n",
    "nmy=12\n",
    "\n",
    "MissingMonths=False\n",
    "\n",
    "first_month=date_time_stamp[0].month\n",
    "last_month=date_time_stamp[-1].month\n",
    "\n",
    "missing_months_beg,missing_months_end=0,0\n",
    "\n",
    "cyear_beg_skip,cyear_end_skip=0,1\n",
    "\n",
    "if(first_month!=1):\n",
    "  missing_months_beg=12-first_month\n",
    "  cyear_beg_skip=1\n",
    "  MissingMonths=True\n",
    "\n",
    "if(last_month!=12):\n",
    "  missing_months_end=12-last_month\n",
    "  cyear_end_skip=2\n",
    "  MissingMonths=True\n",
    "  \n",
    "if(MissingMonths):\n",
    "  print('There are missing months in the set. '+str(missing_months_beg)+' at beginning and '+str(missing_months_end)+' at end.')\n",
    "  print('Currently years with missing months are not used in generating long term monthly climatology.')\n",
    "  print('And missing months will be set to missing in the final time-series.')\n",
    "  ts_beg,ts_end,ts_avg,dt_beg,dt_end,dt_avg=get_timestamp_number(ybeg,yend,1,12,time.units,time.calendar)\n",
    "  year_fraction_monthly_full=fractional_year_from_num2date(ts_avg,time.calendar)\n",
    "  \n",
    "  nino34_monthly_full=ma.masked_all(ydiff*nmy,dtype=float) #ensure missing months are masked out.\n",
    "  last_month_index=ydiff*nmy-last_month\n",
    "  nino34_monthly_full[first_month-1:last_month_index]=nino34_monthly\n",
    "else:\n",
    "  print('All years have 12 months.')\n",
    "  \n",
    "  nino34_monthly_full=nino34_monthly\n",
    "  year_fraction_monthly_full=year_fraction_monthly\n",
    "  \n",
    "nino34_monthly_reshaped=np.reshape(nino34_monthly_full,new_monthly_array_shape(nino34_monthly.shape,ydiff,nmy)) #check this works on large multi-dimensional arrays.\n",
    "\n",
    "'''\n",
    "#print('nino34_monthly_reshaped.shape=',nino34_monthly_reshaped.shape)\n",
    "#print('cyear_beg,end_skip=',cyear_beg_skip,cyear_end_skip)\n",
    "#print('ydiff=',ydiff)\n",
    "#j=nino34_monthly_reshaped[cyear_beg_skip:-cyear_end_skip]\n",
    "#print('j=',j)\n",
    "'''\n",
    "\n",
    "climatology=np.average(nino34_monthly_reshaped[cyear_beg_skip:-cyear_end_skip],axis=0) #average over full years only, this could be an option as could use for all years present.\n",
    "nino34_monthly_climatology=np.expand_dims(climatology,0)\n",
    "nino34_monthly_climatology=np.tile(nino34_monthly_climatology,(ydiff,1))\n",
    "nino34_monthly_climatology_flat=nino34_monthly_climatology.flatten()\n",
    "nino34_monthly_anomaly=nino34_monthly_full-nino34_monthly_climatology_flat\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "zero=np.zeros(len(year_fraction_monthly_full))\n",
    "plt.plot(year_fraction_monthly_full,nino34_monthly_anomaly)\n",
    "plt.plot(year_fraction_monthly_full,zero)\n",
    "\n",
    "plt.title('Monthly anomalies')\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel('Ni$\\~{n}$o34 ($^o$C)')\n",
    "plt.show()\n",
    "\n",
    "#raise Exception('STOP!')\n",
    "print('END')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "hidden": true
       },
       "report_default": {
        "hidden": true
       }
      }
     }
    }
   },
   "outputs": [],
   "source": [
    "from datetime import date\n",
    "\n",
    "# random Person\n",
    "class Person:\n",
    "    def __init__(self, name, age):\n",
    "        self.name = name\n",
    "        self.age = age\n",
    "\n",
    "    @staticmethod\n",
    "    def fromFathersAge(name, fatherAge, fatherPersonAgeDiff):\n",
    "        return Person(name, date.today().year - fatherAge + fatherPersonAgeDiff)\n",
    "\n",
    "    @classmethod\n",
    "    def fromBirthYear(cls, name, birthYear):\n",
    "        return cls(name, date.today().year - birthYear)\n",
    "\n",
    "    def display(self):\n",
    "        print(self.name + \"'s age is: \" + str(self.age))\n",
    "\n",
    "class Man(Person):\n",
    "    sex = 'Male'\n",
    "\n",
    "man = Man.fromBirthYear('John', 1985)\n",
    "print(isinstance(man, Man))\n",
    "\n",
    "man1 = Man.fromFathersAge('John', 1965, 20)\n",
    "print(isinstance(man1, Man))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "man\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "man1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "man.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "man1.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Person:\n",
    "    TITLES = ('Dr', 'Mr', 'Mrs', 'Ms')\n",
    "\n",
    "    def __init__(self, name, surname):\n",
    "        self.name = name\n",
    "        self.surname = surname\n",
    "\n",
    "    def fullname(self): # instance method\n",
    "        # instance object accessible through self\n",
    "        return \"%s %s\" % (self.name, self.surname)\n",
    "\n",
    "    @classmethod\n",
    "    def allowed_titles_starting_with(cls, startswith): # class method\n",
    "        # class or instance object accessible through cls\n",
    "        return [t for t in cls.TITLES if t.startswith(startswith)]\n",
    "\n",
    "    @staticmethod\n",
    "    def allowed_titles_ending_with(endswith): # static method\n",
    "        # no parameter for class or instance object\n",
    "        # we have to use Person directly\n",
    "        return [t for t in Person.TITLES if t.endswith(endswith)]\n",
    "\n",
    "\n",
    "jane = Person(\"Jane\", \"Smith\")\n",
    "\n",
    "print(jane.fullname())\n",
    "\n",
    "print(jane.allowed_titles_starting_with(\"M\"))\n",
    "print(Person.allowed_titles_starting_with(\"M\"))\n",
    "\n",
    "print(jane.allowed_titles_ending_with(\"s\"))\n",
    "print(Person.allowed_titles_ending_with(\"s\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#http://python-textbok.readthedocs.io/en/1.0/Object_Oriented_Programming.html\n",
    "#following fails, possibly need to run python3.\n",
    "\n",
    "class Person:\n",
    "    def __init__(self, name, surname, number):\n",
    "        self.name = name\n",
    "        self.surname = surname\n",
    "        self.number = number\n",
    "\n",
    "\n",
    "class Student(Person):\n",
    "    UNDERGRADUATE, POSTGRADUATE = range(2)\n",
    "\n",
    "    def __init__(self, student_type, *args, **kwargs):\n",
    "        self.student_type = student_type\n",
    "        self.classes = []\n",
    "        super(Student, self).__init__(*args, **kwargs)\n",
    "\n",
    "    def enrol(self, course):\n",
    "        self.classes.append(course)\n",
    "\n",
    "\n",
    "class StaffMember(Person):\n",
    "    PERMANENT, TEMPORARY = range(2)\n",
    "\n",
    "    def __init__(self, employment_type, *args, **kwargs):\n",
    "        self.employment_type = employment_type\n",
    "        super(StaffMember, self).__init__(*args, **kwargs)\n",
    "\n",
    "\n",
    "class Lecturer(StaffMember):\n",
    "    def __init__(self, *args, **kwargs):\n",
    "        self.courses_taught = []\n",
    "        super(Lecturer, self).__init__(*args, **kwargs)\n",
    "\n",
    "    def assign_teaching(self, course):\n",
    "        self.courses_taught.append(course)\n",
    "\n",
    "\n",
    "jane = Student(Student.POSTGRADUATE, \"Jane\", \"Smith\", \"SMTJNX045\")\n",
    "jane.enrol(a_postgrad_course)\n",
    "\n",
    "bob = Lecturer(StaffMember.PERMANENT, \"Bob\", \"Jones\", \"123456789\")\n",
    "bob.assign_teaching(an_undergrad_course)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Song:\n",
    "\n",
    "    def __init__(self, title, artist, album, track_number):\n",
    "        self.title = title\n",
    "        self.artist = artist\n",
    "        self.album = album\n",
    "        self.track_number = track_number\n",
    "\n",
    "        artist.add_song(self)\n",
    "\n",
    "\n",
    "class Album:\n",
    "\n",
    "    def __init__(self, title, artist, year):\n",
    "        self.title = title\n",
    "        self.artist = artist\n",
    "        self.year = year\n",
    "\n",
    "        self.tracks = []\n",
    "\n",
    "        artist.add_album(self)\n",
    "\n",
    "    def add_track(self, title, artist=None):\n",
    "        if artist is None:\n",
    "            artist = self.artist\n",
    "\n",
    "        track_number = len(self.tracks)\n",
    "\n",
    "        song = Song(title, artist, self, track_number)\n",
    "\n",
    "        self.tracks.append(song)\n",
    "\n",
    "\n",
    "class Artist:\n",
    "    def __init__(self, name):\n",
    "        self.name = name\n",
    "\n",
    "        self.albums = []\n",
    "        self.songs = []\n",
    "\n",
    "    def add_album(self, album):\n",
    "        self.albums.append(album)\n",
    "\n",
    "    def add_song(self, song):\n",
    "        self.songs.append(song)\n",
    "\n",
    "\n",
    "class Playlist:\n",
    "    def __init__(self, name):\n",
    "        self.name = name\n",
    "        self.songs = []\n",
    "\n",
    "    def add_song(self, song):\n",
    "        self.songs.append(song)\n",
    "\n",
    "band = Artist(\"Bob's Awesome Band\")\n",
    "album = Album(\"Bob's First Single\", band, 2013)\n",
    "album.add_track(\"A Ballad about Cheese\")\n",
    "album.add_track(\"A Ballad about Cheese (dance remix)\")\n",
    "album.add_track(\"A Third Song to Use Up the Rest of the Space\")\n",
    "\n",
    "playlist = Playlist(\"My Favourite Songs\")\n",
    "\n",
    "for song in album.tracks:\n",
    "    playlist.add_song(song)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "album.artist\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime\n",
    "from dateutil import relativedelta\n",
    "\n",
    "DateTickInterval = 6\n",
    "date1 = datetime.date(2013, 1, 1)\n",
    "date2 = datetime.date(2015, 12, 31)\n",
    "r = relativedelta.relativedelta(date2, date1)\n",
    "NumMonths = r.years*12 + r.months\n",
    "dates = [ date1 + relativedelta.relativedelta(months=n) for n in range(0,NumMonths+1,DateTickInterval)]\n",
    "date_strings = [dt.strftime(\"%d/%m/%Y\") for dt in dates]\n",
    "print date_strings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "extensions": {
   "jupyter_dashboards": {
    "activeView": "report_default",
    "version": 1,
    "views": {
     "grid_default": {
      "cellMargin": 10,
      "defaultCellHeight": 20,
      "maxColumns": 12,
      "name": "grid",
      "type": "grid"
     },
     "report_default": {
      "name": "report",
      "type": "report"
     }
    }
   }
  },
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}