pp_FHIST_calculate_extremes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate and save extremes (both atm and lnd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Settings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Import the necessary python libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "import sys\n",
    "import os\n",
    "from getpass import getuser\n",
    "import string\n",
    "import subprocess\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import netCDF4 as netcdf4\n",
    "import xarray as xr\n",
    "import pandas\n",
    "import regionmask\n",
    "import cartopy.crs as ccrs\n",
    "from IPython.display import display, Math, Latex\n",
    "import warnings\n",
    "from datetime import datetime"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 General Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set directories\n",
    "outdir = '/glade/scratch/ivanderk/'\n",
    "\n",
    "# Define directory where processing is done -- subject to change\n",
    "procdir =  '/glade/work/ivanderk/postprocessing/' \n",
    "\n",
    "# go to processing directory \n",
    "os.chdir(procdir)\n",
    "\n",
    "# ignore all runtime warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 User settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set case name\n",
    "case_res   = 'f.FHIST.f09_f09_mg17.CTL'\n",
    "case_nores = 'f.FHIST.f09_f09_mg17.NORES'\n",
    "\n",
    "# set number of ensemble members\n",
    "n_ens = 5\n",
    "\n",
    "# set individual case names for reference\n",
    "case_res_ind   = 'f.FHIST.f09_f09_mg17.CTL.001'\n",
    "case_nores_ind   = 'f.FHIST.f09_f09_mg17.NORES.001'\n",
    "case   = 'f.FHIST.f09_f09_mg17.CTL.001'\n",
    "\n",
    "# run settings -- change this to terms directly? \n",
    "block = 'atm'  # lnd data\n",
    "               # atm data\n",
    "               # rof data\n",
    "        \n",
    "stream = 'h0'  # h0 output block\n",
    "               # h1 output block\n",
    "               # h2 output block\n",
    "               # xtrm calculated (annual)\n",
    "        \n",
    "# define start and end year\n",
    "spstartyear = '1979'   # spin up start year \n",
    "startyear   = '1979'   # start year, spin up excluded\n",
    "endyear     = '2014'   # last year of the simulation\n",
    "\n",
    "# list of hydrological variables which need to be converted from m/s to mm/day\n",
    "hydrol_vars = ['PRECT','PRECMC', 'PRECC','PRECL','Rx1day']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Functions to open datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# open nc variable as dataset and interpolate lnd variables to atm grid\n",
    "def open_ds(var,case=case,stream=stream, block=block):\n",
    "           \n",
    "    tfreqs = {'h0' : 'month_1'                     , 'h1' : 'day_1'                            , 'h2' : 'month_1'}\n",
    "    tspans = {'h0' : spstartyear+'01-'+endyear+'12', 'h1' : spstartyear+'0101-'+ endyear+'1231', 'h2' : spstartyear+'01-'+endyear+'12'} \n",
    "\n",
    "    model = {'lnd' : 'clm2', 'atm' : 'cam', 'rof' : 'mosart'}\n",
    "\n",
    "    # Define directory where timeseries data is stored\n",
    "    tseriesdir = outdir + 'archive/' + case + '/' + block + '/proc/tseries/' + tfreqs[stream] + '/'\n",
    "\n",
    "    # define filename\n",
    "    fn = case + '.'+ model[block] + '.' + stream + '.' + var + '.' + tspans[stream] +'.nc'\n",
    "\n",
    "    # check if variable timeseries exists and open variable as data array\n",
    "    if not os.path.isfile(tseriesdir + fn):\n",
    "        print(fn + ' does not exists in ')\n",
    "        print(tseriesdir)\n",
    "        return\n",
    "    else: \n",
    "        \n",
    "        # open the dataset\n",
    "        ds = xr.open_dataset(tseriesdir+fn)\n",
    "        \n",
    "        # the lats of the atm and lnd grid differ with about E-7. \n",
    "        # therefore, interpolate land to atm grid to be able to work with exactly the same grids. \n",
    "        if block == 'lnd':\n",
    "            ds_atm = open_ds('TREFHT',block='atm', stream = stream )\n",
    "            ds = ds.interp_like(ds_atm)\n",
    "            \n",
    "            \n",
    "    return ds\n",
    "\n",
    "# function to cut out analysis period out of data-array (1900-2015)\n",
    "def extract_anaperiod(da, stream): \n",
    "    \n",
    "    # number of spin up years\n",
    "    nspinupyears = int(startyear) - int(spstartyear)\n",
    "    \n",
    "    if nspinupyears == 0 :\n",
    "        # no spin up \n",
    "        da = da[:-1,:,:]\n",
    "        \n",
    "    elif stream == 'h1' : # this option still to test \n",
    "        # daily timesteps\n",
    "        # last day of previous year is also saved in variable therefore add one\n",
    "        nspinupdays = (nspinupyears * 365) + 1\n",
    "    \n",
    "        # exclude spin up year and last timestep ()\n",
    "        da = da[nspinupdays:-1,:,:]\n",
    "        \n",
    "    else: \n",
    "        # spin up with monthly timestep\n",
    "        # first month of first year is not saved in variable therefore substract one\n",
    "        nspinupmonths = (nspinupyears * 12) - 1\n",
    "    \n",
    "        # exclude spin up year and last timestep ()\n",
    "        da = da[nspinupmonths:-1,:,:]\n",
    "\n",
    "    return da\n",
    "\n",
    "\n",
    "# open variable as data-array\n",
    "def open_da(var, case=case, stream=stream, block=block):\n",
    "\n",
    "    ds = open_ds(var, case, stream, block=block)\n",
    "\n",
    "    da = ds[var]\n",
    "    \n",
    "    # extract analysis period - not necessary\n",
    "    da = extract_anaperiod(da, stream)\n",
    "    \n",
    "    return da\n",
    "\n",
    "\n",
    "\n",
    "def open_da_delta(var, case, case_ref, stream=stream, block=block): \n",
    "\n",
    "    # Load the two datasets\n",
    "    da_res = open_da(var,case=case, stream=stream, block=block)\n",
    "    da_ctl = open_da(var,case=case_ref, stream=stream, block=block)\n",
    "\n",
    "    # calculate difference and update attributes\n",
    "    da_delta = da_res - da_ctl\n",
    "\n",
    "    da_delta.attrs['long_name'] = '$\\Delta$ '+ da_ctl.long_name\n",
    "    da_delta.attrs['units'] = da_ctl.units\n",
    "    da_delta.name = '$\\Delta$ '+ da_ctl.name\n",
    "\n",
    "    return da_delta\n",
    "\n",
    "\n",
    "\n",
    "# save dataset as nc in postprocessing dir for extremes   \n",
    "def save_da_xtrm(da,var,case=case, block=block):\n",
    "           \n",
    "    tspan = spstartyear+'01-'+endyear+'12' \n",
    "    model = {'lnd' : 'clm2', 'atm' : 'cam', 'rof' : 'mosart'}\n",
    "\n",
    "    savedir = procdir + 'postprocessing/extremes/f09_f09/'\n",
    "\n",
    "    # define filename\n",
    "    fn = case + '.'+ model[block] + '.' + var + '.' + tspan +'.nc'\n",
    "\n",
    "    # check if variable timeseries exists and open variable as data array\n",
    "    if os.path.isfile(savedir + fn):\n",
    "        print(fn + ' already exists')\n",
    "    else: \n",
    "        da.to_dataset().to_netcdf(savedir+fn)\n",
    "    return\n",
    "\n",
    "\n",
    "# open dataset of extremes\n",
    "def open_da_xtrm(var,case=case, block=block):\n",
    "\n",
    "    tspan = spstartyear+'01-'+endyear+'12'\n",
    "    model = {'lnd' : 'clm2', 'atm' : 'cam', 'rof' : 'mosart'}\n",
    "\n",
    "    savedir = outdir + 'postprocessing/extremes/f09_f09/'\n",
    "\n",
    "    # define filename\n",
    "    fn = case + '.'+ model[block] + '.' + var + '.' + tspan +'.nc'\n",
    "\n",
    "    # check if variable timeseries exists and open variable as data array\n",
    "    if not os.path.isfile(savedir + fn):\n",
    "        print(fn + ' does not exist')\n",
    "        return\n",
    "    else: \n",
    "        ds = xr.open_dataset(savedir+fn)\n",
    "        \n",
    "    da = ds[var]\n",
    "        \n",
    "    return da\n",
    "\n",
    "# check if dataset of extremes exists\n",
    "def exist_da_xtrm(var,case=case, block=block):\n",
    "\n",
    "    tspan = spstartyear+'01-'+endyear+'12'\n",
    "    model = {'lnd' : 'clm2', 'atm' : 'cam', 'rof' : 'mosart'}\n",
    "\n",
    "    savedir = outdir + 'postprocessing/extremes/f09_f09/'\n",
    "\n",
    "    # define filename\n",
    "    fn = case + '.'+ model[block] + '.' + var + '.' + tspan +'.nc'\n",
    "\n",
    "    # check if variable timeseries exists and open variable as data array\n",
    "    if not os.path.isfile(savedir + fn): exists = False\n",
    "    else: exists = True\n",
    "        \n",
    "    return exists\n",
    "\n",
    "# remove nc file of extreme (for development purposes)\n",
    "def remove_da_xtrm(var,case=case, block=block):\n",
    "\n",
    "    tspan = spstartyear+'01-'+endyear+'12'\n",
    "    model = {'lnd' : 'clm2', 'atm' : 'cam', 'rof' : 'mosart'}\n",
    "\n",
    "    savedir = outdir + 'postprocessing/extremes/f09_f09/'\n",
    "\n",
    "    # define filename\n",
    "    fn = case + '.'+ model[block] + '.' + var + '.' + tspan +'.nc'\n",
    "\n",
    "    # check if variable timeseries exists and open variable as data array\n",
    "    if os.path.isfile(savedir + fn): \n",
    "        \n",
    "        print('removed file '+fn)\n",
    "        os.system('rm '+savedir + fn)\n",
    "        \n",
    "    return \n",
    "\n",
    "\n",
    "# fucntion to open (and calculate) delta of extreme \n",
    "def open_da_delta_xtrm(var, case, case_ref, block=block): \n",
    "\n",
    "    # Load the two datasets\n",
    "    da_res = open_da_xtrm(var,case=case, block=block)\n",
    "    da_ctl = open_da_xtrm(var,case=case_ref, block=block)\n",
    "\n",
    "    # calculate difference and update attributes\n",
    "    da_delta = da_res - da_ctl\n",
    "\n",
    "    da_delta.attrs['long_name'] = '$\\Delta$ '+ da_ctl.long_name\n",
    "    da_delta.attrs['units'] = da_ctl.units\n",
    "    da_delta.name = '$\\Delta$ '+ da_ctl.name\n",
    "\n",
    "    return da_delta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Functions to do conversions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def conv_m_s_to_mm_day(da_in):\n",
    "\n",
    "    if not da_in.attrs['units'] == 'mm/day':\n",
    "        da_out = da_in * 86400000  \n",
    "        # update attributes and change units\n",
    "        da_out.attrs= da_in.attrs\n",
    "        da_out.attrs['units'] = 'mm/day' \n",
    "    else: \n",
    "        da_out = da_in\n",
    "    return da_out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Functions to calculate extremes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# process TXx: calculate and save annual maximum of maxdaytime temperature \n",
    "def proc_TXx(var_or, case, block=block):\n",
    "    \n",
    "    # define new variable name\n",
    "    var = 'TXx'\n",
    "    \n",
    "    # check if var is already existing\n",
    "    if  exist_da_xtrm(var,case=case, block=block):\n",
    "        print(var +' already exists')\n",
    "    else: # do calculations\n",
    "        \n",
    "        # open da with daily data\n",
    "        da = open_da(var_or,case=case, stream='h1', block=block)\n",
    "\n",
    "        # calculate maximum per year\n",
    "        da_xtrm= da.groupby('time.year').max(keep_attrs=True)\n",
    "\n",
    "        da_xtrm.name = var\n",
    "        da_xtrm.attrs['long_name'] = 'Annual maximum of '+da.long_name\n",
    "        \n",
    "        # save variable into netcdf\n",
    "        save_da_xtrm(da_xtrm,var,case=case, block=block)\n",
    "\n",
    "# process TNn: calculate and save annual maximum of maxdaytime temperature \n",
    "def proc_TNn(var_or, case, block=block):\n",
    "    \n",
    "    # define new variable name\n",
    "    var = 'TNn'\n",
    "    \n",
    "    # check if var is already existing\n",
    "    if  exist_da_xtrm(var,case=case, block=block):\n",
    "        print(var +' already exists')\n",
    "    else: # do calculations\n",
    "        \n",
    "        # open da with daily data\n",
    "        da = open_da(var_or,case=case, stream='h1', block=block)\n",
    "\n",
    "        # calculate minimum per year\n",
    "        da_xtrm = da.groupby('time.year').min(keep_attrs=True)\n",
    "\n",
    "        da_xtrm.name = var\n",
    "        da_xtrm.attrs['long_name'] = 'Annual minimum of '+da.long_name\n",
    "        \n",
    "        # save variable into netcdf\n",
    "        save_da_xtrm(da_xtrm,var,case=case, block=block)\n",
    "\n",
    "        \n",
    "# calculate 99th percentile of max daytime temperatures \n",
    "def proc_TX99(var_or, case, block=block):\n",
    "    \n",
    "    # define new variable name\n",
    "    var = 'TX99'\n",
    "    \n",
    "    # check if var is already existing\n",
    "    if  exist_da_xtrm(var,case=case, block=block):\n",
    "        print(var +' already exists')\n",
    "    else: # do calculations\n",
    "        \n",
    "        # open da with daily data\n",
    "        da = open_da_ens(var_or,case=case, stream='h0', block=block, mode='all')\n",
    "        da_lumped = da.stack(dim=(\"ens_member\", \"time\"))\n",
    "        \n",
    "        # calculate maximum per year\n",
    "        da_xtrm  = da_lumped.quantile(0.99, dim=('dim'))\n",
    "        \n",
    "\n",
    "        da_xtrm.name = var\n",
    "        da_xtrm.attrs['long_name'] = '99th percentile of daily '+da.long_name\n",
    "        da_xtrm.attrs['units'] = 'K'\n",
    "\n",
    "        # save variable into netcdf\n",
    "        save_da_xtrm(da_xtrm,var,case=case, block=block)        \n",
    "\n",
    "        \n",
    "# calculate 1st percentile of min nighttime temperatures \n",
    "def proc_TN01(var_or, case, block=block):\n",
    "    \n",
    "    # define new variable name\n",
    "    var = 'TN01'\n",
    "    \n",
    "    # check if var is already existing\n",
    "    if  exist_da_xtrm(var,case=case, block=block):\n",
    "        print(var +' already exists')\n",
    "    else: # do calculations\n",
    "        \n",
    "        # open da with daily data\n",
    "        da = open_da_ens(var_or,case=case, stream='h0', block=block, mode='all')\n",
    "        da_lumped = da.stack(dim=(\"ens_member\", \"time\"))\n",
    "        \n",
    "        # calculate maximum per year\n",
    "        da_xtrm = da_lumped.quantile(0.01, dim=('dim'))\n",
    "        da_xtrm.name = var\n",
    "        da_xtrm.attrs['long_name'] = '1st percentile of daily '+da.long_name\n",
    "        da_xtrm.attrs['units'] = 'K'\n",
    "        \n",
    "        # save variable into netcdf\n",
    "        save_da_xtrm(da_xtrm,var,case=case, block=block)   \n",
    "\n",
    "def proc_TN10(var_or, case, block=block):\n",
    "    \"\"\"calculcate and save cold days 10pctl of days within period \"\"\"\n",
    "    # define new variable name\n",
    "    var = 'TN10'\n",
    "    \n",
    "    # check if var is already existing\n",
    "    if  exist_da_xtrm(var,case=case, block=block):\n",
    "        print(var +' already exists')\n",
    "    else: # do calculations\n",
    "        \n",
    "        # open da with daily data\n",
    "        da = open_da_ens(var_or,case=case, stream='h0', block=block, mode='all')\n",
    "        da_lumped = da.stack(dim=(\"ens_member\", \"time\"))\n",
    "        \n",
    "        # calculate maximum per year\n",
    "        da_xtrm = da_lumped.quantile(0.1, dim=('dim'))\n",
    "        da_xtrm.name = var\n",
    "        da_xtrm.attrs['long_name'] = '10th percentile of daily'+da.long_name\n",
    "        da_xtrm.attrs['units'] = 'K'\n",
    "        \n",
    "        # save variable into netcdf\n",
    "        save_da_xtrm(da_xtrm,var,case=case, block=block)   \n",
    "\n",
    "def proc_TX90(var_or, case, block=block):\n",
    "    \"\"\"calculcate and save warm days 90th pctl of days within period \"\"\"\n",
    "\n",
    "    # define new variable name\n",
    "    var = 'TX90'\n",
    "    \n",
    "    # check if var is already existing\n",
    "    if  exist_da_xtrm(var,case=case, block=block):\n",
    "        print(var +' already exists')\n",
    "    else: # do calculations\n",
    "        \n",
    "        # open da with daily data\n",
    "        da = open_da_ens(var_or,case=case, stream='h0', block=block, mode='all')\n",
    "        da_lumped = da.stack(dim=(\"ens_member\", \"time\"))\n",
    "        \n",
    "        # calculate maximum per year\n",
    "        da_xtrm  = da_lumped.quantile(0.90, dim=('dim'))\n",
    "        \n",
    "\n",
    "        da_xtrm.name = var\n",
    "        da_xtrm.attrs['long_name'] = '90th percentile of daily'+da.long_name\n",
    "        da_xtrm.attrs['units'] = 'K'\n",
    "\n",
    "        # save variable into netcdf\n",
    "        save_da_xtrm(da_xtrm,var,case=case, block=block)  \n",
    "    \n",
    "def proc_Rx1day(var_or, case, block=block):\n",
    "    \"\"\"process Rx1day: calculate annual maximum 1 day precipitation\"\"\"\n",
    "    # define new variable name\n",
    "    var = 'Rx1day'\n",
    "    \n",
    "    # check if var is already existing\n",
    "    if  exist_da_xtrm(var,case=case, block=block):\n",
    "        print(var +' already exists')\n",
    "    else: # do calculations\n",
    "        \n",
    "        # open da with daily data\n",
    "        da = open_da(var_or,case=case, stream='h1', block=block)\n",
    "\n",
    "        # calculate maximum per year\n",
    "        da_xtrm= da.groupby('time.year').max(keep_attrs=True)\n",
    "\n",
    "        da_xtrm.name = var\n",
    "        da_xtrm.attrs['long_name'] = 'Annual maximum of '+da.long_name\n",
    "        \n",
    "        # save variable into netcdf\n",
    "        save_da_xtrm(da_xtrm,var,case=case, block=block)\n",
    "\n",
    "def proc_R05(var_or, case, block=block):\n",
    "    \"\"\"calculcate and save 5th pctl of monthly precip: drought months  \"\"\"\n",
    "    # define new variable name\n",
    "    var = 'R05'\n",
    "\n",
    "    # check if var is already existing\n",
    "    if  exist_da_xtrm(var,case=case, block=block):\n",
    "        print(var +' already exists')\n",
    "    else: # do calculations\n",
    "\n",
    "        # open da with daily data\n",
    "        da = open_da_ens(var_or,case=case, stream='h0', block=block, mode='all')\n",
    "        da_lumped = da.stack(dim=(\"ens_member\", \"time\"))\n",
    "        # calculate quantile over months\n",
    "\n",
    "        da_xtrm = da_lumped.quantile(0.05, dim=('dim'))\n",
    "        da_xtrm.name = var\n",
    "        da_xtrm.attrs['long_name'] = '5th percentile of monthly'+da.long_name\n",
    "        da_xtrm.attrs['units'] = 'mm/year'\n",
    "\n",
    "        # save variable into netcdf\n",
    "        save_da_xtrm(da_xtrm,var,case=case, block=block)   \n",
    "       \n",
    "def proc_R95(var_or, case, block=block):\n",
    "    \"\"\"calculcate and save 95th pctl of monthly precip: wet months  \"\"\"\n",
    "    # define new variable name\n",
    "    var = 'R95'\n",
    "\n",
    "    # check if var is already existing\n",
    "    if  exist_da_xtrm(var,case=case, block=block):\n",
    "        print(var +' already exists')\n",
    "    else: # do calculations\n",
    "\n",
    "        # open da with daily data\n",
    "        da = open_da_ens(var_or,case=case, stream='h0', block=block, mode='all')\n",
    "        da_lumped = da.stack(dim=(\"ens_member\", \"time\"))\n",
    "        # calculate quantile over months\n",
    "\n",
    "        da_xtrm = da_lumped.quantile(0.95, dim=('dim'))\n",
    "        da_xtrm.name = var\n",
    "        da_xtrm.attrs['long_name'] = '95th percentile of monthly'+da.long_name\n",
    "        da_xtrm.attrs['units'] = 'mm/year'\n",
    "\n",
    "        # save variable into netcdf\n",
    "        save_da_xtrm(da_xtrm,var,case=case, block=block)   \n",
    "        \n",
    "        \n",
    "def proc_colddays(case):\n",
    "    \"\"\"calculcate and save annual percent of cold days (T< 10pctl) of days within period for individual ensemble members\"\"\"\n",
    "    \n",
    "    # define new variable name\n",
    "    var = 'ColdDays_pct'\n",
    "    \n",
    "    # check if var is already existing\n",
    "    if  exist_da_xtrm(var,case=case):\n",
    "        print(var +' already exists')\n",
    "    else: # do calculations\n",
    "        \n",
    "        # open da with daily data\n",
    "        da = open_da('TREFHTMN', case=case, stream='h1', block='atm')\n",
    "        \n",
    "        # calculate percentile over whole period\n",
    "        da_pctl = da.quantile(0.1, dim=('time'))\n",
    "        \n",
    "        ncolddays_annual = (da<da_pctl).groupby('time.year').sum()\n",
    "        ndays_annual = (da<da_pctl).groupby('time.year').count()\n",
    "        \n",
    "        # calculate percent\n",
    "        da_xtrm = ncolddays_annual/ndays_annual *100\n",
    "\n",
    "        da_xtrm.name = var\n",
    "        da_xtrm.attrs['long_name'] = 'Percent of cold days per year (days with TN < TN_10pctl )'\n",
    "        da_xtrm.attrs['units'] = 'K'\n",
    "        \n",
    "        # save variable into netcdf\n",
    "        save_da_xtrm(da_xtrm,var,case=case, block=block)   \n",
    "\n",
    "def proc_warmdays(case):\n",
    "    \"\"\"calculcate and save annual percent of warm days (T> 90pctl) of days within period for individual ensemble members\"\"\"\n",
    "    \n",
    "    # define new variable name\n",
    "    var = 'WarmDays_pct'\n",
    "    \n",
    "    # check if var is already existing\n",
    "    if  exist_da_xtrm(var,case=case):\n",
    "        print(var +' already exists')\n",
    "    else: # do calculations\n",
    "        \n",
    "        # open da with daily data\n",
    "        da = open_da('TREFHTMX', case=case, stream='h1', block='atm')\n",
    "        \n",
    "        # calculate percentile over whole period\n",
    "        da_pctl = da.quantile(0.9, dim=('time'))\n",
    "        \n",
    "        nwarmdays_annual = (da>da_pctl).groupby('time.year').sum()\n",
    "        ndays_annual = (da>da_pctl).groupby('time.year').count()\n",
    "        \n",
    "        # calculate percent\n",
    "        da_xtrm = nwarmdays_annual/ndays_annual *100\n",
    "\n",
    "        da_xtrm.name = var\n",
    "        da_xtrm.attrs['long_name'] = 'Percent of warm days per year (days with TX > TX_90pctl )'\n",
    "        da_xtrm.attrs['units'] = 'K'\n",
    "        \n",
    "        # save variable into netcdf\n",
    "        save_da_xtrm(da_xtrm,var,case=case, block=block)   \n",
    "\n",
    "def proc_TXx_monthly(var_or, case, block=block):\n",
    "    \"\"\" process TXx: calculate and save monthly maximum of daytime temperature and save in monthly folder\"\"\"\n",
    "    \n",
    "    # define new variable name\n",
    "    var = 'TXx_m'\n",
    "    \n",
    "\n",
    "\n",
    "    # open da with daily data\n",
    "    da = open_da(var_or,case=case, stream='h1', block=block)\n",
    "\n",
    "    # calculate maximum per year\n",
    "    da_xtrm = da.resample(time='M').max(keep_attrs=True)\n",
    "\n",
    "    da_xtrm.name = var\n",
    "    da_xtrm.attrs['long_name'] = 'Monthly maximum of '+da.long_name\n",
    "\n",
    "\n",
    "    # save da into netcdf \n",
    "\n",
    "    stream = 'h0'\n",
    "\n",
    "    # save into new file \n",
    "    tfreqs = {'h0' : 'month_1'                     , 'h1' : 'day_1'                            , 'h2' : 'month_1'}\n",
    "    tspans = {'h0' : spstartyear+'01-'+endyear+'12', 'h1' : spstartyear+'0101-'+ endyear+'1231', 'h2' : spstartyear+'01-'+endyear+'12'} \n",
    "\n",
    "    model = {'lnd' : 'clm2', 'atm' : 'cam', 'rof' : 'mosart'}\n",
    "\n",
    "    # Define directory where timeseries data is stored\n",
    "    savedir = outdir + 'archive/' + case + '/' + block + '/proc/tseries/' + tfreqs[stream] + '/'\n",
    "\n",
    "    # define filename\n",
    "    fn = case + '.'+ model[block] + '.' + stream + '.' + var + '.' + tspans[stream] +'.nc'\n",
    "\n",
    "    # check if variable timeseries exists and open variable as data array\n",
    "    #if os.path.isfile(savedir + fn):\n",
    "     #   print(fn + ' already exists')\n",
    "    #else: \n",
    "    da_xtrm.to_dataset().to_netcdf(savedir+fn)\n",
    "\n",
    "# process TXx: calculate and save annual maximum of maxdaytime temperature \n",
    "def proc_TNn_monthly(var_or, case, block=block):\n",
    "        \"\"\" process TNn: calculate and save monthly minimum of nighttime temperature and save in monthly folder\"\"\"\n",
    "\n",
    "        # define new variable name\n",
    "        var = 'TNn_m'\n",
    "\n",
    "        # open da with daily data\n",
    "        da = open_da(var_or,case=case, stream='h1', block=block)\n",
    "\n",
    "        # calculate maximum per year\n",
    "        da_xtrm = da.resample(time='M').min(keep_attrs=True)\n",
    "        da_xtrm.name = var\n",
    "        da_xtrm.attrs['long_name'] = 'Monthly minimum of '+da.long_name\n",
    "\n",
    "        # save da into netcdf\n",
    "        stream = 'h0'\n",
    "\n",
    "        # save into new file \n",
    "        tfreqs = {'h0' : 'month_1'                     , 'h1' : 'day_1'                            , 'h2' : 'month_1'}\n",
    "        tspans = {'h0' : spstartyear+'01-'+endyear+'12', 'h1' : spstartyear+'0101-'+ endyear+'1231', 'h2' : spstartyear+'01-'+endyear+'12'} \n",
    "\n",
    "        model = {'lnd' : 'clm2', 'atm' : 'cam', 'rof' : 'mosart'}\n",
    "\n",
    "        # Define directory where timeseries data is stored\n",
    "        savedir = outdir + 'archive/' + case + '/' + block + '/proc/tseries/' + tfreqs[stream] + '/'\n",
    "\n",
    "        # define filename\n",
    "        fn = case + '.'+ model[block] + '.' + stream + '.' + var + '.' + tspans[stream] +'.nc'\n",
    "\n",
    "        # check if variable timeseries exists and open variable as data array\n",
    "        #if os.path.isfile(savedir + fn):\n",
    "         #   print(fn + ' already exists')\n",
    "        #else: \n",
    "        da_xtrm.to_dataset().to_netcdf(savedir+fn)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Calculate and save extremes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Land variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Atmosphere variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "ename": "PermissionError",
     "evalue": "[Errno 13] Permission denied: b'/glade/work/ivanderk/postprocessing/postprocessing/extremes/f09_f09/f.FHIST.f09_f09_mg17.CTL.001.cam.TXx.197901-201412.nc'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m~/miniconda/envs/ctsm/lib/python3.7/site-packages/xarray/backends/file_manager.py\u001b[0m in \u001b[0;36m_acquire_with_cache_info\u001b[0;34m(self, needs_lock)\u001b[0m\n\u001b[1;32m    197\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 198\u001b[0;31m                 \u001b[0mfile\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_key\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    199\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda/envs/ctsm/lib/python3.7/site-packages/xarray/backends/lru_cache.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m     52\u001b[0m         \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_lock\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 53\u001b[0;31m             \u001b[0mvalue\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     54\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cache\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmove_to_end\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: [<class 'netCDF4._netCDF4.Dataset'>, ('/glade/work/ivanderk/postprocessing/postprocessing/extremes/f09_f09/f.FHIST.f09_f09_mg17.CTL.001.cam.TXx.197901-201412.nc',), 'a', (('clobber', True), ('diskless', False), ('format', 'NETCDF4'), ('persist', False))]",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mPermissionError\u001b[0m                           Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-7-1ba36b45c353>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m     \u001b[0;31m# calculate annual maximum daytime temperature (based on cam variable)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m     \u001b[0mproc_TXx\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'TREFHTMX'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcase_res_mem\u001b[0m\u001b[0;34m,\u001b[0m   \u001b[0mblock\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'atm'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     10\u001b[0m     \u001b[0mproc_TXx\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'TREFHTMX'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcase_nores_mem\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mblock\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'atm'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-6-c104d05baec8>\u001b[0m in \u001b[0;36mproc_TXx\u001b[0;34m(var_or, case, block)\u001b[0m\n\u001b[1;32m     20\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     21\u001b[0m         \u001b[0;31m# save variable into netcdf\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 22\u001b[0;31m         \u001b[0msave_da_xtrm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mda_xtrm\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mvar\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcase\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcase\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mblock\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mblock\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     23\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     24\u001b[0m \u001b[0;31m# process TNn: calculate and save annual maximum of maxdaytime temperature\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-4-1e8d516df977>\u001b[0m in \u001b[0;36msave_da_xtrm\u001b[0;34m(da, var, case, block)\u001b[0m\n\u001b[1;32m    107\u001b[0m         \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfn\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m' already exists'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    108\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 109\u001b[0;31m         \u001b[0mda\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_dataset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_netcdf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msavedir\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mfn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    110\u001b[0m     \u001b[0;32mreturn\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    111\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda/envs/ctsm/lib/python3.7/site-packages/xarray/core/dataset.py\u001b[0m in \u001b[0;36mto_netcdf\u001b[0;34m(self, path, mode, format, group, engine, encoding, unlimited_dims, compute, invalid_netcdf)\u001b[0m\n\u001b[1;32m   1543\u001b[0m             \u001b[0munlimited_dims\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0munlimited_dims\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1544\u001b[0m             \u001b[0mcompute\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcompute\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1545\u001b[0;31m             \u001b[0minvalid_netcdf\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minvalid_netcdf\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1546\u001b[0m         )\n\u001b[1;32m   1547\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda/envs/ctsm/lib/python3.7/site-packages/xarray/backends/api.py\u001b[0m in \u001b[0;36mto_netcdf\u001b[0;34m(dataset, path_or_file, mode, format, group, engine, encoding, unlimited_dims, compute, multifile, invalid_netcdf)\u001b[0m\n\u001b[1;32m   1069\u001b[0m                 \u001b[0;34m\"unrecognized option 'invalid_netcdf' for engine %s\"\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mengine\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1070\u001b[0m             )\n\u001b[0;32m-> 1071\u001b[0;31m     \u001b[0mstore\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mstore_open\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtarget\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mformat\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgroup\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[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1072\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1073\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0munlimited_dims\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda/envs/ctsm/lib/python3.7/site-packages/xarray/backends/netCDF4_.py\u001b[0m in \u001b[0;36mopen\u001b[0;34m(cls, filename, mode, format, group, clobber, diskless, persist, lock, lock_maker, autoclose)\u001b[0m\n\u001b[1;32m    356\u001b[0m             \u001b[0mnetCDF4\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDataset\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwargs\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[1;32m    357\u001b[0m         )\n\u001b[0;32m--> 358\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmanager\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgroup\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mgroup\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlock\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlock\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mautoclose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mautoclose\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    359\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    360\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_acquire\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mneeds_lock\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\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[0;32m~/miniconda/envs/ctsm/lib/python3.7/site-packages/xarray/backends/netCDF4_.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, manager, group, mode, lock, autoclose)\u001b[0m\n\u001b[1;32m    312\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_group\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgroup\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    313\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_mode\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 314\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mds\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata_model\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    315\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_filename\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mds\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfilepath\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    316\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_remote\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mis_remote_uri\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_filename\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda/envs/ctsm/lib/python3.7/site-packages/xarray/backends/netCDF4_.py\u001b[0m in \u001b[0;36mds\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    365\u001b[0m     \u001b[0;34m@\u001b[0m\u001b[0mproperty\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    366\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mds\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\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[0;32m--> 367\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_acquire\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[0m\u001b[1;32m    368\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    369\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mopen_store_variable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvar\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[0;32m~/miniconda/envs/ctsm/lib/python3.7/site-packages/xarray/backends/netCDF4_.py\u001b[0m in \u001b[0;36m_acquire\u001b[0;34m(self, needs_lock)\u001b[0m\n\u001b[1;32m    359\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    360\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_acquire\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mneeds_lock\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\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[0;32m--> 361\u001b[0;31m         \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_manager\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0macquire_context\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mneeds_lock\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mroot\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    362\u001b[0m             \u001b[0mds\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_nc4_require_group\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mroot\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_group\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_mode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    363\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mds\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda/envs/ctsm/lib/python3.7/contextlib.py\u001b[0m in \u001b[0;36m__enter__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    110\u001b[0m         \u001b[0;32mdel\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    111\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 112\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mnext\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgen\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    113\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mStopIteration\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    114\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mRuntimeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"generator didn't yield\"\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda/envs/ctsm/lib/python3.7/site-packages/xarray/backends/file_manager.py\u001b[0m in \u001b[0;36macquire_context\u001b[0;34m(self, needs_lock)\u001b[0m\n\u001b[1;32m    184\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0macquire_context\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mneeds_lock\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\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    185\u001b[0m         \u001b[0;34m\"\"\"Context manager for acquiring a file.\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 186\u001b[0;31m         \u001b[0mfile\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcached\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_acquire_with_cache_info\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mneeds_lock\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    187\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    188\u001b[0m             \u001b[0;32myield\u001b[0m \u001b[0mfile\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda/envs/ctsm/lib/python3.7/site-packages/xarray/backends/file_manager.py\u001b[0m in \u001b[0;36m_acquire_with_cache_info\u001b[0;34m(self, needs_lock)\u001b[0m\n\u001b[1;32m    202\u001b[0m                     \u001b[0mkwargs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\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    203\u001b[0m                     \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"mode\"\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_mode\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 204\u001b[0;31m                 \u001b[0mfile\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_opener\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_args\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[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    205\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_mode\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m\"w\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    206\u001b[0m                     \u001b[0;31m# ensure file doesn't get overriden when opened again\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mnetCDF4/_netCDF4.pyx\u001b[0m in \u001b[0;36mnetCDF4._netCDF4.Dataset.__init__\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mnetCDF4/_netCDF4.pyx\u001b[0m in \u001b[0;36mnetCDF4._netCDF4._ensure_nc_success\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mPermissionError\u001b[0m: [Errno 13] Permission denied: b'/glade/work/ivanderk/postprocessing/postprocessing/extremes/f09_f09/f.FHIST.f09_f09_mg17.CTL.001.cam.TXx.197901-201412.nc'"
     ]
    }
   ],
   "source": [
    "# get string of individual case names\n",
    "case_res_names = [case_res+'.00'+str(i) for i in range(1,n_ens+1)]\n",
    "case_nores_names = [case_nores+'.00'+str(i) for i in range(1,n_ens+1)]\n",
    "\n",
    "# loop over all cases and calculate extremes for all individual cases\n",
    "for case_res_mem,case_nores_mem in zip(case_res_names, case_nores_names):\n",
    "\n",
    "    # calculate annual maximum daytime temperature (based on cam variable)\n",
    "    proc_TXx('TREFHTMX', case_res_mem,   block='atm')\n",
    "    proc_TXx('TREFHTMX', case_nores_mem, block='atm')\n",
    "    \n",
    "    # calculate annual minimum nighttime temperature (based on cam variable)\n",
    "    proc_TNn('TREFHTMN', case_res_mem,   block='atm')\n",
    "    proc_TNn('TREFHTMN', case_nores_mem, block='atm')\n",
    "\n",
    "    # Rx1day annual maximum 1 day precipitation\n",
    "    proc_Rx1day('PRECT', case_res_mem,   block='atm')\n",
    "    proc_Rx1day('PRECT', case_nores_mem, block='atm')\n",
    "    \n",
    "    # calculate annual % of cold and warm days\n",
    "    proc_warmdays(case_res_mem)\n",
    "    proc_warmdays(case_nores_mem)\n",
    "    proc_colddays(case_res_mem)\n",
    "    proc_colddays(case_nores_mem)\n",
    "    \n",
    "    # calculate monthly maximum daytime temperature (based on cam variable)\n",
    "    proc_TXx_monthly('TREFHTMX', case_res_mem,   block='atm')\n",
    "    proc_TXx_monthly('TREFHTMX', case_nores_mem, block='atm')\n",
    "    \n",
    "    # calculate monthly minimum nighttime temperature (based on cam variable)\n",
    "    proc_TNn_monthly('TREFHTMN', case_res_mem,   block='atm')\n",
    "    proc_TNn_monthly('TREFHTMN', case_nores_mem, block='atm')\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TX99 already exists\n",
      "TX99 already exists\n",
      "TN01 already exists\n",
      "TN01 already exists\n",
      "R05 already exists\n",
      "R05 already exists\n",
      "R95 already exists\n",
      "R95 already exists\n"
     ]
    }
   ],
   "source": [
    "from iv_utils import *\n",
    "\n",
    "# calculate 99 pctl of daytime temperature (based on cam variable)\n",
    "proc_TX99('TREFHTMX', case_res,   block='atm')\n",
    "proc_TX99('TREFHTMX', case_nores, block='atm')\n",
    "\n",
    "# calculate 1st pctl of nighttime temperature (based on cam variable)\n",
    "proc_TN01('TREFHTMN', case_res,   block='atm')\n",
    "proc_TN01('TREFHTMN', case_nores, block='atm')\n",
    "\n",
    "proc_TX90('TREFHTMX', case_res,   block='atm')\n",
    "proc_TX90('TREFHTMX', case_nores, block='atm')\n",
    "\n",
    "proc_TN10('TREFHTMN', case_res,   block='atm')\n",
    "proc_TN10('TREFHTMN', case_nores, block='atm')\n",
    "\n",
    "proc_R05('PRECT', case_res, block='atm')\n",
    "proc_R05('PRECT', case_nores, block='atm')\n",
    "\n",
    "proc_R95('PRECT', case_res, block='atm')\n",
    "proc_R95('PRECT', case_nores, block='atm')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Save h1 as h0 time frequency files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function to do conversion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_h1_as_h0(var,case_name, calcsum=False, block=block):\n",
    "    \"\"\"Convert h1 to h0 and save file \"\"\"\n",
    "    # open da with daily data\n",
    "    da = open_da(var,case=case_name, stream='h1', block=block)\n",
    "\n",
    "    if calcsum: \n",
    "         # calculate monthly mean\n",
    "        da_monthly = da.resample(time='1M').sum(keep_attrs=True)\n",
    "    else:\n",
    "        # calculate monthly mean\n",
    "        da_monthly = da.resample(time='1M').mean(keep_attrs=True)\n",
    "\n",
    "    stream_new = 'h0'\n",
    "\n",
    "    # save into new file \n",
    "    tfreqs = {'h0' : 'month_1'                     , 'h1' : 'day_1'                            , 'h2' : 'month_1'}\n",
    "    tspans = {'h0' : spstartyear+'01-'+endyear+'12', 'h1' : spstartyear+'0101-'+ endyear+'1231', 'h2' : spstartyear+'01-'+endyear+'12'} \n",
    "\n",
    "    model = {'lnd' : 'clm2', 'atm' : 'cam', 'rof' : 'mosart'}\n",
    "\n",
    "    # Define directory where timeseries data is stored\n",
    "    savedir = outdir + 'archive/' + case_name + '/' + block + '/proc/tseries/' + tfreqs[stream_new] + '/'\n",
    "\n",
    "    # define filename\n",
    "    fn = case_name + '.'+ model[block] + '.' + stream_new + '.' + var + '.' + tspans[stream_new] +'.nc'\n",
    "\n",
    "    # check if variable timeseries exists and open variable as data array\n",
    "    if os.path.isfile(savedir + fn):\n",
    "        print(fn + ' already exists')\n",
    "    else: \n",
    "        da_monthly.to_dataset().to_netcdf(savedir+fn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f.FHIST.f09_f09_mg17.CTL.001.cam.h0.TREFHTMX.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.001.cam.h0.TREFHTMX.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.001.cam.h0.TREFHTMN.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.001.cam.h0.TREFHTMN.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.001.cam.h0.PRECT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.001.cam.h0.PRECT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.002.cam.h0.TREFHTMX.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.002.cam.h0.TREFHTMX.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.002.cam.h0.TREFHTMN.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.002.cam.h0.TREFHTMN.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.002.cam.h0.PRECT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.002.cam.h0.PRECT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.003.cam.h0.TREFHTMX.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.003.cam.h0.TREFHTMX.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.003.cam.h0.TREFHTMN.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.003.cam.h0.TREFHTMN.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.003.cam.h0.PRECT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.003.cam.h0.PRECT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.004.cam.h0.PRECT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.004.cam.h0.PRECT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.005.cam.h0.PRECT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.005.cam.h0.PRECT.197901-201412.nc already exists\n"
     ]
    }
   ],
   "source": [
    "# loop over all cases\n",
    "del case_res_names,case_nores_names\n",
    "case_res_names = [case_res+'.00'+str(i) for i in range(1,n_ens+1)]\n",
    "case_nores_names = [case_nores+'.00'+str(i) for i in range(1,n_ens+1)]\n",
    "\n",
    "# loop over all cases and calculate extremes for all individual cases\n",
    "for case_res_mem,case_nores_mem in zip(case_res_names, case_nores_names):\n",
    "    save_h1_as_h0('TREFHTMX',case_res_mem)\n",
    "    save_h1_as_h0('TREFHTMX',case_nores_mem)\n",
    "\n",
    "    save_h1_as_h0('TREFHTMN',case_res_mem)\n",
    "    save_h1_as_h0('TREFHTMN',case_nores_mem)\n",
    "\n",
    "    save_h1_as_h0('PRECT',case_res_mem, calcsum=True)\n",
    "    save_h1_as_h0('PRECT',case_nores_mem, calcsum=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function to calculate DTR (and calculate monthly mean too)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def proc_DTR(case_name): \n",
    "    TREFHTMN = open_da('TREFHTMN',case=case_name, stream='h1', block=block)\n",
    "    TREFHTMX = open_da('TREFHTMX',case=case_name, stream='h1', block=block)\n",
    "\n",
    "    DTR = TREFHTMX-TREFHTMN\n",
    "    DTR.attrs['long_name']= 'Diurnal temperature range'\n",
    "    DTR.attrs['units']= 'K'\n",
    "    DTR.name = 'DTR'\n",
    "\n",
    "    DTR = DTR[1:,:,:]\n",
    "\n",
    "    var = 'DTR'\n",
    "    stream = 'h1'\n",
    "\n",
    "    # save into new file \n",
    "    tfreqs = {'h0' : 'month_1'                     , 'h1' : 'day_1'                            , 'h2' : 'month_1'}\n",
    "    tspans = {'h0' : spstartyear+'01-'+endyear+'12', 'h1' : spstartyear+'0101-'+ endyear+'1231', 'h2' : spstartyear+'01-'+endyear+'12'} \n",
    "\n",
    "    model = {'lnd' : 'clm2', 'atm' : 'cam', 'rof' : 'mosart'}\n",
    "\n",
    "    # Define directory where timeseries data is stored\n",
    "    savedir = outdir + 'archive/' + case_name + '/' + block + '/proc/tseries/' + tfreqs[stream] + '/'\n",
    "\n",
    "    # define filename\n",
    "    fn = case_name + '.'+ model[block] + '.' + stream + '.' + var + '.' + tspans[stream] +'.nc'\n",
    "\n",
    "    # check if variable timeseries exists and open variable as data array\n",
    "    if os.path.isfile(savedir + fn):\n",
    "        print(fn + ' already exists')\n",
    "    else: \n",
    "        DTR.to_dataset().to_netcdf(savedir+fn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f.FHIST.f09_f09_mg17.CTL.001.cam.h1.DTR.19790101-20141231.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.001.cam.h1.DTR.19790101-20141231.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.001.cam.h0.DTR.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.001.cam.h0.DTR.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.002.cam.h1.DTR.19790101-20141231.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.002.cam.h1.DTR.19790101-20141231.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.002.cam.h0.DTR.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.002.cam.h0.DTR.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.003.cam.h1.DTR.19790101-20141231.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.003.cam.h1.DTR.19790101-20141231.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.003.cam.h0.DTR.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.003.cam.h0.DTR.197901-201412.nc already exists\n"
     ]
    }
   ],
   "source": [
    "# loop over all cases\n",
    "#case_res_names = [case_res+'.00'+str(i) for i in range(1,n_ens+1)]\n",
    "#case_nores_names = [case_nores+'.00'+str(i) for i in range(1,n_ens+1)]\n",
    "\n",
    "# loop over all cases and calculate extremes for all individual cases\n",
    "for case_res_mem,case_nores_mem in zip(case_res_names, case_nores_names):\n",
    "    proc_DTR(case_res_mem)\n",
    "    proc_DTR(case_nores_mem)\n",
    "    save_h1_as_h0('DTR',case_res_mem)\n",
    "    save_h1_as_h0('DTR',case_nores_mem)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# process albedo\n",
    "def proc_albedo(case_name):\n",
    "    \n",
    "    SWin = open_da('FSDS', block = 'lnd', case = case_name)\n",
    "    SWout= open_da('FSR' , block = 'lnd', case = case_name)\n",
    "\n",
    "    albedo = SWout/SWin\n",
    "    albedo.attrs['long_name']= 'Albedo'\n",
    "    albedo.attrs['units']= '-'\n",
    "    albedo.name = 'albedo'\n",
    "\n",
    "    var = 'albedo'\n",
    "    stream = 'h0'\n",
    "    block = 'lnd'\n",
    "\n",
    "    # save into new file \n",
    "    tfreqs = {'h0' : 'month_1'                     , 'h1' : 'day_1'                            , 'h2' : 'month_1'}\n",
    "    tspans = {'h0' : spstartyear+'01-'+endyear+'12', 'h1' : spstartyear+'0101-'+ endyear+'1231', 'h2' : spstartyear+'01-'+endyear+'12'} \n",
    "\n",
    "    model = {'lnd' : 'clm2', 'atm' : 'cam', 'rof' : 'mosart'}\n",
    "\n",
    "    # Define directory where timeseries data is stored\n",
    "    savedir = outdir + 'archive/' + case_name + '/' + block + '/proc/tseries/' + tfreqs[stream] + '/'\n",
    "\n",
    "    # define filename\n",
    "    fn = case_name + '.'+ model[block] + '.' + stream + '.' + var + '.' + tspans[stream] +'.nc'\n",
    "\n",
    "    # check if variable timeseries exists and open variable as data array\n",
    "    if os.path.isfile(savedir + fn):\n",
    "        print(fn + ' already exists')\n",
    "    else: \n",
    "        albedo.to_dataset().to_netcdf(savedir+fn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f.FHIST.f09_f09_mg17.CTL.001.clm2.h0.albedo.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.001.clm2.h0.albedo.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.002.clm2.h0.albedo.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.002.clm2.h0.albedo.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.003.clm2.h0.albedo.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.003.clm2.h0.albedo.197901-201412.nc already exists\n"
     ]
    }
   ],
   "source": [
    "from iv_utils import *\n",
    "\n",
    "case_res_names = [case_res+'.00'+str(i) for i in range(1,n_ens+1)]\n",
    "case_nores_names = [case_nores+'.00'+str(i) for i in range(1,n_ens+1)]\n",
    "\n",
    "# loop over all cases and calculate extremes for all individual cases\n",
    "for case_res_mem,case_nores_mem in zip(case_res_names, case_nores_names):\n",
    "    proc_albedo(case_res_mem)\n",
    "    proc_albedo(case_nores_mem)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function to calculate Apparent temperature and wet bulb temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def proc_AT(case_name): \n",
    "    \"\"\"  compute Apparent temperature - Version including the effects of temperature, humidity, and wind\n",
    "      # http://www.bom.gov.au/info/thermal_stress/?cid=003bl08\"\"\"\n",
    "    \n",
    "    TREFHT = open_da('TREFHT',case=case_name, stream='h0', block=block)\n",
    "    RHREFHT = open_da('RHREFHT',case=case_name, stream='h0', block=block)\n",
    "    U10 = open_da('U10',case=case_name, stream='h0', block=block)\n",
    "\n",
    "    p\n",
    "\n",
    "    AT = (TREFHT-273.15) + 0.33 * (RHREFHT / 100 * 6.105 * np.exp( 17.27 * (TREFHT-273.15) / ( 237.7 + (TREFHT-273.15) ) ) ) - 0.70 * U10 - 4.00\n",
    "    AT.attrs['long_name']= 'Apparent Temperature'\n",
    "    AT.attrs['units']= '°C'\n",
    "    AT.name = 'AT'\n",
    "\n",
    "    var = 'AT'\n",
    "    stream = 'h0'\n",
    "\n",
    "    # save into new file \n",
    "    tfreqs = {'h0' : 'month_1'                     , 'h1' : 'day_1'                            , 'h2' : 'month_1'}\n",
    "    tspans = {'h0' : spstartyear+'01-'+endyear+'12', 'h1' : spstartyear+'0101-'+ endyear+'1231', 'h2' : spstartyear+'01-'+endyear+'12'} \n",
    "\n",
    "    model = {'lnd' : 'clm2', 'atm' : 'cam', 'rof' : 'mosart'}\n",
    "\n",
    "    # Define directory where timeseries data is stored\n",
    "    savedir = outdir + 'archive/' + case_name + '/' + block + '/proc/tseries/' + tfreqs[stream] + '/'\n",
    "\n",
    "    # define filename\n",
    "    fn = case_name + '.'+ model[block] + '.' + stream + '.' + var + '.' + tspans[stream] +'.nc'\n",
    "\n",
    "    # check if variable timeseries exists and open variable as data array\n",
    "    if os.path.isfile(savedir + fn):\n",
    "        print(fn + ' already exists')\n",
    "    else: \n",
    "        AT.to_dataset().to_netcdf(savedir+fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function to calculate Wet bulb temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def proc_WBGT(case_name): \n",
    "    \"\"\"  Wet Bulb Global Temperature - approximation to the WBGT used by the Bureau of Meteorology\n",
    "      # http://www.bom.gov.au/info/thermal_stress/?cid=003bl08\n",
    "      # !!! uses 24-h average T and RH !!!\"\"\"\n",
    "    \n",
    "    TREFHT = open_da('TREFHT',case=case_name, stream='h1', block=block)\n",
    "    RHREFHT = open_da('RHREFHT',case=case_name, stream='h1', block=block)\n",
    "    TREFHTMX = open_da('TREFHTMX',case=case_name, stream='h1', block=block)\n",
    "\n",
    "\n",
    "    WBGT = 0.567 * (TREFHT-273.15) + 0.393 * (RHREFHT / 100 * 6.105 * np.exp( 17.27 * (TREFHTMX-273.15) / ( 237.7 + (TREFHTMX-273.15) ) )) + 3.94\n",
    "    WBGT.attrs['long_name']= 'Wet Bulb Global Temperature'\n",
    "    WBGT.attrs['units']= '°C'\n",
    "    WBGT.name = 'WBGT'\n",
    "\n",
    "    var = 'WBGT'\n",
    "    stream = 'h1'\n",
    "\n",
    "    # save into new file \n",
    "    tfreqs = {'h0' : 'month_1'                     , 'h1' : 'day_1'                            , 'h2' : 'month_1'}\n",
    "    tspans = {'h0' : spstartyear+'01-'+endyear+'12', 'h1' : spstartyear+'0101-'+ endyear+'1231', 'h2' : spstartyear+'01-'+endyear+'12'} \n",
    "\n",
    "    model = {'lnd' : 'clm2', 'atm' : 'cam', 'rof' : 'mosart'}\n",
    "\n",
    "    # Define directory where timeseries data is stored\n",
    "    savedir = outdir + 'archive/' + case_name + '/' + block + '/proc/tseries/' + tfreqs[stream] + '/'\n",
    "\n",
    "    # define filename\n",
    "    fn = case_name + '.'+ model[block] + '.' + stream + '.' + var + '.' + tspans[stream] +'.nc'\n",
    "\n",
    "    # check if variable timeseries exists and open variable as data array\n",
    "    if os.path.isfile(savedir + fn):\n",
    "        print(fn + ' already exists')\n",
    "    else: \n",
    "        WBGT.to_dataset().to_netcdf(savedir+fn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f.FHIST.f09_f09_mg17.CTL.001.cam.h0.AT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.001.cam.h0.AT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.002.cam.h0.AT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.002.cam.h0.AT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.003.cam.h0.AT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.003.cam.h0.AT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.004.cam.h0.AT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.004.cam.h0.AT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.CTL.005.cam.h0.AT.197901-201412.nc already exists\n",
      "f.FHIST.f09_f09_mg17.NORES.005.cam.h0.AT.197901-201412.nc already exists\n"
     ]
    }
   ],
   "source": [
    "case_res_names = [case_res+'.00'+str(i) for i in range(1,n_ens+1)]\n",
    "case_nores_names = [case_nores+'.00'+str(i) for i in range(1,n_ens+1)]\n",
    "\n",
    "# loop over all cases and calculate extremes for all individual cases\n",
    "for case_res_mem,case_nores_mem in zip(case_res_names, case_nores_names):\n",
    "    proc_AT(case_res_mem)\n",
    "    proc_AT(case_nores_mem)\n",
    "    proc_WBGT(case_res_mem)\n",
    "    proc_WBGT(case_nores_mem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "case_name = 'f.FHIST.f09_f09_mg17.CTL.001'\n",
    "TREFHT = open_da('TREFHT',case=case_name, stream='h0', block=block)\n",
    "RHREFHT = open_da('RHREFHT',case=case_name, stream='h0', block=block)\n",
    "U10 = open_da('U10',case=case_name, stream='h0', block=block)\n",
    "\n",
    "AT = (TREFHT-273.15) + 0.33 * (RHREFHT / 100 * 6.105 * np.exp( 17.27 * (TREFHT-273.15) / ( 237.7 + (TREFHT-273.15) ) ) ) - 0.70 * U10 - 4.00\n",
    "AT.attrs['long_name']= 'Apparent Temperature'\n",
    "AT.attrs['units']= '°C'\n",
    "AT.name = 'AT'\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:miniconda-ctsm]",
   "language": "python",
   "name": "conda-env-miniconda-ctsm-py"
  },
  "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}