# -*- coding: utf-8 -*- 

import numpy as nm
import matplotlib.pyplot as plt # for plotting
import time as tm
import os
from phys import *
import matplotlib.dates as dts
import datetime as dtm
from matplotlib import rc

#Checking if string is a number
def is_number(s):
        try:
                float(s)
                return True
        except ValueError:
                return False

def runningMeanFast(x, N):
	work = nm.convolve(x, nm.ones((N,))/N)[(N-1):]
	work1 = nm.zeros(len(x))
	work1[N-1:] = work[N-1:]
	for i in range(N-1):
		work1[i] = work1[N-1]
	return work1
	

def timeoffset(time,t0disp):
	# Calculating the time offset
	it0disp = len(time)-1
	for i in range(len(time)):
        	if time[i] > t0disp:
                	it0disp = i
                	break
	return it0disp

def month_doy(nmonth,ny):
	#Day of year at the beginning of month nmonth
	ndm = [31,28,31,30,31,30,31,31,30,31,30,31]
	if ny%4 == 0 :#leap year
		ndm[1] = 29
	doy = 0
	if nmonth > 1:
		for i in range(nmonth-1):
			doy += ndm[i]
	#print (doy, nmonth)
	return doy

def readheader(fid,nheader) :
#Reading the header
	for j in range(nheader):
		line = fid.readline()

def plot_timeser(plotid,nvar,legendcommon) :

	font = {'family': 'serif',
        	'weight': 'normal'}
	rc('font', **font)
	rc('text', usetex=True)
	rc('text.latex',unicode=True)
	rc('text.latex',preamble='\usepackage[utf8]{inputenc}')
	rc('text.latex',preamble='\usepackage[russian]{babel}')
	#rc('text.latex',preamble='\usepackage{amsmath}')

	#Numbering of plot items
	num = 1
	tempmod1_ = num; num +=1
	tempmod2_ = num; num +=1
	tempmod3_ = num; num +=1
	senslat_ = num; num +=1
	uspeed_ = num; num +=1
	ustar_ = num; num +=1
	windwork_ = num; num +=1
	mldmod_ = num; num +=1
	thermthickmod_ = num; num +=1
	wumod_ = num; num +=1
	rimldmod_ = num; num +=1
	maxbvfmod_ = num; num +=1
	mincondmod_ = num; num +=1
	tkemldmod_ = num; num +=1
	Wmod_ = num; num +=1
	LNmod_ = num; num +=1
	Rmod_ = num; num +=1
	ch4mod_ = num; num +=1
	o2mod_ = num; num +=1
	co2mod_ = num; num +=1
	hshnmod_ = num; num +=1
	hlsnmod_ = num; num +=1
	ebulmod_ = num; num +=1
	ch4flxmod_ = num; num +=1
	ch4flxMLmod_ = num; num +=1
	ch4prodmod_ = num; num +=1
	rimod_ = num; num +=1
	condmod_ = num; num +=1
	icethickmod_ = num; num +=1
	snowthickmod_ = num; num +=1

        nmod_max = num-1 #The last item in the list above, from the modeled data

	tempobs_ = num; num +=1
	mldobs_ = num; num +=1
	thermthickobs_ = num; num +=1
	maxbvfobs_ = num; num +=1
	co2obs1_ = num; num +=1
	co2obs2_ = num; num +=1
	ch4obs_ = num; num +=1
	o2obs_ = num; num +=1
	ch4flxobs_ = num; num +=1
	co2flxobs_ = num; num +=1
	shflxobs_ = num; num +=1
	ustarobs_ = num; num +=1
	
	# Choosing a plot number to display
	#if not ( nplot in locals() ):
	nplot = input("Choose a variable to plot \n" +	\
	 str(tempmod1_) + " - Temperature of surface, snow, ice and sediments \n" + \
	 str(tempmod2_) + " - Temperature of water at different depths \n" + \
	 str(tempmod3_) + " - Temperature of water at selected points \n" + \
	 str(senslat_) + " - Sensible and latent heat fluxes \n" + \
	 str(uspeed_) + " - x-component of speed \n" + \
	 str(ustar_) + " - Friction velocity \n" + \
	 str(windwork_) + " - Work by wind at the water surface \n" + \
	 str(mldmod_) + " - Mixed-layer depth, m \n" + \
	 str(thermthickmod_) + " - Thermocline thickness, m \n" + \
	 str(wumod_) + " - w_*/u_*, \n" + \
	 str(rimldmod_) + " - bulk Richardson number in the mixed layer, \n" + \
	 str(maxbvfmod_) + " - maximal Brunt-V\"ais\"ala frequency, \n" + \
	 str(mincondmod_) + " - minimal thermal conductance, \n" + \
	 str(tkemldmod_) + " - TKE terms averaged in the mixed layer, \n" + \
	 str(Wmod_) + " - Wedderburn number, \n" + \
	 str(LNmod_) + " - Lake number, \n" + \
	 str(Rmod_) + " - internal Rossby radius, m, \n" + \
	 str(ch4mod_) + " - CH_4 concentration \n" + \
	 str(o2mod_) + " - O_2 concentration \n" + \
	 str(co2mod_) + " - CO_2 concentration \n" + \
	 str(hshnmod_) + " - h_S&h_N \n" + \
	 str(hlsnmod_) + " - thickness of water, ice, snow and deep ice layers \n" + \
	 str(ebulmod_) + " - surface ebullition fluxes above different soil columns \n" + \
	 str(ch4flxmod_) + " - methane fluxes \n" + \
	 str(ch4flxMLmod_) + " - methane budget in the mixed layer \n" + \
	 str(ch4prodmod_) + " - methane production in sediments \n" + \
         str(rimod_) + " - Richardson number \n" + \
	 str(condmod_) + " - heat conductivity \n" + \
	 str(icethickmod_) + " - ice thickness \n" + \
	 str(snowthickmod_) + " - snow thickness \n" + \
	 str(tempobs_) + " - temperature of water measured \n" + \
	 str(mldobs_) + " - mixed-layer depth, measured \n" + \
	 str(thermthickobs_) + " - thermocline thickness, measured \n" + \
	 str(maxbvfobs_) + " - maximal Brunt-V\"ais\"ala frequency, measured \n" + \
	 str(co2obs1_) + " - CO_2 in water measured \n" + \
	 str(co2obs2_) + " - CO_2 in water measured manually \n" + \
	 str(ch4obs_) + " - CH_4 in water measured manually \n" + \
	 str(o2obs_) + " - O_2 in water measured manually \n" + \
	 str(ch4flxobs_) + " - CH_4 fluxes measured \n" + \
	 str(co2flxobs_) + " - CO_2 fluxes measured \n" + \
	 str(shflxobs_) + " - sensible heat flux measured \n" + \
	 str(ustarobs_) + " - friction velocity measured \n")
	
	#path = '../results/Seida1_2007-2008/time_series/' # path to files
	# Number of experiments is limited to 10
	#pathobs = [os.path.expanduser('/home/victor/Files/main/observdata/Mojai2016/')]
	pathobs = [os.path.expanduser('~/Files/main/observdata/Kuivajarvi/campaign2014/')]
	#path_ = ['../results/Kuivajarvi2013_c2/time_series/']#, \
	path_ = ['../results/Kuiva2014septnew/time_series/'] #, \
		 #'../results/Valkea-Kotinen2006_nomomflxpart/time_series/', \
		 #'../results/Valkea-Kotinen2006_nokor/time_series/', \
		 #'../results/Valkea-Kotinen2006_dynpgrad4nokor/time_series/', \
		 #'../results/Valkea-Kotinen2006_dynpgrad4/time_series/', \
		 #'../results/Valkea-Kotinen2006_nomomflxpartdynpgrad4/time_series/']
	         #'../results/BolshoiVilui20153_extwat1/time_series/', \
	         #'../results/BolshoiVilui20153_extwat025/time_series/']#, \
	         #'../results/BolshoiVilui20153_pblpar1/time_series/']
	
	nameexpobs = ['obs']
	nameexp_ = [u'base',u'nomomflxpart',u'nokor',u'dynpgrad4nokor', \
                    u'dynpgrad4',u'nomomflxpartdynpgrad4'] #
	#nameexp_ = [u'$H=$arg$(k_{min})$',u'$H=$arg$(k<10^{-6}$~м$^2$/с$^2$)'] #
	
	#Choose whether to add experiment name as a suffix to legend entries or not
	addsuff = True #False

	# Temperatures - bounds of thermocline
	thermt1 = 14.
	thermt2 = 8.
	
	#Depths of model output
	#z = [0.2,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,6,7,8,10,12] #Kuivajarvi
	z = [0.2,0.5,1.,1.5,1.75,2.,2.25,2.5,2.75,3.,3.25,3.5,4.] #Valkea-Kotinen	

	nlevs_water = 16

        #The time of the first data entry in time series
	ny = 2014 #year
	mdays = 365/12
        nmonth = 9
        nday = 2
        nhour = 0
        nmin = 0
        timezone = None

        # Choose for date-time format of labels on X axis
        time_label_date = True

        time_unit = 1./24. #Time units in input data, expressed in days
        #if nplot >= 1 and nplot <= nmod_max:
        time0_date = dtm.datetime(ny,nmonth,nday,nhour,nmin,0,0,timezone)
        time0_ = dts.date2num(time0_date)
        years  = dts.YearLocator()   # every year
        months = dts.MonthLocator()  # every month
        days   = dts.DayLocator()  # every day
        dateFmt = dts.DateFormatter('%y-%m-%d')

	#Time of begin and end of time series for analysis
	#t1_aver = 4.
	#t2_aver = 9.9
	t1_aver = 0.#month_doy(5,ny)/mdays #+ 17./31.
	t2_aver = 24.*30 #month_doy(11,ny)/mdays #+ 21./31.
	
	# User controls
	if int(nplot) == tempmod1_: # Temperature of surface and soil
		basename = 'T_fluxes  1  1'
		fileout = 'Temp_soilcols_ser'
		nheader = 22 #Number of lines in the header
		ncols = 22 #Total number of columns
		ncol = [8] #,12,13] #,22,23,24,25,26] # ,22,23,24,25] #The number of columns in the file to display
		legend = ['T\_surf\_mod'] #,'T\_icesurf\_mod', 'T\_snowsurf\_mod']#,'Ts1','Ts2','Ts3','Ts4','Ts5'] #,'T_soil_0-1.9m','T_soil_1.9-3.8m','T_soil_3.8-5.7m','T_soil_5.7-6m']
		labelx = 'Time, months'
		labely = 'Temperature, Celsius'
		ctime = 5
		timescale = mdays*24 # 365/12*24 - number of hours in a month
		time0 = month_doy(5,ny)/mdays #+ 17./31. #Initial time
		#time0 = month_doy(7,ny)/mdays + 17./31.
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	if int(nplot) == tempmod2_: # Temperature at different depths
		basename = 'water_temp  1  1'
		fileout  = 'temp_wat'
		nheader  = 6
		ncols	 = 37 #Total number of columns
		ncol	 = [7,9,11,13,15]#,17,19,21,23,25,27,29,31,33,35,37] #The number of columns in the file to display
		#legend	 = ['T_0.2m','T_0.5m','T_1m','T_1.5m','T_2m','T_2.5m', \
		#'T_3m','T_3.5m','T_4m','T_4.5m','T_5m','T_6m','T_7m','T_8m','T_10m','T_12m']
		legend   = ['T_0.2m','T_1.15m','T_2.1m','T_3.05m','T_4m']
		labelx	 = 'Time, months'
		labely	 = 'Depth, m'
		timescale = mdays # 365/12*24 - number of hours in a month
		ctime = 5
		#time0 = month_doy(5,ny)/mdays #Initial time
		time0 = month_doy(7,ny)/mdays + 17./31. #Initial time, Bolshoi Vilui 2015
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	if int(nplot) == tempmod3_: # Temperature at selected points
		basename = 'temploc  1  1' #'temploc  1  1'
		fileout  = 'temploc'
		nheader  = 0
		ncols	 = 3 #Total number of columns
		ncol	 = [2,3] #[2,3]
		legend	 = ['T$_{7m}$, model','T$_{11m}$, model'] #['point 1','point 2']
		labelx	 = 'Time, days'
		labely	 = 'Temperature, C$^{\circ}$'
		timescale = 24. # 365/12*24 - number of hours in a month
		ctime = 1
		time0 = 0. #month_doy(5,ny)/mdays #Initial time
		#time0 = month_doy(7,ny)/mdays + 17./31. #Initial time, Bolshoi Vilui 2015
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
		missingval = -999.
	if int(nplot) == senslat_: # Sensible and latent heat fluxes
		basename = 'T_fluxes  1  1'
		fileout = 'Senslat_ser'
		nheader = 22 #Number of lines in the header
		ncols = 22 #Total number of columns
		ncol = [14]#,15] #The number of columns in the file to display
		legend = ['Sensible heat flux']#,'Latent heat flux'] #,'Ts1','Ts2','Ts3','Ts4','Ts5'] #,'T_soil_0-1.9m','T_soil_1.9-3.8m','T_soil_3.8-5.7m','T_soil_5.7-6m']
		labelx = 'Time, months'
		labely = 'Heat flux, $W/m^2$'
		ctime = 5
		timescale = mdays*24 # 365/12*24 - number of hours in a month
		time0 = month_doy(7,ny)/mdays + 17./31. #Initial time
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
                time_unit = 1./24. #Time units in input data, expressed in days
	if int(nplot) == uspeed_: # x-speed component, m/s
		basename = 'u  1  1'
		fileout  = 'uspeed_wat'
		nheader  = 6
		ncols	 = 127 #Total number of columns
		ncol	 = [7]#,17,19,21,23,25,27,29,31,33,35,37] #The number of columns in the file to display
		#legend	 = ['T_0.2m','T_0.5m','T_1m','T_1.5m','T_2m','T_2.5m', \
		#'T_3m','T_3.5m','T_4m','T_4.5m','T_5m','T_6m','T_7m','T_8m','T_10m','T_12m']
		legend   = ['u_0.0m']
		labelx	 = 'Time, months'
		labely	 = 'Speed, m/s'
		timescale = mdays # 365/12*24 - number of hours in a month
		ctime = 5
		#time0 = month_doy(5,ny)/mdays #Initial time
		time0 = 0. #month_doy(7,ny)/mdays + 17./31. #Initial time, Bolshoi Vilui 2015
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	if int(nplot) == ustar_: # Friction velocity
		basename = 'T_fluxes  1  1'
		fileout = 'Ustar_ser'
		nheader = 22 #Number of lines in the header
		ncols = 22 #Total number of columns
		ncol = [18]#,15] #The number of columns in the file to display
		legend = ['Friction velocity']#,
		labelx = 'Time, months'
		labely = 'Friction velocity, $m/s$'
		ctime = 5
		timescale = mdays*24 # 365/12*24 - number of hours in a month
		time0 = month_doy(7,ny)/mdays + 17./31. #Initial time
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
                time_unit = 1./24. #Time units in input data, expressed in days
	if int(nplot) == windwork_: # Work by wind at the water surface
		basename = 'T_fluxes  1  1'
		fileout = 'Windwork_ser'
		nheader = 22 #Number of lines in the header
		ncols = 22 #Total number of columns
		ncol = [19]#,15] #The number of columns in the file to display
		legend = ['Wind work']#,
		labelx = 'Time, months'
		labely = 'Wind work, $W/m^2$'
		ctime = 5
		timescale = mdays*24 # 365/12*24 - number of hours in a month
		time0 = month_doy(7,ny)/mdays + 17./31. #Initial time
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
                time_unit = 1./24. #Time units in input data, expressed in days
	if int(nplot) == mldmod_: # Mixed-layer depth
		basename = 'conv_series  1  1'
		fileout = 'Mixed_layer_ser'
		nheader = 21 #Number of lines in the header
		ncols = 21 #Total number of columns
		ncol = [12] # ,20,21,22,23] #The number of columns in the file to display
		#ncol = [11,12] # ,20,21,22,23] #The number of columns in the file to display
		#legend = ['MLD from $N^2$ profile','MLD from TKE profile'] #
		legend = [u'Глубина перемешанного слоя по $k-\epsilon$ модели'] #
		#labelx = 'Time, months'
		#labelx = 'Time, hours'
		#labely = 'Depth, m'
                nrunmean = 7*48
		labelx = u'Время, мес'
		labely = u'Глубина, м'
		ctime = 5
		timescale = mdays*24. #mdays*24 # 365/12*24 - number of hours in a month
		time0 = 5. #month_doy(5,ny)/mdays #Initial time
		revert_yaxis = 1
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	if int(nplot) == thermthickmod_: # Thermocline depth
		basename = 'water_temp  1  1'
		fileout  = 'thermthick_mod'
		nheader  = 6
		ncols	 = 37 #Total number of columns
		ncol	 = [7,9,11,13,15,17,19,21,23,25,27,29,31,33,35,37] #The number of columns in the file to display
		legend	 = ['Thermocline thickness']
		labelx	 = 'Time, months'
		labely	 = 'Depth, m'
		timescale = mdays # 365/12*24 - number of hours in a month
		ctime = 5
		time0 = month_doy(5,ny)/mdays #Initial time
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	if int(nplot) == wumod_: # u_*/w_*
		basename = 'conv_series  1  1'
		fileout = 'w*_u*_ser'
		nheader = 18 #Number of lines in the header
		ncols = 18 #Total number of columns
		ncol = [15] # ,20,21,22,23] #The number of columns in the file to display
		legend = ['-w*/u*'] #
		labelx = 'Time, months'
		labely = 'Depth, m'
		ctime = 5
		timescale = mdays*24 # 365/12*24 - number of hours in a month
		time0 = month_doy(5,ny)/mdays #Initial time
		#revert_yaxis = 1
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	if int(nplot) == rimldmod_: # bulk Richardson number in mixed-layer
		basename = 'conv_series  1  1'
		fileout = 'Ri_bulk_ser'
		nheader = 18 #Number of lines in the header
		ncols = 18 #Total number of columns
		ncol = [16] # ,20,21,22,23] #The number of columns in the file to display
		legend = ['$Ri_{bulk}$'] #
		labelx = 'Time, months'
		labely = '$Ri_{bulk}$, n/d'
		ctime = 5
		timescale = mdays*24 # 365/12*24 - number of hours in a month
		time0 = month_doy(5,ny)/mdays #Initial time
		#revert_yaxis = 1
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	if int(nplot) == maxbvfmod_: # maximal Brunt-Vaisala frequency
		basename = 'conv_series  1  1'
		fileout = 'maxN_ser'
		nheader = 18 #Number of lines in the header
		ncols = 18 #Total number of columns
		ncol = [17] # ,20,21,22,23] #The number of columns in the file to display
		legend = ['$\max N^2$'] #
		labelx = 'Time, months'
		labely = '$\max N^2, s^{-1}$'
		ctime = 5
		timescale = mdays*24 # 365/12*24 - number of hours in a month
		time0 = month_doy(5,ny)/mdays #Initial time
		#revert_yaxis = 1
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	if int(nplot) == mincondmod_: # minimal thermal conductance
		basename = 'conv_series  1  1'
		fileout = 'mincond_ser'
		nheader = 18 #Number of lines in the header
		ncols = 18 #Total number of columns
		ncol = [18] # ,20,21,22,23] #The number of columns in the file to display
		legend = ['$\min \lambda_w$'] #
		labelx = 'Time, months'
		labely = '$\min \lambda_w$'
		ctime = 5
		timescale = mdays*24 # 365/12*24 - number of hours in a month
		time0 = month_doy(5,ny)/mdays #Initial time
		#revert_yaxis = 1
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	if int(nplot) == tkemldmod_: # TKE terms integrated in the mixed layer
		basename = 'TKE_integr_conv  1  12013050100'
		fileout = 'TKE_integr_terms_ser'
		nheader = 15 #Number of lines in the header
		ncols = 15 #Total number of columns
		ncol = [7,8,9,10,11,12] # The number of columns in the file to display
		legend = ['B+','B-','S','$\epsilon$','TKE turb trans','$\partial E /\partial t$'] #
		labelx = 'Time, months'
		labely = 'TKE_terms, averaged, $m^2/s^3$'
		ctime = 5
		timescale = mdays*24 # 365/12*24 - number of hours in a month
		time0 = month_doy(5,ny)/mdays #Initial time
		#revert_yaxis = 1
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	if int(nplot) == Wmod_: # Wedderburn number from model
		basename = 'TKE_integr_conv  1  12013050100'
		fileout = 'Wedderburn_number_ser'
		nheader = 15 #Number of lines in the header
		ncols = 15 #Total number of columns
		ncol = [13] # The number of columns in the file to display
		legend = ['Wedderburn number'] #
		labelx = 'Time, months'
		labely = 'W, n/d'
		ctime = 5
		timescale = mdays*24 # 365/12*24 - number of hours in a month
		time0 = month_doy(5,ny)/mdays + 1.#Initial time
		nrunmean = 48 #timesteps
		horizline = 0.5
		#revert_yaxis = 1
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	if int(nplot) == LNmod_: # Lake number from model
		basename = 'TKE_integr_conv  1  12013050100'
		fileout = 'Lake_number_ser'
		nheader = 15 #Number of lines in the header
		ncols = 15 #Total number of columns
		ncol = [14] # The number of columns in the file to display
		legend = ['Lake number'] #
		labelx = 'Time, months'
		labely = '$L_N$, n/d'
		ctime = 5
		timescale = mdays*24 # 365/12*24 - number of hours in a month
		time0 = month_doy(5,ny)/mdays #Initial time
		nrunmean = 48 #timesteps
		horizline = 1.
		#revert_yaxis = 1
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	if int(nplot) == Rmod_: # Internal Rossby radius from the model
		basename = 'TKE_integr_conv  1  12013050100'
		fileout = 'Rossby_radius_ser'
		nheader = 15 #Number of lines in the header
		ncols = 15 #Total number of columns
		ncol = [15] # The number of columns in the file to display
		legend = ['Internal Rossby radius'] #
		labelx = 'Time, months'
		labely = '$L_R$, m'
		ctime = 5
		timescale = mdays*24 # 365/12*24 - number of hours in a month
		time0 = month_doy(5,ny)/mdays + 1.#Initial time
		horizline = 2600. #lake's maximal size
		#revert_yaxis = 1
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	elif int(nplot) == ch4mod_: # Methane
		basename = 'methane_water  1  1'
		fileout = 'Methane_water_ser'
		nheader = 6 #Number of lines in the header
		ncols = 39 #Total number of columns
		ncol = [7]#,11,19,27,31,37] #The number of columns in the file to display
		legend = ['CH4_0.2m']#,'CH4_1m','CH4_3m','CH4_5m','CH4_7m','CH4_12m']
		labelx = 'Time, months'
		labely = 'Methane concentration, $\mu g/l$'
		ctime = 5
		timescale = mdays # 365/12*24 - number of hours in a month
		time0 = month_doy(5,ny)/mdays #Initial time
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	elif int(nplot) == o2mod_: # Oxygen
		basename = 'oxygen_water  1  1'
		fileout = 'Oxygen_water_ser'
		nheader = 6 #Number of lines in the header
		ncols = 39 #Total number of columns
		ncol = [7]#,11,19,27,31,37] #The number of columns in the file to display
		legend = ['O2_0.2m']#,'O2_1m','O2_3m','O2_5m','O2_7m','O2_12m']
		labelx = 'Time, months'
		labely = 'Oxygen concentration, mg/l'
		ctime = 5
		timescale = mdays # 365/12*24 - number of hours in a month
		time0 = month_doy(5,ny)/mdays #Initial time
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	elif int(nplot) == co2mod_: # Carbon dioxide
		basename = 'co2_water  1  1'
		fileout = 'co2_water_ser'
		nheader = 6 #Number of lines in the header
		ncols = 39 #Total number of columns
		ncol = [7]#,11,19,27,31,37] #The number of columns in the file to display
		legend = ['CO2_0.2m']#,'CO2_1m','CO2_3m','CO2_5m','CO2_7m','CO2_12m']
		labelx = 'Time, months'
		labely = 'CO2 concentration, mg/l'
		ctime = 5
		timescale = mdays # 365/12*24 - number of hours in a month
		time0 = month_doy(5,ny)/mdays #Initial time
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	elif int(nplot) == hshnmod_: # water surface level at souther/northern parts
		basename = 'layers  1  1'
		fileout = 'hE-hW_ser'
		nheader = 19 #Number of lines in the header
		ncols = 19 #Total number of columns
		ncol = [7,8,11,12] #The number of columns in the file to display
		legend = ['ML W thickness','ML E thickness','Lower layer W thickness','Lower layer E thickness']
		labelx = 'Time, minutes'
		labely = 'Thickness, m'
		ctime = 5
		timescale = 1./60 # 365/12*24 - number of hours in a month
		time0 = 0 #Initial time
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	elif int(nplot) == hlsnmod_: # thickness of water, ice, snow and deep ice layers
		basename = 'layers  1  1'
		fileout = 'layers_ser'
		nheader = 19 #Number of lines in the header
		ncols = 19 #Total number of columns
		ncol = [6,15,16,17] #The number of columns in the file to display
		legend = ['water thickness','ice thickness','snow thickness','deep ice thickness']
		labelx = 'Time, minutes'
		labely = 'Thickness, m'
		ctime = 5
		timescale = 1./60 # 365/12*24 - number of hours in a month
		time0 = 0 #Initial time
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	elif int(nplot) == ebulmod_: # Ebullition at the surface
		basename = 'methane_series  1  1'
		fileout = 'ch4_ebul_soilcols_ser'
		nheader = 32 #Number of lines in the header
		ncols = 35 #Total number of columns
		ncol = [31,32,33,34,35] #The number of columns in the file to display
		legend = ['Fb_1','Fb_2','Fb_3','Fb_4','Fb_5']
		labelx = 'Time, months'
		labely = 'Surface bubble flux, $mg/(m^2day)$'
		ctime = 5
		timescale = mdays*24 # 365/12*24 - number of hours in a month
		time0 = month_doy(5,ny)/mdays #Initial time
		#valscale = 1e+6 # converting from mol/(m^2*s) to mcmol/(m^2*s)
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	elif int(nplot) == ch4flxmod_: # Methane fluxes
		basename = 'methane_series  1  1'
		fileout = 'methane_flux_ser'
		nheader = 32 #Number of lines in the header
		ncols = 35 #Total number of columns
		ncol = [14,16,17] #The number of columns in the file to display
		legend = ['Methane ebullition flux at the surface', \
		'Methane diffusion flux at the lake bottom','Methane diffusion flux at the lake surface']
		labelx = 'Time, months'
		labely = 'Methane flux, $\mu mol/(m^2*s)$'
		ctime = 5
		timescale = mdays*24 # 365/12*24 - number of hours in a month
		time0 = month_doy(5,ny)/mdays #Initial time
		valscale = [1e+6] # converting from mol/(m^2*s) to mcmol/(m^2*s)
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	elif int(nplot) == ch4flxMLmod_: # Methane fluxes
		basename = 'methane_series  1  1'
		fileout = 'methane_ML_budget'
		nheader = 32 #Number of lines in the header
		ncols = 35 #Total number of columns
		ncol = [21,22,23,24,18,26] #The number of columns in the file to display
		legend = ['Diffusion at the lake surface ', \
	                  'Diffusion at the bottom of mixed layer ', \
	                  'Diffusion+ebullition from mixed-layer sediments ', \
	                  'Ebullition at the bottom of mixed layer ', \
	                  'Ebullition at the lake surface ', \
	                  'Oxidation in the mixed layer']
		labelx = 'Time, months'
		labely = '$CH_4$ flux, norm. by surf. area, $mg/(m^2*day)$'
		ctime = 5
		timescale = mdays*24 # 365/12*24 - number of hours in a month
		time0 = month_doy(5,ny)/mdays #Initial time
		#valscale = 1e+6 # converting from mol/(m^2*s) to mcmol/(m^2*s)
		valscale = [-1,1,1,1,-1,1]
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	elif int(nplot) == ch4prodmod_: # Methane production in sediments
		basename = 'rprod_soil	1  1'
		fileout = 'methane_prod_new_ser'
		nheader = 6 #Number of lines in the header
		ncols = 15 #Total number of columns
		ncol = [7,9,11,13,15] #The number of columns in the file to display
		legend = ['Methane production in soil col 1', \
			  'Methane production in soil col 2', \
			  'Methane production in soil col 3', \
			  'Methane production in soil col 4', \
			  'Methane production in soil col 5' ]
		labelx = 'Time, months'
		labely = 'Methane production rate, $mol/(m^2s)$'
		ctime = 5
		timescale = mdays # 365/12*24 - number of hours in a month
		time0 = month_doy(5,ny)/mdays #Initial time
		path = path_
		nameexp = nameexp_
		npaths = len(path_)
	#if int(nplot) == tempobs_: # Temperature measured
        #        #Kuivajarvi
	#	basename = 'KJdata2013'  #'obs'
	#	fileout = 'Temp_obs'
	#	nheader = 0 #Number of lines in the header
	#	ncols = 21 #Total number of columns
	#	ncol = [2]#,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17] # ,20,21,22,23] #The number of columns in the file to display
	#	legend = ['T_0.2m_obs','T_0.5m','T_1m','T_1.5m','T_2m','T_2.5m', \
	#	 'T_3m','T_3.5m','T_4m','T_4.5m','T_5m','T_6m','T_7m','T_8m', \
	#	 'T_10m','T_12m'] 
	#	labelx = 'Time, months'
	#	labely = 'Temperature, Celsius'
	#	ctime = 1
	#	timescale = mdays # 365/12*24 - number of hours in a month
	#	time0 = 0 #Initial time
	#	path = pathobs
	#	nameexp = nameexpobs
	#	npaths = 1
	#if int(nplot) == tempobs_: # Temperature measured
        #        #BolshoiVilui2015
	#	basename = 'all'  #'obs'
	#	fileout = 'Temp_obs'
	#	nheader = 0 #Number of lines in the header
	#	ncols = 9 #Total number of columns
	#	ncol = [5,6,7,8,9] #,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17] # ,20,21,22,23] #The number of columns in the file to display
	#	legend = ['T_0.2m','T_1.15m','T_2.1m','T_3.05m','T_4m'] #,'T_2.5m', \
	#	# 'T_3m','T_3.5m','T_4m','T_4.5m','T_5m','T_6m','T_7m','T_8m', \
	#	# 'T_10m','T_12m'] 
	#	labelx = 'Time, months'
	#	labely = 'Temperature, Celsius'
	#	ctime = 4
	#	timescale = mdays # 365/12*24 - number of hours in a month
	#	time0 = month_doy(7,ny)/mdays + 17./31. #Initial time
	#	path = pathobs
	#	nameexp = nameexpobs
	#	npaths = 1
	#	delim = ','
	#if int(nplot) == tempobs_: # Temperature measured
        #        #IseoAugust2011
	#	basename = 'TC' #'LDS'  
	#	fileout = 'Temploc_obs'
	#	nheader = 0 #Number of lines in the header
	#	ncols =  7 + 16*2 # LDS: 7 + 21*2 #Total number of columns
        #        ncol = [29] #LDS_14.8m:[37]
	#	legend = ['T\_15.0m'] #['T\_14.80m']
	#	#ncol = nm.arange(9, ncols, 2) #,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17] # ,20,21,22,23] #The number of columns in the file to display
        #        #legend = \
        #        #['T_0.30m',  'T_0.55m',  'T_0.80m',  'T_1.05m',  'T_1.30m',  'T_2.30m',  'T_3.30m', \
        #        # 'T_4.30m',  'T_5.80m',  'T_7.30m',  'T_8.80m',  'T_10.30m', 'T_11.80m', 'T_13.30m', \
        #        # 'T_14.80m', 'T_16.30m', 'T_17.80m', 'T_19.80m', 'T_21.80m', 'T_35.80m', 'T_49.80m'] #LDS
	#	labelx = 'Time, days'
	#	labely = 'Temperature, Celsius'
	#	ctime = 7
	#	timescale = 24. # 365/12*24 - number of hours in a month
	#	time0 = 0. #month_doy(7,ny)/mdays + 17./31. #Initial time
	#	path = pathobs
	#	nameexp = nameexpobs
	#	npaths = 1
	#	delim = ','
	#if int(nplot) == tempobs_: # Temperature measured
        #        #Mojai2016
	#	basename = 'Tsurf_obs_Gella_2016' #'LDS'  
	#	fileout = 'Temploc_obs'
	#	nheader = 0 #Number of lines in the header
	#	ncols =  3 # LDS: 7 + 21*2 #Total number of columns
        #        ncol = [3] #LDS_14.8m:[37]
	#	legend = ['T\_surf'] #['T\_14.80m']
	#	#ncol = nm.arange(9, ncols, 2) #,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17] # ,20,21,22,23] #The number of columns in the file to display
        #        #legend = \
        #        #['T_0.30m',  'T_0.55m',  'T_0.80m',  'T_1.05m',  'T_1.30m',  'T_2.30m',  'T_3.30m', \
        #        # 'T_4.30m',  'T_5.80m',  'T_7.30m',  'T_8.80m',  'T_10.30m', 'T_11.80m', 'T_13.30m', \
        #        # 'T_14.80m', 'T_16.30m', 'T_17.80m', 'T_19.80m', 'T_21.80m', 'T_35.80m', 'T_49.80m'] #LDS
	#	labelx = 'Time, days'
	#	labely = 'Temperature, Celsius'
	#	ctime = 2
	#	timescale = 1. # 365/12*24 - number of hours in a month
	#	time0 = 0. #month_doy(7,ny)/mdays + 17./31. #Initial time
	#	path = pathobs
	#	nameexp = nameexpobs
	#	npaths = 1
	#	delim = ','
	if int(nplot) == tempobs_: # Temperature measured
                #Kuivajarvi2014 campaign September
		basename = 'KuivaMETEO_201409' #'LDS'  
		fileout = 'Temploc_obs'
		nheader = 0 #Number of lines in the header
		ncols = 34 # LDS: 7 + 21*2 #Total number of columns
                ncol = [15,17] #LDS_14.8m:[37]
		legend = ['T$_{7m}$','T$_{10m}$'] #['T\_14.80m']
		#ncol = nm.arange(9, ncols, 2) #,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17] # ,20,21,22,23] #The number of columns in the file to display
                #legend = \
                #['T_0.30m',  'T_0.55m',  'T_0.80m',  'T_1.05m',  'T_1.30m',  'T_2.30m',  'T_3.30m', \
                # 'T_4.30m',  'T_5.80m',  'T_7.30m',  'T_8.80m',  'T_10.30m', 'T_11.80m', 'T_13.30m', \
                # 'T_14.80m', 'T_16.30m', 'T_17.80m', 'T_19.80m', 'T_21.80m', 'T_35.80m', 'T_49.80m'] #LDS
		labelx = 'Time, days'
		labely = 'Temperature, Celsius'
		ctime = 2
		timescale = 1. # 365/12*24 - number of hours in a month
		time0 = 0. #month_doy(7,ny)/mdays + 17./31. #Initial time
		path = pathobs
		nameexp = nameexpobs
		npaths = 1
		#delim = ','
		missingval = -999.
	if int(nplot) == mldobs_: # Mixed-layer depth observed
		basename = 'Valkea-Kotinen_validation_May-Dec_2006_3_mod' #'KJdata2013'  #'obs'
		fileout = 'mld_obs'
		nheader = 0 #Number of lines in the header
		ncols = 24 #Total number of columns
		ncol = [12,13,14,15,16,17,18,19,20,21,22,23,24] # ,20,21,22,23] #The number of columns in the file to display
		legend = [u'Глубина перемешанного слоя']
		labelx = u'Время, мес' #'Time, months'
		labely = u'Глубина, м' #'Mixed-layer depth, m'
		ctime = 1   
		timescale = mdays # 365/12*24 - number of hours in a month
		time0 = 5. #Initial time
		path = pathobs
		nameexp = nameexpobs
		npaths = 1  
		#revert_yaxis = 1
		missingval = -999.
	if int(nplot) == thermthickobs_: # Thermocline thickness observed
		basename = 'KJdata2013'  #'obs'
		fileout = 'thermthick_obs'
		nheader = 0 #Number of lines in the header
		ncols = 21 #Total number of columns
		ncol = [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17] # ,20,21,22,23] #The number of columns in the file to display
		legend = ['Thermocline thickness']
		labelx = 'Time, months'
		labely = 'Thermocline thickness, m'
		ctime = 1
		timescale = mdays # 365/12*24 - number of hours in a month
		time0 = 0 #Initial time
		path = pathobs
		nameexp = nameexpobs
		npaths = 1
		missingval = -999.
		#revert_yaxis = 1
	if int(nplot) == maxbvfobs_: # Maximal Brunt-Vaisala frequency observed
		basename = 'KJdata2013'  #'obs'
		fileout = 'maxN_obs'
		nheader = 0 #Number of lines in the header
		ncols = 21 #Total number of columns
		ncol = [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17] # ,20,21,22,23] #The number of columns in the file to display
		legend = ['Maximal Brunt-Vaisala frequency']
		labelx = 'Time, months'
		labely = 'Maximal Brunt-Vaisala frequency, $s^{-1}$'
		ctime = 1
		timescale = mdays # 365/12*24 - number of hours in a month
		time0 = 0 #Initial time
		path = pathobs
		nameexp = nameexpobs
		npaths = 1
	#if int(nplot) == co2obs1_: # Carbon dioxide measured
	#	basename = 'KJdata2013'
	#	fileout = 'co2_obs'
	#	nheader = 1 #Number of lines in the header
	#	ncols = 21 #Total number of columns
	#	ncol = [18]#,19,20,21]  #The number of columns in the file to display
	#	legend = ['CO2_0.2m']#,'CO2_1.5m','CO2_2.5m','CO2_7m'] 
	#	labelx = 'Time, months'
	#	labely = 'Carbon dioxide, ppm'
	#	ctime = 1
	#	timescale = mdays # 365/12*24 - number of hours in a month
	#	time0 = month_doy(1,ny)/mdays #Initial time
	#	path = pathobs
	#	nameexp = nameexpobs
	#	npaths = 1
	if int(nplot) == co2obs1_: # Carbon dioxide measured
		# Kuivajarvi 2014 campaign september
		basename = 'KuivaMETEO_201409' #'LDS'  
		fileout = 'co2_obs'
		nheader = 0 #Number of lines in the header
		ncols = 34 # LDS: 7 + 21*2 #Total number of columns
                ncol = [30] #LDS_14.8m:[37]
		legend = ['CO$_2$, 7 m'] #['T\_14.80m']
		#ncol = nm.arange(9, ncols, 2) #,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17] # ,20,21,22,23] #The number of columns in the file to display
                #legend = \
                #['T_0.30m',  'T_0.55m',  'T_0.80m',  'T_1.05m',  'T_1.30m',  'T_2.30m',  'T_3.30m', \
                # 'T_4.30m',  'T_5.80m',  'T_7.30m',  'T_8.80m',  'T_10.30m', 'T_11.80m', 'T_13.30m', \
                # 'T_14.80m', 'T_16.30m', 'T_17.80m', 'T_19.80m', 'T_21.80m', 'T_35.80m', 'T_49.80m'] #LDS
		labelx = 'Time, days'
		labely = 'Concentration, ppm'
		ctime = 2
		timescale = 1. # 365/12*24 - number of hours in a month
		time0 = 0. #month_doy(7,ny)/mdays + 17./31. #Initial time
		path = pathobs
		nameexp = nameexpobs
		npaths = 1
		#delim = ','
		missingval = -999.
	if int(nplot) == co2obs2_: # CO_2 measured manually
		basename = 'ch4&co2man2013'
		fileout = 'co2_obs_man'
		nheader = 1 #Number of lines in the header
		ncols = 28 #Total number of columns
		ncol = [6]#,9,12,15,18,21,24,27]	#The number of columns in the file to display
		legend = ['CO2_0.2m']#,'CO2_1m','CO2_3m','CO2_5m','CO2_7m', \
			  #'CO2_9m','CO2_11m','CO2_12m'] 
		labelx = 'Time, months'
		labely = '$CO_2$ concentration, $mg/l$'
		ctime = 3
		timescale = mdays # 365/12*24 - number of hours in a month
		time0 = month_doy(1,ny)/mdays #Initial time
		path = pathobs
		nameexp = nameexpobs
		npaths = 1
	if int(nplot) == ch4obs_: # Methane measured manually
		basename = 'ch4&co2man2013'
		fileout = 'ch4_obs_man'
		nheader = 1 #Number of lines in the header
		ncols = 28 #Total number of columns
		ncol = [5]#,8,11,14,17,20,23,26]	#The number of columns in the file to display
		legend = ['CH4_0.2m']#,'CH4_1m','CH4_3m','CH4_5m','CH4_7m', \
			  #'CH4_9m','CH4_11m','CH4_12m'] 
		labelx = 'Time, months'
		labely = 'Methane concentration, $\mu g/l$'
		ctime = 3
		timescale = mdays # 365/12*24 - number of hours in a month
		time0 = month_doy(1,ny)/mdays #Initial time
		path = pathobs
		nameexp = nameexpobs
		npaths = 1
	if int(nplot) == o2obs_: # Oxygen measured manually
		basename = 'o2man2013'
		fileout = 'o2_obs'
		nheader = 1 #Number of lines in the header
		ncols = 24 #Total number of columns
		ncol = [2]#,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24]  #The number of columns in the file to display
		legend = ['O2_0m','O2_0.5m','O2_1m','O2_1.5m','O2_2m', \
			  'O2_2.5m','O2_3m','O2_3.5m','O2_4m','O2_4.5m','O2_5m', \
			  'O2_5.5m','O2_6m','O2_6.5m','O2_7m','O2_7.5m','O2_8m', \
			  'O2_8.5m','O2_9m','O2_10m','O2_11m','O2_12m','O2_12.5m'] 
		labelx = 'Time, months'
		labely = 'Oxygen concentration, $mg/l$'
		ctime = 1
		timescale = mdays # 365/12*24 - number of hours in a month
		time0 = month_doy(1,ny)/mdays #Initial time
		path = pathobs
		nameexp = nameexpobs
		npaths = 1
	if int(nplot) == ch4flxobs_: # Methane flux measured
		basename = 's13_Kuic_all_trunc'
		fileout =  'ch4flux_obs'
		nheader = 1 #Number of lines in the header
		ncols = 168 #Total number of columns
		ncol = [19]  #The number of columns in the file to display
		legend = ['CH4_flux_surf'] 
		labelx = 'Time, months'
		labely = 'Surface methane flux, $\mu mol/(m^2*s)$'
		ctime = 1
		timescale = mdays # 365/12*24 - number of hours in a month
		time0 = month_doy(1,ny)/mdays #Initial time
		missingval = -999.
		path = pathobs
		nameexp = nameexpobs
		npaths = 1
	if int(nplot) == co2flxobs_: # Carbon dioxide flux measured
		basename = 's13_Kuic_all_trunc'
		fileout =  'co2flux_obs'
		nheader = 1 #Number of lines in the header
		ncols = 168 #Total number of columns
		ncol = [17]  #The number of columns in the file to display
		legend = ['CO2_flux_surf'] 
		labelx = 'Time, months'
		labely = 'Surface $CO_2$ flux, $mol/(m^2*s)$'
		ctime = 1
		timescale = mdays # 365/12*24 - number of hours in a month
		time0 = month_doy(1,ny)/mdays #Initial time
		missingval = -999.
		path = pathobs
		nameexp = nameexpobs
		npaths = 1
	if int(nplot) == shflxobs_: # Sensible heat flux measured
		basename = 'lake_fluxes4LAKE'
		fileout =  'shflux_obs'
		nheader = 1 #Number of lines in the header
		ncols = 2 #Total number of columns
		ncol = [2]  #The number of columns in the file to display
		legend = ['Sensible_heat_flux'] 
		labelx = 'Time, months'
		labely = 'Sensible heat flux, $W/m^2$'
		ctime = 1
		timescale = mdays # 365/12*24 - number of hours in a month
		#time0 = month_doy(1,ny)/mdays #Initial time
		time0 = month_doy(7,ny)/mdays + 17./31. #Initial time, Bolshoi Vilui 2015
		missingval = -999.
		path = pathobs
		nameexp = nameexpobs
		npaths = 1
                delim = ','
	if int(nplot) == ustarobs_: # friction velocity measured
		basename = 'ustar'
		fileout =  'ustar_obs'
		nheader = 1 #Number of lines in the header
		ncols = 4 #Total number of columns
		ncol = [4]  #The number of columns in the file to display
		legend = ['Friction velocity'] 
		labelx = 'Time, months'
		labely = 'Friction velocity, $m/s$'
		ctime = 2
		timescale = mdays # 365/12*24 - number of hours in a month
		#time0 = month_doy(1,ny)/mdays #Initial time
		time0 = month_doy(7,ny)/mdays + 17./31. #Initial time, Bolshoi Vilui 2015
		missingval = -999.
		path = pathobs
		nameexp = nameexpobs
		npaths = 1
                delim = ','
	if int(nplot) == rimod_: # Richardson number
		basename = 'Ri	1  1'
		fileout = 'Ri_timser'
		nheader = 6 #Number of lines in the header
		ncols = 39
		ncol = [11,19,27,31,37]
		legend = ['Ri_1m','Ri_3m','Ri_5m','Ri_7m','Ri_12m']
		labelx = 'Time, months'
		labely = 'Ri'
		timescale = mdays # 365/12*24 - number of hours in a month
		ctime = 5
		time0 = month_doy(5,ny)/mdays #Initial time
		path = path_
		nameexp = nameexp_
		npaths = len(path)
	if int(nplot) == condmod_: # Heat conductance
		basename = 'lamw  1  1'
		fileout = 'lamw_timser'
		nheader = 6 #Number of lines in the header
		ncols = 39
		ncol = [11,19,27,31,37]
		legend = ['HC_1m','HC_3m','HC_5m','HC_7m','HC_12m']
		labelx = 'Time, months'
		labely = 'Heat conductivity, $m^2/s$'
		timescale = mdays # 365/12*24 - number of hours in a month
		ctime = 5
		time0 = month_doy(5,ny)/mdays #Initial time
		path = path_
		nameexp = nameexp_
		npaths = len(path)
	if int(nplot) == icethickmod_: # Ice thickness
		basename = u'layers  1  1'
		fileout = u'icethick_timser'
		nheader = 19 #Number of lines in the header
		ncols = 19
		ncol = [15]
		legend = [u'Толщина слоя льда, м']
		labelx = u'Время, мес' #'Time, months'
		labely = u'Толщина, м' #'Heat conductance, $m^2/s$'
		timescale = 365./12.*24 #number of hours in a month
		ctime = 5
		time0 = 0. #month_doy(5,ny)/mdays #Initial time
		path = path_
		nameexp = nameexp_
		npaths = len(path)
	if int(nplot) == snowthickmod_: # Ice thickness
		basename = u'layers  1  1'
		fileout = u'icethick_timser'
		nheader = 19 #Number of lines in the header
		ncols = 19
		ncol = [16]
		legend = [u'Толщина слоя снега, м']
		labelx = u'Время, мес' #'Time, months'
		labely = u'Толщина, м' #'Heat conductance, $m^2/s$'
		timescale = 365./12.*24 #number of hours in a month
		ctime = 5
		time0 = 0. #month_doy(5,ny)/mdays #Initial time
		path = path_
		nameexp = nameexp_
		npaths = len(path)

	nleg = len(legend)
	legendbase = []*nleg
	legendbase[:] = legend[:]
	
	colorline = ['b','g','r','c','m','y','k']
	typeline = ['-','--','-.',':','.',',','o','v','^','1']
	ncolors = len(colorline)
	nlinetypes = len(typeline)
	
	# Changing default matplotlib settings
	plt.rc('xtick',labelsize=20)
	plt.rc('ytick',labelsize=20)
	plt.rc('lines',linewidth=2)
	
	ncol = nm.array(ncol,dtype='int32')
	ncol[:] = ncol[:] - 1
	
	#Derived variables
	if int(nplot) == thermthickmod_ or \
	   int(nplot) == thermthickobs_ or \
	   int(nplot) == mldobs_ or \
	   int(nplot) == maxbvfobs_:
		nvardisp = 1
	else:#ready variables from files
		nvardisp = len(ncol)
	
	#nfiles = 10 #number of files
	#time = [0]*nfiles
	##Files names suffixes (timestamps)
	#time[0] = '2003010100'
	#time[1] = '2003010200'
	#time[2] = '2003010300'
	#time[3] = '2003010400'
	#time[4] = '2003010500'
	#time[5] = '2003010600'
	#time[6] = '2003010700'
	#time[7] = '2003010800'
	#time[8] = '2003010900'
	#time[9] = '2003011000'
	
	#for i in range(nfiles):
	for ip in range(npaths):
		fid = open(path[ip] + basename + u'.dat','r')
		readheader(fid,nheader)
	
		#Counting the number of rows in a file (assuming the same in all files)
		array = fid.read()
		nrows = array.count('\n')
		print('The number of values is '+str(nrows))
		time = nm.zeros(nrows)
		vals = nm.zeros((nrows,len(ncol))) #[[0]*nvardisp]*nrows
		fid.seek(0) #rewinding file
		readheader(fid,nheader)
	
		#Reading array of depths and values to be plotted
		for j in range(nrows) :
			line = fid.readline()
			if 'delim' in locals():
			        line = line.split(delim)
			else:
                                line = line.split()
			print(line)
			time[j] = float(line[ctime-1])
			for i in range(len(ncol)):
				vals[j,i] = float(line[ncol[i]])
	
		time[:] = time[:]/timescale + time0
                if time_label_date: time2 = time[:]*time_unit + time0_
		for j in range(nvardisp) : 
			vard = [0]*nrows
	
			#Thermocline thickness
			if   int(nplot) == thermthickmod_ or int(nplot) == thermthickobs_: #thermocline thickness
				for i in range(nrows):
					if 'missingval' in locals():
						if any(vals[i,:]) == missingval:
							vard[i] = float('NaN')
						else:
							vard[i] = thermthick(vals[i,:],z,thermt1,thermt2)
					else:
						vard[i] = thermthick(vals[i,:],z,thermt1,thermt2)
	
			#Observed mixed-layer depth calculation
			elif int(nplot) == mldobs_:#observed mixed-layer depth calculation
				for i in range(nrows):
					if 'missingval' in locals():
						if any(vals[i,:]) == missingval:
							vard[i] = float('NaN')
						else:
							vard[i] = mldepth(vals[i,:],z)
					else:
						vard[i] = mldepth(vals[i,:],z)
	
			#Observed maximum buoyancy frequency
			elif int(nplot) == maxbvfobs_:#observed maximum B-V frequency
				for i in range(nrows):
					if 'missingval' in locals():
						if any(vals[i,:]) == missingval:
							vard[i] = float('NaN')
						else:
							vard[i] = maxbvf(vals[i,:],z)
					else:
						vard[i] = maxbvf(vals[i,:],z)
	
			else:
				for i in range(nrows):
					vard[i] = vals[i,j]
                                        #Checking for non-numbers
					if is_number(vard[i]) != True:
                                                print('Non-float: ' + str(i) + ', ' + str(vard[i]))
			# Assigning missing values to NaN
			if 'missingval' in locals():
				for i in range(nrows):
					if vard[i] == missingval:
						vard[i] = float('NaN')
			jcolor = ip+nvar-( (ip+nvar) //ncolors)*ncolors # periodic use of colors in one experiment
			#jcolor = j-( (j) //ncolors)*ncolors # periodic use of colors in one experiment
			jline = j+nvar-( (j+nvar) //nlinetypes)*nlinetypes #ip+nvar
			if 'valscale' in locals():
				if len(valscale) == 1:
					for k in range(len(vard)):
						vard[k] = vard[k]*valscale[0]
				else:
					for k in range(len(vard)):
						vard[k] = vard[k]*valscale[j]
			if 'nrunmean' in locals():
				vard = runningMeanFast(vard,nrunmean)
			if plotid != None:
				# Adding to existent plot new curves
                                if time_label_date:
                                # Constructing time labels in date-time format
                                        time_ = []
                                        for i in range(len(time)):
                                                time_.append(dts.num2date(time2[i]))
					print('time',time_)
					print('vard',vard)
					print('colorline',colorline[jcolor])
					print('typeline',colorline[jcolor])
					print('figure',plotid)
                                        plt.plot(time_,vard,colorline[jcolor]+typeline[jline]) #,figure=plotid) 
                                        #plt.gca().xaxis.set_major_locator(days)
                                        #plt.gca().xaxis.set_major_formatter(dateFmt)
                                        #fig.xaxis.set_minor_locator(months)
                                        #fig.autofmt_xdate()
                                else:
				        plt.plot(time,vard,colorline[jcolor]+typeline[jline])#,figure=plotid)
			else:
				# Creating plot 
                                if time_label_date:
                                # Constructing time labels in date-time format
                                        time_ = []
                                        for i in range(len(time)):
                                                time_.append(dts.num2date(time2[i]))
                                        fig, ax = plt.subplots(1)
                                        plotid = ax.plot(time_,vard,colorline[jcolor]+typeline[ip+nvar]) 
                                        plt.gca().xaxis.set_major_locator(days)
                                        plt.gca().xaxis.set_major_formatter(dateFmt)
                                        #fig.xaxis.set_minor_locator(months)
                                        fig.autofmt_xdate()
                                else:
				        plotid = plt.plot(time,vard,colorline[jcolor]+typeline[ip+nvar])
				if 'horizline' in locals():
					work = nm.ones(len(time))*horizline
					plt.plot(time,work,colorline[jcolor+1]+typeline[ip+nvar+1],figure=plotid)
			# Beginning and end of times series for analysis
			i0 = timeoffset(time,t1_aver) 
			i1 = timeoffset(time,t2_aver) 
			timan = time[i0:i1]
			varan = vard[i0:i1]
			#Calculating mean
			varanm = nm.ma.masked_array(data = varan, mask = nm.isnan(varan)) #removing NaNs
			#print(varan,i0,i1,t1_aver,t2_aver,time[0],time[len(time)-1])
			print('Time series mean ',nm.mean(varanm))
		
		# Adding to legend new entries from other experiments
		if addsuff: 
			legsuff = ',~'+nameexp[ip]
		else:
			legsuff = ''
		if ip > 0:
			for i in range(nleg):
				legend.append(legendbase[i]+legsuff)
		else:
			for i in range(nleg):
				legend[i] = legend[i]+legsuff
	
	
	#Adding current legend to common legend
	legendcommon += legend

	plt.xlabel(labelx,fontsize=26)
	plt.ylabel(labely,fontsize=26)
	plt.legend(legendcommon,fontsize=16,loc='best')
	
	# Revert y axis if needed 
	# (reverts every time, when revert_yaxis is set; be careful when multiple variables
	# are plotted)
	if 'revert_yaxis' in locals():
		plt.gca().invert_yaxis()
	
	#Postfix in filename
	postfix = nameexp[0]
	for i in range(npaths-1):
		postfix = postfix + '&' + nameexp[i+1]
	
	#Saving figure
	#for i in range(npaths):
	#	plt.savefig(path[i] + fileout + '\_' + postfix + '.png')
	
	#plt.draw()

	return (plotid,timan,varan,legendcommon)
