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 [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'][:,:,:,:]
SOGtmaskPath='/ocean/eolson/MEOPAR/northernNO3PaperCalcs/save/SOGtmask.pkl'
(tmaskSOG,_,_,_,_)=pickle.load(open(SOGtmaskPath,'rb'))
In [3]:
idir='/data/eolson/results/MEOPAR/SS36runs/CedarRuns/shortTest1812/'
In [4]:
fP=nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_dian_T*.nc')[0])
fD=nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_dia1_T*.nc')[0])
fM=nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_graz_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', 'PHYSTRNO3', 'NO3RAD', 'SMS_NO3', 'RIV_NO3', 'RIV_NH4', 'REM_NO3', 'REM_PON', 'REM_DON', 'PRD_NO3', 'PRD_NH4', 'REFRN'])
In [6]:
fD.variables.keys()
Out[6]:
dict_keys(['nav_lat_grid_T', 'nav_lon_grid_T', 'bounds_nav_lon_grid_T', 'bounds_nav_lat_grid_T', 'area_grid_T', 'deptht', 'deptht_bounds', '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', 'time_counter', 'time_counter_bounds', 'ATF_NO3', 'ATF_NH4', 'ATF_DON', 'ATF_PON', 'ATF_LIV', 'BFX_PON', 'BFX_DIAT'])
In [7]:
fD.variables['BFX_PON']
Out[7]:
<class 'netCDF4._netCDF4.Variable'>
float32 BFX_PON(time_counter, y_grid_W, x_grid_W)
long_name: Bottom Sinking flux
units: mmol s-1
online_operation: average
interval_operation: 40 s
interval_write: 1 d
cell_methods: time: mean (interval: 40 s)
cell_measures: area: area_grid_W
_FillValue: 1e+20
missing_value: 1e+20
coordinates: time_centered nav_lat_grid_W nav_lon_grid_W
unlimited dimensions: time_counter
current shape = (5, 898, 398)
filling on
In [8]:
np.min(fD.variables['BFX_PON']),np.max(fD.variables['BFX_PON'])
Out[8]:
(1e+20, 1e+20)
In [50]:
np.min(fD.variables['ATF_LIV']),np.max(fD.variables['ATF_LIV'])
Out[50]:
(1e+20, 1e+20)
In [53]:
np.min(fD.variables['ATF_PON']),np.max(fD.variables['ATF_NH4'])
Out[53]:
(1e+20, 1e+20)
In [52]:
np.min(fD.variables['ATF_NO3']),np.max(fD.variables['ATF_NO3'])
Out[52]:
(-3.0665686, 1e+20)
In [9]:
np.min(fD.variables['BFX_PON']),np.max(fD.variables['BFX_DIAT'])
Out[9]:
(1e+20, 1e+20)
In [10]:
plt.pcolormesh(fD.variables['BFX_PON'][4,:,:])
Out[10]:
<matplotlib.collections.QuadMesh at 0x7f9ed51d4370>
<__array_function__ internals>:5: UserWarning: Warning: converting a masked element to nan. /home/eolson/anaconda3/envs/py38/lib/python3.8/site-packages/numpy/ma/core.py:722: UserWarning: Warning: converting a masked element to nan. data = np.array(a, copy=False, subok=subok)
In [11]:
fM.variables.keys()
Out[11]:
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', 'MRU_MYRI', 'MSZ_DIAT', 'MSZ_MYRI', 'MSZ_PHY', 'MSZ_MICZ', 'MSZ_PON', 'MIZ_DIAT', 'MIZ_MYRI', 'MIZ_PHY', 'MIZ_MICZ', 'MIZ_PON'])
In [12]:
np.min(tmaskSOG*fM.variables['MIZ_DIAT']),np.max(tmaskSOG*fM.variables['MIZ_DIAT'])
Out[12]:
(-18.051245, 0.0)
In [13]:
plt.pcolormesh(fM.variables['MIZ_DIAT'][4,0,:,:])
Out[13]:
<matplotlib.collections.QuadMesh at 0x7f9ed50829a0>
In [14]:
#allSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['PHYSTRNO3'][:,:,:,:],3),2),1)
afiltSum=np.sum(np.sum(np.sum(tmaskSOG*fD.variables['ATF_NO3'][:,:,:,:],3),2),1)
radSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['NO3RAD'][:,:,:,:],3),2),1)
#no3Sum=np.sum(np.sum(np.sum(tmaskSOG*e1t*e2t*fP.variables['NO3_E3TSNAP'][:,:,:,:],3),2),1)
smsSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['SMS_NO3'][:,:,:,:],3),2),1)
nitrSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['REM_NO3'][:,:,:,:],3),2),1)
PPSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['PRD_NO3'][:,:,:,:],3),2),1)
PPSumNH=np.sum(np.sum(np.sum(tmaskSOG*(fP.variables['PRD_NH4'][:,:,:,:]),3),2),1)
rivSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['RIV_NO3'][:,:,:,:],3),2),1)
physSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['PHYSTRNO3'][:,:,:,:],3),2),1)
refrSum=np.sum(np.sum(np.sum(tmaskSOG*fP.variables['REFRN'][:,:,:,:],3),2),1)
In [15]:
np.min(rivSum),np.max(rivSum)
Out[15]:
(606.4505, 636.7251)
In [16]:
np.min
Out[16]:
<function numpy.amin(a, axis=None, out=None, keepdims=<no value>, initial=<no value>, where=<no value>)>
In [17]:
np.min(refrSum),np.max(refrSum)
Out[17]:
(-8.153931, -6.6515636)
In [18]:
np.max(radSum),np.min(radSum)
Out[18]:
(0.059010472, 0.0)
In [19]:
np.max(afiltSum),np.min(afiltSum)
Out[19]:
(-78895.67, -133606.23)
In [20]:
PPSum
Out[20]:
masked_array(data=[-1473524.25, -1839533.0, -1936257.125, -1770761.125,
-1434424.0],
mask=[False, False, False, False, False],
fill_value=1e+20,
dtype=float32)
In [22]:
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),PPSumNH,'c--',label='all')
ax[0].plot(range(0,5),refrSum,'y--',label='all')
#ax[0].plot(range(0,5),afiltSum,'b--',label='afilt')
ax[0].legend()
ax[1].plot(range(0,5),physSum,'k--',label='phys')
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),smsSum,'b--',label='sms')
#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].plot(range(0,5),rivSum,'c--',label='riv')
ax[2].legend()
Out[22]:
<matplotlib.legend.Legend at 0x7f9ed4e38f40>
In [54]:
with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_Malaspina_U*.nc')[0]) as 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)
print(np.min(f.variables['ATX_NO3'][:,:,:,0]),np.max(f.variables['ATX_NO3'][:,:,:,0]))
print(np.min(f.variables['ATX_NH4'][:,:,:,0]),np.max(f.variables['ATX_NH4'][:,:,:,0]))
print(np.min(f.variables['ATX_LIV'][:,:,:,0]),np.max(f.variables['ATX_LIV'][:,:,:,0]))
-24293.174 25672.338 -- -- -- --
In [24]:
with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_Haro_V*.nc')[0]) as f:
harVA=np.sum(np.sum(f.variables['ATY_NO3'][:,:,0,:],2),1)
harVD=np.sum(np.sum(f.variables['DTY_NO3'][:,:,0,:],2),1)
In [25]:
with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_SJC_V*.nc')[0]) as 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)
In [26]:
with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_Rosario_V*.nc')[0]) as 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)
In [27]:
with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_Sutil_V*.nc')[0]) as 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)
In [28]:
with nc.Dataset(glob.glob(idir+'/SalishSea_1d_*_Discovery_V*.nc')[0]) as 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)
In [29]:
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+rivSum,'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[29]:
<matplotlib.legend.Legend at 0x7f9ed4d32220>
In [30]:
adv[0],dif[0],rivSum[0],adv[0]+dif[0],adv[0]+dif[0]+rivSum[0],physSum[0]
Out[30]:
(1445421.1, 471.7424, 606.4505, 1445892.9, 1446499.4, 1493579.0)
In [31]:
## phys difference
adv+dif+rivSum-physSum
Out[31]:
masked_array(data=[-47079.625, -49406.15625, -49654.5546875,
-50259.703125, -50446.30078125],
mask=[False, False, False, False, False],
fill_value=1e+20,
dtype=float32)
In [34]:
fk=nc.Dataset(idir+'SalishSea_1d_20150501_20150505_layer6m_W_20150501-20150505.nc')
In [35]:
fk.variables.keys()
Out[35]:
dict_keys(['nav_lat', 'nav_lon', 'bounds_nav_lon', 'bounds_nav_lat', 'area', 'layer6m_W', 'layer6m_W_bounds', 'time_centered', 'time_centered_bounds', 'time_counter', 'time_counter_bounds', 'VMIXNO3', 'VMIXNH4', 'VMIXPON', 'VMIXDON', 'VMIXLIV', 'ATZ_NO3', 'ATZ_NH4', 'ATZ_PON', 'ATZ_DON', 'ATZ_LIV'])
In [36]:
fk.variables['VMIXNO3']
Out[36]:
<class 'netCDF4._netCDF4.Variable'>
float32 VMIXNO3(time_counter, layer6m_W, y, x)
long_name: upward flux due to vert mixing with horiz comp
units: mmol N/m2/s
online_operation: average
interval_operation: 40 s
interval_write: 1 d
cell_methods: time: mean (interval: 40 s)
cell_measures: area: area
_FillValue: 1e+20
missing_value: 1e+20
coordinates: time_centered nav_lat nav_lon
unlimited dimensions: time_counter
current shape = (5, 1, 898, 398)
filling on
In [40]:
plt.pcolormesh(fk.variables['ATZ_NO3'][4,0,:,:])
Out[40]:
<matplotlib.collections.QuadMesh at 0x7f9ed4fbb940>
In [38]:
plt.pcolormesh(fk.variables['ATZ_PON'][4,0,:,:])
Out[38]:
<matplotlib.collections.QuadMesh at 0x7f9ed4eb8670>
In [ ]:
In [ ]:
In [ ]:
In [42]:
fk2=nc.Dataset(idir+'SalishSea_1d_20150501_20150505_layer6m_W_2_20150501-20150505.nc')
In [43]:
fk2.variables.keys()
Out[43]:
dict_keys(['nav_lat', 'nav_lon', 'bounds_nav_lon', 'bounds_nav_lat', 'area', 'layer6m_W', 'layer6m_W_bounds', 'time_centered', 'time_centered_bounds', 'time_counter', 'time_counter_bounds', 'WLDFNO3', 'WLDFNH4', 'WLDFPON', 'WLDFDON', 'WLDFLIV', 'BSZ_DIAT', 'BSZ_PON', 'VMIXNO3H'])
In [47]:
np.min(fk2.variables['WLDFNH4']),np.max(fk2.variables['WLDFNH4'])
Out[47]:
(1e+20, 1e+20)
In [49]:
plt.pcolormesh(fk2.variables['WLDFNO3'][4,0,:,:])
Out[49]:
<matplotlib.collections.QuadMesh at 0x7f9ec5373970>
In [ ]: