# Plot time series plot of Figure 13 in the paper 

import numpy as np
import numpy.ma as ma
#import Nio,Ngl,os
import matplotlib.pyplot as plt
plt.switch_backend('agg')

from netCDF4 import Dataset

# ---------------- Data source in rhea, the two files save the same data ------------------------
input_file11 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/tseries/20180215.DECKv1b_H1.ne30_oEC.edison.clm2.aavg.h0.yrly.RAIN.nc'
input_file12 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/tseries/20180215.DECKv1b_H1.ne30_oEC.edison.clm2.aavg.h0.yrly.SNOW.nc'
input_file13 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/tseries/20180215.DECKv1b_H1.ne30_oEC.edison.clm2.aavg.h0.yrly.QRUNOFF.nc'
input_file14 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/tseries/20180215.DECKv1b_H1.ne30_oEC.edison.clm2.aavg.h0.yrly.QSOIL.nc'
input_file2 = '/lustre/atlas1/cli115/world-shared/4ue/racmo23_GRN_monthly/tseries/racmo23_GRN_monthly.smb.ann_tseries_aavg.nc'
# --------------------------------------------------------------

# read input_file11 
# time from X-X, both time and acab have 156 data
ncid11 = Dataset(input_file11, mode='r')
time = ncid11.variables['mcdate'][:]
rain = ncid11.variables['RAIN'][:] #m/yr
ncid11.set_auto_scale(True) # put the scale-factor in the time-series plot

#time = (time-time[1])/365 + 1851 # time from 1851-2006
timeB = [x for x in range(1, 1+len(time))] 
timeB = [x + 1960 for x in timeB]
#print "timeB:", timeB

# read input_file12 
# time from X-X, both time and acab have 156 data
ncid12 = Dataset(input_file12, mode='r')
time = ncid12.variables['time'][:]
snow = ncid12.variables['SNOW'][:] #m/yr
ncid12.set_auto_scale(True) # put the scale-factor in the time-series plot

# read input_file13 
# time from X-X, both time and acab have 156 data
ncid13 = Dataset(input_file13, mode='r')
qrun = ncid13.variables['QRUNOFF'][:] #m/yr
ncid13.set_auto_scale(True) # put the scale-factor in the time-series plot

# read input_file14 
# time from X-X, both time and acab have 156 data
ncid14 = Dataset(input_file14, mode='r')
qsoil = ncid14.variables['QSOIL'][:] #m/yr
ncid14.set_auto_scale(True) # put the scale-factor in the time-series plot

# read input_file2 
ncid2 = Dataset(input_file2, mode='r')
time2 = ncid2.variables['time'][:]
smb = ncid2.variables['smb'][:]  # mmWE
#ncid2.set_auto_scale(True) # put the scale-factor in the time-series plot

#smb = smb #unit transfer
time2 = (time2+6.5)/12 + 1957 # add months to make the full year, divide my months to get years, add to get 1958
#print "time2:", time2

qice = (rain + snow - qrun - qsoil)*3.1567e7 # convert to kg/m^2/yr

maxqice=(np.max(qice))
minqice=(np.min(qice))
maxtime=(np.max(timeB))
mintime=(np.min(timeB))
maxsmb=(np.max(smb))
minsmb=(np.min(smb))
maxtime2=(np.max(time2))
mintime2=(np.min(time2))

print "Min QICE E3SM:", minqice
print "Max QICE E3SM:", maxqice
#print "Max time:", maxtime
#print "Min time:", mintime
print "Min smb RACMO:", minsmb
print "Max smb RACMO:", maxsmb
#print "Max time2:", maxtime2
#print "Min time2:", mintime2

##------- Time Series PLOT --------
fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(111)

ax.plot(timeB,qice,'r-')
ax.plot(time2,smb,'b-')

ax.set_ylim([-100, 400])
ax.set_xlim([1955, 2015])
ax.xaxis.set_ticks(np.arange(1950,2010,10))
#ax.set_xticklabels(['1950','1960','1970','1980','1990','2000','2010'])
ax.set_ylabel('SMB (kg/m2/yr)',fontsize=14)
ax.set_xlabel('(years)',fontsize=14)
##ax.text(-12,5,'T-2m, JJA',fontsize=18)
##ax.legend([s1,s2],['>20%','>99%'],loc='lower right')

plt.savefig('/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/plots/TimeSeries_smblike.ps')
#plt.savefig('TimeSeries_SMB.ps')
plt.close()

##------ Box plot of statistics ---------
## If I have 3 sets of data, put them a matrix with 3 columns,
## then boxplot(matrix) in one figure data=[a1,a2,a3]

fig = plt.figure(figsize=(4.8,4))

#--- the whisker shows [Q1-1.5IQR, Q3+1.5IQR] without outliners
#ax = fig.add_subplot(121)
#ax.boxplot(qice,0,'')
#ax.set_ylim([0,0])
#ax.set_ylabel('qice(kg/m2/yr)',fontsize=14)
#ax.set_xticklabels('')

#--- the whisker shows [Q1-1.5IQR, Q3+1.5IQR] with outliners
#--- where the Q1 is 25% percentile and Q3 is 75% percential, IQR=Q3-Q1
#--- the outlier is data point located outside the whiskers.
ax = fig.add_subplot(121)
ax.boxplot(qice,0,'gd')
ax.set_ylabel('E3SM SMB(km/m2/yr)',fontsize=14)
ax.set_ylim([0,500])
ax.set_xticklabels('')
ax = fig.add_subplot(122)
ax.boxplot(smb,0,'gd')
ax.set_ylabel('RACMO SMB(km/m2/yr)',fontsize=14)
ax.set_ylim([0,500])
ax.set_xticklabels('')

plt.savefig('/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/plots/BoxPlot_smblike.ps')
#plt.savefig('BoxPlot_SMB.ps')
plt.close()