#!/usr/bin/env python
# coding: utf-8

# In[1]:


import numpy as np
import xarray as xr
import pandas as pd
from salishsea_tools import viz_tools, places, visualisations
from matplotlib import pyplot as plt, dates
from datetime import datetime, timedelta
from calendar import month_name
from scipy.io import loadmat
from tqdm.notebook import tqdm
from salishsea_tools import nc_tools
from dask.diagnostics import ProgressBar
import cmocean

get_ipython().run_line_magic('matplotlib', 'inline')


# #### Note that all biological rates (e.g., grazing) are scaled (multiplied by 1.111) to correct for a time splitting/Roberts-Asselin filter error in hte model as reported in Olson et al. 2020
from IPython.display import HTML

HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>

<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
# In[2]:


plt.rcParams.update({'font.size': 12, 'axes.titlesize': 'medium'})


# # Load files from monthly averages

# ## Z1 grazing on Diatoms

# In[3]:


#years, months, data
monthly_array_Z1diatoms_depthint_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), ['GRMICZDIAT']
# Temporary list dict
data = {}
# Permanent aggregate dict
#aggregates = {var: {} for var in variables}
#monthlydat = {var: {} for var in variables}

# Loop through years
for year in [2007,2008,2009,2010,2011,2012,2015,2016,2017,2018,2019,2020]:
    # 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 + '_dia2_T.nc') as ds:
            q = np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data
            q2 = q[0,:,:]
            monthly_array_Z1diatoms_depthint_slice[year-2007,month-1,:,:] = q2 #year2015 is index 0 along 1st dimension
            for var in ['GRMICZDIAT']:
                data[var].append(np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data)
    
    
# Loop through years for wrap files
for year in [2013,2014]:
    # 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_wrap/SalishSea_1m_{datestr}_{datestr}'
        
         # Load grazing variables
        with xr.open_dataset(prefix + '_dia2_T.nc') as ds:
            q = np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data
            q2 = q[0,:,:]
            monthly_array_Z1diatoms_depthint_slice[year-2007,month-1,:,:] = q2 #year2015 is index 0 along 1st dimension
            for var in ['GRMICZDIAT']:
                data[var].append(np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data)
    
    
    
    
    
    
    
    # Concatenate months
    #for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)
# # Calculate 5 year mean and anomalies
# for var in variables:
#     aggregates[var][‘mean’] = np.concatenate([aggregates[var][year][None, ...] for year in years]).mean(axis=0)
#     for year in years: aggregates[var][year] = aggregates[var][year] - aggregates[var][‘mean’]









# In[4]:


monthly_array_Z1diatoms_depthint_slice=np.where(np.isnan(monthly_array_Z1diatoms_depthint_slice), 0, monthly_array_Z1diatoms_depthint_slice)
monthly_array_Z1diatoms_depthint_slicemean = \
np.nanmean(np.nanmean(monthly_array_Z1diatoms_depthint_slice, axis = 2),axis = 2)
print(np.shape(monthly_array_Z1diatoms_depthint_slicemean))


# ### Select 4 warmest and 4 coldest years; leave NPGO "neutral" years out

# In[5]:


#2008, 2010, 2011, 2012
NPGO_C_Z1diat=np.array([[monthly_array_Z1diatoms_depthint_slicemean[1,:]],[monthly_array_Z1diatoms_depthint_slicemean[3,:]],\
        [monthly_array_Z1diatoms_depthint_slicemean[4,:]],[monthly_array_Z1diatoms_depthint_slicemean[5,:]]])


# In[6]:


NPGO_C_Z1diat_mean=NPGO_C_Z1diat.mean(axis=0).flatten()*86400*1.111


# In[7]:


NPGO_C_Z1diat_std=NPGO_C_Z1diat.std(axis=0).flatten()*86400*1.111


# In[8]:


#2015, 2018, 2019, 2020
NPGO_W_Z1diat=np.array([[monthly_array_Z1diatoms_depthint_slicemean[8,:]],[monthly_array_Z1diatoms_depthint_slicemean[11,:]],\
        [monthly_array_Z1diatoms_depthint_slicemean[12,:]],[monthly_array_Z1diatoms_depthint_slicemean[13,:]]])


# In[9]:


NPGO_W_Z1diat_mean=NPGO_W_Z1diat.mean(axis=0).flatten()*86400*1.111


# In[10]:


NPGO_W_Z1diat_std=NPGO_W_Z1diat.std(axis=0).flatten()*86400*1.111


# In[11]:


NPGO_C_Z1diat_mean.max()


# In[12]:


NPGO_W_Z1diat_mean.max()


# # Load files from monthly averages

# ## Z1 grazing on Nanoflagellates

# In[13]:


#years, months, data
monthly_array_Z1flag_depthint_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), ['GRMICZPHY']
# Temporary list dict
data = {}
# Permanent aggregate dict
#aggregates = {var: {} for var in variables}
#monthlydat = {var: {} for var in variables}

# Loop through years
for year in [2007,2008,2009,2010,2011,2012,2015,2016,2017,2018,2019,2020]:
    # 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 + '_dia2_T.nc') as ds:
            q = np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data
            q2 = q[0,:,:]
            monthly_array_Z1flag_depthint_slice[year-2007,month-1,:,:] = q2 #year2015 is index 0 along 1st dimension
            for var in ['GRMICZPHY']:
                data[var].append(np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data)
    
# Loop through years for wrap files
for year in [2013,2014]:
    # 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_wrap/SalishSea_1m_{datestr}_{datestr}'
        
        # Load grazing variables
        with xr.open_dataset(prefix + '_dia2_T.nc') as ds:
            q = np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data
            q2 = q[0,:,:]
            monthly_array_Z1flag_depthint_slice[year-2007,month-1,:,:] = q2 #year2015 is index 0 along 1st dimension
            for var in ['GRMICZPHY']:
                data[var].append(np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data)
    
    
    
    
    # Concatenate months
    #for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)
# # Calculate 5 year mean and anomalies
# for var in variables:
#     aggregates[var][‘mean’] = np.concatenate([aggregates[var][year][None, ...] for year in years]).mean(axis=0)
#     for year in years: aggregates[var][year] = aggregates[var][year] - aggregates[var][‘mean’]







# In[14]:


monthly_array_Z1flag_depthint_slice[monthly_array_Z1flag_depthint_slice == 0 ] = np.nan
monthly_array_Z1flag_depthint_slicemean = \
np.nanmean(np.nanmean(monthly_array_Z1flag_depthint_slice, axis = 2),axis = 2)
print(np.shape(monthly_array_Z1flag_depthint_slicemean))


# ### Select 4 warmest and 4 coldest years; leave NPGO "neutral" years out

# In[15]:


#2008, 2010, 2011, 2012
NPGO_C_Z1flag=np.array([[monthly_array_Z1flag_depthint_slicemean[1,:]],[monthly_array_Z1flag_depthint_slicemean[3,:]],\
        [monthly_array_Z1flag_depthint_slicemean[4,:]],[monthly_array_Z1flag_depthint_slicemean[5,:]]])


# In[16]:


NPGO_C_Z1flag_mean=NPGO_C_Z1flag.mean(axis=0).flatten()*86400*1.111


# In[17]:


NPGO_C_Z1flag_std=NPGO_C_Z1flag.std(axis=0).flatten()*86400*1.111


# In[18]:


#2015, 2018, 2019, 2020
NPGO_W_Z1flag=np.array([[monthly_array_Z1flag_depthint_slicemean[8,:]],[monthly_array_Z1flag_depthint_slicemean[11,:]],\
        [monthly_array_Z1flag_depthint_slicemean[12,:]],[monthly_array_Z1flag_depthint_slicemean[13,:]]])


# In[19]:


NPGO_W_Z1flag_mean=NPGO_W_Z1flag.mean(axis=0).flatten()*86400*1.111


# In[20]:


NPGO_W_Z1flag_std=NPGO_W_Z1flag.std(axis=0).flatten()*86400*1.111


# In[21]:


NPGO_C_Z1flag_mean.max()


# In[22]:


NPGO_W_Z1flag_mean.max()


# #### Z1 Grazing on both diatoms and nanoflagelles for Cold Years

# In[23]:


NPGO_C_Z1Both=NPGO_C_Z1diat_mean+NPGO_C_Z1flag_mean


# In[24]:


NPGO_C_Z1Both


# In[25]:


#### Z1 Grazing on both diatoms and nanoflagelles for Warm Years


# In[26]:


NPGO_W_Z1Both=NPGO_W_Z1diat_mean+NPGO_W_Z1flag_mean


# In[27]:


NPGO_W_Z1Both


# In[28]:


NPGO_C_Z1Both.mean()


# In[29]:


NPGO_W_Z1Both.mean()


# ### Load monthly average files for Z1 biomass

# In[30]:


#years, months, data
monthly_array_microzooplankton_depthint_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), ['microzooplankton']
# Temporary list dict
data = {}
# Permanent aggregate dict
aggregates = {var: {} for var in variables}
monthlydat = {var: {} for var in variables}

# Loop through years
for year in [2007,2008,2009,2010,2011,2012,2016,2017,2018,2019,2020]:
    # 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 = np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data
            q2 = q[0,:,:]
            monthly_array_microzooplankton_depthint_slice[year-2007,month-1,:,:] = q2 #year2015 is index 0 along 1st dimension
            for var in ['microzooplankton']:
                data[var].append(np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data)
    # Concatenate months
    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)

# Loop through years for wrap files
for year in [2013,2014,2015]:
    # 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_wrap/SalishSea_1m_{datestr}_{datestr}'
        # Load grazing variables
        with xr.open_dataset(prefix + '_ptrc_T.nc') as ds:
            q = np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data
            q2 = q[0,:,:]
            monthly_array_microzooplankton_depthint_slice[year-2007,month-1,:,:] = q2 #year2015 is index 0 along 1st dimension
            for var in ['microzooplankton']:
                data[var].append(np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data)
    # Concatenate months
    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)      
        
        
# # Calculate 5 year mean and anomalies
# for var in variables:
#     aggregates[var][‘mean’] = np.concatenate([aggregates[var][year][None, ...] for year in years]).mean(axis=0)
#     for year in years: aggregates[var][year] = aggregates[var][year] - aggregates[var][‘mean’]




# In[31]:


monthly_array_microzooplankton_depthint_slice[monthly_array_microzooplankton_depthint_slice == 0 ] = np.nan
monthly_array_microzooplankton_depthint_slicemean = \
np.nanmean(np.nanmean(monthly_array_microzooplankton_depthint_slice, axis = 2),axis = 2)
print(np.shape(monthly_array_microzooplankton_depthint_slicemean))


# ### Select 4 warmest and 4 coldest years; leave NPGO "neutral" years out

# In[32]:


#2008, 2010, 2011, 2012
NPGO_C_micro=np.array([[monthly_array_microzooplankton_depthint_slicemean[1,:]],[monthly_array_microzooplankton_depthint_slicemean[3,:]],\
        [monthly_array_microzooplankton_depthint_slicemean[4,:]],[monthly_array_microzooplankton_depthint_slicemean[5,:]]])


# In[33]:


NPGO_C_micro_mean=NPGO_C_micro.mean(axis=0).flatten()*5.7*12/1000


# In[34]:


NPGO_C_micro_std=NPGO_C_micro.std(axis=0).flatten()*5.7*12/1000


# In[35]:


#2015, 2018, 2019, 2020
NPGO_W_micro=np.array([[monthly_array_microzooplankton_depthint_slicemean[8,:]],[monthly_array_microzooplankton_depthint_slicemean[11,:]],\
        [monthly_array_microzooplankton_depthint_slicemean[12,:]],[monthly_array_microzooplankton_depthint_slicemean[13,:]]])


# In[36]:


NPGO_W_micro_mean=NPGO_W_micro.mean(axis=0).flatten()*5.7*12/1000


# In[37]:


NPGO_W_micro_std=NPGO_W_micro.std(axis=0).flatten()*5.7*12/1000


# In[ ]:





# In[38]:


## Plot Figure 4a

fig, ax = plt.subplots(figsize=(7,5))
bbox = {'boxstyle': 'round', 'facecolor': 'w', 'alpha': 0.9}
cmap = plt.get_cmap('tab10')
palette = [cmap(0), cmap(0.2), 'k', cmap(0.1), cmap(0.3)]
xticks=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov',"Dec"]

ax2=ax.twinx()

ax2.plot(xticks, NPGO_C_micro_mean,label='Z1',color='b',linestyle='--',linewidth=3)


ax.bar(xticks,NPGO_C_Z1diat_mean,color='grey',edgecolor='k',hatch='//',alpha=0.7,label='Diatoms')
ax.bar(xticks,NPGO_C_Z1flag_mean,color='grey',edgecolor='k',label='Nanoflagellates',alpha=0.5,bottom=NPGO_C_Z1diat_mean)



ax2.set_ylim(0,5)
ax2.set_ylabel('Biomass (g C m$^{-2}$)',fontsize=14)

ax.legend(frameon=False,loc=2)
ax2.legend(frameon=False,loc=1)
ax.set_ylim(0,12)
ax.set_ylabel('Grazing (mmol N m$^{-2}$ d$^{-1}$)',fontsize=14)
ax.set_xlabel('',fontsize=14)
#ax.tick_params(axis='x', labelrotation=0)

ax.text(-.5, 12.5, '(a)', fontsize=15,  color='k')

ax.set_title('Cold Years',fontsize=18)

#plt.savefig('Fig6a.png', bbox_inches='tight',dpi=1000,transparent=False)


# In[39]:


## Plot Figure 4b

fig, ax = plt.subplots(figsize=(7,5))
bbox = {'boxstyle': 'round', 'facecolor': 'w', 'alpha': 0.9}
cmap = plt.get_cmap('tab10')
palette = [cmap(0), cmap(0.2), 'k', cmap(0.1), cmap(0.3)]
xticks=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov',"Dec"]

ax2=ax.twinx()

ax2.plot(xticks, NPGO_W_micro_mean,label='Z1',color='r',linestyle='--',linewidth=3)


ax.bar(xticks,NPGO_W_Z1diat_mean,color='grey',edgecolor='k',hatch='//',alpha=0.7,label='Diatoms')
ax.bar(xticks,NPGO_W_Z1flag_mean,color='grey',edgecolor='k',label='Nanoflagellates',alpha=0.5,bottom=NPGO_W_Z1diat_mean)



ax2.set_ylim(0,5)
ax2.set_ylabel('Biomass (g C m$^{-2}$)',fontsize=14)

ax2.legend(frameon=False,loc=1)
ax.set_ylim(0,12)
ax.set_ylabel('Grazing (mmol N m$^{-2}$ d$^{-1}$)',fontsize=14)
ax.set_xlabel('',fontsize=14)

#ax.tick_params(axis='x', labelrotation=0)
ax.text(-.5, 12.5, '(b)', fontsize=15,  color='k')

ax.set_title('Warm Years',fontsize=18)

#plt.savefig('Fig6b.png', bbox_inches='tight',dpi=1000,transparent=False)


# In[ ]:





# In[ ]:





# ## Z2 grazing on Diatoms

# In[40]:


#years, months, data
monthly_array_Z2diatoms_depthint_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), ['GRMESZDIAT']
# Temporary list dict
data = {}
# Permanent aggregate dict
#aggregates = {var: {} for var in variables}
#monthlydat = {var: {} for var in variables}

# Loop through years
for year in [2007,2008,2009,2010,2011,2012,2015,2016,2017,2018,2019,2020]:
    # 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 + '_dia2_T.nc') as ds:
            q = np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data
            q2 = q[0,:,:]
            monthly_array_Z2diatoms_depthint_slice[year-2007,month-1,:,:] = q2 #year2015 is index 0 along 1st dimension
            for var in ['GRMESZDIAT']:
                data[var].append(np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data)
    
# Loop through years for wrap files
for year in [2013,2014]:
    # 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_wrap/SalishSea_1m_{datestr}_{datestr}'
        
        # Load grazing variables
        with xr.open_dataset(prefix + '_dia2_T.nc') as ds:
            q = np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data
            q2 = q[0,:,:]
            monthly_array_Z2diatoms_depthint_slice[year-2007,month-1,:,:] = q2 #year2015 is index 0 along 1st dimension
            for var in ['GRMESZDIAT']:
                data[var].append(np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data)

    
    
    
    
    # Concatenate months
    #for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)
# # Calculate 5 year mean and anomalies
# for var in variables:
#     aggregates[var][‘mean’] = np.concatenate([aggregates[var][year][None, ...] for year in years]).mean(axis=0)
#     for year in years: aggregates[var][year] = aggregates[var][year] - aggregates[var][‘mean’]






# In[41]:


monthly_array_Z2diatoms_depthint_slice[monthly_array_Z2diatoms_depthint_slice == 0 ] = np.nan
monthly_array_Z2diatoms_depthint_slicemean = \
np.nanmean(np.nanmean(monthly_array_Z2diatoms_depthint_slice, axis = 2),axis = 2)
print(np.shape(monthly_array_Z2diatoms_depthint_slicemean))


# ### Select 4 warmest and 4 coldest years; leave NPGO "neutral" years out

# In[42]:


#2008, 2010, 2011, 2012
NPGO_C_Z2diat=np.array([[monthly_array_Z2diatoms_depthint_slicemean[1,:]],[monthly_array_Z2diatoms_depthint_slicemean[3,:]],\
        [monthly_array_Z2diatoms_depthint_slicemean[4,:]],[monthly_array_Z2diatoms_depthint_slicemean[5,:]]])


# In[43]:


NPGO_C_Z2diat_mean=NPGO_C_Z2diat.mean(axis=0).flatten()*86400*1.111


# In[44]:


NPGO_C_Z2diat_std=NPGO_C_Z2diat.std(axis=0).flatten()*86400*1.111


# In[45]:


#2015, 2018, 2019, 2020
NPGO_W_Z2diat=np.array([[monthly_array_Z2diatoms_depthint_slicemean[8,:]],[monthly_array_Z2diatoms_depthint_slicemean[11,:]],\
        [monthly_array_Z2diatoms_depthint_slicemean[12,:]],[monthly_array_Z2diatoms_depthint_slicemean[13,:]]])


# In[46]:


NPGO_W_Z2diat_mean=NPGO_W_Z2diat.mean(axis=0).flatten()*86400*1.111


# In[47]:


NPGO_W_Z2diat_std=NPGO_W_Z2diat.std(axis=0).flatten()*86400*1.111


# In[48]:


NPGO_C_Z2diat_mean.max()


# In[49]:


NPGO_W_Z2diat_mean.max()


# ## Z2 grazing on Nanoflagellates

# In[50]:


#years, months, data
monthly_array_Z2flag_depthint_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), ['GRMESZPHY']
# Temporary list dict
data = {}
# Permanent aggregate dict
#aggregates = {var: {} for var in variables}
#monthlydat = {var: {} for var in variables}

# Loop through years
for year in [2007,2008,2009,2010,2011,2012,2015,2016,2017,2018,2019,2020]:
    # 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 + '_dia2_T.nc') as ds:
            q = np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data
            q2 = q[0,:,:]
            monthly_array_Z2flag_depthint_slice[year-2007,month-1,:,:] = q2 #year2015 is index 0 along 1st dimension
            for var in ['GRMESZPHY']:
                data[var].append(np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data)
    
# Loop through years for wrap files
for year in [2013,2014]:
    # 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_wrap/SalishSea_1m_{datestr}_{datestr}'
        
        # Load grazing variables
        with xr.open_dataset(prefix + '_dia2_T.nc') as ds:
            q = np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data
            q2 = q[0,:,:]
            monthly_array_Z2flag_depthint_slice[year-2007,month-1,:,:] = q2 #year2015 is index 0 along 1st dimension
            for var in ['GRMESZPHY']:
                data[var].append(np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data)
    
    
    # Concatenate months
    #for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)
# # Calculate 5 year mean and anomalies
# for var in variables:
#     aggregates[var][‘mean’] = np.concatenate([aggregates[var][year][None, ...] for year in years]).mean(axis=0)
#     for year in years: aggregates[var][year] = aggregates[var][year] - aggregates[var][‘mean’]





# In[51]:


monthly_array_Z2flag_depthint_slice[monthly_array_Z2flag_depthint_slice == 0 ] = np.nan
monthly_array_Z2flag_depthint_slicemean = \
np.nanmean(np.nanmean(monthly_array_Z2flag_depthint_slice, axis = 2),axis = 2)
print(np.shape(monthly_array_Z2flag_depthint_slicemean))


# ### Select 4 warmest and 4 coldest years; leave NPGO "neutral" years out

# In[52]:


#2008, 2010, 2011, 2012
NPGO_C_Z2flag=np.array([[monthly_array_Z2flag_depthint_slicemean[1,:]],[monthly_array_Z2flag_depthint_slicemean[3,:]],\
        [monthly_array_Z2flag_depthint_slicemean[4,:]],[monthly_array_Z2flag_depthint_slicemean[5,:]]])


# In[53]:


NPGO_C_Z2flag_mean=NPGO_C_Z2flag.mean(axis=0).flatten()*86400*1.111


# In[54]:


NPGO_C_Z2flag_std=NPGO_C_Z2flag.std(axis=0).flatten()*86400*1.111


# In[55]:


#2015, 2018, 2019, 2020
NPGO_W_Z2flag=np.array([[monthly_array_Z2flag_depthint_slicemean[8,:]],[monthly_array_Z2flag_depthint_slicemean[11,:]],\
        [monthly_array_Z2flag_depthint_slicemean[12,:]],[monthly_array_Z2flag_depthint_slicemean[13,:]]])


# In[56]:


NPGO_W_Z2flag_mean=NPGO_W_Z2flag.mean(axis=0).flatten()*86400*1.111


# In[57]:


NPGO_W_Z2flag_std=NPGO_W_Z2flag.std(axis=0).flatten()*86400*1.111


# In[58]:


NPGO_C_Z2flag_mean.max()


# In[59]:


NPGO_W_Z2flag_mean.max()


# #### Z2 Grazing on both diatoms and nanoflagelles for Cold Years

# In[60]:


NPGO_C_Z2Both=NPGO_C_Z2diat_mean+NPGO_C_Z2flag_mean


# In[61]:


NPGO_C_Z2Both


# #### Z2 Grazing on both diatoms and nanoflagelles for Warm Years

# In[62]:


NPGO_W_Z2Both=NPGO_W_Z2diat_mean+NPGO_W_Z2flag_mean


# In[63]:


NPGO_W_Z2Both


# In[64]:


NPGO_C_Z2Both.mean()


# In[65]:


NPGO_W_Z2Both.mean()


# ### Load monthly average files for Z2 biomass

# In[66]:


#years, months, data
monthly_array_mesozooplankton_depthint_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), ['mesozooplankton']
# Temporary list dict
data = {}
# Permanent aggregate dict
aggregates = {var: {} for var in variables}
monthlydat = {var: {} for var in variables}

# Loop through years
for year in [2007,2008,2009,2010,2011,2012,2016,2017,2018,2019,2020]:
    # 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 = np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data
            q2 = q[0,:,:]
            monthly_array_mesozooplankton_depthint_slice[year-2007,month-1,:,:] = q2 #year2015 is index 0 along 1st dimension
            for var in ['mesozooplankton']:
                data[var].append(np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data)
    # Concatenate months
    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)

# Loop through years for wrap files
for year in [2013,2014,2015]:
    # 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_wrap/SalishSea_1m_{datestr}_{datestr}'
        # Load grazing variables
        with xr.open_dataset(prefix + '_ptrc_T.nc') as ds:
            q = np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data
            q2 = q[0,:,:]
            monthly_array_mesozooplankton_depthint_slice[year-2007,month-1,:,:] = q2 #year2015 is index 0 along 1st dimension
            for var in ['mesozooplankton']:
                data[var].append(np.ma.masked_where(tmask == 0, ds[var].isel(deptht=slice(None, 27), **slc).values * e3t).sum(axis=1).data)
    # Concatenate months
    for var in variables: aggregates[var][year] = np.concatenate(data[var]).mean(axis=0)      
        
        
# # Calculate 5 year mean and anomalies
# for var in variables:
#     aggregates[var][‘mean’] = np.concatenate([aggregates[var][year][None, ...] for year in years]).mean(axis=0)
#     for year in years: aggregates[var][year] = aggregates[var][year] - aggregates[var][‘mean’]




# In[67]:


monthly_array_mesozooplankton_depthint_slice[monthly_array_mesozooplankton_depthint_slice == 0 ] = np.nan
monthly_array_mesozooplankton_depthint_slicemean = \
np.nanmean(np.nanmean(monthly_array_mesozooplankton_depthint_slice, axis = 2),axis = 2)
print(np.shape(monthly_array_mesozooplankton_depthint_slicemean))


# ### Select 4 warmest and 4 coldest years; leave NPGO "neutral" years out

# In[68]:


#2008, 2010, 2011, 2012
NPGO_C_meso=np.array([[monthly_array_mesozooplankton_depthint_slicemean[1,:]],[monthly_array_mesozooplankton_depthint_slicemean[3,:]],\
        [monthly_array_mesozooplankton_depthint_slicemean[4,:]],[monthly_array_mesozooplankton_depthint_slicemean[5,:]]])


# In[69]:


NPGO_C_meso_mean=NPGO_C_meso.mean(axis=0).flatten()*5.7*12/1000


# In[70]:


NPGO_C_meso_std=NPGO_C_meso.std(axis=0).flatten()*5.7*12/1000


# In[71]:


#2015, 2018, 2019, 2020
NPGO_W_meso=np.array([[monthly_array_mesozooplankton_depthint_slicemean[8,:]],[monthly_array_mesozooplankton_depthint_slicemean[11,:]],\
        [monthly_array_mesozooplankton_depthint_slicemean[12,:]],[monthly_array_mesozooplankton_depthint_slicemean[13,:]]])


# In[72]:


NPGO_W_meso_mean=NPGO_W_meso.mean(axis=0).flatten()*5.7*12/1000


# In[73]:


NPGO_W_meso_std=NPGO_W_meso.std(axis=0).flatten()*5.7*12/1000


# In[ ]:





# In[74]:


### Plot Figure 6C

fig, ax = plt.subplots(figsize=(7,5))
bbox = {'boxstyle': 'round', 'facecolor': 'w', 'alpha': 0.9}
cmap = plt.get_cmap('tab10')
palette = [cmap(0), cmap(0.2), 'k', cmap(0.1), cmap(0.3)]
xticks=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov',"Dec"]

ax2=ax.twinx()

ax2.plot(xticks, NPGO_C_meso_mean,label='Z2',color='b',linewidth=3)

ax.bar(xticks,NPGO_C_Z2diat_mean,color='grey',edgecolor='k',hatch='//',alpha=0.7,label='Diatoms')
ax.bar(xticks,NPGO_C_Z2flag_mean,color='grey',edgecolor='k',label='Nanoflagellates',alpha=0.5,bottom=NPGO_C_Z2diat_mean)



ax2.set_ylim(0,5)
ax2.set_ylabel('Biomass (g C m$^{-2}$)',fontsize=14)

#ax.legend(frameon=False,loc=2)
ax2.legend(frameon=False,loc=1)
ax.set_ylim(0,12)
ax.set_ylabel('Grazing (mmol N m$^{-2}$ d$^{-1}$)',fontsize=14)
ax.set_xlabel('',fontsize=14)
#ax.tick_params(axis='x', labelrotation=0)

ax.text(-.5, 12.5, '(c)', fontsize=15,  color='k')

ax.set_title('',fontsize=18)

#plt.savefig('Fig6c.png', bbox_inches='tight',dpi=1000,transparent=False)


# In[75]:


### Plot Figure 6d


fig, ax = plt.subplots(figsize=(7,5))
bbox = {'boxstyle': 'round', 'facecolor': 'w', 'alpha': 0.9}
cmap = plt.get_cmap('tab10')
palette = [cmap(0), cmap(0.2), 'k', cmap(0.1), cmap(0.3)]
xticks=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov',"Dec"]

ax2=ax.twinx()

ax2.plot(xticks, NPGO_W_meso_mean,label='Z2',color='r',linestyle='-',linewidth=3)


ax.bar(xticks,NPGO_W_Z2diat_mean,color='grey',edgecolor='k',hatch='//',alpha=0.7,label='Diatoms')
ax.bar(xticks,NPGO_W_Z2flag_mean,color='grey',edgecolor='k',label='Nanoflagellates',alpha=0.5,bottom=NPGO_W_Z2diat_mean)



ax2.set_ylim(0,5)
ax2.set_ylabel('Biomass (g C m$^{-2}$)',fontsize=14)

ax2.legend(frameon=False,loc=1)
ax.set_ylim(0,12)
ax.set_ylabel('Grazing (mmol N m$^{-2}$ d$^{-1}$)',fontsize=14)
ax.set_xlabel('',fontsize=14)

#ax.tick_params(axis='x', labelrotation=0)
ax.text(-.5, 12.5, '(d)', fontsize=15,  color='k')
ax.set_title('',fontsize=18)

#plt.savefig('Fig6d.png', bbox_inches='tight',dpi=1000,transparent=False)


# ### Plot all on one figure

# In[116]:


fig, ax = plt.subplots(2,2,figsize=(13,10))
bbox = {'boxstyle': 'round', 'facecolor': 'w', 'alpha': 0.9}
cmap = plt.get_cmap('tab10')
palette = [cmap(0), cmap(0.2), 'k', cmap(0.1), cmap(0.3)]
xticks=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov',"Dec"]

fig.tight_layout(pad=4.5)


### PLOT FOR Z1 COLD YEARS

ax2 = ax[0,0].twinx()

ax2.errorbar(xticks, NPGO_C_micro_mean,label='Z1',color='b',linestyle='--',linewidth=3)
ax2.fill_between(xticks, NPGO_C_micro_mean-NPGO_C_micro_std, NPGO_C_micro_mean+NPGO_C_micro_std,
    alpha=0.5, edgecolor='royalblue', facecolor='royalblue',
    linewidth=0)

ax[0,0].bar(xticks,NPGO_C_Z1diat_mean,yerr=NPGO_C_Z1diat_std,color='grey',capsize=3,edgecolor='k',hatch='//',alpha=0.7,label='Diatoms')
ax[0,0].bar(xticks,NPGO_C_Z1flag_mean,yerr=NPGO_C_Z1flag_std,color='grey',ecolor='dimgrey',capsize=3,edgecolor='k',label='Nanoflagellates',alpha=0.5,bottom=NPGO_C_Z1diat_mean)



ax2.set_ylim(0,5)
ax2.set_ylabel('Biomass (g C m$^{-2}$)',fontsize=14)

ax[0,0].legend(frameon=False,loc=2)
ax2.legend(frameon=False,loc=1)
ax[0,0].set_ylim(0,12)
ax[0,0].set_ylabel('Grazing (mmol N m$^{-2}$ d$^{-1}$)',fontsize=14)
ax[0,0].set_xlabel('',fontsize=14)
#ax.tick_params(axis='x', labelrotation=0)

ax[0,0].text(-.5, 12.5, '(a)', fontsize=15,  color='k')

ax[0,0].set_title('Cold Years',fontsize=18)


### PLOT FOR Z1 WARM YEARS
ax3=ax[0,1].twinx()

ax3.errorbar(xticks, NPGO_W_micro_mean,label='Z1',color='r',linestyle='--',linewidth=3)
ax3.fill_between(xticks, NPGO_W_micro_mean-NPGO_W_micro_std, NPGO_W_micro_mean+NPGO_W_micro_std,
    alpha=0.5, edgecolor='tomato', facecolor='tomato',
    linewidth=0)

ax[0,1].bar(xticks,NPGO_W_Z1diat_mean,yerr=NPGO_W_Z1diat_std,color='grey',capsize=3,edgecolor='k',hatch='//',alpha=0.7,label='Diatoms')
ax[0,1].bar(xticks,NPGO_W_Z1flag_mean,yerr=NPGO_W_Z1flag_std,color='grey',ecolor='dimgrey',capsize=3,edgecolor='k',label='Nanoflagellates',alpha=0.5,bottom=NPGO_W_Z1diat_mean)



ax3.set_ylim(0,5)
ax3.set_ylabel('Biomass (g C m$^{-2}$)',fontsize=14)

ax3.legend(frameon=False,loc=1)
ax[0,1].set_ylim(0,12)
ax[0,1].set_ylabel('Grazing (mmol N m$^{-2}$ d$^{-1}$)',fontsize=14)
ax[0,1].set_xlabel('',fontsize=14)

#ax.tick_params(axis='x', labelrotation=0)
ax[0,1].text(-.5, 12.5, '(b)', fontsize=15,  color='k')

ax[0,1].set_title('Warm Years',fontsize=18)



### PLOT FOR Z2 COLD YEARS
ax4=ax[1,0].twinx()

ax4.errorbar(xticks, NPGO_C_meso_mean,label='Z2',color='b',linewidth=3)
ax4.fill_between(xticks, NPGO_C_meso_mean-NPGO_C_meso_std, NPGO_C_meso_mean+NPGO_C_meso_std,
    alpha=0.5, edgecolor='royalblue', facecolor='royalblue',
    linewidth=0)
ax[1,0].bar(xticks,NPGO_C_Z2diat_mean,yerr=NPGO_C_Z2diat_std,color='grey',edgecolor='k',capsize=3,hatch='//',alpha=0.7,label='Diatoms')
ax[1,0].bar(xticks,NPGO_C_Z2flag_mean,yerr=NPGO_C_Z2flag_std,color='grey',edgecolor='k',ecolor='dimgrey',capsize=3,label='Nanoflagellates',alpha=0.5,bottom=NPGO_C_Z2diat_mean)



ax4.set_ylim(0,5)
ax4.set_ylabel('Biomass (g C m$^{-2}$)',fontsize=14)

#ax.legend(frameon=False,loc=2)
ax4.legend(frameon=False,loc=1)
ax[1,0].set_ylim(0,12)
ax[1,0].set_ylabel('Grazing (mmol N m$^{-2}$ d$^{-1}$)',fontsize=14)
ax[1,0].set_xlabel('',fontsize=14)
#ax.tick_params(axis='x', labelrotation=0)

ax[1,0].text(-.5, 12.5, '(c)', fontsize=15,  color='k')

ax[1,0].set_title('',fontsize=18)


### PLOT FOR Z2 WARM YEARS
ax5=ax[1,1].twinx()

ax5.errorbar(xticks, NPGO_W_meso_mean,label='Z2',color='r',linestyle='-',linewidth=3)
ax5.fill_between(xticks, NPGO_W_meso_mean-NPGO_W_meso_std, NPGO_W_meso_mean+NPGO_W_meso_std,
    alpha=0.5, edgecolor='tomato', facecolor='tomato',
    linewidth=0)

ax[1,1].bar(xticks,NPGO_W_Z2diat_mean,yerr=NPGO_W_Z2diat_std,color='grey',edgecolor='k',capsize=3,hatch='//',alpha=0.7,label='Diatoms')
ax[1,1].bar(xticks,NPGO_W_Z2flag_mean,yerr=NPGO_W_Z2flag_std,color='grey',ecolor='dimgrey',capsize=3,edgecolor='k',label='Nanoflagellates',alpha=0.5,bottom=NPGO_W_Z2diat_mean)



ax5.set_ylim(0,5)
ax5.set_ylabel('Biomass (g C m$^{-2}$)',fontsize=14)

ax5.legend(frameon=False,loc=1)
ax[1,1].set_ylim(0,12)
ax[1,1].set_ylabel('Grazing (mmol N m$^{-2}$ d$^{-1}$)',fontsize=14)
ax[1,1].set_xlabel('',fontsize=14)

#ax.tick_params(axis='x', labelrotation=0)
ax[1,1].text(-.5, 12.5, '(d)', fontsize=15,  color='k')
ax[1,1].set_title('',fontsize=18)

#plt.savefig('Figure10_ZooplanktonGrazing_revised.png', bbox_inches='tight',dpi=1000,transparent=False)


# In[ ]: