# First, calculate river concentrations from EC data at Hope and Gravesend
# NOTE: am assuming all forms of N are reported in grams of N, so that mole conversion is 
#  by division by 14.006720 (mol. weight N), rather than weights of full molecules. This 
#  is consistent with NO23 values in Susan's notebook.
#  This is consistent with report in 'precdetectlimits.pdf', although it is rather old.
#  However, assume Si in grams of SiO2 (mo. weight = 60.08). This is consistent with 
#  Si data in Susan's notebook (http://nbviewer.jupyter.org/url/bitbucket.org/salishsea/tools/
#     raw/tip/I_ForcingFiles/Rivers/FraserRiverNutrients.ipynb?at=default&fileviewer=file-
#     view-default)
# Also appears to result in rough consistency with Debby's river Si and NO3 values
mwN=14.006720
#mwSiO3=76.083820
mwSiO2=60.08
#
#  nh4 muM  (N) = nh4_cst from mean of data at Hope/Gravesend
#  don muM  (N) = don_cst from mean of data at Hope/Gravesend
#     - only 1 data point available. Should use it or set to zero? Sutton (2013) suggests DON may be unreactive on timescales
#       of exchange of water in SoG, but Clair (2013) reports that organic nitrogen dominated total nitrogen in 
#       Canadian rivers.
#       -> for now, using available data point
#  no3 muM  (N) -> daily from smoothed data
#  sil muM (Si) -> daily from smoothed data
#  dia muM  (N) = 0.001
dia = 0.0 #0.001
#  phy muM  (N) = 0.001
phy = 0.0 #0.001
#  mes muM  (N) = 0.001
mes = 0.0 #0.001
#  zoo muM  (N) = 0.001
zoo = 0.0 #0.001
#  pon muM  (N) = 0.0
pon = 0.0
#  bsi muM (Si) = 0.0
don=0.0
bsi = 0.0
#  oxy ? empty
oxy = 160
# other rivers const vals:
const_si=55.69
const_no=6.26
const_nh=4.09
const_tur=1.0 # default river turbidity (NTU)

# load EC Rivers database
dbpath='/ocean/eolson/MEOPAR/obs/ECRivers/ECRiversDB'
engine = create_engine('sqlite:///'+dbpath+'.sqlite')
Base = automap_base()
Base.prepare(engine, reflect=True)
Profs=Base.classes.profiles
session = create_session(bind = engine, autocommit = False, autoflush = True)
# relevant data types: 
#('NITROGEN DISSOLVED NITRATE',)
#('NITROGEN DISSOLVED NO3 & NO2',)
#('NITROGEN DISSOLVED ORGANIC (CALCD.)',)
#('NITROGEN NITRITE',)
#('NITROGEN TOTAL',)
#('NITROGEN TOTAL DISSOLVED',)
#('NITROGEN TOTAL ORGANIC (CALCD.)',)
#('SALINITY',)
#('SILICA DISSOLVED',)
#('SILICA EXTRACTABLE',)
#('SILICA REACTIVE',)
#('SILICON DISSOLVED',)
#('SILICON EXTRACTABLE',)
#('SILICON TOTAL',)
#('TURBIDITY',)

#define function to plot results of query as repeating fxn of year-day if output rows contain Year, Month, Day, value:
def plotYMDV(qout,pcolor,ptitle):
    Y=[]
    M=[]
    D=[]
    Amm=[]
    date=[]
    for Y0, M0, D0, Amm0 in qout:
        Y.append(Y0)
        M.append(M0)
        D.append(D0)
        Amm.append(Amm0)
        date.append(dt.date(Y0,M0,D0))
    Y=np.array(Y)
    M=np.array(M)
    D=np.array(D)
    Amm=np.array(Amm)
    date=np.array(date)
    YD=0.0*Y
    for i in range(0,len(Y)):
        YD[i]=date[i].timetuple().tm_yday
    plt.plot(YD,Amm,'.',color=pcolor)
    plt.plot(YD+365.0,Amm,'.',color=pcolor)
    plt.title(ptitle)
    plt.plot((366,366),(0,np.max(Amm)),'k--')
    return(YD,Amm)
    
# definte function to smooth data with gaussian filter, plot result, and return dict of days, values
def gsmooth(YD,Amm,L):
    allt=np.arange(1,367)
    fil=np.empty(np.size(allt))
    s=L/2.355
    sdict={}
    for t in allt:
        diff=[min(abs(x-t),abs(x-t+365), abs(x-t-365)) for x in YD]
        weight=[np.exp(-.5*x**2/s**2) if x <= 3*L else 0.0 for x in diff]
        #weight=[np.exp(-.5*x**2/s**2) for x in diff]
        fil[t-1]=np.sum(weight*Amm)/np.sum(weight)
        sdict[t]=fil[t-1]
    plt.plot(allt,fil,'k-')
    plt.plot(allt+365,fil,'k-')
    diff_ex=[min(abs(x-t),abs(x-t+365), abs(x-t-365)) for x in allt]
    weight_ex=[np.exp(-.5*x**2/s**2) if x <= 3*L else 0.0 for x in diff_ex]
    plt.plot(allt,weight_ex/np.sum(weight_ex)*np.max(2.0*Amm),'k:')
    plt.plot(allt+365,weight_ex/np.sum(weight_ex)*np.max(2.0*Amm),'k:')
    return(sdict)
    #plt.plot(YD,weight,'*')

# create constants file
f=nc.Dataset('/results/forcing/rivers/R201702DFraCElse_y2016m01d23.nc') # example for dims
year=2016
fnameC='/data/eolson/MEOPAR/NEMO-forcing-new/rivers-climatology/bio/R201812_bioConst'
new=nc.Dataset(fnameC+'.nc','w')
#Copy dimensions
for dname, the_dim in f.dimensions.items():
    #print (dname, len(the_dim) if not the_dim.isunlimited() else None)
    new.createDimension(dname, len(the_dim) if not the_dim.isunlimited() else None)

#new_tc=new.createVariable('time_counter',float,('time_counter'),zlib=True)
#new_tc[:]=f.variables['time_counter'][:]

new_run=new.createVariable('no3',float,('time_counter', 'y', 'x'),zlib=True)
new_run[:,:,:]=const_no # other rivers
new_run[:,400:520,380:]=np.mean([no3dict[yearday]/mwN*1000.0 for yearday in range(1,367)])

new_run=new.createVariable('nh4',float,('time_counter', 'y', 'x'),zlib=True)
new_run[:,:,:]=const_nh # other river value
new_run[:,400:520,380:]=nh4_cst/mwN*1000.0 # convert to muM; Fraser value

new_run=new.createVariable('sil',float,('time_counter', 'y', 'x'),zlib=True)
new_run[:]=const_si # other river value
new_run[:,400:520,380:]=np.mean([sidict[yearday]/mwSiO2*1000.0 for yearday in range(1,367)])

new_run=new.createVariable('dia',float,('time_counter', 'y', 'x'),zlib=True)
new_run[:]=dia

new_run=new.createVariable('phy',float,('time_counter', 'y', 'x'),zlib=True)
new_run[:]=phy

new_run=new.createVariable('mes',float,('time_counter', 'y', 'x'),zlib=True)
new_run[:]=mes

new_run=new.createVariable('zoo',float,('time_counter', 'y', 'x'),zlib=True)
new_run[:]=zoo

new_run=new.createVariable('don',float,('time_counter', 'y', 'x'),zlib=True)
new_run[:]=don

new_run=new.createVariable('pon',float,('time_counter', 'y', 'x'),zlib=True)
new_run[:]=pon

new_run=new.createVariable('bsi',float,('time_counter', 'y', 'x'),zlib=True)
new_run[:]=bsi

new_run=new.createVariable('tur',float,('time_counter', 'y', 'x'),zlib=True)
new_run[:]=const_tur

new.close()

session.close()
engine.dispose()

f.close()
mesh.close()