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
%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)
SOGtmaskPath='/ocean/eolson/MEOPAR/northernNO3PaperCalcs/save/SOGtmask.pkl'
(tmaskSOG,_,_,_,_)=pickle.load(open(SOGtmaskPath,'rb'))
In [3]:
idir='/data/eolson/results/MEOPAR/SS36runs/CedarRuns/shortTestATFRAD/'
In [4]:
fP=nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_dian_T*.nc')[0])
In [5]:
fP.variables.keys()
Out[5]:
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', 'PPDIATNO3V', 'PPPHYNO3V', 'PPMRUBNO3V', 'RIVNO3', 'BIOTRNO3', 'PHYSTRNO3', 'ALLTRNO3', 'AFILTNO3', 'NITR', 'NO3RADB', 'NO3RADN', 'NO3_E3T', 'e3t'])
In [6]:
np.min(fP.variables['PPDIATNO3V'][0,0,:,:])
Out[6]:
0.0
In [7]:
fig,ax=plt.subplots(1,5,figsize=(18,5))
fig.subplots_adjust(wspace=.5)
m0=ax[0].pcolormesh(fP.variables['PPDIATNO3V'][0,2,:,:])
plt.colorbar(m0,ax=ax[0])
m1=ax[1].pcolormesh(fP.variables['PPPHYNO3V'][0,2,:,:])
plt.colorbar(m1,ax=ax[1])
m2=ax[2].pcolormesh(fP.variables['PPMRUBNO3V'][0,2,:,:])
plt.colorbar(m2,ax=ax[2])
riv=np.ma.masked_where(tmask==0,fP.variables['RIVNO3'][0,:,:,:])
m3=ax[3].pcolormesh(riv[0,:,:])
plt.colorbar(m3,ax=ax[3])
m4=ax[4].pcolormesh(fP.variables['NITR'][0,2,:,:])
plt.colorbar(m4,ax=ax[4])
test=np.argwhere(riv[0,:,:]>0)
print(np.min(riv),np.max(riv))
ax[3].plot(test[:,1],test[:,0],'r*')
0.0 9931.943
Out[7]:
[<matplotlib.lines.Line2D at 0x7f9b7e89f910>]
In [8]:
ii=(test[:,1]>300)&(test[:,0]>490)&(test[:,0]<510)
test[ii,:]
Out[8]:
array([[500, 394]])
In [9]:
fig,ax=plt.subplots(1,4,figsize=(12,5))
fig.subplots_adjust(wspace=.5)
m0=ax[0].pcolormesh(fP.variables['AFILTNO3'][0,2,:,:])
plt.colorbar(m0,ax=ax[0])
m1=ax[1].pcolormesh(fP.variables['NO3_E3T'][0,2,:,:])
plt.colorbar(m1,ax=ax[1])
m2=ax[2].pcolormesh(fP.variables['PHYSTRNO3'][0,2,:,:])
plt.colorbar(m2,ax=ax[2])
bio=np.ma.masked_where(tmask==0,fP.variables['BIOTRNO3'][0,:,:,:])
m3=ax[3].pcolormesh(bio[2,:,:])
plt.colorbar(m3,ax=ax[3])
Out[9]:
<matplotlib.colorbar.Colorbar at 0x7f9b7dddf4f0>
In [10]:
np.shape(tmaskSOG)
Out[10]:
(1, 40, 898, 398)
In [11]:
no3sum=np.sum(np.sum(np.sum(tmaskSOG*e1t*e2t*fP.variables['NO3_E3T'][:,:,:,:],3),2),1)
physSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['PHYSTRNO3'][:,:,:,:],3),2),1)
bioSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['BIOTRNO3'][:,:,:,:],3),2),1)
allSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['ALLTRNO3'][:,:,:,:],3),2),1)
rBSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['NO3RADB'][:,:,:,:],3),2),1)
rNSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['NO3RADN'][:,:,:,:],3),2),1)
rivSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['RIVNO3'][:,:,:,:],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)
afiltSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['AFILTNO3'][:,:,:,:],3),2),1)
In [12]:
no3diff=(no3sum[1:]-no3sum[:-1])/(24*3600)
In [13]:
24*3600/2160
Out[13]:
40.0
In [14]:
rivSum/1000 #mmol/2
Out[14]:
masked_array(data=[48.51604080200195, 49.706382751464844,
50.44744873046875, 50.93800735473633,
50.799522399902344],
mask=[False, False, False, False, False],
fill_value=1e+20,
dtype=float32)
In [15]:
afiltSum
Out[15]:
masked_array(data=[-72447.953125, -110218.71875, -121727.296875,
-108173.546875, -77526.484375],
mask=[False, False, False, False, False],
fill_value=1e+20,
dtype=float32)
In [16]:
np.max(afiltSum-physSum)
Out[16]:
-407797.44
In [17]:
fig,ax=plt.subplots(1,3,figsize=(16,5))
ax[0].plot(np.arange(0.5,4,1),no3diff,'k-',label='no3diff')
ax[0].plot(range(0,5),physSum,'g--',label='phys')
ax[0].plot(range(0,5),bioSum,'c--',label='bio')
ax[0].plot(range(0,5),physSum+bioSum,'r--',label='phys+bio')
#ax[0].plot(range(0,5),afiltSum,'b--',label='afilt')
ax[0].legend()
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),physSum,'k-',label='phys')
ax[2].plot(range(0,5),rivSum,'b--',label='rivers')
ax[2].legend()
Out[17]:
<matplotlib.legend.Legend at 0x7f9b657c3850>
In [18]:
fig,ax=plt.subplots(1,3,figsize=(16,5))
ax[0].plot(np.arange(0.5,4,1),no3diff,'k-',label='no3diff')
ax[0].plot(range(0,5),allSum,'g--',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),rNSum+rBSum,'c--',label='rad')
ax[1].plot(range(0,5),physSum+bioSum+afiltSum,'r--',label='phys+bio+afilt')
#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),physSum,'k-',label='phys')
#ax[2].plot(range(0,5),rivSum,'b--',label='rivers')
#ax[2].legend()
Out[18]:
<matplotlib.legend.Legend at 0x7f9b656cb4f0>
In [21]:
## max and min of RADB and RADN (negative correction)
print(np.min(rNSum),np.max(rNSum))
print(np.min(rBSum),np.max(rBSum))
0.0 0.0030240756 0.0 0.0
In [19]:
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 [20]:
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 [21]:
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 [22]:
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 [23]:
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 [24]:
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 [25]:
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,'b-',label='lateral')
ax[2].plot(range(0,5),adv+dif+rivSum,'r--',label='rivers+lateral')
#ax[2].plot(np.arange(0.5,4,1),no3diff,'k--',label='no3diff')
ax[2].legend()
Out[25]:
<matplotlib.legend.Legend at 0x7fd3ba390b50>
In [ ]: