Notebook

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
%matplotlib inline

In [31]:

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'][:,:,:,:]
SOGtmaskPath='/ocean/eolson/MEOPAR/northernNO3PaperCalcs/save/SOGtmask.pkl'
(tmaskSOG,_,_,_,_)=pickle.load(open(SOGtmaskPath,'rb'))

In [3]:

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

In [4]:

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

In [5]:

fP2=nc.Dataset(glob.glob('/data/eolson/results/MEOPAR/SS36runs/CedarRuns/shortTestSNAPe3tn/SalishSea_1d_*_dian_T*.nc')[0])

In [6]:

fP.variables.keys()

Out[6]:

dict_keys(['nav_lat', 'nav_lon', 'bounds_nav_lon', 'bounds_nav_lat', 'area', 'deptht', 'deptht_bounds', 'time_centered', 'time_centered_bounds', 'time_counter', 'time_counter_bounds', 'time_instant', 'time_instant_bounds', 'RIVNO3', 'PHYSTRNO3', 'ALLTRNO3', 'NO3_E3TSNAP', 'PPDIATNO3V', 'PPPHYNO3V', 'PPMRUBNO3V', 'NITR', 'NO3RAD', 'ATF_NO3', 'SMS_NO3', 'RDB_NO3', 'RDN_NO3', 'REM_NO3'])

In [7]:

fP2.variables.keys()

Out[7]:

dict_keys(['nav_lat', 'nav_lon', 'bounds_nav_lon', 'bounds_nav_lat', 'area', 'deptht', 'deptht_bounds', 'time_centered', 'time_centered_bounds', 'time_counter', 'time_counter_bounds', 'time_instant', 'time_instant_bounds', 'RIVNO3', 'BIOTRNO3', 'PHYSTRNO3', 'ALLTRNO3', 'AFILTNO3', 'NO3_E3T', 'NO3SNAP', 'NO3_E3TSNAP', 'PPDIATNO3V', 'PPPHYNO3V', 'PPMRUBNO3V', 'NITR', 'NO3RADB', 'NO3RADN'])

In [8]:

allSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['ALLTRNO3'][:,:,:,:],3),2),1)
afiltSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['ATF_NO3'][:,:,:,:],3),2),1)
radSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['RDB_NO3'][:,:,:,:],3),2),1)+\
       np.sum(np.sum(np.sum(tmaskSOG*fP.variables['RDN_NO3'][:,:,:,:],3),2),1)
no3Sum=np.sum(np.sum(np.sum(tmaskSOG*e1t*e2t*fP.variables['NO3_E3TSNAP'][:,:,:,:],3),2),1)
nitrSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['NITR'][:,:,:,:],3),2),1)
PPSum=np.sum(np.sum(np.sum(tmaskSOG*(fP.variables['PPDIATNO3V'][:,:,:,:]+\
                                             fP.variables['PPPHYNO3V'][:,:,:,:]+\
                                             fP.variables['PPMRUBNO3V'][:,:,:,:]),3),2),1)
PPDSum=np.sum(np.sum(np.sum(tmaskSOG*(fP.variables['PPDIATNO3V'][:,:,:,:]),3),2),1)
PPPSum=np.sum(np.sum(np.sum(tmaskSOG*(fP.variables['PPPHYNO3V'][:,:,:,:]),3),2),1)
rivSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['RIVNO3'][:,:,:,:],3),2),1)-nitrSum+PPSum
physSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['PHYSTRNO3'][:,:,:,:],3),2),1)-rivSum-nitrSum+PPSum

In [9]:

physSum2=np.sum(np.sum(np.sum(tmaskSOG*fP2.variables['PHYSTRNO3'][:,:,:,:],3),2),1)
allSum2=np.sum(np.sum(np.sum(tmaskSOG*fP2.variables['ALLTRNO3'][:,:,:,:],3),2),1)
afiltSum2=np.sum(np.sum(np.sum(tmaskSOG*fP2.variables['AFILTNO3'][:,:,:,:],3),2),1)
radSum2=np.sum(np.sum(np.sum(tmaskSOG*fP2.variables['NO3RADB'][:,:,:,:],3),2),1)+\
       np.sum(np.sum(np.sum(tmaskSOG*fP2.variables['NO3RADN'][:,:,:,:],3),2),1)
no3Sum2=np.sum(np.sum(np.sum(tmaskSOG*e1t*e2t*fP2.variables['NO3_E3TSNAP'][:,:,:,:],3),2),1)
nitrSum2=np.sum(np.sum(np.sum(tmaskSOG*fP2.variables['NITR'][:,:,:,:],3),2),1)
PPSum2=np.sum(np.sum(np.sum(tmaskSOG*(fP2.variables['PPDIATNO3V'][:,:,:,:]+\
                                             fP2.variables['PPPHYNO3V'][:,:,:,:]+\
                                             fP2.variables['PPMRUBNO3V'][:,:,:,:]),3),2),1)
PPDSum2=np.sum(np.sum(np.sum(tmaskSOG*(fP2.variables['PPDIATNO3V'][:,:,:,:]),3),2),1)
PPPSum2=np.sum(np.sum(np.sum(tmaskSOG*(fP2.variables['PPPHYNO3V'][:,:,:,:]),3),2),1)
rivSum2=np.sum(np.sum(np.sum(tmaskSOG*fP2.variables['RIVNO3'][:,:,:,:],3),2),1)

In [10]:

no3diff=(no3Sum[1:]-no3Sum[:-1])/(24*3600)
no3diff2=(no3Sum2[1:]-no3Sum2[:-1])/(24*3600)
#no3diff2=np.concatenate((no3diff,[(np.sum(np.sum(np.sum(tmaskSOG*e1t*e2t*fT.variables['NO3_E3T'][-1,:,:,:])))-no3sum[-1])/(23*3600)]))

In [11]:

afiltSum

Out[11]:

masked_array(data=[-67.60630798339844, -22.575477600097656,
                   9.601781845092773, 56.61271667480469,
                   -1.2884745597839355],
             mask=[False, False, False, False, False],
       fill_value=1e+20,
            dtype=float32)

In [12]:

np.max(radSum),np.min(radSum),np.max(radSum2),np.min(radSum2)

Out[12]:

(0.0, 0.0, 0.0030240756, 0.0)

In [ ]:

In [13]:

fig,ax=plt.subplots(1,3,figsize=(16,5))
ax[0].plot(range(1,5),no3diff,'k-',label='no3diff')
ax[0].plot(range(0,5),allSum,'g--',label='all')
ax[0].plot(range(0,5),radSum,'r--',label='all')
#ax[0].plot(range(0,5),afiltSum,'b--',label='afilt')
ax[0].legend()
ax[1].plot(range(0,5),allSum,'g-',label='all')
ax[1].plot(range(0,5),physSum,'k--',label='phys+bio+afilt')
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(0,5),nitrSum-PPSum,'r--',label='NITR-PP')
ax[2].legend()

Out[13]:

<matplotlib.legend.Legend at 0x7f7115875f10>

In [14]:

fig,ax=plt.subplots(1,3,figsize=(16,5))
ax[0].plot(range(1,5),no3diff,'k-',label='no3diff')
ax[0].plot(range(1,5),no3diff2,'r--',label='no3diff2')
ax[0].plot(range(0,5),allSum,'g-',label='all')
ax[0].plot(range(0,5),allSum2,'c--',label='all2')
ax[0].plot(0,allSum[0]*1.1,'k*')
ax[0].legend()
ax[0].set_ylim(-1200000,1200000)
ax[1].plot(range(0,5),physSum,'k-',label='phys')
ax[1].plot(range(0,5),physSum2,'r-',label='phys2')
ax[1].plot(range(0,5),nitrSum-PPSum,'m-',label='bio')
ax[1].plot(range(0,5),nitrSum2-PPSum2,'c--',label='bio2')
#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[1].set_ylim(-1200000,1600000)
ax[2].plot(range(0,5),nitrSum-PPSum,'r--',label='NITR-PP')
ax[2].plot(range(0,5),rivSum,'y-',label='riv')
ax[2].plot(range(0,5),rivSum2,'g--',label='riv2')
ax[2].legend()
ax[2].set_ylim(-800000,100000)

Out[14]:

(-800000, 100000)

In [15]:

fig,ax=plt.subplots(1,1,figsize=(5,5))
ax.plot(range(0,5),PPSum,'k-',label='PPTot')
ax.plot(range(0,5),PPSum2,'r--',label='PPTot2')
ax.plot(range(0,5),PPDSum,'b-',label='PPDiat')
ax.plot(range(0,5),PPDSum2,'c--',label='PPDiat2')
ax.plot(range(0,5),PPPSum,'y-',label='PPPhy')
ax.plot(range(0,5),PPPSum2,'g--',label='PPPhy2')
ax.legend()

Out[15]:

<matplotlib.legend.Legend at 0x7f71156cab20>

In [16]:

(nitrSum2-PPSum2)/(nitrSum-PPSum)

Out[16]:

masked_array(data=[0.989433228969574, 0.9959542155265808,
                   1.0051114559173584, 1.0044801235198975,
                   0.9778334498405457],
             mask=[False, False, False, False, False],
       fill_value=1e+20,
            dtype=float32)

In [ ]:

In [17]:

with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_Malaspina_U*.nc')[0]) as f:
    malUA=np.sum(np.sum(f.variables['NO3TVDX'][:,:,:,0],2),1)
    malUD=np.sum(np.sum(f.variables['ULDFNO3'][:,:,:,0],2),1)

In [18]:

with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_Haro_V*.nc')[0]) as f:
    harVA=np.sum(np.sum(f.variables['NO3TVDY'][:,:,0,:],2),1)
    harVD=np.sum(np.sum(f.variables['VLDFNO3'][:,:,0,:],2),1)

In [19]:

with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_SJC_V*.nc')[0]) as f:
    sjcVA=np.sum(np.sum(f.variables['NO3TVDY'][:,:,0,:],2),1)
    sjcVD=np.sum(np.sum(f.variables['VLDFNO3'][:,:,0,:],2),1)

In [20]:

with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_Rosario_V*.nc')[0]) as f:
    rosVA=np.sum(np.sum(f.variables['NO3TVDY'][:,:,0,:],2),1)
    rosVD=np.sum(np.sum(f.variables['VLDFNO3'][:,:,0,:],2),1)

In [21]:

with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_Sutil_V*.nc')[0]) as f:
    sutVA=np.sum(np.sum(f.variables['NO3TVDY'][:,:,0,:],2),1)
    sutVD=np.sum(np.sum(f.variables['VLDFNO3'][:,:,0,:],2),1)

In [22]:

with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_Discovery_V*.nc')[0]) as f:
    disVA=np.sum(np.sum(f.variables['NO3TVDY'][:,:,0,:],2),1)
    disVD=np.sum(np.sum(f.variables['VLDFNO3'][:,:,0,:],2),1)

In [23]:

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),physSum,'k-',label='phys')
ax[2].plot(range(0,5),adv+dif,'c--',label='lateral')
ax[2].plot(range(0,5),adv+dif+rivSum,'r--',label='rivers+lateral')
ax[2].plot(0,adv[0]+dif[0]+rivSum[0],'k*')
#ax[2].plot(np.arange(0.5,4,1),no3diff,'k--',label='no3diff')
ax[2].legend()

Out[23]:

<matplotlib.legend.Legend at 0x7f711552bf40>

In [24]:

adv[0],dif[0],rivSum[0],adv[0]+dif[0],adv[0]+dif[0]+rivSum[0],physSum[0]

Out[24]:

(1445586.5, 472.48126, 48516.0, 1446059.0, 1494575.0, 1444547.5)

np.max(np.abs(harVA-harTVDY))/np.max(np.abs(harTVDY))*100np.max(np.abs(harVD-harYM))/np.max(np.abs(harYM))*100

In [25]:

## phys difference
adv+dif+rivSum-physSum

Out[25]:

masked_array(data=[50027.5, 49391.359375, 50610.1015625, 50935.828125,
                   50657.890625],
             mask=[False, False, False, False, False],
       fill_value=1e+20,
            dtype=float32)

In [28]:

#old file, day 1:
f0=nc.Dataset('/results/SalishSea/nowcast-green.201812/01may15/SalishSea_1d_20150501_20150501_ptrc_T.nc')
# old file, day 5:
f4=nc.Dataset('/results/SalishSea/nowcast-green.201812/05may15/SalishSea_1d_20150505_20150505_ptrc_T.nc')
# new file:
f=nc.Dataset(idir+'SalishSea_1d_20150501_20150505_ptrc_T_20150501-20150505.nc')

In [27]:

fig,ax=plt.subplots(1,5,figsize=(16,5))
iax=ax[0]
fig.subplots_adjust(wspace=.4)
m0=iax.pcolormesh(np.ma.masked_where(tmask[0,:,:]==0,f.variables['NO3'][4,0,:,:]/f4.variables['nitrate'][0,0,:,:]),
                cmap=cmocean.cm.balance,vmin=0,vmax=2)
fig.colorbar(m0,ax=iax)
iax.set_title('leap-euler\n NO3')

iax=ax[1]
fig.subplots_adjust(wspace=.4)
m0=iax.pcolormesh(np.ma.masked_where(tmask[0,:,:]==0,f.variables['DIAT'][4,0,:,:]/f4.variables['diatoms'][0,0,:,:]),
                cmap=cmocean.cm.balance,vmin=0,vmax=2)
fig.colorbar(m0,ax=iax)
iax.set_title('leap-euler\n diatoms')

iax=ax[2]
fig.subplots_adjust(wspace=.4)
m0=iax.pcolormesh(np.ma.masked_where(tmask[0,:,:]==0,f.variables['PON'][4,0,:,:]/f4.variables['particulate_organic_nitrogen'][0,0,:,:]),
                cmap=cmocean.cm.balance,vmin=0,vmax=2)
fig.colorbar(m0,ax=iax)
iax.set_title('leap-euler\n PON')

iax=ax[3]
fig.subplots_adjust(wspace=.4)
m0=iax.pcolormesh(np.ma.masked_where(tmask[0,:,:]==0,f.variables['Si'][4,0,:,:]/f4.variables['silicon'][0,0,:,:]),
                cmap=cmocean.cm.balance,vmin=0,vmax=2)
fig.colorbar(m0,ax=iax)
iax.set_title('leap-euler\n Si')

iax=ax[4]
fig.subplots_adjust(wspace=.4)
m0=iax.pcolormesh(np.ma.masked_where(tmask[0,:,:]==0,f.variables['bSi'][4,0,:,:]/f4.variables['biogenic_silicon'][0,0,:,:]),
                cmap=cmocean.cm.balance,vmin=0,vmax=2)
fig.colorbar(m0,ax=iax)
iax.set_title('leap-euler\n bSi')

Out[27]:

Text(0.5, 1.0, 'leap-euler\n bSi')

In [37]:

fig,ax=plt.subplots(1,5,figsize=(16,5))
for jj in range(455,460): #range(450,475):
    for ii in range(265,270): #range(260,275):
        ax[0].plot(f0.variables['diatoms'][0,:,jj,ii],-1*f0.variables['deptht'][:],'k-',alpha=.1)
        ax[4].plot(f4.variables['diatoms'][0,:,jj,ii],-1*f0.variables['deptht'][:],'k-',alpha=.1)
for jj in range(455,460): #range(450,475):
    for ii in range(265,270): #range(260,275):
        ax[0].plot(f.variables['DIAT'][0,:,jj,ii],-1*f.variables['deptht'][:])
        ax[1].plot(f.variables['DIAT'][1,:,jj,ii],-1*f.variables['deptht'][:])
        ax[2].plot(f.variables['DIAT'][2,:,jj,ii],-1*f.variables['deptht'][:])
        ax[3].plot(f.variables['DIAT'][3,:,jj,ii],-1*f.variables['deptht'][:])
        ax[4].plot(f.variables['DIAT'][4,:,jj,ii],-1*f.variables['deptht'][:])
for iax in ax:
    iax.set_ylim(-100,0)
print('Input sinking rate: 0.5 -- 1.2 m/d')
z0=np.mean(np.mean(np.nansum(f.variables['DIAT'][0,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(f.variables['deptht'][:],(40,1,1)),0)/np.sum(f.variables['DIAT'][0,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z1=np.mean(np.mean(np.nansum(f.variables['DIAT'][1,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(f.variables['deptht'][:],(40,1,1)),0)/np.sum(f.variables['DIAT'][1,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z2=np.mean(np.mean(np.nansum(f.variables['DIAT'][2,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(f.variables['deptht'][:],(40,1,1)),0)/np.sum(f.variables['DIAT'][2,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z3=np.mean(np.mean(np.nansum(f.variables['DIAT'][3,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(f.variables['deptht'][:],(40,1,1)),0)/np.sum(f.variables['DIAT'][3,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z4=np.mean(np.mean(np.nansum(f.variables['DIAT'][4,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(f.variables['deptht'][:],(40,1,1)),0)/np.sum(f.variables['DIAT'][4,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
print(z1,z2,z3,z4)
print('Diagnosed sinking rate: ',(z4-z0)/4,'m/d')

Input sinking rate: 0.5 -- 1.2 m/d
26.127551393759376 21.70290536010795 36.49435152888682 57.49751129395327
Diagnosed sinking rate:  5.808659174567762 m/d

In [ ]: