#!/usr/bin/env python # coding: utf-8 # In[1]: import netCDF4 as nc import datetime as dt import subprocess import requests import matplotlib.pyplot as plt import cmocean import numpy as np import os import glob import dateutil as dutil from salishsea_tools import viz_tools get_ipython().run_line_magic('matplotlib', 'inline') # In[2]: with nc.Dataset('/ocean/eolson/MEOPAR/NEMO-forcing/grid/mesh_mask201702_noLPE.nc') as fm: tmask=np.copy(fm.variables['tmask']) umask=np.copy(fm.variables['umask']) vmask=np.copy(fm.variables['vmask']) navlon=np.copy(fm.variables['nav_lon']) navlat=np.copy(fm.variables['nav_lat']) dept=np.copy(fm.variables['gdept_1d']) e3t_0=np.copy(fm.variables['e3t_0']) e3u_0=np.copy(fm.variables['e3u_0']) e3v_0=np.copy(fm.variables['e3v_0']) e1t=np.copy(fm.variables['e1t']) e2t=np.copy(fm.variables['e2t']) e1v=np.copy(fm.variables['e1v']) e2u=np.copy(fm.variables['e2u']) A=fm.variables['e1t'][0,:,:]*fm.variables['e2t'][0,:,:]*tmask[0,0,:,:] # In[3]: t0=dt.datetime(2014,12,1) # 1st start date of run #te=dt.datetime(2016,12,1)# last start date of runfnum=18 stm=np.shape(tmask) SiN=2.0 #nlen=36*2 nlen=53 dlist=[t0+dt.timedelta(days=ii*10) for ii in range(0,nlen)] # In[4]: #sdir0='/results/SalishSea/nowcast-green/' sdir1='/results/SalishSea/hindcast/' #sdir3='/data/eolson/MEOPAR/SS36runs/CedarRuns/spring2015_HCMZ/' # In[5]: tmaskC=np.copy(tmask) tmaskC[:,:,370:490,:12]=0 tmaskC[:,:,887:,30:70]=0 # In[6]: tlist=dlist SiGlobalTot=dict() SiTot=dict() BSiTot=dict() DiatTot=dict() changeSiGlobalTot=dict() for idir in (sdir1,): fformat1='%d%b%y/' if idir.startswith('/data/eolson/MEOPAR/SS36runs/CedarRuns/'): fformatT='SalishSea_1h_*_ptrc_T_%Y%m%d-*.nc' fformatP='SalishSea_1h_*_ptrc_T_%Y%m%d-*.nc' #elif idir==sdir0: # fformatT='SalishSea_1h_%Y%m%d_%Y%m%d_ptrc_T.nc' # fformatP='SalishSea_1h_%Y%m%d_%Y%m%d_grid_T.nc' elif idir==sdir1: fformatT='SalishSea_1h_%Y%m%d_%Y%m%d_ptrc_T.nc' fformatP='SalishSea_1h_%Y%m%d_%Y%m%d_carp_T.nc' sumSi=np.zeros((len(tlist),stm[1])) sumBSi=np.zeros((len(tlist),stm[1])) sumDiat=np.zeros((len(tlist),stm[1])) ind=-1 for idt0 in tlist: ind=ind+1 cdir=idt0.strftime(fformat1).lower() iffT=idt0.strftime(fformatT) iffP=idt0.strftime(fformatP) if idir.startswith('/data/eolson/MEOPAR/SS36runs/CedarRuns/'): sffT=idir+iffT sffP=idir+iffP elif idir.startswith('/results/'): sffT=idir+cdir+iffT sffP=idir+cdir+iffP f=nc.Dataset(glob.glob(sffT)[0]) print(sffT) fP=nc.Dataset(glob.glob(sffP)[0]) #if idir==sdir0: # e3t=np.expand_dims((1+fP.variables['sossheig'][0,:,:]/np.sum(e3t_0*tmask,1)),0)*e3t_0 if idir==sdir1: e3t=fP.variables['e3t'][:2,:,:,:] Vol=A*np.ones(np.shape(e3t)) sumSi[ind,:]=1e-3*np.sum(np.sum(tmaskC[0,:,:,:]*Vol[0,:,:,:]*f.variables['silicon'][0,:,:,:],2),1) #mmol/m3*m3*10^-3=mol sumBSi[ind,:]=1e-3*np.sum(np.sum(tmaskC[0,:,:,:]*Vol[0,:,:,:]*f.variables['biogenic_silicon'][0,:,:,:],2),1) #mmol/m3*m3*10^-3=mol sumDiat[ind,:]=SiN*1e-3*np.sum(np.sum(tmaskC[0,:,:,:]*Vol[0,:,:,:]*f.variables['diatoms'][0,:,:,:],2),1) #mmol/m3*m3*10^-3=mol f.close() fP.close() SiGlobalTot[idir]=sumSi+sumBSi+sumDiat # In[7]: siplot=SiGlobalTot[idir]-SiGlobalTot[idir][0,:] # In[8]: cm1=cmocean.cm.balance #plt.pcolormesh([t0+dt.timedelta(10*ii) for ii in range(0,30)],dept[0,:],np.transpose(siplot[:30,:]),cmap=cm1,vmin=-1e9,vmax=1e9) plt.pcolormesh([t0+dt.timedelta(10*ii) for ii in range(0,nlen)],dept[0,:],np.transpose(siplot),cmap=cm1,vmin=-1e8,vmax=1e8) plt.ylim(400,0) plt.colorbar() # ## repeat for N # In[9]: tlist=dlist NGlobalTot=dict() VolTot=dict() NO3Tot=dict() NH4Tot=dict() PONTot=dict() DONTot=dict() DiatTot=dict() MyriTot=dict() NanoTot=dict() MiZoTot=dict() changeNGlobalTot=dict() for idir in (sdir1,): fformat1='%d%b%y/' if idir.startswith('/data/eolson/MEOPAR/SS36runs/CedarRuns/'): fformatT='SalishSea_1h_*_ptrc_T_%Y%m%d-*.nc' fformatP='SalishSea_1h_*_ptrc_T_%Y%m%d-*.nc' #elif idir==sdir0: # fformatT='SalishSea_1h_%Y%m%d_%Y%m%d_ptrc_T.nc' # fformatP='SalishSea_1h_%Y%m%d_%Y%m%d_grid_T.nc' elif idir==sdir1: fformatT='SalishSea_1h_%Y%m%d_%Y%m%d_ptrc_T.nc' fformatP='SalishSea_1h_%Y%m%d_%Y%m%d_carp_T.nc' sumNO3=np.zeros((len(tlist),stm[1])) sumVol=np.zeros((len(tlist),stm[1])) sumNH4=np.zeros((len(tlist),stm[1])) sumPON=np.zeros((len(tlist),stm[1])) sumDON=np.zeros((len(tlist),stm[1])) sumDiat=np.zeros((len(tlist),stm[1])) sumMyri=np.zeros((len(tlist),stm[1])) sumNano=np.zeros((len(tlist),stm[1])) sumMiZo=np.zeros((len(tlist),stm[1])) ind=-1 for idt0 in tlist: ind=ind+1 cdir=idt0.strftime(fformat1).lower() iffT=idt0.strftime(fformatT) iffP=idt0.strftime(fformatP) if idir.startswith('/data/eolson/MEOPAR/SS36runs/CedarRuns/'): sffT=idir+iffT sffP=idir+iffP elif idir.startswith('/results/'): sffT=idir+cdir+iffT sffP=idir+cdir+iffP f=nc.Dataset(glob.glob(sffT)[0]) print(sffT) fP=nc.Dataset(glob.glob(sffP)[0]) #if idir==sdir0: # e3t=np.expand_dims((1+fP.variables['sossheig'][0,:,:]/np.sum(e3t_0*tmask,1)),0)*e3t_0 if idir==sdir1: e3t=fP.variables['e3t'][:2,:,:,:] Vol=A*np.ones(np.shape(e3t)) sumVol[ind,:]=1e-3*np.sum(np.sum(tmaskC[0,:,:,:]*Vol[0,:,:,:],2),1) #mmol/m3*m3*10^-3=mol sumNO3[ind,:]=1e-3*np.sum(np.sum(tmaskC[0,:,:,:]*Vol[0,:,:,:]*f.variables['nitrate'][0,:,:,:],2),1) #mmol/m3*m3*10^-3=mol sumNH4[ind,:]=1e-3*np.sum(np.sum(tmaskC[0,:,:,:]*Vol[0,:,:,:]*f.variables['ammonium'][0,:,:,:],2),1) #mmol/m3*m3*10^-3=mol sumPON[ind,:]=1e-3*np.sum(np.sum(tmaskC[0,:,:,:]*Vol[0,:,:,:]*f.variables['particulate_organic_nitrogen'][0,:,:,:],2),1) #mmol/m3*m3*10^-3=mol sumDON[ind,:]=1e-3*np.sum(np.sum(tmaskC[0,:,:,:]*Vol[0,:,:,:]*f.variables['dissolved_organic_nitrogen'][0,:,:,:],2),1) #mmol/m3*m3*10^-3=mol sumDiat[ind,:]=1e-3*np.sum(np.sum(tmaskC[0,:,:,:]*Vol[0,:,:,:]*f.variables['diatoms'][0,:,:,:],2),1) #mmol/m3*m3*10^-3=mol sumMyri[ind,:]=1e-3*np.sum(np.sum(tmaskC[0,:,:,:]*Vol[0,:,:,:]*f.variables['ciliates'][0,:,:,:],2),1) #mmol/m3*m3*10^-3=mol sumMiZo[ind,:]=1e-3*np.sum(np.sum(tmaskC[0,:,:,:]*Vol[0,:,:,:]*f.variables['microzooplankton'][0,:,:,:],2),1) #mmol/m3*m3*10^-3=mol f.close() fP.close() NGlobalTot[idir]=sumNO3+sumNH4+sumPON+sumDON+sumDiat+sumMyri+sumNano+sumMiZo # In[10]: Nplot=NGlobalTot[idir]-NGlobalTot[idir][0,:] # In[11]: #plt.pcolormesh([t0+dt.timedelta(10*ii) for ii in range(0,30)],dept[0,:],np.transpose(siplot[:30,:]),cmap=cm1,vmin=-1e9,vmax=1e9) plt.pcolormesh([t0+dt.timedelta(10*ii) for ii in range(0,nlen)],dept[0,:],np.transpose(Nplot),cmap=cm1,vmin=-1e8,vmax=1e8) plt.ylim(400,0) plt.colorbar() # In[ ]: # In[ ]: # In[ ]: