#!/usr/bin/env python # coding: utf-8 # ## Data for Figure 8 # #### Based on Elise Olson's https://github.com/SalishSeaCast/analysis-elise-2/blob/master/notebooks/modelEqs/NewRateComparison.ipynb notebook # In[1]: import numpy as np import netCDF4 as nc import matplotlib.pyplot as plt from salishsea_tools import bio_tools as bt, places import xarray as xr import os import glob get_ipython().run_line_magic('matplotlib', 'inline') # In[2]: nml=bt.load_nml_bio(resDir='/ocean/eolson/MEOPAR/NEMO-3.6-code/NEMOGCM/CONFIG/SalishSeaCast/EXP00/', nmlname='nampisprod',bioRefName='namelist_smelt_cfg_HC201905equiv',bioCfgName='namelist_smelt_cfg_HC201905equiv') # In[3]: nml # In[4]: def phyto_Tdep_Factor(TT, zz_rate_maxtemp, zz_rate_temprange): if hasattr(TT,'__len__'): # assume 1-d array or similar and return array return np.array([phyto_Tdep_Factor(el,zz_rate_maxtemp, zz_rate_temprange) for el in TT]) else: return np.exp(0.07 * (TT - 20)) * min(max((zz_rate_maxtemp - TT), 0.0),zz_rate_temprange) / (zz_rate_temprange + 1e-10) # In[5]: def calc_T_Factors(TT,nampisprod): Tdep_Diat=phyto_Tdep_Factor(TT,nampisprod['zz_rate_maxtemp_diat'],nampisprod['zz_rate_temprange_diat']) Tdep_Myri=phyto_Tdep_Factor(TT,nampisprod['zz_rate_maxtemp_myri'],nampisprod['zz_rate_temprange_myri']) Tdep_Nano=phyto_Tdep_Factor(TT,nampisprod['zz_rate_maxtemp_nano'],nampisprod['zz_rate_temprange_nano']) return Tdep_Diat, Tdep_Myri, Tdep_Nano # ### Import Temperature Files # In[6]: # Original CY monthly_array_temp_orig_slice = np.zeros([14,12,50,50]) # Load monthly averages mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc') slc = {'y': slice(450,500), 'x': slice(250,300)} e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')] years, variables = range(2007, 2021), ['votemper'] # Temporary list dict data = {} # Permanent aggregate dict aggregates = {var: {} for var in variables} monthlydat = {var: {} for var in variables} ### 2008 using higher temperature threshold # Add experiment year for year in [2008]: # Initialize lists for var in variables: data[var] = [] # Load monthly averages for month in range(1, 13): datestr = f'{year}{month:02d}' prefix = f'/data/sallen/results/MEOPAR/v201905r/SalishSea_1m_{datestr}_{datestr}' with xr.open_dataset(prefix + '_grid_T.nc') as ds: q = ds.votemper.isel(deptht=0, **slc).values q2 = q[0,:,:] monthly_array_temp_orig_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension for var in ['votemper']: data[var].append(ds.votemper.isel(deptht=0, **slc).values) # Concatenate months for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) # In[7]: monthly_array_temp_orig_slice[monthly_array_temp_orig_slice == 0 ] = np.nan monthly_array_temp_orig_slicemean = \ np.nanmean(np.nanmean(monthly_array_temp_orig_slice, axis = 2),axis = 2) print(np.shape(monthly_array_temp_orig_slicemean)) # In[8]: ## Experiment 10 monthly_array_temp_exp_slice = np.zeros([14,12,50,50]) # Load monthly averages mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc') slc = {'y': slice(450,500), 'x': slice(250,300)} e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')] years, variables = range(2007, 2021), ['votemper'] # Temporary list dict data = {} # Permanent aggregate dict aggregates = {var: {} for var in variables} monthlydat = {var: {} for var in variables} ### 2008 using higher temperature threshold # Add experiment year for year in [2008]: # Initialize lists for var in variables: data[var] = [] # Load monthly averages for month in range(1, 7): datestr = f'{year}{month:02d}' prefix = f'/data/sallen/results/MEOPAR/Karyn/01jan08_river19/SalishSea_1m_{datestr}_{datestr}' with xr.open_dataset(prefix + '_grid_T.nc') as ds: q = ds.votemper.isel(deptht=0, **slc).values q2 = q[0,:,:] monthly_array_temp_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension for var in ['votemper']: data[var].append(ds.votemper.isel(deptht=0, **slc).values) # Concatenate months for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) for year in [2008]: # Initialize lists for var in variables: data[var] = [] # Load monthly averages for month in range(7, 13): datestr = f'{year}{month:02d}' prefix = f'/data/sallen/results/MEOPAR/Karyn/01jul08_river19/SalishSea_1m_{datestr}_{datestr}' with xr.open_dataset(prefix + '_grid_T.nc') as ds: q = ds.votemper.isel(deptht=0, **slc).values q2 = q[0,:,:] monthly_array_temp_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension for var in ['votemper']: data[var].append(ds.votemper.isel(deptht=0, **slc).values) # Concatenate months for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) # In[9]: monthly_array_temp_exp_slice[monthly_array_temp_exp_slice == 0 ] = np.nan monthly_array_temp_exp_slicemean = \ np.nanmean(np.nanmean(monthly_array_temp_exp_slice, axis = 2),axis = 2) print(np.shape(monthly_array_temp_exp_slicemean)) # In[10]: monthly_array_temp_exp_slicemean[1,:] ## Experiment 2 SST # In[11]: monthly_array_temp_orig_slicemean[1,:] ## Original SST # In[ ]: # ### Temperature response factors # # In[12]: TdepDiatOrig,__,TdepNanoOrig=calc_T_Factors(monthly_array_temp_orig_slicemean[1,:],nml) TdepDiatExp,__,TdepNanoExp=calc_T_Factors(monthly_array_temp_exp_slicemean[1,:],nml) # #### Open Nutrient files # In[13]: ## Original CY monthly_array_nitrate_orig_slice = np.zeros([14,12,50,50]) # Load monthly averages mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc') slc = {'y': slice(450,500), 'x': slice(250,300)} e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')] years, variables = range(2007, 2021), ['nitrate'] # Temporary list dict data = {} # Permanent aggregate dict aggregates = {var: {} for var in variables} monthlydat = {var: {} for var in variables} # Add experiment year for year in [2008]: # Initialize lists for var in variables: data[var] = [] # Load monthly averages for month in range(1, 13): datestr = f'{year}{month:02d}' prefix = f'/data/sallen/results/MEOPAR/v201905r/SalishSea_1m_{datestr}_{datestr}' # Load grazing variables with xr.open_dataset(prefix + '_ptrc_T.nc') as ds: q = ds.nitrate.isel(deptht=0, **slc).values q2 = q[0,:,:] monthly_array_nitrate_orig_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension for var in ['nitrate']: data[var].append(ds.nitrate.isel(deptht=0, **slc).values) # Concatenate months for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) # In[14]: monthly_array_nitrate_orig_slice[monthly_array_nitrate_orig_slice == 0 ] = np.nan monthly_array_nitrate_orig_slicemean = \ np.nanmean(np.nanmean(monthly_array_nitrate_orig_slice, axis = 2),axis = 2) print(np.shape(monthly_array_nitrate_orig_slicemean)) # In[15]: ## Experiment 10 monthly_array_nitrate_exp_slice = np.zeros([14,12,50,50]) # Load monthly averages mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc') slc = {'y': slice(450,500), 'x': slice(250,300)} e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')] years, variables = range(2007, 2021), ['nitrate'] # Temporary list dict data = {} # Permanent aggregate dict aggregates = {var: {} for var in variables} monthlydat = {var: {} for var in variables} # Add experiment year for year in [2008]: # Initialize lists for var in variables: data[var] = [] # Load monthly averages for month in range(1, 7): datestr = f'{year}{month:02d}' prefix = f'/data/sallen/results/MEOPAR/Karyn/01jan08_river19/SalishSea_1m_{datestr}_{datestr}' # Load grazing variables with xr.open_dataset(prefix + '_ptrc_T.nc') as ds: q = ds.nitrate.isel(deptht=0, **slc).values q2 = q[0,:,:] monthly_array_nitrate_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension for var in ['nitrate']: data[var].append(ds.nitrate.isel(deptht=0, **slc).values) # Concatenate months for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) # Add experiment year for year in [2008]: # Initialize lists for var in variables: data[var] = [] # Load monthly averages for month in range(7, 13): datestr = f'{year}{month:02d}' prefix = f'/data/sallen/results/MEOPAR/Karyn/01jul08_river19/SalishSea_1m_{datestr}_{datestr}' # Load grazing variables with xr.open_dataset(prefix + '_ptrc_T.nc') as ds: q = ds.nitrate.isel(deptht=0, **slc).values q2 = q[0,:,:] monthly_array_nitrate_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension for var in ['nitrate']: data[var].append(ds.nitrate.isel(deptht=0, **slc).values) # Concatenate months for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) # In[16]: monthly_array_nitrate_exp_slice[monthly_array_nitrate_exp_slice == 0 ] = np.nan monthly_array_nitrate_exp_slicemean = \ np.nanmean(np.nanmean(monthly_array_nitrate_exp_slice, axis = 2),axis = 2) print(np.shape(monthly_array_nitrate_exp_slicemean)) # In[17]: ## Original CY monthly_array_silicon_orig_slice = np.zeros([14,12,50,50]) # Load monthly averages mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc') slc = {'y': slice(450,500), 'x': slice(250,300)} e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')] years, variables = range(2007, 2021), ['silicon'] # Temporary list dict data = {} # Permanent aggregate dict aggregates = {var: {} for var in variables} monthlydat = {var: {} for var in variables} # Add experiment year for year in [2008]: # Initialize lists for var in variables: data[var] = [] # Load monthly averages for month in range(1, 13): datestr = f'{year}{month:02d}' prefix = f'/data/sallen/results/MEOPAR/v201905r/SalishSea_1m_{datestr}_{datestr}' # Load grazing variables with xr.open_dataset(prefix + '_ptrc_T.nc') as ds: q = ds.silicon.isel(deptht=0, **slc).values q2 = q[0,:,:] monthly_array_silicon_orig_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension for var in ['silicon']: data[var].append(ds.silicon.isel(deptht=0, **slc).values) # Concatenate months for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) # In[18]: monthly_array_silicon_orig_slice[monthly_array_silicon_orig_slice == 0 ] = np.nan monthly_array_silicon_orig_slicemean = \ np.nanmean(np.nanmean(monthly_array_silicon_orig_slice, axis = 2),axis = 2) print(np.shape(monthly_array_silicon_orig_slicemean)) # In[19]: ## Experiment 10 monthly_array_silicon_exp_slice = np.zeros([14,12,50,50]) # Load monthly averages mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc') slc = {'y': slice(450,500), 'x': slice(250,300)} e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')] years, variables = range(2007, 2021), ['silicon'] # Temporary list dict data = {} # Permanent aggregate dict aggregates = {var: {} for var in variables} monthlydat = {var: {} for var in variables} # Add experiment year for year in [2008]: # Initialize lists for var in variables: data[var] = [] # Load monthly averages for month in range(1, 7): datestr = f'{year}{month:02d}' prefix = f'/data/sallen/results/MEOPAR/Karyn/01jan08_river19/SalishSea_1m_{datestr}_{datestr}' # Load grazing variables with xr.open_dataset(prefix + '_ptrc_T.nc') as ds: q = ds.silicon.isel(deptht=0, **slc).values q2 = q[0,:,:] monthly_array_silicon_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension for var in ['silicon']: data[var].append(ds.silicon.isel(deptht=0, **slc).values) # Concatenate months for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) # Add experiment year for year in [2008]: # Initialize lists for var in variables: data[var] = [] # Load monthly averages for month in range(7, 13): datestr = f'{year}{month:02d}' prefix = f'/data/sallen/results/MEOPAR/Karyn/01jul08_river19/SalishSea_1m_{datestr}_{datestr}' # Load grazing variables with xr.open_dataset(prefix + '_ptrc_T.nc') as ds: q = ds.silicon.isel(deptht=0, **slc).values q2 = q[0,:,:] monthly_array_silicon_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension for var in ['silicon']: data[var].append(ds.silicon.isel(deptht=0, **slc).values) # Concatenate months for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) # In[20]: monthly_array_silicon_exp_slice[monthly_array_silicon_exp_slice == 0 ] = np.nan monthly_array_silicon_exp_slicemean = \ np.nanmean(np.nanmean(monthly_array_silicon_exp_slice, axis = 2),axis = 2) print(np.shape(monthly_array_silicon_exp_slicemean)) # In[21]: ## Original CY monthly_array_ammonium_orig_slice = np.zeros([14,12,50,50]) # Load monthly averages mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc') slc = {'y': slice(450,500), 'x': slice(250,300)} e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')] years, variables = range(2007, 2021), ['ammonium'] # Temporary list dict data = {} # Permanent aggregate dict aggregates = {var: {} for var in variables} monthlydat = {var: {} for var in variables} # Add experiment year for year in [2008]: # Initialize lists for var in variables: data[var] = [] # Load monthly averages for month in range(1, 13): datestr = f'{year}{month:02d}' prefix = f'/data/sallen/results/MEOPAR/v201905r/SalishSea_1m_{datestr}_{datestr}' # Load grazing variables with xr.open_dataset(prefix + '_ptrc_T.nc') as ds: q = ds.ammonium.isel(deptht=0, **slc).values q2 = q[0,:,:] monthly_array_ammonium_orig_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension for var in ['ammonium']: data[var].append(ds.ammonium.isel(deptht=0, **slc).values) # Concatenate months for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) # In[22]: monthly_array_ammonium_orig_slice[monthly_array_ammonium_orig_slice == 0 ] = np.nan monthly_array_ammonium_orig_slicemean = \ np.nanmean(np.nanmean(monthly_array_ammonium_orig_slice, axis = 2),axis = 2) print(np.shape(monthly_array_ammonium_orig_slicemean)) # In[23]: ## Experiment 10 monthly_array_ammonium_exp_slice = np.zeros([14,12,50,50]) # Load monthly averages mask = xr.open_dataset('/data/eolson/results/MEOPAR/NEMO-forcing-new/grid/mesh_mask201702.nc') slc = {'y': slice(450,500), 'x': slice(250,300)} e3t, tmask = [mask[var].isel(z=slice(None, 27),**slc).values for var in ('e3t_0', 'tmask')] years, variables = range(2007, 2021), ['ammonium'] # Temporary list dict data = {} # Permanent aggregate dict aggregates = {var: {} for var in variables} monthlydat = {var: {} for var in variables} # Add experiment year for year in [2008]: # Initialize lists for var in variables: data[var] = [] # Load monthly averages for month in range(1, 7): datestr = f'{year}{month:02d}' prefix = f'/data/sallen/results/MEOPAR/Karyn/01jan08_river19/SalishSea_1m_{datestr}_{datestr}' # Load grazing variables with xr.open_dataset(prefix + '_ptrc_T.nc') as ds: q = ds.ammonium.isel(deptht=0, **slc).values q2 = q[0,:,:] monthly_array_ammonium_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension for var in ['ammonium']: data[var].append(ds.ammonium.isel(deptht=0, **slc).values) # Concatenate months for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) # Add experiment year for year in [2008]: # Initialize lists for var in variables: data[var] = [] # Load monthly averages for month in range(7, 13): datestr = f'{year}{month:02d}' prefix = f'/data/sallen/results/MEOPAR/Karyn/01jul08_river19/SalishSea_1m_{datestr}_{datestr}' # Load grazing variables with xr.open_dataset(prefix + '_ptrc_T.nc') as ds: q = ds.ammonium.isel(deptht=0, **slc).values q2 = q[0,:,:] monthly_array_ammonium_exp_slice[year-2007,month-1,:,:] = q2 #year2007 is index 0 along 1st dimension for var in ['ammonium']: data[var].append(ds.ammonium.isel(deptht=0, **slc).values) # Concatenate months for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0) # In[24]: monthly_array_ammonium_exp_slice[monthly_array_ammonium_exp_slice == 0 ] = np.nan monthly_array_ammonium_exp_slicemean = \ np.nanmean(np.nanmean(monthly_array_ammonium_exp_slice, axis = 2),axis = 2) print(np.shape(monthly_array_ammonium_exp_slicemean)) # In[25]: # for now just set light to constant and ignore 'limiter' and 'limval' DiatLimOrig, __, NanoLimOrig = bt.calc_p_limiters(10*np.ones(np.shape(monthly_array_nitrate_orig_slicemean[1,:])), NO=monthly_array_nitrate_orig_slicemean[1,:], NH=monthly_array_ammonium_orig_slicemean[1,:], Si=monthly_array_silicon_orig_slicemean[1,:], tmask=np.ones(np.shape(monthly_array_nitrate_orig_slicemean[1,:])), nampisprod=nml) # In[26]: # for now just set light to constant and ignore 'limiter' and 'limval' DiatLimExp, __, NanoLimExp = bt.calc_p_limiters(10*np.ones(np.shape(monthly_array_nitrate_exp_slicemean[1,:])), NO=monthly_array_nitrate_exp_slicemean[1,:], NH=monthly_array_ammonium_exp_slicemean[1,:], Si=monthly_array_silicon_exp_slicemean[1,:], tmask=np.ones(np.shape(monthly_array_nitrate_exp_slicemean[1,:])), nampisprod=nml) # In[27]: DiatLimOrig # In[28]: NutLimDiatOrig=np.where(DiatLimOrig['SiLim']