Ice_Crash_ACCESS-OM2-01.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CICE crash investigation in ACCESS-OM2_01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Author: Andrew Kiss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from glob import glob\n",
    "import os,sys\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import netCDF4 as nc\n",
    "from tqdm import tqdm_notebook\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "from calendar import month_abbr\n",
    "from netCDF4 import Dataset\n",
    "import cmocean as cm\n",
    "import string"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# figdir = 'figures/access-om2-01/Ice_Daily_ACCESS-OM2-01/'\n",
    "figdir = '/g/data/v45/aek156/figures/access-om2-01/Ice_Crash_ACCESS-OM2-01/'\n",
    "# figdir = 'scalewind'\n",
    "if not os.path.exists(figdir):\n",
    "    os.makedirs (figdir)\n",
    "\n",
    "def savefigure(fname):\n",
    "#     plt.savefig(os.path.join(figdir, fname),dpi=100, bbox_inches=\"tight\")  # comment out to disable saving\n",
    "    plt.savefig(os.path.join(figdir, fname),dpi=50, bbox_inches=\"tight\")  # comment out to disable saving\n",
    "    return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Investigate extreme winds in `/g/data1/ua8/JRA55-do/RYF/v1-3/RYF.u_10.1984_1985.nc`, `/g/data1/ua8/JRA55-do/RYF/v1-3/RYF.v_10.1984_1985.nc`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "ncFile = nc.Dataset('/g/data1/ua8/JRA55-do/RYF/v1-3/RYF.u_10.1984_1985.nc', mode='r')\n",
    "u10 = ncFile.variables['uas_10m'][...]\n",
    "\n",
    "ncFile = nc.Dataset('/g/data1/ua8/JRA55-do/RYF/v1-3/RYF.v_10.1984_1985.nc', mode='r')\n",
    "v10 = ncFile.variables['vas_10m'][...]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-45.070335, 40.71067, -43.399582, 40.83132]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[np.min(u10), np.max(u10), np.min(v10), np.max(v10)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These values are large but not preposterous."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# model data paths:\n",
    "basepath = '/g/data3/hh5/tmp/cosima/'\n",
    "model = 'access-om2-01'\n",
    "# expt = '01deg_jra55v13_ryf8485_spinup6'\n",
    "expt = '01deg_jra55v13_ryf8485_spinup6_newtopog'\n",
    "expt = '01deg_jra55v13_ryf8485_spinup6_newtopog_scalewind'\n",
    "expt = '01deg_jra55v13_ryf9091'\n",
    "# expt = '01deg_jra55v13_ryf8485_spinup6_dragio_0'\n",
    "# expt = '01deg_jra55v13_ryf8485_spinup6_tripole_bug'\n",
    "\n",
    "# model = 'access-om2-025'\n",
    "# expt = '025deg_jra55v13_ryf8485_KDS75'\n",
    "\n",
    "# model = 'access-om2'\n",
    "# expt = '1deg_jra55_ryf8485_gfdl50'\n",
    "# # expt = '1deg_jra55_ryf04_kds50'\n",
    "\n",
    "# # Paul's mom01v5 MOM-SIS runs\n",
    "# # /short/v45/pas561/mom/archive/kds75_wp2/output473\n",
    "# # it will be moved to /g/data3/hh5/tmp/cosima/mom01v5/kds75_wp2/ at some point.\n",
    "# model = 'mom01v5'\n",
    "# expt = 'kds75_wp2'\n",
    "\n",
    "DataDir = os.path.join(basepath, model)\n",
    "expdir = os.path.join(DataDir, expt)\n",
    "# date = '0004-07-05' # nicest view\n",
    "# date = '0004-07-12'\n",
    "# fn01 = '/g/data3/hh5/tmp/cosima/' + model + '/' +  expt + '/output023/ice/OUTPUT/iceh.' + date + '.nc'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get model grid data:\n",
    "gridFileList = glob(os.path.join(expdir, 'output*/ocean/ocean_grid.nc'))\n",
    "# gridFileList = glob('/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf8485_spinup6_dragio_0/output*/ocean/ocean_grid.nc')\n",
    "# gridFileList = ['/short/v45/pas561/mom/archive/kds75_wp2/output473/ocean_grid.nc']\n",
    "# # Paul's mom01v5 MOM-SIS runs\n",
    "# gridFileList = ['/g/data3/hh5/tmp/cosima/mom01v5/kds75_wp2/output473/ocean_grid.nc']\n",
    "\n",
    "\n",
    "# gridFileList = glob('/g/data3/hh5/tmp/cosima/access-om2-025/025deg_jra55v13_ryf8485_KDS75/output017/ocean/ocean_grid.nc')\n",
    "\n",
    "gridFileList.sort()\n",
    "ncFile = nc.Dataset(gridFileList[0], mode='r')\n",
    "xt_ocean = ncFile.variables['xt_ocean'][...]\n",
    "yt_ocean = ncFile.variables['yt_ocean'][...]\n",
    "lon_t = ncFile.variables['geolon_t'][...]\n",
    "lat_t = ncFile.variables['geolat_t'][...]\n",
    "area_t = ncFile.variables['area_t'][...]\n",
    "# NLAT_half = int(np.shape(area_t)[0]/2)\n",
    "ht = ncFile.variables['ht'][...]\n",
    "land_mask = np.copy(ht)\n",
    "land_mask[np.where(ht>30)] = 0\n",
    "land_mask[np.where(ht<=30)] = 1\n",
    "ncFile.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# icehFileList = glob(os.path.join(expdir, 'output*/ice/OUTPUT/iceh.0017-08*.nc'))\n",
    "# icehFileList = glob(os.path.join(expdir, 'output173/ice/OUTPUT/iceh.0017-08*.nc'))\n",
    "# icehFileList = glob(os.path.join(expdir, 'output*/ice/OUTPUT/iceh.0026-*-*.nc'))\n",
    "\n",
    "# icehFileList = glob(os.path.join(expdir, 'output*/ice/OUTPUT/iceh.0027-12-*.nc'))\n",
    "\n",
    "# icehFileList = glob(os.path.join(expdir, 'output*/ice/OUTPUT/iceh.0029-01-0*.nc'))\n",
    "\n",
    "# icehFileList = glob(os.path.join(expdir, 'output429/ice/OUTPUT/iceh.1938-08-[123]*.nc'))\n",
    "\n",
    "# icehFileList = glob(os.path.join(expdir, 'output*/ice/OUTPUT/iceh.1939-08-12.nc'))\n",
    "\n",
    "icehFileList = glob(os.path.join(expdir, 'output*/ice/OUTPUT/iceh.193?-0[78].nc'))\n",
    "\n",
    "\n",
    "\n",
    "# icehFileList = glob(os.path.join(expdir, 'output*/ice/OUTPUT/iceh.0026-12-3*.nc'))\n",
    "\n",
    "# icehFileList = glob(os.path.join(expdir, 'output*/ice/OUTPUT/iceh.0027-01-0*.nc'))\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# icehFileList = glob(os.path.join(expdir, 'output173//ocean/ocean_month.nc'))\n",
    "\n",
    "\n",
    "# icehFileList = glob(os.path.join(expdir, 'output*/ocean/o2i.nc'))\n",
    "\n",
    "# icehFileList = glob(os.path.join(expdir, 'output*/ocean/ocean_month.nc'))\n",
    "\n",
    "# icehFileList = glob(os.path.join(expdir, 'output*/ice/OUTPUT/iceh.????-07.nc'))\n",
    "# icehFileList = glob('/g/data3/hh5/tmp/cosima/access-om2-025/025deg_jra55v13_ryf8485_KDS75/output*/ice/OUTPUT/iceh*07.nc')\n",
    "\n",
    "icehFileList.sort()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091/output118/ice/OUTPUT/iceh.1930-07.nc',\n",
       " '/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091/output118/ice/OUTPUT/iceh.1930-08.nc',\n",
       " '/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091/output122/ice/OUTPUT/iceh.1931-07.nc',\n",
       " '/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091/output122/ice/OUTPUT/iceh.1931-08.nc',\n",
       " '/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091/output126/ice/OUTPUT/iceh.1932-07.nc',\n",
       " '/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091/output126/ice/OUTPUT/iceh.1932-08.nc',\n",
       " '/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091/output130/ice/OUTPUT/iceh.1933-07.nc',\n",
       " '/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091/output130/ice/OUTPUT/iceh.1933-08.nc',\n",
       " '/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091/output134/ice/OUTPUT/iceh.1934-07.nc',\n",
       " '/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091/output134/ice/OUTPUT/iceh.1934-08.nc',\n",
       " '/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091/output138/ice/OUTPUT/iceh.1935-07.nc',\n",
       " '/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091/output138/ice/OUTPUT/iceh.1935-08.nc',\n",
       " '/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091/output142/ice/OUTPUT/iceh.1936-07.nc',\n",
       " '/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091/output142/ice/OUTPUT/iceh.1936-08.nc',\n",
       " '/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091/output146/ice/OUTPUT/iceh.1937-07.nc',\n",
       " '/g/data3/hh5/tmp/cosima/access-om2-01/01deg_jra55v13_ryf9091/output146/ice/OUTPUT/iceh.1937-08.nc']"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "icehFileList"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='variables', max=1, style=ProgressStyle(description_width='ini…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='day', max=16, style=ProgressStyle(description_width='initial'…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Skipping 01deg_jra55v13_ryf9091_tripole_bug_shear_1930-07.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_Arctic_closeup_shear_1930-07.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_tripole_bug_shear_1930-08.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_Arctic_closeup_shear_1930-08.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_tripole_bug_shear_1931-07.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_Arctic_closeup_shear_1931-07.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_tripole_bug_shear_1931-08.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_Arctic_closeup_shear_1931-08.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_tripole_bug_shear_1932-07.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_Arctic_closeup_shear_1932-07.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_tripole_bug_shear_1932-08.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_Arctic_closeup_shear_1932-08.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_tripole_bug_shear_1933-07.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_Arctic_closeup_shear_1933-07.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_tripole_bug_shear_1933-08.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_Arctic_closeup_shear_1933-08.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_tripole_bug_shear_1934-07.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_Arctic_closeup_shear_1934-07.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_tripole_bug_shear_1934-08.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_Arctic_closeup_shear_1934-08.png (file exists)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-19.04/lib/python3.6/site-packages/ipykernel_launcher.py:123: RuntimeWarning: divide by zero encountered in log10\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='views', max=2, style=ProgressStyle(description_width='initial…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='views', max=2, style=ProgressStyle(description_width='initial…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Skipping 01deg_jra55v13_ryf9091_tripole_bug_shear_1936-07.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_Arctic_closeup_shear_1936-07.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_tripole_bug_shear_1936-08.png (file exists)\n",
      "Skipping 01deg_jra55v13_ryf9091_Arctic_closeup_shear_1936-08.png (file exists)\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='views', max=2, style=ProgressStyle(description_width='initial…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='views', max=2, style=ProgressStyle(description_width='initial…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r"
     ]
    }
   ],
   "source": [
    "varnames = ['hi', 'aice', 'uvel', 'vvel', 'divu', 'shear', 'sig1', 'sig2', \n",
    "            'strairx', 'strairy', 'uocn', 'vocn', 'sst', 'sss', 'icespeed', 'oceanspeed',\n",
    "            'strocnx', 'strocny']\n",
    "\n",
    "varnames = ['hi', 'aice','uvel', 'vvel', 'shear', 'uocn', 'vocn', 'strairx', 'strairy', 'divu']\n",
    "\n",
    "# varnames = ['uvel', 'vvel', 'strairx', 'strairy']\n",
    "\n",
    "varnames = ['aice']\n",
    "\n",
    "# varnames = ['uvel']\n",
    "\n",
    "varnames = ['strocnx', 'strocny']\n",
    "\n",
    "# varnames = ['shear','uocn','uvel']\n",
    "\n",
    "# varnames = ['ssu_i', 'ssv_i'] # from ocean/o2i.nc\n",
    "\n",
    "# varnames = ['tau_x',  'tau_y'] # from ocean/ocean_month.nc\n",
    "\n",
    "\n",
    "# # for Paul's mom01v5 MOM-SIS runs\n",
    "# icehFileList = glob(os.path.join(expdir, 'output473/ice__*.nc'))\n",
    "# icehFileList.sort()\n",
    "# varnames = ['UI', 'VI', 'SIGI', 'SIGII']\n",
    "# # varnames = ['UI', 'VI', 'shear_stress']\n",
    "\n",
    "# varnames = ['UI', 'VI']\n",
    "\n",
    "varnames = ['oceanspeed']\n",
    "\n",
    "varnames = ['shear', 'hi']\n",
    "\n",
    "varnames = ['strairy', 'strairx', 'aice', 'uvel', 'vvel']\n",
    "\n",
    "varnames = ['strairy', 'strairx']\n",
    "\n",
    "varnames = ['shear']\n",
    "\n",
    "\n",
    "views = ['Aug_crash', 'May_crash', 'tripole_bug', 'Arctic_closeup', 'Arctic', 'Antarctic', 'EAC', 'Agulhas', 'Kuroshio', 'Gulf_Stream', 'ACC']\n",
    "\n",
    "views = ['EAC', 'Agulhas', 'Kuroshio', 'Gulf_Stream', 'ACC']\n",
    "\n",
    "views = ['Arctic', 'Antarctic']\n",
    "\n",
    "views = ['Aug_crash','Arctic']\n",
    "\n",
    "views = ['Aug_crash']\n",
    "\n",
    "views = ['tripole_bug', 'Arctic_closeup']\n",
    "\n",
    "\n",
    "\n",
    "# views = ['ACC']\n",
    "\n",
    "for varname in tqdm_notebook(varnames, leave=False, desc='variables'):\n",
    "    fnum = 0  # enumerate doesn't work with tqdm_notebook...\n",
    "    for fpath in tqdm_notebook(icehFileList, leave=False, desc='day'):\n",
    "#             print(fpath)\n",
    "        fnum += fnum + 1\n",
    "        fname = os.path.basename(fpath)\n",
    "        try:\n",
    "            date = fname.split('.')[1]\n",
    "            month = int(date.split('-')[1])\n",
    "        except:\n",
    "            date = str(fnum)\n",
    "            month = ''\n",
    "            \n",
    "        viewstodo = []\n",
    "        outnames = []\n",
    "        for view in views:\n",
    "            outname = expt + '_' + view + '_' + varname + '_' + date + '.png'\n",
    "            if os.path.exists(os.path.join(figdir, outname)):\n",
    "                print('Skipping ' + outname + ' (file exists)')\n",
    "            else:\n",
    "                viewstodo.append(view)\n",
    "                outnames.append(outname)\n",
    "        if len(viewstodo)>0:\n",
    "            try:\n",
    "                del var\n",
    "            except:\n",
    "                pass\n",
    "            try:\n",
    "                del varmask\n",
    "            except:\n",
    "                pass\n",
    "            ncFile = Dataset(fpath, mode='r')\n",
    "            if varname.endswith('speed'):\n",
    "                if varname == 'icespeed':\n",
    "                    try:\n",
    "                        u = ncFile.variables['uvel'][:][0,:,:]\n",
    "                        v = ncFile.variables['vvel'][:][0,:,:]\n",
    "                    except:\n",
    "                        u = ncFile.variables['uvel_m'][:][0,:,:]\n",
    "                        v = ncFile.variables['vvel_m'][:][0,:,:]\n",
    "                elif varname == 'oceanspeed':\n",
    "                    try:\n",
    "                        u = ncFile.variables['uocn'][:][0,:,:]\n",
    "                        v = ncFile.variables['vocn'][:][0,:,:]\n",
    "                    except:\n",
    "                        u = ncFile.variables['uocn_m'][:][0,:,:]\n",
    "                        v = ncFile.variables['vocn_m'][:][0,:,:]\n",
    "                else:\n",
    "                    print('Error: ' + varname + ' unknown')\n",
    "                    raise \n",
    "                var = np.ma.sqrt(u**2 + v**2)  # u,v are co-located on B grid\n",
    "#                 elif varname=='shear_stress':\n",
    "#                     # NB: this is the maximum shearing stress, not the shearing strain\n",
    "#                     # https://en.wikipedia.org/wiki/Cauchy_stress_tensor#Maximum_and_minimum_shear_stresses\n",
    "#                     sigi  = ncFile.variables['SIGI' ][:][0,:,:]\n",
    "#                     sigii = ncFile.variables['SIGII'][:][0,:,:]\n",
    "#                     var = abs(sigi-sigii)/2 # BUG: INCORRECT! SIGI, SIGII are invariants, not stresses - see notebook p588\n",
    "            else:\n",
    "                try:\n",
    "                    var = ncFile.variables[varname][:][0,:,:]\n",
    "                except:\n",
    "                    var = ncFile.variables[varname + '_m'][:][0,:,:]\n",
    "            ncFile.close()\n",
    "\n",
    "            varmask = np.ma.getmask(var)\n",
    "            if varname.startswith('shear'):\n",
    "                var = np.where(var>0, np.log10(var), -2)\n",
    "\n",
    "#             varmask = np.ma.getmask(var)\n",
    "#             np.ma.masked_where(varmask,var)\n",
    "\n",
    "            var = np.ma.masked_where(varmask,var)\n",
    "\n",
    "\n",
    "       \n",
    "            vnum = -1  # enumerate doesn't work with tqdm_notebook...\n",
    "            for view in tqdm_notebook(viewstodo, leave=False, desc='views'):\n",
    "                vnum += 1\n",
    "                outname = outnames[vnum]\n",
    "#                 print('doing ' + outname)\n",
    "                fact = 1  # scaling factor for var\n",
    "                parlabels = [False,False,False,False]\n",
    "\n",
    "                if varname == 'aice':\n",
    "                    desc = 'Ice concentration'\n",
    "                    levels = np.arange(0,1.01,.01)\n",
    "                    ticks = [0,.2,.4,.6,.8,1]\n",
    "                    ext = 'neither'\n",
    "                    cmp = cm.cm.ice\n",
    "                elif varname == 'hi':\n",
    "                    desc = 'Ice thickness (m)'\n",
    "                    levels = np.arange(0,6.01,.06)\n",
    "                    ticks = range(7)\n",
    "                    ext = 'max'\n",
    "                    cmp = cm.cm.ice\n",
    "                elif varname in ['uvel', 'UI']:\n",
    "                    desc = 'Ice x velocity (m/s)'\n",
    "                    levels = np.arange(-1,1.01,.01)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.2)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname in ['vvel', 'VI']:\n",
    "                    desc = 'Ice y velocity (m/s)'\n",
    "                    levels = np.arange(-1,1.01,.01)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.2)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname == 'uocn':\n",
    "                    desc = 'Ocean x velocity (m/s)'\n",
    "                    levels = np.arange(-1,1.01,.01)\n",
    "    #                     levels = np.arange(-.25,.2501,.01)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.2)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname == 'vocn':\n",
    "                    desc = 'Ocean v velocity (m/s)'\n",
    "                    levels = np.arange(-1,1.01,.01)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.2)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname == 'ssu_i':\n",
    "                    desc = 'ssu_i (m/s)'\n",
    "                    levels = np.arange(-1,1.01,.01)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.2)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname == 'ssv_i':\n",
    "                    desc = 'ssv_i (m/s)'\n",
    "                    levels = np.arange(-1,1.01,.01)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.2)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname == 'icespeed':\n",
    "                    desc = 'Ice speed (m/s)'\n",
    "                    levels = np.arange(0,1.01,.01)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.2)\n",
    "                    ext = 'max'\n",
    "                    cmp = cm.cm.thermal\n",
    "                elif varname == 'oceanspeed':\n",
    "                    desc = 'Ocean surface speed (m/s)'\n",
    "                    levels = np.arange(0,1.01,.01)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.2)\n",
    "                    ext = 'max'\n",
    "                    cmp = cm.cm.thermal\n",
    "                elif varname == 'divu':\n",
    "                    desc = 'Ice divergence rate (%/day)'\n",
    "                    levels = np.arange(-1000,1000.01,20)\n",
    "                    ticks = np.rint(np.arange(levels[0],levels[-1]+.001,250.0))\n",
    "    #                     var = np.copysign(np.log10(np.abs(var)), var)\n",
    "    #                     desc = 'Log of ice divergence rate (%/day)'\n",
    "    #                     levels = np.arange(-3,3.01,.05)\n",
    "    #                     ticks = np.rint(np.arange(levels[0],levels[-1],1.0))\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname.startswith('shear'):\n",
    "    # #                     desc = 'Ice strain rate (shear) (%/day)'\n",
    "    # #                     levels = np.arange(0,300.01,10)\n",
    "    # #                     ticks = np.rint(np.arange(levels[0],levels[-1],30.0))\n",
    "    # #                     ext = 'max'\n",
    "    # #                     noice = (var <= 0).astype('int')\n",
    "    #                     varcopy = var\n",
    "    # # #                     var = np.copysign(np.log10(np.abs(var)), var)\n",
    "    # #                     var = np.log10(var)\n",
    "    # #                     var = np.where(varcopy<=0, -2, var)\n",
    "    # #                     var = np.where(varcopy==np.nan, np.nan, var)\n",
    "#                     if vnum==0:\n",
    "#                         var = np.where(var>0, np.log10(var), -2)\n",
    "#     #                     var = np.ma.masked_where(np.ma.getmask(varcopy),var)\n",
    "                    desc = 'Log of ice shear rate (%/day)'\n",
    "                    levels = np.arange(-2,3.01,.05)\n",
    "                    ticks = np.rint(np.arange(levels[0],levels[-1]+.001,1.0))\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.thermal\n",
    "                elif varname in ['strairx', 'tau_x']:\n",
    "                    desc = 'Wind stress, x direction (Pa)'\n",
    "#                     levels = np.arange(-0.1,0.101,.002)\n",
    "#                     ticks = np.arange(levels[0],levels[-1]+.001,0.1)\n",
    "                    levels = np.arange(-0.5,0.501,.02)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.01,0.1)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname in ['strairy', 'tau_y']:\n",
    "                    desc = 'Wind stress, y direction (Pa)'\n",
    "#                     levels = np.arange(-0.1,0.101,.002)\n",
    "#                     ticks = np.arange(levels[0],levels[-1]+.001,0.1)\n",
    "                    levels = np.arange(-0.5,0.501,.02)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.01,0.1)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname == 'strocnx':\n",
    "                    desc = '- Ocean-ice stress, x direction (Pa)'\n",
    "                    levels = np.arange(-0.1,0.101,.002)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.1)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                    fact = -1\n",
    "                elif varname == 'strocny':\n",
    "                    desc = '- Ocean-ice stress, y direction (Pa)'\n",
    "                    levels = np.arange(-0.1,0.101,.002)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.1)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                    fact = -1\n",
    "                elif varname == 'sst':\n",
    "                    desc = 'SST (C)'\n",
    "                    levels = np.arange(-3,10.1,.25)\n",
    "    #                     ticks = np.arange(levels[0],levels[-1]+.001,2)\n",
    "                    ticks = np.arange(-4,100,2)\n",
    "                    ext = 'both'\n",
    "    #                     cmp = plt.cm.gist_ncar\n",
    "                    cmp = cm.cm.thermal\n",
    "                elif varname.startswith('SIG'):\n",
    "                    desc = 'Stress invariant ' + varname + ' (units??)'\n",
    "                    levels = np.linspace(-1e6,1e6,50)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+1e4,1e5)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                else:\n",
    "                    desc = '***UNKNOWN***'\n",
    "                    levels = np.arange(0,1.01,.1)\n",
    "                    ticks = [0,.2,.4,.6,.8,1]\n",
    "                    ext = 'neither'\n",
    "                    cmp = cm.cm.balance\n",
    "\n",
    "\n",
    "\n",
    "#                 var = \n",
    "#                 np.ma.masked_where(varmask,var)\n",
    "\n",
    "\n",
    "                font = {'size':24}\n",
    "                tick_font=24\n",
    "\n",
    "                fig = plt.figure(figsize=(30,30))\n",
    "                if view == 'Arctic':\n",
    "                    m = Basemap(projection ='npstere',boundinglat=45,lon_0=-10,resolution='l')\n",
    "                elif view == 'Arctic_closeup':\n",
    "                    m = Basemap(projection ='npstere',boundinglat=70,lon_0=-10,resolution='l')\n",
    "                elif view == 'Aug_crash':\n",
    "                    m = Basemap(projection='stere',llcrnrlat=67,urcrnrlat=75,\n",
    "                                llcrnrlon=-110,urcrnrlon=-90,lat_0=75,lon_0=-100,resolution='l')\n",
    "                elif view == 'May_crash':\n",
    "                    m = Basemap(projection='stere',llcrnrlat=66,urcrnrlat=85,\n",
    "                                llcrnrlon=70,urcrnrlon=110,lat_0=75,lon_0=80,resolution='l')\n",
    "                elif view == 'tripole_bug':\n",
    "                    m = Basemap(projection='stere',llcrnrlat=72,urcrnrlat=85,\n",
    "                                llcrnrlon=70,urcrnrlon=110,lat_0=75,lon_0=80,resolution='l')\n",
    "                elif view == 'EAC':\n",
    "                    m = Basemap(projection='cyl',\n",
    "                                llcrnrlon=142-360,urcrnrlon=175-360,\n",
    "                                llcrnrlat=-48,urcrnrlat=-20,\n",
    "                                resolution='l')\n",
    "                    parlabels = [True,False,False,True]\n",
    "                elif view == 'Agulhas':\n",
    "                    m = Basemap(projection='cyl',\n",
    "                                llcrnrlon=0,urcrnrlon=50,\n",
    "                                llcrnrlat=-48,urcrnrlat=-15,\n",
    "                                resolution='l')\n",
    "                    parlabels = [True,False,False,True]\n",
    "                elif view == 'Kuroshio':\n",
    "                    m = Basemap(projection='cyl',\n",
    "                                llcrnrlon=125-360,urcrnrlon=165-360,\n",
    "                                llcrnrlat=25,urcrnrlat=50,\n",
    "                                resolution='l')\n",
    "                    parlabels = [True,False,False,True]\n",
    "                elif view == 'Gulf_Stream':\n",
    "                    m = Basemap(projection='cyl',\n",
    "                                llcrnrlon=-98,urcrnrlon=-50,\n",
    "                                llcrnrlat=15,urcrnrlat=45,\n",
    "                                resolution='l')\n",
    "                    parlabels = [True,False,False,True]\n",
    "                elif view == 'ACC':\n",
    "                    m = Basemap(projection ='spstere',boundinglat=-30,lon_0=170,resolution='l')\n",
    "                    parlabels = [False,False,False,False]\n",
    "                else:\n",
    "                    m = Basemap(projection ='spstere',boundinglat=-50,lon_0=170,resolution='l')\n",
    "\n",
    "                # m = Basemap(projection ='ortho',resolution='l',lat_0=90,lon_0=45,\n",
    "                #            llcrnrx=-m.urcrnrx/2.,llcrnry=-m.urcrnry/2.,urcrnrx=m.urcrnrx/2.,urcrnry=m.urcrnry/2.)\n",
    "                # #             llcrnrlat=60,urcrnrlat=60,\n",
    "                # #             llcrnrlon=-45,urcrnrlon=45)\n",
    "                # #             llcrnrlon=loncorners[0],urcrnrlon=loncorners[2]) ,boundinglat=48\n",
    "                x,y = m(*(lon_t,lat_t))\n",
    "                ctr = m.contourf(x,y,fact*var,levels=levels,cmap=cmp, extend=ext) #, vmin=0, vmax=1.1)\n",
    "    #                 ctr.cmap.set_under(color=cmp(np.min(levels)), alpha=None)\n",
    "    #                 ctr.cmap.set_over(color=cmp(np.max(levels)), alpha=None)\n",
    "                ctr.cmap.set_under(color=cmp(0), alpha=None)\n",
    "                ctr.cmap.set_over(color=cmp(.999), alpha=None)\n",
    "    #                 if view.startswith('Arctic'):\n",
    "    #                     xobs,yobs = m(*(obs_lon_NH,obs_lat_NH))\n",
    "    #                     m.contour(xobs,yobs,CN_obs_NH[month-1],[0.3],colors='r')\n",
    "    #                 else:\n",
    "    #                     xobs,yobs = m(*(obs_lon_SH,obs_lat_SH))\n",
    "    #                     m.contour(xobs,yobs,CN_obs_SH[month-1],[0.3],colors='r')\n",
    "\n",
    "                plt.title(desc + ', ' + date,font)\n",
    "                cbar = m.colorbar(ctr, location = 'bottom', pad = \"6%\")\n",
    "                cbar.set_label(desc,size=tick_font)\n",
    "                cbar.set_ticks(ticks)\n",
    "                cbar_labels=plt.getp(cbar.ax.axes,'xticklabels')\n",
    "                plt.setp(cbar_labels,fontsize=tick_font)\n",
    "                m.drawmapboundary(fill_color='gray') # background color - for non-ocean areas\n",
    "                # fill continents, set lake color same as ocean color.\n",
    "                # m.fillcontinents(color='gray',lake_color='gray')\n",
    "                # draw parallels and meridians.\n",
    "                # label parallels on right and top\n",
    "                # meridians on bottom and left\n",
    "                parallels = np.arange(-80.,81,10.)\n",
    "                # labels = [left,right,top,bottom]\n",
    "    #                 m.drawparallels(parallels,color='white')#,labels=[False,True,True,False],size=tick_font)\n",
    "                m.drawparallels(parallels,color='white',labels=parlabels,size=tick_font)\n",
    "\n",
    "    #                 meridians = np.arange(0.,351.,20.)\n",
    "                meridians = np.arange(0.,351.,10.)\n",
    "                m.drawmeridians(meridians,labels=[True,False,False,True],size=tick_font,color='white')\n",
    "                savefigure(outname)\n",
    "                plt.close()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "varnames = ['hi', 'aice', 'uvel', 'vvel', 'divu', 'shear', 'sig1', 'sig2', \n",
    "            'strairx', 'strairy', 'uocn', 'vocn', 'sst', 'sss', 'icespeed', 'oceanspeed',\n",
    "            'strocnx', 'strocny']\n",
    "\n",
    "varnames = ['hi', 'aice','uvel', 'vvel', 'shear', 'uocn', 'vocn', 'strairx', 'strairy', 'divu']\n",
    "\n",
    "# varnames = ['uvel', 'vvel', 'strairx', 'strairy']\n",
    "\n",
    "varnames = ['aice']\n",
    "\n",
    "# varnames = ['uvel']\n",
    "\n",
    "varnames = ['strocnx', 'strocny']\n",
    "\n",
    "# varnames = ['shear','uocn','uvel']\n",
    "\n",
    "# varnames = ['ssu_i', 'ssv_i'] # from ocean/o2i.nc\n",
    "\n",
    "# varnames = ['tau_x',  'tau_y'] # from ocean/ocean_month.nc\n",
    "\n",
    "\n",
    "# # for Paul's mom01v5 MOM-SIS runs\n",
    "# icehFileList = glob(os.path.join(expdir, 'output473/ice__*.nc'))\n",
    "# icehFileList.sort()\n",
    "# varnames = ['UI', 'VI', 'SIGI', 'SIGII']\n",
    "# # varnames = ['UI', 'VI', 'shear_stress']\n",
    "\n",
    "# varnames = ['UI', 'VI']\n",
    "\n",
    "varnames = ['oceanspeed']\n",
    "\n",
    "\n",
    "views = ['Aug_crash', 'May_crash', 'tripole_bug', 'Arctic_closeup', 'Arctic', 'Antarctic', 'EAC', 'Agulhas', 'Kuroshio', 'Gulf_Stream', 'ACC']\n",
    "\n",
    "views = ['tripole_bug', 'Arctic_closeup']\n",
    "\n",
    "views = ['EAC', 'Agulhas', 'Kuroshio', 'Gulf_Stream', 'ACC']\n",
    "\n",
    "# views = ['ACC']\n",
    "\n",
    "\n",
    "# TODO: swap loop order to reuse data for multiple views (need to change skipping logic)\n",
    "\n",
    "for varname in tqdm_notebook(varnames, leave=False, desc='variables'):\n",
    "    for view in tqdm_notebook(views, leave=False, desc='views'):\n",
    "        fnum = 0  # enumerate doesn't work with tqdm_notebook...\n",
    "        for fpath in tqdm_notebook(icehFileList, leave=False, desc='day'):\n",
    "#             print(fpath)\n",
    "            fnum = fnum + 1\n",
    "            fname = os.path.basename(fpath)\n",
    "            try:\n",
    "                date = fname.split('.')[1]\n",
    "                month = int(date.split('-')[1])\n",
    "            except:\n",
    "                date = str(fnum)\n",
    "                month = ''\n",
    "            outname = expt + '_' + view + '_' + varname + '_' + date + '.png'\n",
    "\n",
    "#             if False:\n",
    "            if os.path.exists(os.path.join(figdir, outname)):\n",
    "                print('Skipping ' + outname + ' (file exists)')\n",
    "            else:\n",
    "                try:\n",
    "                    del var\n",
    "                except:\n",
    "                    pass\n",
    "                try:\n",
    "                    del varmask\n",
    "                except:\n",
    "                    pass\n",
    "                ncFile = Dataset(fpath, mode='r')\n",
    "                if varname.endswith('speed'):\n",
    "                    if varname == 'icespeed':\n",
    "                        try:\n",
    "                            u = ncFile.variables['uvel'][:][0,:,:]\n",
    "                            v = ncFile.variables['vvel'][:][0,:,:]\n",
    "                        except:\n",
    "                            u = ncFile.variables['uvel_m'][:][0,:,:]\n",
    "                            v = ncFile.variables['vvel_m'][:][0,:,:]\n",
    "                    elif varname == 'oceanspeed':\n",
    "                        try:\n",
    "                            u = ncFile.variables['uocn'][:][0,:,:]\n",
    "                            v = ncFile.variables['vocn'][:][0,:,:]\n",
    "                        except:\n",
    "                            u = ncFile.variables['uocn_m'][:][0,:,:]\n",
    "                            v = ncFile.variables['vocn_m'][:][0,:,:]\n",
    "                    else:\n",
    "                        print('Error: ' + varname + ' unknown')\n",
    "                        raise \n",
    "                    var = np.ma.sqrt(u**2 + v**2)  # u,v are co-located on B grid\n",
    "#                 elif varname=='shear_stress':\n",
    "#                     # NB: this is the maximum shearing stress, not the shearing strain\n",
    "#                     # https://en.wikipedia.org/wiki/Cauchy_stress_tensor#Maximum_and_minimum_shear_stresses\n",
    "#                     sigi  = ncFile.variables['SIGI' ][:][0,:,:]\n",
    "#                     sigii = ncFile.variables['SIGII'][:][0,:,:]\n",
    "#                     var = abs(sigi-sigii)/2 # BUG: INCORRECT! SIGI, SIGII are invariants, not stresses - see notebook p588\n",
    "                else:\n",
    "                    try:\n",
    "                        var = ncFile.variables[varname][:][0,:,:]\n",
    "                    except:\n",
    "                        var = ncFile.variables[varname + '_m'][:][0,:,:]\n",
    "                varmask = np.ma.getmask(var)\n",
    "                ncFile.close()\n",
    "                fact = 1\n",
    "                parlabels = [False,False,False,False]\n",
    "\n",
    "                if varname == 'aice':\n",
    "                    desc = 'Ice concentration'\n",
    "                    levels = np.arange(0,1.01,.01)\n",
    "                    ticks = [0,.2,.4,.6,.8,1]\n",
    "                    ext = 'neither'\n",
    "                    cmp = cm.cm.ice\n",
    "                elif varname == 'hi':\n",
    "                    desc = 'Ice thickness (m)'\n",
    "                    levels = np.arange(0,6.01,.06)\n",
    "                    ticks = range(7)\n",
    "                    ext = 'max'\n",
    "                    cmp = cm.cm.ice\n",
    "                elif varname in ['uvel', 'UI']:\n",
    "                    desc = 'Ice x velocity (m/s)'\n",
    "                    levels = np.arange(-1,1.01,.01)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.2)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname in ['vvel', 'VI']:\n",
    "                    desc = 'Ice y velocity (m/s)'\n",
    "                    levels = np.arange(-1,1.01,.01)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.2)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname == 'uocn':\n",
    "                    desc = 'Ocean x velocity (m/s)'\n",
    "                    levels = np.arange(-1,1.01,.01)\n",
    "#                     levels = np.arange(-.25,.2501,.01)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.2)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname == 'vocn':\n",
    "                    desc = 'Ocean v velocity (m/s)'\n",
    "                    levels = np.arange(-1,1.01,.01)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.2)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname == 'ssu_i':\n",
    "                    desc = 'ssu_i (m/s)'\n",
    "                    levels = np.arange(-1,1.01,.01)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.2)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname == 'ssv_i':\n",
    "                    desc = 'ssv_i (m/s)'\n",
    "                    levels = np.arange(-1,1.01,.01)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.2)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname == 'icespeed':\n",
    "                    desc = 'Ice speed (m/s)'\n",
    "                    levels = np.arange(0,1.01,.01)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.2)\n",
    "                    ext = 'max'\n",
    "                    cmp = cm.cm.thermal\n",
    "                elif varname == 'oceanspeed':\n",
    "                    desc = 'Ocean surface speed (m/s)'\n",
    "                    levels = np.arange(0,1.01,.01)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.2)\n",
    "                    ext = 'max'\n",
    "                    cmp = cm.cm.thermal\n",
    "                elif varname == 'divu':\n",
    "                    desc = 'Ice divergence rate (%/day)'\n",
    "                    levels = np.arange(-1000,1000.01,20)\n",
    "                    ticks = np.rint(np.arange(levels[0],levels[-1]+.001,250.0))\n",
    "#                     var = np.copysign(np.log10(np.abs(var)), var)\n",
    "#                     desc = 'Log of ice divergence rate (%/day)'\n",
    "#                     levels = np.arange(-3,3.01,.05)\n",
    "#                     ticks = np.rint(np.arange(levels[0],levels[-1],1.0))\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname.startswith('shear'):\n",
    "# #                     desc = 'Ice strain rate (shear) (%/day)'\n",
    "# #                     levels = np.arange(0,300.01,10)\n",
    "# #                     ticks = np.rint(np.arange(levels[0],levels[-1],30.0))\n",
    "# #                     ext = 'max'\n",
    "# #                     noice = (var <= 0).astype('int')\n",
    "#                     varcopy = var\n",
    "# # #                     var = np.copysign(np.log10(np.abs(var)), var)\n",
    "# #                     var = np.log10(var)\n",
    "# #                     var = np.where(varcopy<=0, -2, var)\n",
    "# #                     var = np.where(varcopy==np.nan, np.nan, var)\n",
    "                    var = np.where(var>0, np.log10(var), -2)\n",
    "#                     var = np.ma.masked_where(np.ma.getmask(varcopy),var)\n",
    "                    desc = 'Log of ice shear rate (%/day)'\n",
    "                    levels = np.arange(-2,3.01,.05)\n",
    "                    ticks = np.rint(np.arange(levels[0],levels[-1]+.001,1.0))\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.thermal\n",
    "                elif varname in ['strairx', 'tau_x']:\n",
    "                    desc = 'Wind stress, x direction (Pa)'\n",
    "                    levels = np.arange(-0.1,0.101,.002)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.1)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname in ['strairy', 'tau_y']:\n",
    "                    desc = 'Wind stress, y direction (Pa)'\n",
    "                    levels = np.arange(-0.1,0.101,.002)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.1)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                elif varname == 'strocnx':\n",
    "                    desc = '- Ocean-ice stress, x direction (Pa)'\n",
    "                    levels = np.arange(-0.1,0.101,.002)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.1)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                    fact = -1\n",
    "                elif varname == 'strocny':\n",
    "                    desc = '- Ocean-ice stress, y direction (Pa)'\n",
    "                    levels = np.arange(-0.1,0.101,.002)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+.001,0.1)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                    fact = -1\n",
    "                elif varname == 'sst':\n",
    "                    desc = 'SST (C)'\n",
    "                    levels = np.arange(-3,10.1,.25)\n",
    "#                     ticks = np.arange(levels[0],levels[-1]+.001,2)\n",
    "                    ticks = np.arange(-4,100,2)\n",
    "                    ext = 'both'\n",
    "#                     cmp = plt.cm.gist_ncar\n",
    "                    cmp = cm.cm.thermal\n",
    "                elif varname.startswith('SIG'):\n",
    "                    desc = 'Stress invariant ' + varname + ' (units??)'\n",
    "                    levels = np.linspace(-1e6,1e6,50)\n",
    "                    ticks = np.arange(levels[0],levels[-1]+1e4,1e5)\n",
    "                    ext = 'both'\n",
    "                    cmp = cm.cm.balance\n",
    "                else:\n",
    "                    desc = '***UNKNOWN***'\n",
    "                    levels = np.arange(0,1.01,.1)\n",
    "                    ticks = [0,.2,.4,.6,.8,1]\n",
    "                    ext = 'neither'\n",
    "                    cmp = cm.cm.balance\n",
    "\n",
    "\n",
    "                    \n",
    "                var = np.ma.masked_where(varmask,var)\n",
    "\n",
    "\n",
    "                font = {'size':24}\n",
    "                tick_font=24\n",
    "\n",
    "                fig = plt.figure(figsize=(30,30))\n",
    "                if view == 'Arctic':\n",
    "                    m = Basemap(projection ='npstere',boundinglat=45,lon_0=-10,resolution='l')\n",
    "                elif view == 'Arctic_closeup':\n",
    "                    m = Basemap(projection ='npstere',boundinglat=70,lon_0=-10,resolution='l')\n",
    "                elif view == 'Aug_crash':\n",
    "                    m = Basemap(projection='stere',llcrnrlat=67,urcrnrlat=75,\n",
    "                                llcrnrlon=-110,urcrnrlon=-90,lat_0=75,lon_0=-100,resolution='l')\n",
    "                elif view == 'May_crash':\n",
    "                    m = Basemap(projection='stere',llcrnrlat=66,urcrnrlat=85,\n",
    "                                llcrnrlon=70,urcrnrlon=110,lat_0=75,lon_0=80,resolution='l')\n",
    "                elif view == 'tripole_bug':\n",
    "                    m = Basemap(projection='stere',llcrnrlat=72,urcrnrlat=85,\n",
    "                                llcrnrlon=70,urcrnrlon=110,lat_0=75,lon_0=80,resolution='l')\n",
    "                elif view == 'EAC':\n",
    "                    m = Basemap(projection='cyl',\n",
    "                                llcrnrlon=142-360,urcrnrlon=175-360,\n",
    "                                llcrnrlat=-48,urcrnrlat=-20,\n",
    "                                resolution='l')\n",
    "                    parlabels = [True,False,False,True]\n",
    "                elif view == 'Agulhas':\n",
    "                    m = Basemap(projection='cyl',\n",
    "                                llcrnrlon=0,urcrnrlon=50,\n",
    "                                llcrnrlat=-48,urcrnrlat=-15,\n",
    "                                resolution='l')\n",
    "                    parlabels = [True,False,False,True]\n",
    "                elif view == 'Kuroshio':\n",
    "                    m = Basemap(projection='cyl',\n",
    "                                llcrnrlon=125-360,urcrnrlon=165-360,\n",
    "                                llcrnrlat=25,urcrnrlat=50,\n",
    "                                resolution='l')\n",
    "                    parlabels = [True,False,False,True]\n",
    "                elif view == 'Gulf_Stream':\n",
    "                    m = Basemap(projection='cyl',\n",
    "                                llcrnrlon=-98,urcrnrlon=-50,\n",
    "                                llcrnrlat=15,urcrnrlat=45,\n",
    "                                resolution='l')\n",
    "                    parlabels = [True,False,False,True]\n",
    "                elif view == 'ACC':\n",
    "                    m = Basemap(projection ='spstere',boundinglat=-30,lon_0=170,resolution='l')\n",
    "                    parlabels = [False,False,False,False]\n",
    "                else:\n",
    "                    m = Basemap(projection ='spstere',boundinglat=-50,lon_0=170,resolution='l')\n",
    "\n",
    "                # m = Basemap(projection ='ortho',resolution='l',lat_0=90,lon_0=45,\n",
    "                #            llcrnrx=-m.urcrnrx/2.,llcrnry=-m.urcrnry/2.,urcrnrx=m.urcrnrx/2.,urcrnry=m.urcrnry/2.)\n",
    "                # #             llcrnrlat=60,urcrnrlat=60,\n",
    "                # #             llcrnrlon=-45,urcrnrlon=45)\n",
    "                # #             llcrnrlon=loncorners[0],urcrnrlon=loncorners[2]) ,boundinglat=48\n",
    "                x,y = m(*(lon_t,lat_t))\n",
    "                ctr = m.contourf(x,y,fact*var,levels=levels,cmap=cmp, extend=ext) #, vmin=0, vmax=1.1)\n",
    "#                 ctr.cmap.set_under(color=cmp(np.min(levels)), alpha=None)\n",
    "#                 ctr.cmap.set_over(color=cmp(np.max(levels)), alpha=None)\n",
    "                ctr.cmap.set_under(color=cmp(0), alpha=None)\n",
    "                ctr.cmap.set_over(color=cmp(.999), alpha=None)\n",
    "#                 if view.startswith('Arctic'):\n",
    "#                     xobs,yobs = m(*(obs_lon_NH,obs_lat_NH))\n",
    "#                     m.contour(xobs,yobs,CN_obs_NH[month-1],[0.3],colors='r')\n",
    "#                 else:\n",
    "#                     xobs,yobs = m(*(obs_lon_SH,obs_lat_SH))\n",
    "#                     m.contour(xobs,yobs,CN_obs_SH[month-1],[0.3],colors='r')\n",
    "\n",
    "                plt.title(desc + ', ' + date,font)\n",
    "                cbar = m.colorbar(ctr, location = 'bottom', pad = \"6%\")\n",
    "                cbar.set_label(desc,size=tick_font)\n",
    "                cbar.set_ticks(ticks)\n",
    "                cbar_labels=plt.getp(cbar.ax.axes,'xticklabels')\n",
    "                plt.setp(cbar_labels,fontsize=tick_font)\n",
    "                m.drawmapboundary(fill_color='gray') # background color - for non-ocean areas\n",
    "                # fill continents, set lake color same as ocean color.\n",
    "                # m.fillcontinents(color='gray',lake_color='gray')\n",
    "                # draw parallels and meridians.\n",
    "                # label parallels on right and top\n",
    "                # meridians on bottom and left\n",
    "                parallels = np.arange(-80.,81,10.)\n",
    "                # labels = [left,right,top,bottom]\n",
    "#                 m.drawparallels(parallels,color='white')#,labels=[False,True,True,False],size=tick_font)\n",
    "                m.drawparallels(parallels,color='white',labels=parlabels,size=tick_font)\n",
    "\n",
    "#                 meridians = np.arange(0.,351.,20.)\n",
    "                meridians = np.arange(0.,351.,10.)\n",
    "                m.drawmeridians(meridians,labels=[True,False,False,True],size=tick_font,color='white')\n",
    "                savefigure(outname)\n",
    "                plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}