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 [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]:

temp=np.copy(tmask[35,:,:])
temp[450:475,260:275]=2
plt.pcolormesh(temp)
plt.colorbar()

Out[3]:

<matplotlib.colorbar.Colorbar at 0x7ff7fa23cc10>

In [4]:

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

In [5]:

fP=nc.Dataset(glob.glob(idir+'/SalishSea_1h_*ptrc_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', 'NO3', 'NH4', 'PHY', 'DIAT', 'Si', 'bSi', 'PON'])

In [7]:

np.shape(e3t_0)

Out[7]:

(1, 40, 898, 398)

In [8]:

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(fP.variables['DIAT'][0,:,jj,ii],-1*fP.variables['deptht'][:])
        ax[1].plot(fP.variables['DIAT'][1*24,:,jj,ii],-1*fP.variables['deptht'][:])
        ax[2].plot(fP.variables['DIAT'][2*24,:,jj,ii],-1*fP.variables['deptht'][:])
        ax[3].plot(fP.variables['DIAT'][3*24,:,jj,ii],-1*fP.variables['deptht'][:])
        ax[4].plot(fP.variables['DIAT'][4*24,:,jj,ii],-1*fP.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(fP.variables['DIAT'][0,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['DIAT'][0*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z1=np.mean(np.mean(np.nansum(fP.variables['DIAT'][1*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['DIAT'][1*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z2=np.mean(np.mean(np.nansum(fP.variables['DIAT'][2*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['DIAT'][2*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z3=np.mean(np.mean(np.nansum(fP.variables['DIAT'][3*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['DIAT'][3*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z4=np.mean(np.mean(np.nansum(fP.variables['DIAT'][4*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['DIAT'][4*24,:,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
6.552980618263962 10.393471727376022 11.6143544998826 13.394177148332306
Diagnosed sinking rate:  1.789894259933398 m/d

In [9]:

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(fP.variables['bSi'][0,:,jj,ii],-1*fP.variables['deptht'][:])
        ax[1].plot(fP.variables['bSi'][1*24,:,jj,ii],-1*fP.variables['deptht'][:])
        ax[2].plot(fP.variables['bSi'][2*24,:,jj,ii],-1*fP.variables['deptht'][:])
        ax[3].plot(fP.variables['bSi'][3*24,:,jj,ii],-1*fP.variables['deptht'][:])
        ax[4].plot(fP.variables['bSi'][4*24,:,jj,ii],-1*fP.variables['deptht'][:])
for iax in ax:
    iax.set_ylim(-450,0)
print('bSi sinking rate: 11.23m/d')
z0=np.mean(np.mean(np.nansum(fP.variables['bSi'][0,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['bSi'][0*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z1=np.mean(np.mean(np.nansum(fP.variables['bSi'][1*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['bSi'][1*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z2=np.mean(np.mean(np.nansum(fP.variables['bSi'][2*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['bSi'][2*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z3=np.mean(np.mean(np.nansum(fP.variables['bSi'][3*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['bSi'][3*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z4=np.mean(np.mean(np.nansum(fP.variables['bSi'][4*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['bSi'][4*24,:,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')

bSi sinking rate: 11.23m/d
17.516393791055954 37.02603073822324 46.582902057490436 65.82292345888409
Diagnosed sinking rate:  14.871882742617911 m/d

In [10]:

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(fP.variables['PON'][0,:,jj,ii],-1*fP.variables['deptht'][:])
        ax[1].plot(fP.variables['PON'][1*24,:,jj,ii],-1*fP.variables['deptht'][:])
        ax[2].plot(fP.variables['PON'][2*24,:,jj,ii],-1*fP.variables['deptht'][:])
        ax[3].plot(fP.variables['PON'][3*24,:,jj,ii],-1*fP.variables['deptht'][:])
        ax[4].plot(fP.variables['PON'][4*24,:,jj,ii],-1*fP.variables['deptht'][:])
for iax in ax:
    iax.set_ylim(-450,0)
print('PON sinking rate: 8.64 m/d')
z0=np.mean(np.mean(np.nansum(fP.variables['PON'][0,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['PON'][0*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z1=np.mean(np.mean(np.nansum(fP.variables['PON'][1*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['PON'][1*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z2=np.mean(np.mean(np.nansum(fP.variables['PON'][2*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['PON'][2*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z3=np.mean(np.mean(np.nansum(fP.variables['PON'][3*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['PON'][3*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z4=np.mean(np.mean(np.nansum(fP.variables['PON'][4*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['PON'][4*24,:,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')

PON sinking rate: 8.64 m/d
14.450184783677994 29.323760222528353 35.655545388018474 53.70322165840983
Diagnosed sinking rate:  11.848672222833674 m/d

In [11]:

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(fP.variables['NO3'][0,:,jj,ii],-1*fP.variables['deptht'][:])
        ax[1].plot(fP.variables['NO3'][1*24,:,jj,ii],-1*fP.variables['deptht'][:])
        ax[2].plot(fP.variables['NO3'][2*24,:,jj,ii],-1*fP.variables['deptht'][:])
        ax[3].plot(fP.variables['NO3'][3*24,:,jj,ii],-1*fP.variables['deptht'][:])
        ax[4].plot(fP.variables['NO3'][4*24,:,jj,ii],-1*fP.variables['deptht'][:])
for iax in ax:
    iax.set_ylim(-450,0)
print('NO3 sinking rate: 0 m/d')
z0=np.mean(np.mean(np.nansum(fP.variables['NO3'][0,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['NO3'][0*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z1=np.mean(np.mean(np.nansum(fP.variables['NO3'][1*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['NO3'][1*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z2=np.mean(np.mean(np.nansum(fP.variables['NO3'][2*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['NO3'][2*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z3=np.mean(np.mean(np.nansum(fP.variables['NO3'][3*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['NO3'][3*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270],0)))
z4=np.mean(np.mean(np.nansum(fP.variables['NO3'][4*24,:,455:460,265:270]*e3t_0[0,:,455:460,265:270]*np.reshape(fP.variables['deptht'][:],(40,1,1)),0)/np.sum(fP.variables['NO3'][4*24,:,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')

NO3 sinking rate: 0 m/d
5.186553668517444 8.814400683700363 9.29954783606704 10.110866587230795
Diagnosed sinking rate:  0.9718796187345566 m/d

In [ ]: