{ "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_3771895/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 3\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", "### 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, 7):\n", " datestr = f'{year}{month:02d}'\n", " prefix = f'/data/sallen/results/MEOPAR/Karyn/01jan19_tsc/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, 13):\n", " datestr = f'{year}{month:02d}'\n", " prefix = f'/data/sallen/results/MEOPAR/Karyn/01jul19_tsc/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)" ] }, { "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_3771895/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.67567996, 5.16967481, 7.29056174, 10.17487312, 14.0136752 ,\n", " 17.62500597, 20.65200325, 20.13192256, 16.657524 , 11.60955676,\n", " 8.68565119, 6.48617353])" ] }, "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_3771895/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 3\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_NSi/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, 13):\n", " datestr = f'{year}{month:02d}'\n", " prefix = f'/data/sallen/results/MEOPAR/Karyn/01jul19_NSi/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" ] }, { "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_3771895/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_3771895/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 3\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_tsc/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, 13):\n", " datestr = f'{year}{month:02d}'\n", " prefix = f'/data/sallen/results/MEOPAR/Karyn/01jul19_tsc/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" ] }, { "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_3771895/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_3771895/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 3\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_tsc/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, 13):\n", " datestr = f'{year}{month:02d}'\n", " prefix = f'/data/sallen/results/MEOPAR/Karyn/01jul19_tsc/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" ] }, { "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_3771895/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']