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

# In[1]:


import netCDF4 as nc
from matplotlib import pyplot as plt
import numpy as np
import glob
import pickle
from salishsea_tools import evaltools as et
import datetime as dt
import os
import re
import cmocean
get_ipython().run_line_magic('matplotlib', 'inline')


# In[2]:


with nc.Dataset('/ocean/eolson/MEOPAR/NEMO-forcing/grid/mesh_mask201702_noLPE.nc') as mesh:
    tmask=mesh.variables['tmask'][0,:,:,:]
    e1t=np.expand_dims(mesh.variables['e1t'][:,:,:],1)
    e2t=np.expand_dims(mesh.variables['e2t'][:,:,:],1)
    e3t_0=mesh.variables['e3t_0'][:,:,:,:]
    vol_0=e1t*e2t*e3t_0*mesh.variables['tmask'][:,:,:,:]
cvolr=np.where(vol_0>0,1/vol_0,0)
SOGtmaskPath='/ocean/eolson/MEOPAR/northernNO3PaperCalcs/save/SOGtmask.pkl'
(tmaskSOG,_,_,_,_)=pickle.load(open(SOGtmaskPath,'rb'))


# In[3]:


np.shape(cvolr)


# In[4]:


idir='/data/eolson/results/MEOPAR/SS36runs/CedarRuns/shortTest1_11/'


# In[5]:


fP=nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_dian_T*.nc')[0])
fD=nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_dia1_T*.nc')[0])
fD2=nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_dia2_T*.nc')[0])


# In[6]:


len(np.shape(fP.variables['PRD_NH4']))


# In[7]:


len(np.shape(fD.variables['BFX_PON']))


# In[9]:


ignorelist=('nav_lat','nav_lon', 'bounds_nav_lon', 'bounds_nav_lat', 'area', 
            'deptht', 'deptht_bounds', 'time_centered', 'time_centered_bounds', 
            'time_counter', 'time_counter_bounds','nav_lat_grid_T', 'nav_lon_grid_T', 
            'bounds_nav_lon_grid_T', 'bounds_nav_lat_grid_T', 'area_grid_T', 
            'nav_lat_grid_W', 'nav_lon_grid_W', 'bounds_nav_lon_grid_W', 
            'bounds_nav_lat_grid_W', 'area_grid_W', 'time_centered', 'time_centered_bounds',
            'depthu', 'depthu_bounds','depthv', 'depthv_bounds','depthw', 'depthw_bounds',
            'layer6m_W', 'layer6m_W_bounds')
def checkall(ff):
    fkeys=ff.variables.keys()
    print('fP.variables.keys():',fkeys)
    print('Min/Max:')
    for var in fkeys:
        if var not in ignorelist:
            if len(np.shape(ff.variables[var]))==4:
                print(var,':',np.min(np.ma.masked_where(tmask[:,:,:]==0,ff.variables[var][-1,:,:,:])),
                      np.max(np.ma.masked_where(tmask[:,:,:]==0,ff.variables[var][-1,:,:,:])))
            elif len(np.shape(ff.variables[var]))==3:
                print(var,':',np.min(np.ma.masked_where(tmask[0,:,:]==0,ff.variables[var][-1,:,:])),
                      np.max(np.ma.masked_where(tmask[0,:,:]==0,fD.variables[var][-1,:,:])))
            else: 
                print('unknown shape: ',var,len(np.shape(ff.variables[var])))
def checkallSlice(ff):
    fkeys=ff.variables.keys()
    print('fP.variables.keys():',fkeys)
    print('Min/Max:')
    for var in fkeys:
        if var not in ignorelist:
            print(var,':',np.min(ff.variables[var][-1,:,:,:]),
                      np.max(ff.variables[var][-1,:,:,:]))


# In[10]:


checkall(fP)


# In[11]:


checkall(fD)


# In[12]:


checkall(fD2)


# In[13]:


plt.pcolormesh(fD.variables['ATF_NH4'][0,0,:,:])


# In[14]:


afiltSum=np.sum(np.sum(np.sum(tmaskSOG*(fD.variables['ATF_NO3'][:,:,:,:]+fD.variables['ATF_NH4']+\
            fD.variables['ATF_PON']+fD.variables['ATF_DON']+fD.variables['ATF_LIV']),3),2),1)


# In[15]:


afiltSumNO3=np.sum(np.sum(np.sum(tmaskSOG*fD.variables['ATF_NO3'][:,:,:,:],3),2),1)
afiltSumNH4=np.sum(np.sum(np.sum(tmaskSOG*fD.variables['ATF_NH4'],3),2),1)
afiltSumPON=np.sum(np.sum(np.sum(tmaskSOG*fD.variables['ATF_PON'],3),2),1)
afiltSumDON=np.sum(np.sum(np.sum(tmaskSOG*fD.variables['ATF_DON'],3),2),1)
afiltSumLiv=np.sum(np.sum(np.sum(tmaskSOG*fD.variables['ATF_LIV'],3),2),1)


# In[16]:


smsSumNO3=np.sum(np.sum(np.sum(tmaskSOG*fD2.variables['SMS_NO3'],3),2),1)
smsSumNH4=np.sum(np.sum(np.sum(tmaskSOG*fD2.variables['SMS_NH4'],3),2),1)
smsSumPON=np.sum(np.sum(np.sum(tmaskSOG*fD2.variables['SMS_PON'],3),2),1)
smsSumDON=np.sum(np.sum(np.sum(tmaskSOG*fD2.variables['SMS_DON'],3),2),1)
smsSumLiv=np.sum(np.sum(np.sum(tmaskSOG*fD2.variables['SMS_LIV'],3),2),1)


# In[17]:


phsSumNO3=np.sum(np.sum(np.sum(tmaskSOG*fD2.variables['PHS_NO3'],3),2),1)
phsSumNH4=np.sum(np.sum(np.sum(tmaskSOG*fD2.variables['PHS_NH4'],3),2),1)
phsSumPON=np.sum(np.sum(np.sum(tmaskSOG*fD2.variables['PHS_PON'],3),2),1)
phsSumDON=np.sum(np.sum(np.sum(tmaskSOG*fD2.variables['PHS_DON'],3),2),1)
phsSumLiv=np.sum(np.sum(np.sum(tmaskSOG*fD2.variables['PHS_LIV'],3),2),1)


# In[18]:


np.shape(e1t)


# In[19]:


bfxSumPON=np.sum(np.sum(tmaskSOG[:,0,:,:]*e1t[:,0,:,:]*e2t[:,0,:,:]*fD.variables['BFX_PON'],2),1)
bfxSumDIAT=np.sum(np.sum(tmaskSOG[:,0,:,:]*e1t[:,0,:,:]*e2t[:,0,:,:]*fD.variables['BFX_DIAT'],2),1)


# In[20]:


smsSumNO3


# In[21]:


smsSumNH4


# In[22]:


smsSumPON


# In[23]:


smsSumDON


# In[24]:


afiltSumNO3


# In[25]:


afiltSumNH4


# In[26]:


afiltSumPON


# In[27]:


afiltSumDON


# In[28]:


afiltSumLiv


# In[29]:


nitrSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['REM_NO3'][:,:,:,:],3),2),1)


# In[30]:


PPNO3Sum=-1*np.sum(np.sum(np.sum(tmaskSOG*(fP.variables['PRD_NO3'][:,:,:,:]),3),2),1)
PPNH4Sum=-1*np.sum(np.sum(np.sum(tmaskSOG*fP.variables['PRD_NH4'][:,:,:,:],3),2),1)


# In[31]:


REFRNSum=np.sum(np.sum(np.sum(tmaskSOG*(fP.variables['REFRN'][:,:,:,:]),3),2),1)


# In[32]:


remSum=-1*np.sum(np.sum(np.sum(tmaskSOG*(fP.variables['REM_DON'][:,:,:,:]+fP.variables['REM_PON'][:,:,:,:]),3),2),1)


# In[33]:


rivSum=np.sum(np.sum(np.sum(tmaskSOG*(fP.variables['RIV_NO3'][:,:,:,:]+fP.variables['RIV_NH4'][:,:,:,:]),3),2),1)


# In[34]:


rivSumNO3=np.sum(np.sum(np.sum(tmaskSOG*(fP.variables['RIV_NO3'][:,:,:,:]),3),2),1)


# In[35]:


rivSumNH4=np.sum(np.sum(np.sum(tmaskSOG*(fP.variables['RIV_NH4'][:,:,:,:]),3),2),1)


# In[36]:


fS=nc.Dataset(idir+'SalishSea_1d_20150501_20150505_snp_T_20150501-20150505.nc')


# In[37]:


checkall(fS)


# In[38]:


fS.variables.keys()


# In[39]:


no3Sum=np.sum(np.sum(np.sum(tmaskSOG*e1t*e2t*fS.variables['NO3_E3TSNAP'][:,:,:,:],3),2),1)
nh4Sum=np.sum(np.sum(np.sum(tmaskSOG*e1t*e2t*fS.variables['NH4_E3TSNAP'][:,:,:,:],3),2),1)
ponSum=np.sum(np.sum(np.sum(tmaskSOG*e1t*e2t*fS.variables['PON_E3TSNAP'][:,:,:,:],3),2),1)
donSum=np.sum(np.sum(np.sum(tmaskSOG*e1t*e2t*fS.variables['DON_E3TSNAP'][:,:,:,:],3),2),1)
livSum=np.sum(np.sum(np.sum(tmaskSOG*e1t*e2t*fS.variables['LIV_E3TSNAP'][:,:,:,:],3),2),1)


# In[40]:


no3Sum


# In[41]:


livSum


# In[42]:


no3diff=(no3Sum[1:]-no3Sum[:-1])/(24*3600)
nh4diff=(nh4Sum[1:]-nh4Sum[:-1])/(24*3600)
pondiff=(ponSum[1:]-ponSum[:-1])/(24*3600)
dondiff=(donSum[1:]-donSum[:-1])/(24*3600)
livdiff=(livSum[1:]-livSum[:-1])/(24*3600)


# In[43]:


plt.pcolormesh(fS.variables['LIV_E3TSNAP'][0,0,:,:])
plt.colorbar()


# In[44]:


phsSumNO3


# In[45]:


smsSumNO3


# In[46]:


fig,ax=plt.subplots(1,3,figsize=(16,5))
#ax[0].plot(range(0,5),afiltSum,'k-',label='afilt')
ax[0].plot(range(0,5),PPNO3Sum,'r--',label='PPNO3')
ax[0].plot(range(0,5),PPNH4Sum,'m--',label='PPNH4')
ax[0].plot(range(0,5),rivSum,'b--',label='riv')
ax[0].legend()
ax[1].plot(range(0,5),remSum,'g-',label='REM->NH4')
ax[1].plot(range(0,5),REFRNSum,'k--',label='REFR')
ax[0].plot(range(0,5),nitrSum,'g--',label='nitr')
#ax[1].plot(range(0,5),afiltSum,'y--',label='afilt')
#ax[1].plot(range(0,5),nitrSum-PPSum,'m--',label='bio')
#ax[1].plot(range(0,5),physSum+nitrSum-PPSum+rivSum,'c--',label='phys+bio+riv')
#ax[1].plot(range(0,5),bioSum,'k-',label='bio')
#ax[1].plot(range(0,5),nitrSum-PPSum,'r--',label='NITR-PP')
ax[1].legend()
ax[2].plot(range(1,5),no3diff,'k-',label='NO3 diff')
ax[2].plot(range(0,5),smsSumNO3,'m--',label='NO3 sms')
ax[2].plot(range(0,5),phsSumNO3+smsSumNO3,'r--',label='NO3 sms+phys')
#ax[2].plot(range(1,5),livdiff,'k-',label='liv diff')
#ax[2].plot(range(0,5),phsSumLiv+smsSumLiv,'r--',label='liv sms+phys')
#ax[2].legend()


# In[47]:


fig,ax=plt.subplots(1,3,figsize=(16,5))
ax[0].plot(range(0,5),smsSumNO3,'r-',label='sms')
ax[0].plot(range(0,5),nitrSum-PPNO3Sum,'b--',label='nitr-ppno3')
ax[0].plot(range(1,5),no3diff,'k-',label='diff')
ax[0].plot(range(0,5),phsSumNO3,'c-',label='phys')
ax[0].plot(range(0,5),phsSumNO3+smsSumNO3+afiltSumNO3+rivSumNO3,'m--',label='phs+sms+aft+riv')
ax[0].legend()
ax[0].set_title('NO3')
ax[1].plot(range(0,5),smsSumNH4,'r-',label='sms')
ax[1].plot(range(1,5),nh4diff,'k-',label='diff')
ax[1].plot(range(0,5),phsSumNH4,'c-',label='phys')
ax[1].plot(range(0,5),phsSumNH4+smsSumNH4+afiltSumNH4+rivSumNH4,'m--',label='phs+sms+aft+riv')
ax[1].legend()
ax[1].set_title('NH4')
ax[2].plot(range(0,5),smsSumDON,'r-',label='sms')
ax[2].plot(range(1,5),dondiff,'k-',label='diff')
ax[2].plot(range(0,5),phsSumDON,'c-',label='phys')
ax[2].plot(range(0,5),phsSumDON+smsSumDON+afiltSumDON,'m--',label='phs+sms+aft+riv')
ax[2].legend()
ax[2].set_title('DON')


# In[48]:


bfxSumPON


# In[49]:


afiltSumPON


# In[50]:


fig,ax=plt.subplots(1,3,figsize=(16,5))
ax[0].plot(range(1,5),pondiff,'k-',label='diff')
ax[0].plot(range(0,5),smsSumPON,'r-',label='sms')
ax[0].plot(range(0,5),bfxSumPON,'g-',label='sinking')
ax[0].plot(range(0,5),afiltSumPON,'y-',label='afilt')
ax[0].plot(range(0,5),phsSumPON,'c-',label='phys')
ax[0].plot(range(0,5),phsSumPON+smsSumPON+afiltSumPON+bfxSumPON,'m--',label='phs+sms+aft+sink')
ax[0].legend()
ax[0].set_title('PON')
ax[1].plot(range(1,5),livdiff,'k-',label='diff')
ax[1].plot(range(0,5),smsSumLiv,'r-',label='sms')
ax[1].plot(range(0,5),bfxSumDIAT,'g-',label='sinking')
ax[1].plot(range(0,5),afiltSumLiv,'y-',label='afilt')
ax[1].plot(range(0,5),phsSumLiv,'c-',label='phys')
ax[1].plot(range(0,5),phsSumLiv+smsSumLiv+afiltSumLiv+bfxSumDIAT,'m--',label='phs+sms+aft+sink')
ax[1].legend()
ax[1].set_title('Liv')


# In[ ]:





# In[51]:


with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_Malaspina_U*.nc')[0]) as f:
    checkallSlice(f)
    malUA=np.sum(np.sum(f.variables['ATX_NO3'][:,:,:,0],2),1)
    malUD=np.sum(np.sum(f.variables['DTX_NO3'][:,:,:,0],2),1)
    malUANH4=np.sum(np.sum(f.variables['ATX_NH4'][:,:,:,0],2),1)
    malUDNH4=np.sum(np.sum(f.variables['DTX_NH4'][:,:,:,0],2),1)
    malUADON=np.sum(np.sum(f.variables['ATX_DON'][:,:,:,0],2),1)
    malUDDON=np.sum(np.sum(f.variables['DTX_DON'][:,:,:,0],2),1)
    malUAPON=np.sum(np.sum(f.variables['ATX_PON'][:,:,:,0],2),1)
    malUDPON=np.sum(np.sum(f.variables['DTX_PON'][:,:,:,0],2),1)
    malUALIV=np.sum(np.sum(f.variables['ATX_LIV'][:,:,:,0],2),1)
    malUDLIV=np.sum(np.sum(f.variables['DTX_LIV'][:,:,:,0],2),1)


# In[52]:


with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_Haro_V*.nc')[0]) as f:
    checkallSlice(f)
    print(f.variables.keys())
    harVA=np.sum(np.sum(f.variables['ATY_NO3'][:,:,0,:],2),1)
    harVD=np.sum(np.sum(f.variables['DTY_NO3'][:,:,0,:],2),1)
    harVANH4=np.sum(np.sum(f.variables['ATY_NH4'][:,:,0,:],2),1)
    harVDNH4=np.sum(np.sum(f.variables['DTY_NH4'][:,:,0,:],2),1)
    harVADON=np.sum(np.sum(f.variables['ATY_DON'][:,:,0,:],2),1)
    harVDDON=np.sum(np.sum(f.variables['DTY_DON'][:,:,0,:],2),1)
    harVAPON=np.sum(np.sum(f.variables['ATY_PON'][:,:,0,:],2),1)
    harVDPON=np.sum(np.sum(f.variables['DTY_PON'][:,:,0,:],2),1)
    harVALIV=np.sum(np.sum(f.variables['ATY_LIV'][:,:,0,:],2),1)
    harVDLIV=np.sum(np.sum(f.variables['DTY_LIV'][:,:,0,:],2),1)


# In[53]:


harVANH4


# In[54]:


with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_SJC_V*.nc')[0]) as f:
    checkallSlice(f)
    sjcVA=np.sum(np.sum(f.variables['ATY_NO3'][:,:,0,:],2),1)
    sjcVD=np.sum(np.sum(f.variables['DTY_NO3'][:,:,0,:],2),1)
    sjcVANH4=np.sum(np.sum(f.variables['ATY_NH4'][:,:,0,:],2),1)
    sjcVDNH4=np.sum(np.sum(f.variables['DTY_NH4'][:,:,0,:],2),1)
    sjcVADON=np.sum(np.sum(f.variables['ATY_DON'][:,:,0,:],2),1)
    sjcVDDON=np.sum(np.sum(f.variables['DTY_DON'][:,:,0,:],2),1)
    sjcVAPON=np.sum(np.sum(f.variables['ATY_PON'][:,:,0,:],2),1)
    sjcVDPON=np.sum(np.sum(f.variables['DTY_PON'][:,:,0,:],2),1)
    sjcVALIV=np.sum(np.sum(f.variables['ATY_LIV'][:,:,0,:],2),1)
    sjcVDLIV=np.sum(np.sum(f.variables['DTY_LIV'][:,:,0,:],2),1)


# In[55]:


with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_Rosario_V*.nc')[0]) as f:
    checkallSlice(f)
    rosVA=np.sum(np.sum(f.variables['ATY_NO3'][:,:,0,:],2),1)
    rosVD=np.sum(np.sum(f.variables['DTY_NO3'][:,:,0,:],2),1)
    rosVANH4=np.sum(np.sum(f.variables['ATY_NH4'][:,:,0,:],2),1)
    rosVDNH4=np.sum(np.sum(f.variables['DTY_NH4'][:,:,0,:],2),1)
    rosVADON=np.sum(np.sum(f.variables['ATY_DON'][:,:,0,:],2),1)
    rosVDDON=np.sum(np.sum(f.variables['DTY_DON'][:,:,0,:],2),1)
    rosVAPON=np.sum(np.sum(f.variables['ATY_PON'][:,:,0,:],2),1)
    rosVDPON=np.sum(np.sum(f.variables['DTY_PON'][:,:,0,:],2),1)
    rosVALIV=np.sum(np.sum(f.variables['ATY_LIV'][:,:,0,:],2),1)
    rosVDLIV=np.sum(np.sum(f.variables['DTY_LIV'][:,:,0,:],2),1)


# In[56]:


with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_Sutil_V*.nc')[0]) as f:
    checkallSlice(f)
    sutVA=np.sum(np.sum(f.variables['ATY_NO3'][:,:,0,:],2),1)
    sutVD=np.sum(np.sum(f.variables['DTY_NO3'][:,:,0,:],2),1)
    sutVANH4=np.sum(np.sum(f.variables['ATY_NH4'][:,:,0,:],2),1)
    sutVDNH4=np.sum(np.sum(f.variables['DTY_NH4'][:,:,0,:],2),1)
    sutVADON=np.sum(np.sum(f.variables['ATY_DON'][:,:,0,:],2),1)
    sutVDDON=np.sum(np.sum(f.variables['DTY_DON'][:,:,0,:],2),1)
    sutVAPON=np.sum(np.sum(f.variables['ATY_PON'][:,:,0,:],2),1)
    sutVDPON=np.sum(np.sum(f.variables['DTY_PON'][:,:,0,:],2),1)
    sutVALIV=np.sum(np.sum(f.variables['ATY_LIV'][:,:,0,:],2),1)
    sutVDLIV=np.sum(np.sum(f.variables['DTY_LIV'][:,:,0,:],2),1)


# In[57]:


with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_Discovery_V*.nc')[0]) as f:
    checkallSlice(f)
    disVA=np.sum(np.sum(f.variables['ATY_NO3'][:,:,0,:],2),1)
    disVD=np.sum(np.sum(f.variables['DTY_NO3'][:,:,0,:],2),1)
    disVANH4=np.sum(np.sum(f.variables['ATY_NH4'][:,:,0,:],2),1)
    disVDNH4=np.sum(np.sum(f.variables['DTY_NH4'][:,:,0,:],2),1)
    disVADON=np.sum(np.sum(f.variables['ATY_DON'][:,:,0,:],2),1)
    disVDDON=np.sum(np.sum(f.variables['DTY_DON'][:,:,0,:],2),1)
    disVAPON=np.sum(np.sum(f.variables['ATY_PON'][:,:,0,:],2),1)
    disVDPON=np.sum(np.sum(f.variables['DTY_PON'][:,:,0,:],2),1)
    disVALIV=np.sum(np.sum(f.variables['ATY_LIV'][:,:,0,:],2),1)
    disVDLIV=np.sum(np.sum(f.variables['DTY_LIV'][:,:,0,:],2),1)


# In[58]:


fig,ax=plt.subplots(1,3,figsize=(16,6))
ax[0].plot(range(0,5),-1*malUA,'r-',label='Malaspina')
ax[0].plot(range(0,5),-1*sutVA,'m-',label='Sutil')
ax[0].plot(range(0,5),-1*disVA,'y-',label='Discovery')
ax[0].plot(range(0,5),harVA,'b-',label='Haro')
ax[0].plot(range(0,5),sjcVA,'c-',label='SJC')
ax[0].plot(range(0,5),rosVA,'g-',label='Rosario')
adv=harVA+sjcVA+rosVA-disVA-sutVA-malUA
ax[0].plot(range(0,5),adv,'k-',label='sum')
ax[0].legend()
ax[0].set_title('Advection')
ax[1].plot(range(0,5),-1*malUD,'r-',label='Malaspina')
ax[1].plot(range(0,5),-1*sutVD,'m-',label='Sutil')
ax[1].plot(range(0,5),-1*disVD,'y-',label='Discovery')
ax[1].plot(range(0,5),harVD,'b-',label='Haro iso')
ax[1].plot(range(0,5),sjcVD,'c-',label='SJC')
ax[1].plot(range(0,5),rosVD,'g-',label='Rosario')
dif=harVD+sjcVD+rosVD-disVD-sutVD-malUD
ax[1].plot(range(0,5),dif,'k-',label='sum')
ax[1].legend()
ax[1].set_title('Diffusion')
ax[2].plot(range(0,5),phsSumNO3,'k-',label='phys')
ax[2].plot(range(0,5),adv+dif,'c--',label='lateral')
ax[2].plot(range(0,5),adv+dif+rivSumNO3,'r--',label='rivers+lateral')
ax[2].plot(0,adv[0]+dif[0]+rivSumNO3[0],'k*')
ax[2].legend()


# In[59]:


fig,ax=plt.subplots(1,3,figsize=(16,6))
ax[0].plot(range(0,5),-1*malUANH4,'r-',label='Malaspina')
ax[0].plot(range(0,5),-1*sutVANH4,'m-',label='Sutil')
ax[0].plot(range(0,5),-1*disVANH4,'y-',label='Discovery')
ax[0].plot(range(0,5),harVANH4,'b-',label='Haro')
ax[0].plot(range(0,5),sjcVANH4,'c-',label='SJC')
ax[0].plot(range(0,5),rosVANH4,'g-',label='Rosario')
adv=harVANH4+sjcVANH4+rosVANH4-disVANH4-sutVANH4-malUANH4
#ax[0].plot(range(0,5),adv,'k-',label='sum')
ax[0].legend()
ax[0].set_title('Advection')
ax[1].plot(range(0,5),-1*malUDNH4,'r-',label='Malaspina')
ax[1].plot(range(0,5),-1*sutVDNH4,'m-',label='Sutil')
ax[1].plot(range(0,5),-1*disVDNH4,'y-',label='Discovery')
ax[1].plot(range(0,5),harVDNH4,'b-',label='Haro iso')
ax[1].plot(range(0,5),sjcVDNH4,'c-',label='SJC')
ax[1].plot(range(0,5),rosVDNH4,'g-',label='Rosario')
dif=harVDNH4+sjcVDNH4+rosVDNH4-disVDNH4-sutVDNH4-malUDNH4
#ax[1].plot(range(0,5),dif,'k-',label='sum')
ax[1].legend()
ax[1].set_title('Diffusion')
ax[2].plot(range(0,5),phsSumNH4,'k-',label='phys')
ax[2].plot(range(0,5),adv+dif,'c--',label='lateral')
ax[2].plot(range(0,5),adv+dif+rivSumNH4,'r--',label='rivers+lateral')
ax[2].plot(0,adv[0]+dif[0]+rivSumNH4[0],'k*')
ax[2].legend()


# In[60]:


fig,ax=plt.subplots(1,3,figsize=(16,6))
ax[0].plot(range(0,5),-1*malUAPON,'r-',label='Malaspina')
ax[0].plot(range(0,5),-1*sutVAPON,'m-',label='Sutil')
ax[0].plot(range(0,5),-1*disVAPON,'y-',label='Discovery')
ax[0].plot(range(0,5),harVAPON,'b-',label='Haro')
ax[0].plot(range(0,5),sjcVAPON,'c-',label='SJC')
ax[0].plot(range(0,5),rosVAPON,'g-',label='Rosario')
adv=harVAPON+sjcVAPON+rosVAPON-disVAPON-sutVAPON-malUAPON
#ax[0].plot(range(0,5),adv,'k-',label='sum')
ax[0].legend()
ax[0].set_title('Advection')
ax[1].plot(range(0,5),-1*malUDPON,'r-',label='Malaspina')
ax[1].plot(range(0,5),-1*sutVDPON,'m-',label='Sutil')
ax[1].plot(range(0,5),-1*disVDPON,'y-',label='Discovery')
ax[1].plot(range(0,5),harVDPON,'b-',label='Haro iso')
ax[1].plot(range(0,5),sjcVDPON,'c-',label='SJC')
ax[1].plot(range(0,5),rosVDPON,'g-',label='Rosario')
dif=harVDPON+sjcVDPON+rosVDPON-disVDPON-sutVDPON-malUDPON
#ax[1].plot(range(0,5),dif,'k-',label='sum')
ax[1].legend()
ax[1].set_title('Diffusion')
ax[2].plot(range(0,5),phsSumPON,'k-',label='phys')
ax[2].plot(range(0,5),adv+dif,'c--',label='lateral')
ax[2].plot(range(0,5),adv+dif,'r--',label='rivers+lateral')
ax[2].plot(0,adv[0]+dif[0],'k*')
ax[2].legend()


# In[61]:


fig,ax=plt.subplots(1,3,figsize=(16,6))
ax[0].plot(range(0,5),-1*malUADON,'r-',label='Malaspina')
ax[0].plot(range(0,5),-1*sutVADON,'m-',label='Sutil')
ax[0].plot(range(0,5),-1*disVADON,'y-',label='Discovery')
ax[0].plot(range(0,5),harVADON,'b-',label='Haro')
ax[0].plot(range(0,5),sjcVADON,'c-',label='SJC')
ax[0].plot(range(0,5),rosVADON,'g-',label='Rosario')
adv=harVADON+sjcVADON+rosVADON-disVADON-sutVADON-malUADON
#ax[0].plot(range(0,5),adv,'k-',label='sum')
ax[0].legend()
ax[0].set_title('Advection')
ax[1].plot(range(0,5),-1*malUDDON,'r-',label='Malaspina')
ax[1].plot(range(0,5),-1*sutVDDON,'m-',label='Sutil')
ax[1].plot(range(0,5),-1*disVDDON,'y-',label='Discovery')
ax[1].plot(range(0,5),harVDDON,'b-',label='Haro iso')
ax[1].plot(range(0,5),sjcVDDON,'c-',label='SJC')
ax[1].plot(range(0,5),rosVDDON,'g-',label='Rosario')
dif=harVDDON+sjcVDDON+rosVDDON-disVDDON-sutVDDON-malUDDON
#ax[1].plot(range(0,5),dif,'k-',label='sum')
ax[1].legend()
ax[1].set_title('Diffusion')
ax[2].plot(range(0,5),phsSumDON,'k-',label='phys')
ax[2].plot(range(0,5),adv+dif,'c--',label='lateral')
ax[2].plot(range(0,5),adv+dif,'r--',label='rivers+lateral')
ax[2].plot(0,adv[0]+dif[0],'k*')
ax[2].legend()


# In[62]:


fig,ax=plt.subplots(1,3,figsize=(16,6))
ax[0].plot(range(0,5),-1*malUALIV,'r-',label='Malaspina')
ax[0].plot(range(0,5),-1*sutVALIV,'m-',label='Sutil')
ax[0].plot(range(0,5),-1*disVALIV,'y-',label='Discovery')
ax[0].plot(range(0,5),harVALIV,'b-',label='Haro')
ax[0].plot(range(0,5),sjcVALIV,'c-',label='SJC')
ax[0].plot(range(0,5),rosVALIV,'g-',label='Rosario')
adv=harVALIV+sjcVALIV+rosVALIV-disVALIV-sutVALIV-malUALIV
#ax[0].plot(range(0,5),adv,'k-',label='sum')
ax[0].legend()
ax[0].set_title('Advection')
ax[1].plot(range(0,5),-1*malUDLIV,'r-',label='Malaspina')
ax[1].plot(range(0,5),-1*sutVDLIV,'m-',label='Sutil')
ax[1].plot(range(0,5),-1*disVDLIV,'y-',label='Discovery')
ax[1].plot(range(0,5),harVDLIV,'b-',label='Haro iso')
ax[1].plot(range(0,5),sjcVDLIV,'c-',label='SJC')
ax[1].plot(range(0,5),rosVDLIV,'g-',label='Rosario')
dif=harVDLIV+sjcVDLIV+rosVDLIV-disVDLIV-sutVDLIV-malUDLIV
#ax[1].plot(range(0,5),dif,'k-',label='sum')
ax[1].legend()
ax[1].set_title('Diffusion')
ax[2].plot(range(0,5),phsSumLiv,'k-',label='phys')
ax[2].plot(range(0,5),adv+dif,'c--',label='lateral')
ax[2].plot(range(0,5),adv+dif,'r--',label='rivers+lateral')
ax[2].plot(0,adv[0]+dif[0],'k*')
ax[2].legend()


# In[63]:


fW=nc.Dataset(idir+'SalishSea_1d_20150501_20150505_layer6m_W_20150501-20150505.nc')
fW2=nc.Dataset(idir+'SalishSea_1d_20150501_20150505_layer6m_W_2_20150501-20150505.nc')


# In[64]:


checkallSlice(fW)


# In[65]:


checkallSlice(fW2)


# In[66]:


fr=nc.Dataset(idir+'SalishSea_1d_20150501_20150505_rad_T_20150501-20150505.nc')


# In[67]:


checkallSlice(fr)


# In[ ]: