fun.MEFavg011 <- function(PLOT=T)
	{

## Function designed to go through raw 1-s trace gas data from Manitou Forest and :
## (1) make 1 min. averages (leaving out the first 10-20 sec.)
######## 
#### For Manitou Exp forest data - 2009

### Input text file with names of input/output files in it.
### Will need to update "c:/data/Manitou09/Datafiles09.txt" to analyze different months.
### will create monthly files (intially 1 min. averages)
#datafile <- scan("c:/data/Manitou/Manitou09/Datafiles09.txt",what="character")
datafile <- scan("c:/data/Manitou/Manitou011/Datafiles011.txt",what="character")

## Need list of raw data files:
infile <- datafile[36]		### list of filenames of 1 sec data.
datafile.1min <- datafile[52]	### 1 min. avg'd data to fill in with.
### Output data
outfile <- datafile[37]
outfile2 <- datafile[1]
########################
month <- datafile[30]
cat("\n Month:",month)
### read in list of filenames
datafile <- scan(infile,skip=1,sep=",",list(index=0,filenames="",starttime=""))

filenames <- datafile$filenames
Findex <- datafile$index
t.starts <- datafile$starttime
t.startA <- strptime(t.starts, "%m/%d/%Y %H:%M")
t.starts <- as.POSIXct(t.startA)

## read in the names of data files
#### when making up these files - get rid of last 10 columns!!!  EXCEPT Std(NOstate)!!!!!!!!  Keep that one!
browser()
data.1min <- matrix(scan(datafile.1min,skip=1,sep=",",list(datetime="",junk=0,Vso2=0,Vo3=0,Vno=0,Vco2=0,Vh2o=0,c1=0,c2=0,c3=0,c4=0,c5=0,Zero=0,Span=0,NOstate=0,Lev=0,Fstd=0,Ftot=0,sjunk=0,sVso2=0,sVo3=0,sVno=0,sVco2=0,sVh2o=0,sNOstate=0,sLevel=0)),byrow=T,ncol=26)
# parse the data
date1m <- as.ts(data.1min[[1]])
date1m <- strptime(date1m, "%m/%d/%Y %H:%M")
date1m <- fun.Ltimetag(date1m)
#### now have numeric dates - matrix of yy,jday,hrx100,min,sec
## get decimal day:
dday1m <- date1m[,2] + (date1m[,3]/2400) + (date1m[,4]/1440)
dday1m <- signif(dday1m,d=7)
date1m <- cbind(date1m,dday1m)		## add decimal day to date.
#####################################################################################
### other 5 min. data:
Vso2.1m <- as.ts(data.1min[[3]])
Vo3.1m <- as.ts(data.1min[[4]])
Vno.1m <- as.ts(data.1min[[5]])
Vco2.1m <- as.ts(data.1min[[6]])
Vh2o.1m <- as.ts(data.1min[[7]])
Zero.1m <- as.ts(data.1min[[13]])
Span.1m <- as.ts(data.1min[[14]])
NOstate.1m <- as.ts(data.1min[[15]])
Level.1m <- as.ts(data.1min[[16]])
Fstd.1m <- as.ts(data.1min[[17]])
Ftot.1m <- as.ts(data.1min[[18]])

sso2.1m <- as.ts(data.1min[[20]])
so3.1m <- as.ts(data.1min[[21]])
sno.1m <- as.ts(data.1min[[22]])
sco2.1m <- as.ts(data.1min[[23]])
sh2o.1m <- as.ts(data.1min[[24]])
sNOstate.1m <- as.ts(data.1min[[25]])
sLevel.1m <- as.ts(data.1min[[26]])
Flag.1m <- rep(2,length=length(sso2.1m))
NOFlag.1m <- ifelse(sNOstate.1m == 0,1,0)	
######### Put back together as matrix//same format as the final output matrix from teh 1 sec data.
DATA.1min <- cbind(date1m,Vso2.1m,Vo3.1m,Vno.1m,Vco2.1m,Vh2o.1m,Zero.1m,Span.1m,NOstate.1m,Level.1m,Fstd.1m,Ftot.1m,sso2.1m,so3.1m,sno.1m,sco2.1m,sh2o.1m,sNOstate.1m,sLevel.1m,NOFlag.1m,Flag.1m)
cat("\n Correct 1 min avg'd data?")
#browser()


i <- 1			### will denote file number
j <- 1			### j will lines in output file matrix.
Lfile <- length(Findex)
err1ct <- 0		### some error counting.
err2ct <- 0

## generate output matrix
Outlength <-  (60*Lfile) + 100			### 60 min/file, 100 extra rows for missing data.		
M <- matrix(NA,ncol=25,nrow=Outlength)	
	
### one separate for the date (since it is text)
#D <- matrix(NA,ncol=1,nrow=(60*Lfile))

#browser()

while(i <= Lfile) {

cat("\n Index:",Findex[i],"   Loop #:",i)
cat("\n\n")
filename.i <- as.numeric(filenames[i])
p <- is.na(filename.i)

#if(i == 1)	browser()
#if(i == 4)	browser()

if(p != 0)	{

### before 11/09
#d <- matrix(scan(filenames[i],sep="\t",fill=T,list(date="",msec=0,junk=0,so2=0,o3=0,no=0,co2=0,h2o=0,c1=0,c2=0,c3=0,c4=0,c5=0,Z=0,S=0,Sno=0,Level=0,Fst=0,Ftot=0,junk=0)),byrow=T,ncol=20)

### new program after 11/09
#d <- matrix(scan(filenames[i],sep="\t",fill=T,list(date="",msec=0,junk=0,so2=0,o3=0,no=0,co2=0,h2o=0,c1=0,c2=0,c3=0,c4=0,c5=0,Z=0,S=0,Sno=0,Level=0,Fst=0,Ftot=0,junk=0,junk2=0)),byrow=T,ncol=21)

### on June 22, 2010 - added another column to output files:
#if(month == "jun" & i >= 496)	{
#d <- matrix(scan(filenames[i],sep="\t",fill=T,list(date="",msec=0,junk=0,so2=0,o3=0,no=0,co2=0,h2o=0,c1=0,c2=0,c3=0,c4=0,c5=0,Z=0,S=0,Sno=0,Level=0,Fst=0,Ftot=0,junk=0,Pli=0,junk2=0)),byrow=T,ncol=22)
#	}

#if(i==496)	browser()
#if(i==497)	browser()
#if(i==498)	browser()
### After June 2010:
d <- matrix(scan(filenames[i],sep="\t",fill=T,list(date="", msec=0,junk=0,so2=0, o3=0,no=0,co2=0,h2o=0,c1=0,c2=0,c3=0,c4=0,c5=0,Z=0,S=0,Sno=0,Level=0,Fst=0,Ftot=0,junk=0,Pli=0,junk2=0)),byrow=T,ncol=22)

### xx is just a blank column
### parse the data and make into 10 Hz time series....
date <- as.ts(d[[1]])
so2 <- as.ts(d[[4]])
o3 <- as.ts(d[[5]])
no <- as.ts(d[[6]])
co2 <- as.ts(d[[7]])
h2o <- as.ts(d[[8]])
#####################
Zero <- as.ts(d[[14]])
Span <- as.ts(d[[15]])
NOstate <- as.ts(d[[16]])
Level <- as.ts(d[[17]])
Fstd <- as.ts(d[[18]])
Ftot <- as.ts(d[[19]])

###### Get time tags into something workable ###############
############################################################
date1 <- strptime(date, "%m/%d/%Y %H:%M:%S")
date2 <- fun.Ltimetag(date1)
#### now have numeric dates - matrix of yy,jday,hrx100,min,sec
## get decimal day:
dday <- date2[,2] + (date2[,3]/2400) + (date2[,4]/1440) + (date2[,5]/86400)
dday <- signif(dday,d=7)
date2 <- cbind(date2,dday)		## add decimal day to date.
## also create an index - sequential numbers to reference
Ldata <- length(date2[,1])
 
##############################################################
## also create an index - sequential numbers to reference
index <- seq(from=1,to=Ldata,by=1)
############################################

#### PUt into matrix for analysis
data <- cbind(so2,o3,no,co2,h2o,Zero,Span,NOstate,Level,Fstd,Ftot)

cat("\n File # = ",i)


#################################################
### now we want a flag for every new minute ####
Fmin <- diff(date2[,4])		
Fmin <-  c(1,Fmin)		## offset by 1, now have "1" each new min.
				## also have a "1" as the first point.
### get rid of any other numbers (i.e., -59)
Fmin <- ifelse(Fmin < -0.2 | Fmin > 1.2, 0,Fmin)
Fmin <- as.logical(Fmin)

###### now get index numbers of new minutes (60 values):
Mindex <- index[Fmin]
lastmin <- date2[1,4]

k <- 1			### will denote lines in Mindex
LMindex <- length(Mindex)
TTAGflag <- T
	
### another while loop
#if(i == 1)	browser()
while(k <= LMindex)	{
	
	### get minute's worth of data (from one Mindex to the next)
	t.start <- Mindex[k]
	######## get the end of the minute or the end of the file
	if(k < length(Mindex))	{	t.end <- Mindex[(k+1)]		}
	else			{	t.end <- length(data[,1])	}

	deltaT <- t.end - t.start
### need at least 5 data points 
#if(i==1)	browser()
if(deltaT > 5)	{

		data.m <- data[t.start:t.end,]
		date.m <- date2[t.end,]		### use end of time period for time tag!!!!!
		########## test to see if sec > 50 (usually 58 or 59)
		currsec <- date.m[5]
		currmin <- date.m[4]
		currhr <- date.m[3]
		currday <- date.m[2]
		t1 <- date.m
		if(currsec > 50)	{
			date.m <- fun.Addminute(t1)
			err1ct <- err1ct + 1		## add an error
				}
		
		######## now check for identical consecutive time tags
		######## use minutes as our test (since we have corrected for the 59 sec thing)
		currmin <- date.m[4]
		if(currmin == lastmin)		{
			cat("\n Same time tag!", i,k,j)
			### How many seconds of data do we have?
			deltaS <- date2[t.end,5] - date2[t.start,5]
			if(deltaS > 30)		{
				### add minute to time tag and keep the data
				date.m <- fun.Addminute(t1)
				cat("\n 30-60 sec file")
					}
			else		{
				TTAGflag <- F	#### don't keep the data
				err2ct <- err2ct + 1
					}
						}	### end of if(currmin == lastmin)
	#if(i == 177 & k == 49)		browser()

	### Get median Level and stateNO  (should be pretty constant!)
	Level.m <- median(data.m[,9],na.rm=T)
	NOstate.m <- median(data.m[,8],na.rm=T)
	Zero.m <- median(data.m[,6],na.rm=T)
	Span.m <- median(data.m[,7],na.rm=T)
	status.m <- c(Zero.m,Span.m,NOstate.m,Level.m)
	### get subset of data.m - go 15 pts (secs) from beginning - then avg. 40 pts (sec)
	### this gives 5 sec at the end of a minute to throw away.
	### check the length!!!!!!!!!
	L.m <- length(data.m[,1])
		if(L.m >= 55)	{
			## L.m should be ~ 60 (60 sec/min)
			data.m <- data.m[15:55,]		## throw away 1st 15 sec, last 5 sec.
				}
		if(L.m > 10 & L.m < 55)		{
			### just throw away first 5 sec (these are really bad ones)
			data.m <- data.m[5:L.m,]
				}
	######## if(L.m < 10)  -- just keep all of the data.

	########## Now get averages and standard deviations for columns
	C.avg <- apply(data.m,2,mean,na.rm=T)
	C.avg <- signif(C.avg,d=6)
	C.std <- (apply(data.m,2,var,na.rm=T))^.5
	C.std <- signif(C.std,d=6)
	##########
	#### make sure that Level and stateNO are constant through-out !!!!!!!!!!
	################### Std dev = 0 ################################
	if(C.std[9] == 0)	Lflag <- 1		### all have same NOstatus!!)
	else			Lflag <- 0		### bad somewhere - flag - will go into output matrix as VFair

	##3 now collect the terms you want
	if(TTAGflag == T)	{ 

		if(j > length(M[,1]))	{
			newrow <- rep(NA,25)
			M <- rbind(M,newrow)
					}

		### if TTAGflag == F - just skip putting data in matrix
		M[j,1:6] <- date.m		### time tag at the end of the minute
		M[j,7:11] <- C.avg[1:5]
		M[j,12:15] <- status.m
		M[j,16:17] <- C.avg[10:11]	
		M[j,18:22] <- C.std[1:5]		### STd deviations of trace gas voltages
		M[j,23:24] <- C.std[8:9]			### std dev. of NOstate (tells if we are all on 1 NO status) and Level
		M[j,25] <- Lflag
		j <- j + 1
				}
	############################ will leave concentrations as NAs
	#if(j == 11000)	browser()
	#if(j==10)	browser()
	#cat("\n j =",j)
	#if(j == 1190)	browser()
	###### increment and move on
			}		## end of if(deltaT > 10)

	k <- k + 1
	lastmin <- currmin
	TTAGflag <- T
	}		## end of while(k <= LMindex)



		}		## end of if(filename.i !=0)
else		{
		j <- j + 1
		Flength <- length(M[1,])
		M[j,] <- rep(NA,length=Flength)
		}			### end of else


	#if(i == 365)	browser()
	i <- i + 1	

	}		## end of while (i <= Lfile)

## output error counts
cat("\n 59sec errors:",err1ct)
cat("\n Number of files:",Lfile)
cat("\n Double ttag errors - no data:",err2ct)
	
cat("\n Finish with output file")
cat("\n i = ",i,"j=",j)
#browser()
## cut off extra rows at end.			
if(j < length(M[,1]))		M <- M[1:j,]	
	
############### should now have our 1 min. avg'd output file: M


if(PLOT == T)	{
##### Make some plots to see what's going on:
par(mfrow=c(2,1))
plot(M[,6],M[,7],type="l",ylab="so2")
plot(M[,6],M[,8],type="l",ylab="O3")
browser()
plot(M[,6],M[,10],type="l",ylab="CO2")
plot(M[,6],M[,11],type="l",ylab="H2O")
browser()
## for NO - break into types (by stateNO)
NOmode <- M[,14] == 1
NOxmode <- M[,14] == 3
NOymode <- M[,14] == 5
plot(M[NOmode,6],M[NOmode,9],type="l",ylab="NO")
plot(M[NOymode,6],M[NOymode,9],type="l",ylab="NOy")
lines(M[NOxmode,6],M[NOxmode,9],col=2,lty=2)
cat("\n end of plots - ready for output")
#browser()

	}	## end of if (plot==T)



### Put date column (D)into M-matrix
### Put in 1 more flag at end - to denote data derived from 1 sec data
Flag1s <- rep(1,length(M[,1]))
M <- cbind(M,Flag1s)

fun.asci(M,outfile)

#browser()
### Now we will have a loop to insert Labview avg'd data into what we have just averaged. 
### (means of filling gaps when Labview was averaging but not saving the 1 sec data!!)
#########################################################################################


Ldata <- length(DATA.1min[,1])
i <- 1		## index of original calculated 1 min data
k <- 1		## index of Labview avg'd data		
j <- 1		## substitution counter

#N <- matrix(NA,ncol=26,nrow=length(DATA.1min[,1]))
cat("\n Start gap-filling loop.  Ldata=", Ldata)
#### Loop to substitute in Labview avg'd data.
while(k <= Ldata)	{

test2 <- 0
#if(k == 30)	browser()
## get dec. day time tag for both the calc avg's and Labview 1 min. avg's
ttag.i <- M[i,6]
cat("\n i = ",i,"Time tag=",ttag.i)
p <- is.na(ttag.i)

	if(p == 0)	{
		d.ttag <- DATA.1min[k,6] - ttag.i
		## if they are the same - then put computed data in the row.
		test <- is.na(d.ttag)
		if(test==1)	browser()

		if(abs(d.ttag) < 0.0005)	DATA.1min[k,] <- M[i,]
		else 	{
		### 2 cases - common case - d.ttag will be negative- so 1 sec data is later than 5 min.
		### Case 1:
			if(d.ttag < 0)	{
			### need to increment k until they match times again.
			cat("\n Time tag diff: k=",k)
				while(abs(d.ttag) > 0.0006)	{
					k <- k + 1
					test2 <- test2 + 1
					d.ttag <- DATA.1min[k,6] - ttag.i
						cat("\n delta ttag", d.ttag)
						test <- is.na(d.ttag)
						if(test==1)	browser()
						j <- j + 1
						### if ttag "jumps" over the matching time.
						if(d.ttag > 0)	 d.ttag <- 0
						if(test2 > 1000)		cat("\n Substitution:",test2)
						}
					#browser()
					#test2 <- 0
					DATA.1min[k,] <- M[i,]		## once matched again - put computed data in row.
							}		## end of if(d.ttag < 0)
			### Case 2: 1 sec data is before the 5 min data (mostly at beginning of file!)
			if(d.ttag > 0.0006)	{
				if((k-1) == 0)	{
				### special case - at the beginning -  just add the row at the beginning		
					cat("\n Time tag diff - file beginning: k=", k)
					j <- j + 1
					#browser()
					DATA.1min <- rbind(M[i,],DATA.1min)
					
							}	## end of if((k-1) -- 0)
				else		{
				### put this row between rows(k - 1) and k
					cat("\n Time tag diff - Insert row: k=", k)
					j <- j + 1
					#browser()
					LK <- length(DATA.1min[,1])
					DATA.1min <- rbind(DATA.1min[1:(k-1),],M[i,],DATA.1min[k:LK,])
						}	
					}	## end of if(d.ttag > 0.0006)

			}		## end of else
		## go to the next sample.
		k <- k + 1
		}	## end of if(ttag.i is not an NA, p==0)
else		{
	cat("\n Time tag error",k,i)
		}
#cat("\n Stop for testing")
#browser()
	i <- i + 1

}		## end of substitution while loop


cat("\n Number of Substitutions:",j)
cat("\n Total data lines:",length(DATA.1min[,1]))
cat("\n end of substitution - ready for output again")

if(PLOT == T)	{
##### Make some plots to see what's going on:
par(mfrow=c(2,1))
plot(DATA.1min[,6],DATA.1min[,7],type="l",ylab="so2")
lines(M[,6],M[,7],col=2, lty=2)
plot(DATA.1min[,6],DATA.1min[,8],type="l",ylab="O3")
lines(M[,6],M[,8],col=2,lty=2)
browser()
plot(DATA.1min[,6],DATA.1min[,10],type="l",ylab="CO2")
lines(M[,6],M[,10],col=2, lty=2)

plot(DATA.1min[,6],DATA.1min[,11],type="l",ylab="H2O")
lines(M[,6],M[,11],col=2, lty=2)
#browser()
## for NO - break into types (by stateNO)
NOmode <- M[,14] == 1
NOxmode <- M[,14] == 3
NOymode <- M[,14] == 5
plot(M[NOmode,6],M[NOmode,9],type="l",ylab="NO")
plot(M[NOymode,6],M[NOymode,9],type="l",ylab="NOy")
lines(M[NOxmode,6],M[NOxmode,9],col=2,lty=2)
browser()

	}	## end of if (plot==T)




browser()

fun.asci(DATA.1min,outfile2)




	}

########## end of function ############################################################33

######### OUTPUT FILE ##########################################
#### REAd in data #######
#### works for 2009 data #############
#M <- matrix(scan(infile,sep=",",fill=T,skip=2,list(date="",FS=0,Vso2=0,Vo3=0,Vno=0,Vco2=0,Vh2o=0,Cso2=0,Co3c=0,Cno=0,Cco2=0,Ch2o=0,Zero=0,Span=0,NOstate=0,Level=0,Vfstd=0,Vftot=0,Vfair=0,sdVso2=0,sdVo3=0,sdVno=0,sdVco2=0,sdVh2o=0,sdSO2=0,sdO3=0,sdNO=0,sdCO2=0,sdH2O=0,sdZ=0,sdS=0,sdNOs=0,sdL=0,sdFstd=0,sdFtot=0)),byrow=T,ncol=35)
## parse data:

########################################################################################
##############3 subroutines ########################

fun.despike <- function(x, name, n)
	{

# Subroutine which takes a time series, x, with a given "name" 
# and filters out any obvious spikes in the data. 
# Can be used for sonic data(name = "u","v",etc), CO2 data (name =
# "co2","h2o","t.ec",etc), Kr Hygrometer data (name = "h2o.kr"), etc.
# The filter used is a essentially n*standard deviation.  Anything
# values that fall outside of n*sigma will be replaced with the
# median value.

# Define the filter - Our cut-off :
		y <- n*(var(x,na.rm=T)^0.5)

# Look at positive spikes first :
	test.1 <- x > (median(x,na.rm=T) + y)
		z2 <- sum(test.1,na.rm=T)
	if(z2 > 0)	{
		cat(" :\n replacing ", z2)
		cat(" (+) spikes in ",name)
		x <- ifelse(x > median(x,na.rm=T)+y, median(x,na.rm=T), x)	}

# Look at negative spikes :
	test.1 <- x < (median(x,na.rm=T) - y)
		z3 <- sum(test.1,na.rm=T)
	if(z3 > 0)	{
		cat(" :\n replacing ", z3)
		cat(" (-) spikes in ",name)
		x <- ifelse(x < median(x,na.rm=T)-y, median(x,na.rm=T), x)	}

# Generate a flag which tells when greater than 5% of the time series
# has spikes in it.
	per.spike <- (z2 + z3)/ length(x)
	if(per.spike>0.04)	{
		SPIKE <- F	}
	else	{
		SPIKE <- T	}

# return list of filtered x and SPIKE flag
list(x=x, SPIKE=SPIKE)

}

########################################################################

fun.fastplot <- function(t,x, Label)
{

# subroutine to plot profiles - can change xax and yax

## t = time tag matrix - has six columns with six time tags
## x = matrix with 6 levels (columns 1-6, levels 1-6)
## Label = characters for ylabel

p <- is.na(x[,1])
ind <- seq(from=1,to=length(p),by=1)
samples <- ind[p==0]
I.start <- min(samples,na.rm=T)
I.end <- max(samples,na.rm=T)
xax <- c(t[I.start,1],t[I.end,1])
y.min <- min(x,na.rm=T)
y.max <- max(x,na.rm=T)
yax <- c(y.min,y.max)

#browser()

plot(t[,1],x[,1],type="l",xax,ylab=Label,ylim=yax)
lines(t[,2],x[,2],col=2)
lines(t[,3],x[,3],col=3)
lines(t[,4],x[,4],col=4)
lines(t[,5],x[,5],col=5)
lines(t[,6],x[,6],col=6)

#### create legend
leg.names <- c("z=2m","z=10.9m","z=16.7m","z=23.9m","z=30.3m","z=39.7m")
#colors <- c("black","red","green","blue","cyano","magenta")
#legend(locator(1),leg.names,col=c(1,2,3,4,5,6),fill=T)
#legend(5,1,leg.names,col=c(1,2,3,4,5,6),fill=T)
text(locator(1),"black=2m, red=10.9m, grn=16.7m, blue=23.9m, cyano=30.3m, mag.=39.7m")


}

##################################################

fun.asci <- function(x,data.file)
{

## function to output an asci file

y <- matrix(NA,nrow=length(x[1,]),ncol=length(x[,1]))
i 	<- 1
while(i <= length(x[1,])) 
{	
	y[i,] 	<- x[,i]
	i 	<- i+1
}

### Output data:
write(y,data.file,ncol=length(x[1,]))


}



######################################################################
fun.Camptimetag <- function(day,hrmin,sec)
#
# Converts Campbell-type time tags into decimal day time tags:
# note - we assume that we know the year!!!!
{

## get hours and minutes separated
hr <- floor(hrmin/100)			## give hour of day
min <- (hrmin/100) - floor(hrmin/100)	## min/100
min <- 100*min

## put everything in days and add it up:
min.sec <- min + sec/60
hr.min <- hr + min.sec/60
ttag <- day + hr.min/24		
ttag <- signif(ttag,d=8)
# output ttag
	ttag

}


###############################################################

	
fun.linearextrap <- function(index,x,xx,yy)
	{
## subroutine to take a time series (x) and and index (0->length)
## and a set of index pts/calibration points (xx and yy) and create a 
## time series of calibration constants.
### does linear extrapolation between different xx points.
### xx == vector of indices
### yy == corresponding signals (for zero's) or gain coeff. (for spans)

cat("\nTest linear extrap")
#browser()

### do line fit between each pt.
Lz <- length(xx) - 1
i <- 1

while(i <= Lz)	{
	
j <- i + 1
x1 <- xx[i:j]
y1 <- yy[i:j]
#browser()
Lfit <- lm.fit <- lm(y1 ~ x1)
Int <- Lfit$coeff[1]
m <- Lfit$coeff[2]
x[xx[i]:xx[j]] <- Int + (m*index[xx[i]:xx[j]])
i <- i + 1

}		###end of loop

#### then use last zero points to fill out the zero time series
Lz <- length(yy)
Last.zero <- yy[Lz]
Li <- length(index)
Lend <- xx[Lz]
x[Lend:Li] <- yy[Lz]



x


}

#####################################################################
fun.skewness <- function(x,OUT=F)
	{
## function designed to take a time series of x and 
## calculate the mean, median, std. dev., skewness and
## kurtosis.

# get rid of NAs for the rest :
xx <- is.na(x)
if(sum(xx,na.rm=T) < length(xx))	{

N <- length(x) - sum(xx,na.rm=T)
xbar <- mean(x,na.rm=T)
xmed <- median(x,na.rm=T)
x.min <- min(x,na.rm=T)
x.max <- max(x,na.rm=T)

sdx <- var(x,na.rm=T)^.5

Skx <- mean((x - mean(x))^3) / (sdx^3)

Krx <- mean((x - mean(x))^4) / (sdx^4)

if(OUT==T)	{
	cat("\n N = ",N)
	cat("\n xbar = ",xbar)
	cat("\n xmed = ",xmed)
	cat("\n sdx = ",sdx)
	cat("\n x.min = ",x.min)
	cat("\n x.max = ",x.max)
		}
	}	## end of if(sum < length)
else	{
	N <- 0
	xbar <- NA
	xmed <- NA
	sdx <- NA
	Skx <- NA
	Krx <- NA
		}
	
list(N=N,xbar=xbar,xmed=xmed,sdx=sdx,Skx=Skx,Krx=Krx,x.min=x.min,x.max=x.max)

#invisible()

# end function
}
##############################################################

######################################################



fun.Ltimetag <- function(t1)
	{
# function to convert a POSIXct object to a numeric string of the type:
#	Julian day,  hour of day, min of day,  sec
# output as a vector of numbers:

year <- strftime(t1, "%Y")
year <- (as.numeric(year))
jday <- strftime(t1,"%j")
jday <- as.numeric(jday)
hday <- strftime(t1,"%H")
hday <- (as.numeric(hday))*100
mday <- strftime(t1, "%M")
mday <- as.numeric(mday)
Sday <- strftime(t1, "%S")
Sday <- as.numeric(Sday)

## put in to vector
ttag <- cbind(year,jday,hday,mday,Sday)

}
 
##################################################

fun.Addminute <- function(t1)
		{
### function to "round" a time tag to the next minute
### t1 = vector of yy,jday,hrx100,min,sec,DOY

### Assume we want to round up to the next minute:
	t1[5] <- 0
	t1[4] <- t1[4] + 1
	cat("\n 59 sec alarm!")
		if(t1[4] == 60)	{
				t1[4] <- t1[4] - 60
				t1[3] <- t1[3] + 100
				cat("\n hour alarm!")
			if(t1[3] == 2400)	{
				t1[3] <- t1[3] - 24
				t1[2] <- t1[2] + 1
				cat("\n day alarm!")
	### until Dec. 31 - won't worry about changing years.
								}
							}

## output t1
t1
######## end of function
						}		

######################################################