#!/usr/bin/env python # coding: utf-8 # 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 get_ipython().run_line_magic('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_b/' # 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() # In[6]: fD.variables.keys() # In[7]: fD.variables['BFX_PON'] # In[8]: np.min(fD.variables['BFX_PON']),np.max(fD.variables['BFX_PON']) # In[9]: np.min(fD.variables['ATF_LIV']),np.max(fD.variables['ATF_LIV']) # In[10]: np.min(fD.variables['ATF_PON']),np.max(fD.variables['ATF_NH4']) # In[11]: np.min(fD.variables['ATF_NO3']),np.max(fD.variables['ATF_NO3']) # In[12]: np.min(fD.variables['BFX_PON']),np.max(fD.variables['BFX_DIAT']) # In[13]: plt.pcolormesh(fD.variables['BFX_PON'][4,:,:]) # In[14]: fM.variables.keys() # In[15]: np.min(tmaskSOG*fM.variables['MIZ_DIAT']),np.max(tmaskSOG*fM.variables['MIZ_DIAT']) # In[16]: plt.pcolormesh(fM.variables['MIZ_DIAT'][4,0,:,:]) # In[17]: #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[18]: np.min(rivSum),np.max(rivSum) # In[19]: np.min # In[20]: np.min(refrSum),np.max(refrSum) # In[21]: np.max(radSum),np.min(radSum) # In[22]: np.max(afiltSum),np.min(afiltSum) # In[23]: PPSum # In[24]: 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() # In[25]: 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])) # In[26]: 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[27]: 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[28]: 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[29]: 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[30]: 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[31]: 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() # In[32]: adv[0],dif[0],rivSum[0],adv[0]+dif[0],adv[0]+dif[0]+rivSum[0],physSum[0] np.max(np.abs(harVA-harTVDY))/np.max(np.abs(harTVDY))*100np.max(np.abs(harVD-harYM))/np.max(np.abs(harYM))*100 # In[33]: ## phys difference adv+dif+rivSum-physSum # In[34]: fk=nc.Dataset(idir+'SalishSea_1d_20150501_20150505_layer6m_W_20150501-20150505.nc') # In[35]: fk.variables.keys() # In[36]: fk.variables['VMIXNO3'] # In[37]: plt.pcolormesh(fk.variables['ATZ_NO3'][4,0,:,:]) # In[38]: plt.pcolormesh(fk.variables['ATZ_PON'][4,0,:,:]) # In[ ]: # In[ ]: # In[ ]: # In[39]: fk2=nc.Dataset(idir+'SalishSea_1d_20150501_20150505_layer6m_W_2_20150501-20150505.nc') # In[40]: fk2.variables.keys() # In[41]: np.min(fk2.variables['WLDFNH4']),np.max(fk2.variables['WLDFNH4']) # In[42]: plt.pcolormesh(fk2.variables['WLDFNO3'][4,0,:,:]) # In[ ]: