{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data for Figure 8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Based on Elise Olson's https://github.com/SalishSeaCast/analysis-elise-2/blob/master/notebooks/modelEqs/NewRateComparison.ipynb notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import netCDF4 as nc\n",
    "import matplotlib.pyplot as plt\n",
    "from salishsea_tools import bio_tools as bt, places\n",
    "import xarray as xr\n",
    "import os\n",
    "import glob\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "nml=bt.load_nml_bio(resDir='/ocean/eolson/MEOPAR/NEMO-3.6-code/NEMOGCM/CONFIG/SalishSeaCast/EXP00/',\n",
    "                 nmlname='nampisprod',bioRefName='namelist_smelt_cfg_HC201905equiv',bioCfgName='namelist_smelt_cfg_HC201905equiv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Namelist([('zz_rate_r_diat', 6.0495e-05),\n",
       "          ('zz_rate_r_myri', 2.22e-05),\n",
       "          ('zz_rate_r_nano', 2.109e-05),\n",
       "          ('zz_rate_maxtemp_diat', 26.0),\n",
       "          ('zz_rate_maxtemp_myri', 31.0),\n",
       "          ('zz_rate_maxtemp_nano', 31.0),\n",
       "          ('zz_rate_temprange_diat', 14.0),\n",
       "          ('zz_rate_temprange_myri', 13.0),\n",
       "          ('zz_rate_temprange_nano', 13.0),\n",
       "          ('zz_rate_iopt_diat', 45.0),\n",
       "          ('zz_rate_iopt_myri', 37.0),\n",
       "          ('zz_rate_iopt_nano', 10.0),\n",
       "          ('zz_rate_gamma_diat', 0.0),\n",
       "          ('zz_rate_gamma_myri', 0.0),\n",
       "          ('zz_rate_gamma_nano', 0.0),\n",
       "          ('zz_rate_k_si_diat', 2.2),\n",
       "          ('zz_rate_k_si_myri', 0.0),\n",
       "          ('zz_rate_k_si_nano', 0.0),\n",
       "          ('zz_rate_kapa_diat', 1.0),\n",
       "          ('zz_rate_kapa_myri', 0.5),\n",
       "          ('zz_rate_kapa_nano', 0.3),\n",
       "          ('zz_rate_k_diat', 2.0),\n",
       "          ('zz_rate_k_myri', 0.5),\n",
       "          ('zz_rate_k_nano', 0.2),\n",
       "          ('zz_rate_si_ratio_diat', 1.8),\n",
       "          ('zz_rate_si_ratio_myri', 0.0),\n",
       "          ('zz_rate_si_ratio_nano', 0.0)])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def phyto_Tdep_Factor(TT, zz_rate_maxtemp, zz_rate_temprange):\n",
    "    if hasattr(TT,'__len__'): # assume 1-d array or similar and return array\n",
    "        return np.array([phyto_Tdep_Factor(el,zz_rate_maxtemp, zz_rate_temprange) for el in TT])\n",
    "    else:\n",
    "        return np.exp(0.07 * (TT - 20)) * min(max((zz_rate_maxtemp - TT), 0.0),zz_rate_temprange) / (zz_rate_temprange + 1e-10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_T_Factors(TT,nampisprod):\n",
    "    Tdep_Diat=phyto_Tdep_Factor(TT,nampisprod['zz_rate_maxtemp_diat'],nampisprod['zz_rate_temprange_diat'])\n",
    "    Tdep_Myri=phyto_Tdep_Factor(TT,nampisprod['zz_rate_maxtemp_myri'],nampisprod['zz_rate_temprange_myri'])\n",
    "    Tdep_Nano=phyto_Tdep_Factor(TT,nampisprod['zz_rate_maxtemp_nano'],nampisprod['zz_rate_temprange_nano'])\n",
    "    return Tdep_Diat, Tdep_Myri, Tdep_Nano"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Import Temperature Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Original WY\n",
    "\n",
    "\n",
    "monthly_array_temp_orig_slice = np.zeros([14,12,50,50])\n",
    "# Load monthly averages\n",
    "mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc')\n",
    "slc = {'y': slice(450,500), 'x': slice(250,300)}\n",
    "e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')]\n",
    "years, variables = range(2007, 2021), ['votemper']\n",
    "# Temporary list dict\n",
    "data = {}\n",
    "# Permanent aggregate dict\n",
    "aggregates = {var: {} for var in variables}\n",
    "monthlydat = {var: {} for var in variables}\n",
    "\n",
    "### 2008 using higher temperature threshold       \n",
    "# Add experiment year\n",
    "for year in [2019]:\n",
    "    # Initialize lists\n",
    "    for var in variables: data[var] = []\n",
    "    # Load monthly averages\n",
    "    for month in range(1, 13):\n",
    "        datestr = f'{year}{month:02d}'\n",
    "        prefix = f'/data/sallen/results/MEOPAR/v201905r/SalishSea_1m_{datestr}_{datestr}'\n",
    "        \n",
    "        with xr.open_dataset(prefix + '_grid_T.nc') as ds:\n",
    "            q = ds.votemper.isel(deptht=0, **slc).values\n",
    "            q2 = q[0,:,:]\n",
    "            monthly_array_temp_orig_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension\n",
    "            for var in ['votemper']:\n",
    "                data[var].append(ds.votemper.isel(deptht=0, **slc).values)\n",
    "    # Concatenate months\n",
    "    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(14, 12)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_3795255/1734163470.py:3: RuntimeWarning: Mean of empty slice\n",
      "  np.nanmean(np.nanmean(monthly_array_temp_orig_slice, axis = 2),axis = 2)\n"
     ]
    }
   ],
   "source": [
    "monthly_array_temp_orig_slice[monthly_array_temp_orig_slice == 0 ] = np.nan\n",
    "monthly_array_temp_orig_slicemean = \\\n",
    "np.nanmean(np.nanmean(monthly_array_temp_orig_slice, axis = 2),axis = 2)\n",
    "print(np.shape(monthly_array_temp_orig_slicemean))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Experiment 9\n",
    "\n",
    "\n",
    "monthly_array_temp_exp_slice = np.zeros([14,12,50,50])\n",
    "# Load monthly averages\n",
    "mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc')\n",
    "slc = {'y': slice(450,500), 'x': slice(250,300)}\n",
    "e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')]\n",
    "years, variables = range(2007, 2021), ['votemper']\n",
    "# Temporary list dict\n",
    "data = {}\n",
    "# Permanent aggregate dict\n",
    "aggregates = {var: {} for var in variables}\n",
    "monthlydat = {var: {} for var in variables}\n",
    "\n",
    "for year in [2019]:\n",
    "    # Initialize lists\n",
    "    for var in variables: data[var] = []\n",
    "    # Load monthly averages\n",
    "    for month in range(1, 7):\n",
    "        datestr = f'{year}{month:02d}'\n",
    "        prefix = f'/data/sallen/results/MEOPAR/Karyn/01jan19_river08/SalishSea_1m_{datestr}_{datestr}'\n",
    "        \n",
    "        with xr.open_dataset(prefix + '_grid_T.nc') as ds:\n",
    "            q = ds.votemper.isel(deptht=0, **slc).values\n",
    "            q2 = q[0,:,:]\n",
    "            monthly_array_temp_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension\n",
    "            for var in ['votemper']:\n",
    "                data[var].append(ds.votemper.isel(deptht=0, **slc).values)\n",
    "    # Concatenate months\n",
    "    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)\n",
    "\n",
    "for year in [2019]:\n",
    "    # Initialize lists\n",
    "    for var in variables: data[var] = []\n",
    "    # Load monthly averages\n",
    "    for month in range(7, 10):\n",
    "        datestr = f'{year}{month:02d}'\n",
    "        prefix = f'/data/sallen/results/MEOPAR/Karyn/01jul19_river08/SalishSea_1m_{datestr}_{datestr}'\n",
    "        \n",
    "        with xr.open_dataset(prefix + '_grid_T.nc') as ds:\n",
    "            q = ds.votemper.isel(deptht=0, **slc).values\n",
    "            q2 = q[0,:,:]\n",
    "            monthly_array_temp_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension\n",
    "            for var in ['votemper']:\n",
    "                data[var].append(ds.votemper.isel(deptht=0, **slc).values)\n",
    "    # Concatenate months\n",
    "    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)\n",
    "        \n",
    "\n",
    "for year in [2019]:\n",
    "    # Initialize lists\n",
    "    for var in variables: data[var] = []\n",
    "    # Load monthly averages\n",
    "    for month in range(10, 13):\n",
    "        datestr = f'{year}{month:02d}'\n",
    "        prefix = f'/data/sallen/results/MEOPAR/Karyn/11oct19_river08/SalishSea_1m_{datestr}_{datestr}'\n",
    "        \n",
    "        with xr.open_dataset(prefix + '_grid_T.nc') as ds:\n",
    "            q = ds.votemper.isel(deptht=0, **slc).values\n",
    "            q2 = q[0,:,:]\n",
    "            monthly_array_temp_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension\n",
    "            for var in ['votemper']:\n",
    "                data[var].append(ds.votemper.isel(deptht=0, **slc).values)\n",
    "    # Concatenate months\n",
    "    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(14, 12)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_3795255/2338418733.py:3: RuntimeWarning: Mean of empty slice\n",
      "  np.nanmean(np.nanmean(monthly_array_temp_exp_slice, axis = 2),axis = 2)\n"
     ]
    }
   ],
   "source": [
    "monthly_array_temp_exp_slice[monthly_array_temp_exp_slice == 0 ] = np.nan\n",
    "monthly_array_temp_exp_slicemean = \\\n",
    "np.nanmean(np.nanmean(monthly_array_temp_exp_slice, axis = 2),axis = 2)\n",
    "print(np.shape(monthly_array_temp_exp_slicemean))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 6.68492084,  5.18552462,  7.25557932, 10.16054886, 14.0224531 ,\n",
       "       17.57203394, 20.68068024, 20.043708  , 16.67098008, 11.61331057,\n",
       "        8.74718019,  6.62090044])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "monthly_array_temp_exp_slicemean[12,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 6.66574269,  5.18524643,  7.24007046, 10.16174981, 14.03568474,\n",
       "       17.7114818 , 20.61275981, 20.03991294, 16.73050814, 11.58581398,\n",
       "        8.60023705,  6.47922814])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "monthly_array_temp_orig_slicemean[12,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Temperature response factors\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "TdepDiatOrig,__,TdepNanoOrig=calc_T_Factors(monthly_array_temp_orig_slicemean[12,:],nml)\n",
    "TdepDiatExp,__,TdepNanoExp=calc_T_Factors(monthly_array_temp_exp_slicemean[12,:],nml)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Open Nitrate files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Original WY\n",
    "\n",
    "monthly_array_nitrate_orig_slice = np.zeros([14,12,50,50])\n",
    "# Load monthly averages\n",
    "mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc')\n",
    "slc = {'y': slice(450,500), 'x': slice(250,300)} \n",
    "e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')]\n",
    "years, variables = range(2007, 2021), ['nitrate']\n",
    "# Temporary list dict\n",
    "data = {}\n",
    "# Permanent aggregate dict\n",
    "aggregates = {var: {} for var in variables}\n",
    "monthlydat = {var: {} for var in variables}\n",
    "        \n",
    "# Add experiment year\n",
    "for year in [2019]:\n",
    "    # Initialize lists\n",
    "    for var in variables: data[var] = []\n",
    "    # Load monthly averages\n",
    "    for month in range(1, 13):\n",
    "        datestr = f'{year}{month:02d}'\n",
    "        prefix = f'/data/sallen/results/MEOPAR/v201905r/SalishSea_1m_{datestr}_{datestr}'\n",
    "        \n",
    "        \n",
    "        # Load grazing variables\n",
    "        with xr.open_dataset(prefix + '_ptrc_T.nc') as ds:\n",
    "            q = ds.nitrate.isel(deptht=0, **slc).values\n",
    "            q2 = q[0,:,:]\n",
    "            monthly_array_nitrate_orig_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension\n",
    "            for var in ['nitrate']:\n",
    "                data[var].append(ds.nitrate.isel(deptht=0, **slc).values)\n",
    "                \n",
    "    # Concatenate months\n",
    "    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)        \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(14, 12)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_3795255/3312634990.py:3: RuntimeWarning: Mean of empty slice\n",
      "  np.nanmean(np.nanmean(monthly_array_nitrate_orig_slice, axis = 2),axis = 2)\n"
     ]
    }
   ],
   "source": [
    "monthly_array_nitrate_orig_slice[monthly_array_nitrate_orig_slice == 0 ] = np.nan\n",
    "monthly_array_nitrate_orig_slicemean = \\\n",
    "np.nanmean(np.nanmean(monthly_array_nitrate_orig_slice, axis = 2),axis = 2)\n",
    "print(np.shape(monthly_array_nitrate_orig_slicemean))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Experiment 9\n",
    "\n",
    "monthly_array_nitrate_exp_slice = np.zeros([14,12,50,50])\n",
    "# Load monthly averages\n",
    "mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc')\n",
    "slc = {'y': slice(450,500), 'x': slice(250,300)} \n",
    "e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')]\n",
    "years, variables = range(2007, 2021), ['nitrate']\n",
    "# Temporary list dict\n",
    "data = {}\n",
    "# Permanent aggregate dict\n",
    "aggregates = {var: {} for var in variables}\n",
    "monthlydat = {var: {} for var in variables}\n",
    "        \n",
    "# Add experiment year\n",
    "for year in [2019]:\n",
    "    # Initialize lists\n",
    "    for var in variables: data[var] = []\n",
    "    # Load monthly averages\n",
    "    for month in range(1, 7):\n",
    "        datestr = f'{year}{month:02d}'\n",
    "        prefix = f'/data/sallen/results/MEOPAR/Karyn/01jan19_river08/SalishSea_1m_{datestr}_{datestr}'\n",
    "        \n",
    "        \n",
    "        # Load grazing variables\n",
    "        with xr.open_dataset(prefix + '_ptrc_T.nc') as ds:\n",
    "            q = ds.nitrate.isel(deptht=0, **slc).values\n",
    "            q2 = q[0,:,:]\n",
    "            monthly_array_nitrate_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension\n",
    "            for var in ['nitrate']:\n",
    "                data[var].append(ds.nitrate.isel(deptht=0, **slc).values)\n",
    "                \n",
    "    # Concatenate months\n",
    "    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)   \n",
    "        \n",
    "\n",
    "# Add experiment year\n",
    "for year in [2019]:\n",
    "    # Initialize lists\n",
    "    for var in variables: data[var] = []\n",
    "    # Load monthly averages\n",
    "    for month in range(7, 10):\n",
    "        datestr = f'{year}{month:02d}'\n",
    "        prefix = f'/data/sallen/results/MEOPAR/Karyn/01jul19_river08/SalishSea_1m_{datestr}_{datestr}'\n",
    "        \n",
    "        \n",
    "        # Load grazing variables\n",
    "        with xr.open_dataset(prefix + '_ptrc_T.nc') as ds:\n",
    "            q = ds.nitrate.isel(deptht=0, **slc).values\n",
    "            q2 = q[0,:,:]\n",
    "            monthly_array_nitrate_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension\n",
    "            for var in ['nitrate']:\n",
    "                data[var].append(ds.nitrate.isel(deptht=0, **slc).values)\n",
    "                \n",
    "    # Concatenate months\n",
    "    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) \n",
    "\n",
    "        \n",
    "\n",
    "\n",
    "for year in [2019]:\n",
    "    # Initialize lists\n",
    "    for var in variables: data[var] = []\n",
    "    # Load monthly averages\n",
    "    for month in range(10, 13):\n",
    "        datestr = f'{year}{month:02d}'\n",
    "        prefix = f'/data/sallen/results/MEOPAR/Karyn/11oct19_river08/SalishSea_1m_{datestr}_{datestr}'\n",
    "        \n",
    "        \n",
    "        # Load grazing variables\n",
    "        with xr.open_dataset(prefix + '_ptrc_T.nc') as ds:\n",
    "            q = ds.nitrate.isel(deptht=0, **slc).values\n",
    "            q2 = q[0,:,:]\n",
    "            monthly_array_nitrate_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension\n",
    "            for var in ['nitrate']:\n",
    "                data[var].append(ds.nitrate.isel(deptht=0, **slc).values)\n",
    "                \n",
    "    # Concatenate months\n",
    "    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(14, 12)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_3795255/1010059431.py:3: RuntimeWarning: Mean of empty slice\n",
      "  np.nanmean(np.nanmean(monthly_array_nitrate_exp_slice, axis = 2),axis = 2)\n"
     ]
    }
   ],
   "source": [
    "monthly_array_nitrate_exp_slice[monthly_array_nitrate_exp_slice == 0 ] = np.nan\n",
    "monthly_array_nitrate_exp_slicemean = \\\n",
    "np.nanmean(np.nanmean(monthly_array_nitrate_exp_slice, axis = 2),axis = 2)\n",
    "print(np.shape(monthly_array_nitrate_exp_slicemean))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Original WY\n",
    "\n",
    "\n",
    "monthly_array_silicon_orig_slice = np.zeros([14,12,50,50])\n",
    "# Load monthly averages\n",
    "mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc')\n",
    "slc = {'y': slice(450,500), 'x': slice(250,300)} \n",
    "e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')]\n",
    "years, variables = range(2007, 2021), ['silicon']\n",
    "# Temporary list dict\n",
    "data = {}\n",
    "# Permanent aggregate dict\n",
    "aggregates = {var: {} for var in variables}\n",
    "monthlydat = {var: {} for var in variables}\n",
    "\n",
    "        \n",
    "# Add experiment year\n",
    "for year in [2019]:\n",
    "    # Initialize lists\n",
    "    for var in variables: data[var] = []\n",
    "    # Load monthly averages\n",
    "    for month in range(1, 13):\n",
    "        datestr = f'{year}{month:02d}'\n",
    "        prefix = f'/data/sallen/results/MEOPAR/v201905r/SalishSea_1m_{datestr}_{datestr}'\n",
    "        \n",
    "        \n",
    "        # Load grazing variables\n",
    "        with xr.open_dataset(prefix + '_ptrc_T.nc') as ds:\n",
    "            q = ds.silicon.isel(deptht=0, **slc).values\n",
    "            q2 = q[0,:,:]\n",
    "            monthly_array_silicon_orig_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension\n",
    "            for var in ['silicon']:\n",
    "                data[var].append(ds.silicon.isel(deptht=0, **slc).values)\n",
    "    \n",
    "    # Concatenate months\n",
    "    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)        \n",
    "\n",
    "       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(14, 12)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_3795255/241793216.py:3: RuntimeWarning: Mean of empty slice\n",
      "  np.nanmean(np.nanmean(monthly_array_silicon_orig_slice, axis = 2),axis = 2)\n"
     ]
    }
   ],
   "source": [
    "monthly_array_silicon_orig_slice[monthly_array_silicon_orig_slice == 0 ] = np.nan\n",
    "monthly_array_silicon_orig_slicemean = \\\n",
    "np.nanmean(np.nanmean(monthly_array_silicon_orig_slice, axis = 2),axis = 2)\n",
    "print(np.shape(monthly_array_silicon_orig_slicemean))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Experiment 9\n",
    "\n",
    "monthly_array_silicon_exp_slice = np.zeros([14,12,50,50])\n",
    "# Load monthly averages\n",
    "mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc')\n",
    "slc = {'y': slice(450,500), 'x': slice(250,300)} \n",
    "e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')]\n",
    "years, variables = range(2007, 2021), ['silicon']\n",
    "# Temporary list dict\n",
    "data = {}\n",
    "# Permanent aggregate dict\n",
    "aggregates = {var: {} for var in variables}\n",
    "monthlydat = {var: {} for var in variables}\n",
    "        \n",
    "# Add experiment year\n",
    "for year in [2019]:\n",
    "    # Initialize lists\n",
    "    for var in variables: data[var] = []\n",
    "    # Load monthly averages\n",
    "    for month in range(1, 7):\n",
    "        datestr = f'{year}{month:02d}'\n",
    "        prefix = f'/data/sallen/results/MEOPAR/Karyn/01jan19_river08/SalishSea_1m_{datestr}_{datestr}'\n",
    "        \n",
    "        \n",
    "        # Load grazing variables\n",
    "        with xr.open_dataset(prefix + '_ptrc_T.nc') as ds:\n",
    "            q = ds.silicon.isel(deptht=0, **slc).values\n",
    "            q2 = q[0,:,:]\n",
    "            monthly_array_silicon_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension\n",
    "            for var in ['silicon']:\n",
    "                data[var].append(ds.silicon.isel(deptht=0, **slc).values)\n",
    "                \n",
    "    # Concatenate months\n",
    "    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)   \n",
    "        \n",
    "\n",
    "# Add experiment year\n",
    "for year in [2019]:\n",
    "    # Initialize lists\n",
    "    for var in variables: data[var] = []\n",
    "    # Load monthly averages\n",
    "    for month in range(7, 10):\n",
    "        datestr = f'{year}{month:02d}'\n",
    "        prefix = f'/data/sallen/results/MEOPAR/Karyn/01jul19_river08/SalishSea_1m_{datestr}_{datestr}'\n",
    "        \n",
    "        \n",
    "        # Load grazing variables\n",
    "        with xr.open_dataset(prefix + '_ptrc_T.nc') as ds:\n",
    "            q = ds.silicon.isel(deptht=0, **slc).values\n",
    "            q2 = q[0,:,:]\n",
    "            monthly_array_silicon_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension\n",
    "            for var in ['silicon']:\n",
    "                data[var].append(ds.silicon.isel(deptht=0, **slc).values)\n",
    "                \n",
    "    # Concatenate months\n",
    "    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) \n",
    "\n",
    "\n",
    "for year in [2019]:\n",
    "    # Initialize lists\n",
    "    for var in variables: data[var] = []\n",
    "    # Load monthly averages\n",
    "    for month in range(10, 13):\n",
    "        datestr = f'{year}{month:02d}'\n",
    "        prefix = f'/data/sallen/results/MEOPAR/Karyn/11oct19_river08/SalishSea_1m_{datestr}_{datestr}'\n",
    "        \n",
    "        \n",
    "        # Load grazing variables\n",
    "        with xr.open_dataset(prefix + '_ptrc_T.nc') as ds:\n",
    "            q = ds.silicon.isel(deptht=0, **slc).values\n",
    "            q2 = q[0,:,:]\n",
    "            monthly_array_silicon_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension\n",
    "            for var in ['silicon']:\n",
    "                data[var].append(ds.silicon.isel(deptht=0, **slc).values)\n",
    "                \n",
    "    # Concatenate months\n",
    "    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(14, 12)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_3795255/4227105976.py:3: RuntimeWarning: Mean of empty slice\n",
      "  np.nanmean(np.nanmean(monthly_array_silicon_exp_slice, axis = 2),axis = 2)\n"
     ]
    }
   ],
   "source": [
    "monthly_array_silicon_exp_slice[monthly_array_silicon_exp_slice == 0 ] = np.nan\n",
    "monthly_array_silicon_exp_slicemean = \\\n",
    "np.nanmean(np.nanmean(monthly_array_silicon_exp_slice, axis = 2),axis = 2)\n",
    "print(np.shape(monthly_array_silicon_exp_slicemean))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "## Original WY\n",
    "\n",
    "monthly_array_ammonium_orig_slice = np.zeros([14,12,50,50])\n",
    "# Load monthly averages\n",
    "mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc')\n",
    "slc = {'y': slice(450,500), 'x': slice(250,300)} \n",
    "e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')]\n",
    "years, variables = range(2007, 2021), ['ammonium']\n",
    "# Temporary list dict\n",
    "data = {}\n",
    "# Permanent aggregate dict\n",
    "aggregates = {var: {} for var in variables}\n",
    "monthlydat = {var: {} for var in variables}\n",
    "\n",
    "        \n",
    "# Add experiment year\n",
    "for year in [2019]:\n",
    "    # Initialize lists\n",
    "    for var in variables: data[var] = []\n",
    "    # Load monthly averages\n",
    "    for month in range(1, 13):\n",
    "        datestr = f'{year}{month:02d}'\n",
    "        prefix = f'/data/sallen/results/MEOPAR/v201905r/SalishSea_1m_{datestr}_{datestr}'\n",
    "        \n",
    "        \n",
    "        # Load grazing variables\n",
    "        with xr.open_dataset(prefix + '_ptrc_T.nc') as ds:\n",
    "            q = ds.ammonium.isel(deptht=0, **slc).values\n",
    "            q2 = q[0,:,:]\n",
    "            monthly_array_ammonium_orig_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension\n",
    "            for var in ['ammonium']:\n",
    "                data[var].append(ds.ammonium.isel(deptht=0, **slc).values)\n",
    "    \n",
    "    # Concatenate months\n",
    "    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)        \n",
    "\n",
    "       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(14, 12)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_3795255/3976818894.py:3: RuntimeWarning: Mean of empty slice\n",
      "  np.nanmean(np.nanmean(monthly_array_ammonium_orig_slice, axis = 2),axis = 2)\n"
     ]
    }
   ],
   "source": [
    "monthly_array_ammonium_orig_slice[monthly_array_ammonium_orig_slice == 0 ] = np.nan\n",
    "monthly_array_ammonium_orig_slicemean = \\\n",
    "np.nanmean(np.nanmean(monthly_array_ammonium_orig_slice, axis = 2),axis = 2)\n",
    "print(np.shape(monthly_array_ammonium_orig_slicemean))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Experiment 9\n",
    "\n",
    "monthly_array_ammonium_exp_slice = np.zeros([14,12,50,50])\n",
    "# Load monthly averages\n",
    "mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc')\n",
    "slc = {'y': slice(450,500), 'x': slice(250,300)} \n",
    "e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')]\n",
    "years, variables = range(2007, 2021), ['ammonium']\n",
    "# Temporary list dict\n",
    "data = {}\n",
    "# Permanent aggregate dict\n",
    "aggregates = {var: {} for var in variables}\n",
    "monthlydat = {var: {} for var in variables}\n",
    "        \n",
    "# Add experiment year\n",
    "for year in [2019]:\n",
    "    # Initialize lists\n",
    "    for var in variables: data[var] = []\n",
    "    # Load monthly averages\n",
    "    for month in range(1, 7):\n",
    "        datestr = f'{year}{month:02d}'\n",
    "        prefix = f'/data/sallen/results/MEOPAR/Karyn/01jan19_river08/SalishSea_1m_{datestr}_{datestr}'\n",
    "        \n",
    "        \n",
    "        # Load grazing variables\n",
    "        with xr.open_dataset(prefix + '_ptrc_T.nc') as ds:\n",
    "            q = ds.ammonium.isel(deptht=0, **slc).values\n",
    "            q2 = q[0,:,:]\n",
    "            monthly_array_ammonium_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension\n",
    "            for var in ['ammonium']:\n",
    "                data[var].append(ds.ammonium.isel(deptht=0, **slc).values)\n",
    "                \n",
    "    # Concatenate months\n",
    "    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)   \n",
    "        \n",
    "\n",
    "# Add experiment year\n",
    "for year in [2019]:\n",
    "    # Initialize lists\n",
    "    for var in variables: data[var] = []\n",
    "    # Load monthly averages\n",
    "    for month in range(7, 10):\n",
    "        datestr = f'{year}{month:02d}'\n",
    "        prefix = f'/data/sallen/results/MEOPAR/Karyn/01jul19_river08/SalishSea_1m_{datestr}_{datestr}'\n",
    "        \n",
    "        \n",
    "        # Load grazing variables\n",
    "        with xr.open_dataset(prefix + '_ptrc_T.nc') as ds:\n",
    "            q = ds.ammonium.isel(deptht=0, **slc).values\n",
    "            q2 = q[0,:,:]\n",
    "            monthly_array_ammonium_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension\n",
    "            for var in ['ammonium']:\n",
    "                data[var].append(ds.ammonium.isel(deptht=0, **slc).values)\n",
    "                \n",
    "    # Concatenate months\n",
    "    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) \n",
    "\n",
    "for year in [2019]:\n",
    "    # Initialize lists\n",
    "    for var in variables: data[var] = []\n",
    "    # Load monthly averages\n",
    "    for month in range(10, 13):\n",
    "        datestr = f'{year}{month:02d}'\n",
    "        prefix = f'/data/sallen/results/MEOPAR/Karyn/11oct19_river08/SalishSea_1m_{datestr}_{datestr}'\n",
    "        \n",
    "        \n",
    "        # Load grazing variables\n",
    "        with xr.open_dataset(prefix + '_ptrc_T.nc') as ds:\n",
    "            q = ds.ammonium.isel(deptht=0, **slc).values\n",
    "            q2 = q[0,:,:]\n",
    "            monthly_array_ammonium_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension\n",
    "            for var in ['ammonium']:\n",
    "                data[var].append(ds.ammonium.isel(deptht=0, **slc).values)\n",
    "                \n",
    "    # Concatenate months\n",
    "    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(14, 12)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_3795255/3603250393.py:3: RuntimeWarning: Mean of empty slice\n",
      "  np.nanmean(np.nanmean(monthly_array_ammonium_exp_slice, axis = 2),axis = 2)\n"
     ]
    }
   ],
   "source": [
    "monthly_array_ammonium_exp_slice[monthly_array_ammonium_exp_slice == 0 ] = np.nan\n",
    "monthly_array_ammonium_exp_slicemean = \\\n",
    "np.nanmean(np.nanmean(monthly_array_ammonium_exp_slice, axis = 2),axis = 2)\n",
    "print(np.shape(monthly_array_ammonium_exp_slicemean))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for now just set light to constant and ignore 'limiter' and 'limval'\n",
    "DiatLimOrig, __, NanoLimOrig = bt.calc_p_limiters(10*np.ones(np.shape(monthly_array_nitrate_orig_slicemean[12,:])),\n",
    "                                               NO=monthly_array_nitrate_orig_slicemean[12,:],\n",
    "                                               NH=monthly_array_ammonium_orig_slicemean[12,:],\n",
    "                                               Si=monthly_array_silicon_orig_slicemean[12,:],\n",
    "                                               tmask=np.ones(np.shape(monthly_array_nitrate_orig_slicemean[12,:])),\n",
    "                                               nampisprod=nml)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for now just set light to constant and ignore 'limiter' and 'limval'\n",
    "DiatLimExp, __, NanoLimExp = bt.calc_p_limiters(10*np.ones(np.shape(monthly_array_nitrate_exp_slicemean[12,:])),\n",
    "                                               NO=monthly_array_nitrate_exp_slicemean[12,:],\n",
    "                                               NH=monthly_array_ammonium_exp_slicemean[12,:],\n",
    "                                               Si=monthly_array_silicon_exp_slicemean[12,:],\n",
    "                                               tmask=np.ones(np.shape(monthly_array_nitrate_exp_slicemean[12,:])),\n",
    "                                               nampisprod=nml)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'ILim': array([0.51559692, 0.51559692, 0.51559692, 0.51559692, 0.51559692,\n",
       "        0.51559692, 0.51559692, 0.51559692, 0.51559692, 0.51559692,\n",
       "        0.51559692, 0.51559692]),\n",
       " 'NLim': array([0.92177655, 0.91938623, 0.88324047, 0.71279276, 0.65085919,\n",
       "        0.39827465, 0.38265472, 0.40354223, 0.77676333, 0.89482032,\n",
       "        0.9016511 , 0.9135385 ]),\n",
       " 'SiLim': array([0.95638142, 0.95736831, 0.94794749, 0.87408713, 0.83030771,\n",
       "        0.91917568, 0.93623009, 0.93732263, 0.94506319, 0.95022164,\n",
       "        0.95311615, 0.95570312]),\n",
       " 'limiter': array([0., 0., 0., 0., 0., 2., 2., 2., 0., 0., 0., 0.]),\n",
       " 'limval': array([0.51559692, 0.51559692, 0.51559692, 0.51559692, 0.51559692,\n",
       "        2.39827465, 2.38265472, 2.40354223, 0.51559692, 0.51559692,\n",
       "        0.51559692, 0.51559692])}"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DiatLimOrig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Overall nutrient response factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "NutLimDiatOrig=np.where(DiatLimOrig['SiLim']<DiatLimOrig['NLim'],DiatLimOrig['SiLim'],DiatLimOrig['NLim'])\n",
    "NutLimNanoOrig=NanoLimOrig['NLim']\n",
    "NutLimDiatExp=np.where(DiatLimExp['SiLim']<DiatLimExp['NLim'],DiatLimExp['SiLim'],DiatLimExp['NLim'])\n",
    "NutLimNanoExp=NanoLimExp['NLim']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now calculate diatom to flagellate overall growth ratios (multiply by max growth) when temperature and nutrient responses are considered in isolation:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### * these rates would never occur in the model, even under repleat/optimal conditions, as the maximum of the nutrient and temperature responses are not equal to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "mumaxDiat=nml['zz_rate_r_diat']\n",
    "mumaxNano=nml['zz_rate_r_nano']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get values to use in summary plots in Figure8_TempNutDependence.ipynb notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.71312929, 2.70705032, 2.63109434, 2.25416726, 2.06368532,\n",
       "       1.61223576, 1.44798715, 1.44531944, 2.39809586, 2.65546763,\n",
       "       2.67459232, 2.69711616])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mumaxDiat/mumaxNano*NutLimDiatExp/NutLimNanoExp #'WY with CY Rivers Nutrient dependence only'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.86842105, 2.86842105, 2.86842105, 2.86842105, 2.45404626,\n",
       "       1.72678252, 1.37297695, 1.4480068 , 1.91139694, 2.86842105,\n",
       "       2.86842105, 2.86842105])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mumaxDiat/mumaxNano*TdepDiatExp/TdepNanoExp #'WY with CY Rivers Temperature dependence only'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.71312929, 2.70705032, 2.63109434, 2.25416726, 1.76556341,\n",
       "       0.97056202, 0.69308269, 0.72961129, 1.59799172, 2.65546763,\n",
       "       2.67459232, 2.69711616])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mumaxDiat/mumaxNano*TdepDiatExp/TdepNanoExp*NutLimDiatExp/NutLimNanoExp #'WY with CY Rivers both')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (py39)",
   "language": "python",
   "name": "py39"
  },
  "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.9.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}