#!/python3
# read and plot time series of total KE  of GLBb0.08
# from log files

# to avoid attribute error due to version in proplot ...
from importlib.metadata import version
from packaging.version import parse

import sys
prfx='/discover/nobackup/projects/gmao/cal_ocn/abozec1/PYTHON/'
sys.path.append(prfx)

import os
import proplot as plot
import netCDF4 as nc
import numpy as np
import xarray as xr
import matplotlib.dates as mdates
import datetime


print('Starting here')
## Relative winds path
expt='10.6'
io='/discover/nobackup/projects/gmao/cal_ocn/abozec1/MOM6-expt/GLBb0.08/expt_'+expt+'/'
file='KE_energy_'+expt+'_2010-2017.log'

# open a first time to get the size
f=open(io+file,'r')
# Loop over lines and extract variables of interest
t=0
for line in f:
    columns = line.split()
    t=t+1
f.close()
print('t:',t)

date0=np.zeros([t]) ; KE0=np.zeros([t])
pt0=np.zeros([t])   ; sal0=np.zeros([t])
# populate the array
f=open(io+file,'r')
# Loop over lines and extract variables of interest
t=0
for line in f:
    columns = line.split()
    #date0[t]=columns[2]
    #KE0[t]=float(columns[6][:-1])*1e4
    pt0[t]=float(columns[14])
    sal0[t]=float(columns[12][:-1])
    t=t+1
f.close()
time0=np.arange(t)/365+2010

## Absolute winds path
expt='10.7'
io='/discover/nobackup/projects/gmao/cal_ocn/abozec1/MOM6-expt/GLBb0.08/expt_'+expt+'/'
file='KE_energy_'+expt+'_2010-2017.log'

# open a first time to get the size
f=open(io+file,'r')
# Loop over lines and extract variables of interest
t=0
for line in f:
    columns = line.split()
    t=t+1
f.close()
print('t:',t)

date1=np.zeros([t]) ; KE1=np.zeros([t])
pt1=np.zeros([t])   ; sal1=np.zeros([t])
# populate the array
f=open(io+file,'r')
# Loop over lines and extract variables of interest
t=0
for line in f:
    columns = line.split()
    #date1[t]=columns[2]
    #KE1[t]=float(columns[6][:-1])*1e4
    pt1[t]=float(columns[14])
    sal1[t]=float(columns[12][:-1])
    t=t+1
f.close()
time1=np.arange(t)/365+2010


## plot
fig, ax = plot.subplots(ncols=1,nrows=2,refaspect=3, refwidth='15cm')
ax[0].plot(time0[:],pt0[:], lw=2,legend='ur',labels=['potT rel. winds'])
ax[0].plot(time1[:],pt1[:], lw=2,legend='ur',labels=['potT abs. winds'])
ax[0].format(title='Pot. Temperature ', xlabel='Time (days)', \
          ylabel='Pot. T (C)',ylim=(3.5,4.3),xlim=(2010,2018))


ax[1].plot(time0[:],sal0[:], lw=2,legend='lr',labels=['sal rel. winds'])
ax[1].plot(time1[:],sal1[:], lw=2,legend='lr',labels=['sal abs. winds'])
ax[1].format(title='Salinity ', xlabel='Time (days)', \
          ylabel='Salinity (psu)',ylim=(34.71, 34.77),xlim=(2010,2018))



ps_dir=prfx+'/PNG/GLBb0.08/'
file_ps='time_TS_mom6_glb8_106-107_2010-2017.png'
print(file_ps)
fig.savefig(ps_dir+file_ps,dpi=150,\
            facecolor='w', edgecolor='w',transparent=False) ## .pdf,.eps
plot.show()