fun.MEFtracegas011 <- function(PLOT=T,GRAD=T,RAWPLOT=F,YEAR=2011, Z2=0,NewZ=1,UseNewZ=0)
	{

## Function designed to go through trace gas data from Manitou Forest and :
## (1) calculate approx. concentrations (using lab cals. // or field cals.)
## (2) put in calibrations (if available)
## (3) flag bad data -create flag files
## (3) Make some plots the raw and filtered data.
######## (will write second program for looking at profiles)

## originally written for CHATs
## modified for AMAZE 2008
## modified for Beachon-SRM08
### works from Oct. 2009 -> present.
## PLOT - shows data plots
## RAWPLOT == T, plots raw data, == F, plots filtered data
## Z2  == 0, no zero substitution.  Z2 = 1, put in a zero for NOx and NOy based on Sept. 2010 zero air tests.
##		For use with data before Sept. 9, 2010 (no zero air offsets).
## NewZ = 0, Old Data format, no zero air offsets.  == 1:  New cal file format - includes zero air offsets.
## UseNewZ == 0, Pre-rxn NO offsets,  ==1, uses zero air offsets for NOx and NOy.
##   (note - for UseNewZ to be == 1, need NewZ to also == 1!!!)

### 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 (5 min. averages)
#datafile <- scan("c:/data/Manitou/Manitou09/Datafiles09.txt",what="character")
datafile <- scan("c:/data/Manitou/Manitou011/Datafiles011.txt",what="character")

## Input files:
infile <- datafile[53]
Calfile <- datafile[21]
inflags <- datafile[16]
MONTH <- datafile[30]
### Output data
outfile <- datafile[2]
outflags <- datafile[17]
########################


######## assign a number to MONTH
if(MONTH == "jan")	month <- 1
if(MONTH == "feb")	month <- 2
if(MONTH == "mar")	month <- 3
if(MONTH == "apr")	month <- 4
if(MONTH == "may")	month <- 5
if(MONTH == "jun")	month <- 6
if(MONTH == "jul")	month <- 7
if(MONTH == "aug")	month <- 8
if(MONTH == "sepA")	month <- 9
if(MONTH == "sepB")	month <- 9
if(MONTH == "oct")	month <- 10
if(MONTH == "nov")	month <- 11
if(MONTH == "dec")	month <- 12
cat("\n Month of data:",MONTH, "   ",month)
cat("\n")

#### REAd in data #######
#### works for 2009 data #############
d <- matrix(scan(infile,sep=",",fill=T,skip=1,list(yr=0,jday=0,hr=0,mn=0,sec=0,doy=0,Vso2=0,Vo3=0,Vno=0,Vco2=0,Vh2o=0,Zero=0,Span=0,NOstate=0,Level=0,Vfstd=0,Vftot=0,sdVso2=0,sdVo3=0,sdVno=0,sdVco2=0,sdVh2o=0,sdNOs=0,sdL=0,Lflag=0,Flag1s=0)),byrow=T,ncol=26)
## parse data:
#### Raw voltages
yr <- as.ts(d[[1]])
jday <- as.ts(d[[2]])
hr <- as.ts(d[[3]])
mn <- as.ts(d[[4]])
sc <- as.ts(d[[5]])
dday <- as.ts(d[[6]])
ttag <- cbind(yr,jday,hr,mn,sc,dday)
########################
so2 <- as.ts(d[[7]])
o3 <- as.ts(d[[8]])
no <- as.ts(d[[9]])
co2 <- as.ts(d[[10]])
h2o <- as.ts(d[[11]])
### STatus ########
Zero <- as.ts(d[[12]])
Span <- as.ts(d[[13]])
statNO <- as.ts(d[[14]])
Level <- as.ts(d[[15]])
### Flows (voltages) ###3
VFstd <- as.ts(d[[16]])
VFtot <- as.ts(d[[17]])
#### Standard deviations ####
sdSO2 <- as.ts(d[[18]])
sdO3 <- as.ts(d[[19]])
sdNO <- as.ts(d[[20]])
sdCO2 <- as.ts(d[[21]])
sdH2O <- as.ts(d[[22]])
Lflag <- as.ts(d[[25]])

### could be a few values of "60" in minutes. replace with "0"
mn <- ifelse(mn > 59.5,0,mn)
### could be a few values of "2376" in hr. replace with "0"
hr <- ifelse(hr > 2350,0,hr)
###########################################
#### Read in Flags ######
Flags <- matrix(scan(inflags, sep=",", fill=T,skip=1, list(date=0, Fso2=0,Fo3=0,Fno=0,Fnox=0,Fnoy=0, Fco2=0,Fh2o=0)),byrow=T,ncol=8)

### Read in Calibration files ##### 
if(NewZ == 0)	{ 
calfile <- matrix(scan(Calfile, sep=" ", fill=T, list(year=0,day=0,hr=0,min=0,sec=0,dday=0,index=0,Zo3=0,Zso2=0,Zno=0,Znox=0,Znoy=0,Zco2=0,Zh2o=0,So3=0,Sso2=0,Sno=0,Snox=0,Snoy=0,Sco2=0,Sh2o=0)),byrow=T,ncol=21)
		}
else		{
calfile <- matrix(scan(Calfile, sep=" ", fill=T, list(year=0,day=0,hr=0,min=0,sec=0,dday=0,index=0,Zo3=0,Zso2=0,Zno=0,Znox=0,Znoy=0,Zco2=0,Zh2o=0,So3=0,Sso2=0,Sno=0,Snox=0,Snoy=0,Sco2=0,Sh2o=0,Z2.no=0,Z2.nox=0,Z2.noy=0)),byrow=T,ncol=24)
		}

################ Flags #######################
Fdate <- as.ts(Flags[[1]])
Fso2 <- as.ts(Flags[[2]])
Fo3 <- as.ts(Flags[[3]])
Fco2 <- as.ts(Flags[[7]])
Fh2o <- as.ts(Flags[[8]])
Fno <- as.ts(Flags[[4]])
Fnox <- as.ts(Flags[[5]])
Fnoy <- as.ts(Flags[[6]])
########## Flag Codes:
## 1 == Okay, sampling
## 0 == bad data (for whatever reason)
###### REMEMBER - to replace autocalibration time periods with "no data" and flag!!!!

################## Calibrations (zero/gains) ##################
Cdday <- as.ts(calfile[[6]])
Cindex <- as.ts(calfile[[7]])
Z.o3 <- as.ts(calfile[[8]])
Z.so2 <- as.ts(calfile[[9]])
Z.no <- as.ts(calfile[[10]])
Z.nox <- as.ts(calfile[[11]])
Z.noy <- as.ts(calfile[[12]])
Z.co2 <- as.ts(calfile[[13]])
Z.h2o <- as.ts(calfile[[14]])
G.o3 <- as.ts(calfile[[15]])
G.so2 <- as.ts(calfile[[16]])
G.no <- as.ts(calfile[[17]])
### assume G.no = G.nox = G.noy to start with.
G.nox <- G.no
G.noy <- as.ts(calfile[[19]])
G.co2 <- as.ts(calfile[[20]])
G.h2o <- as.ts(calfile[[21]])
if(NewZ == 1)	{
	Z2.no <- as.ts(calfile[[22]])
	Z2.nox <- as.ts(calfile[[23]])
	Z2.noy <- as.ts(calfile[[24]])
		}
### make sure that Fday and day are the same size
Lday <- length(dday)
Lfday <- length(Fdate)
Lcday <- length(Cdday)
########### if not equal - no more calculations:

cat("\n check on data!")
browser()

############ CAN GET ALL OF THIS FROM INPUT FILES !!!!!!!! ##############3
if(Lday == Lfday & Lday == Lcday)	{
###########################################################################
### get time tags
#date1 <- strptime(date1, "%m/%d/%Y %H:%M")
#date1 <- fun.Ltimetag(date1)
#### now have numeric dates - matrix of yy,jday,hrx100,min,sec
## get decimal day:
#dday <- date1[,2] + (date1[,3]/2400) + (date1[,4]/1440)
#dday <- signif(dday,d=7)
#date1 <- cbind(date1,dday)		## add decimal day to date.
## also create an index - sequential numbers to reference
date1 <- cbind(yr,jday,hr,mn,sc,dday)		## matrix with 6 columns.
Ldata <- length(dday)
index <- seq(from=1,to=Ldata,by=1)
############################################

cat("\n Look at standard deviations and think about a flag!!")

par(mfrow=c(2,1))
## Plot the raw voltages one at a time.
plot(dday,so2,type="l",ylab="SO2")
plot(dday,o3,type="l",ylab="O3")
browser()

plot(dday,co2,type="l",ylab="CO2")
plot(dday,h2o,type="l",ylab="H2O")
browser()

plot(dday,no,type="l",ylab="NO")
#lines(dday,Z.no,lty=2,col=2)
cat("\n End of voltage plots")
browser()

### get rid of NAs in data -set them to something crazy that will get flagged out
p <- is.na(no)
no <- ifelse(p == T,-100,no)
p <- is.na(o3)
o3 <- ifelse(p == T,-100,o3)
p <- is.na(so2)
so2 <- ifelse(p == T,-700,so2)
p <- is.na(co2)
co2 <- ifelse(p == T,-100,co2)
p <- is.na(h2o)
h2o <- ifelse(p == T,-100,h2o)

##### Now filter data for easily found bad points:
###### put into flags first ###########################
Fo3 <- ifelse(o3 < 0,0,Fo3)
Fso2 <- ifelse(so2 < -500, 0,Fso2)
Fno <- ifelse(no < 300,0,Fno)
Fnox <- ifelse(no < 300,0,Fnox)
Fnoy <- ifelse(no < 300,0,Fnoy)
Fco2 <- ifelse(co2 < 0,0,Fco2)
Fh2o <- ifelse(h2o < 0,0,Fh2o)

Fo3 <- ifelse(o3 > 1100,0,Fo3)		## > 215 ppbv O3
Fso2 <- ifelse(so2 > 2390,0,Fso2)	## 
Fno <- ifelse(no > 2390,0,Fno)
Fnox <- ifelse(no > 2390,0,Fnox)
Fnoy <- ifelse(no > 2390,0,Fnoy)
Fco2 <- ifelse(co2 > 2400,0,Fco2)
Fh2o <- ifelse(h2o > 2400,0,Fh2o)
######## maybe create a flag based on the standard deviations !!! (rapidly changing data).
##########################################################################################

############ Now Filter the data using the Flag vectors ############################
####################################################################
o3 <- ifelse(Fo3 == 0,NA,o3)
so2 <- ifelse(Fso2 == 0,NA,so2)
no <- ifelse(Fno == 0,NA,no)
co2 <- ifelse(Fco2 == 0,NA,co2)
h2o <- ifelse(Fh2o == 0,NA,h2o)
############ GET RID OF ZERO AND SPAN'S ################
### DEfault
Z1 <- Zero > 0.3 & Zero < 1.4		## zero air offsets // works no matter what!
no <- ifelse(Z1==T ,NA,no)
co2 <- ifelse(Z1==T ,NA,co2)
h2o <- ifelse(Z1==T ,NA,h2o)
so2 <- ifelse(Z1 == T,NA,so2)
o3 <- ifelse(Z1==T,NA,o3)
Z2 <- Zero > 1.4			## NO pre-rxn offset (new cal. method)
no <- ifelse(Z2==T ,NA,no)

## Spans are a bit crazier.
######################### need to correct here!!! 
##################################################
S2 <- Span > 0.4 & Span < 1.4
S <- S2 == T & Level == 3
no <- ifelse(S==T ,NA,no)
so2 <- ifelse(S == T,NA,so2)

co2 <- ifelse(Span > 1.5,NA,co2)
h2o <- ifelse(Span > 1.5,NA,h2o)

## Special Case:  Jan-Mar,2011 was running the wrong code
if(month <= 3 & YEAR==2011)	{
	so2 <- ifelse(Span > 1.5,NA,so2)
	no <- ifelse(Span > 1.5,NA,no)
	cat("\n Get rid of NO and SO2 spans!")
				}

### finally - get rid of NAs in Level and statNO
p <- is.na(Level)
Level <- ifelse(p == T,0,Level)
p <- is.na(statNO)
statNO <- ifelse(p == T,0,statNO)

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

### Now have our time series to work with
##########################################

################### PUT IN CALIBRATIONS (FROM CALFILE) #############################
###############
### Ozone #####
###############		
o3.ppb <- (o3 - Z.o3)/G.o3		### units = ppbv

###############
### SO2   #####
###############
so2.ppb <- (so2 - Z.so2)/G.so2		### units = ppbv		

###############################################
##### CO2 and H2O ########
###############################################
co2.ppm <- (co2 - Z.co2)/G.co2		### units = ppmv
h2o.mmol <- (h2o - Z.h2o)/G.h2o		### units = mmol/mol

########################################################################
### NO   #####
##############
### now have time series of NO,Nox and NOy zeros and gains.
### have to parse by NO status: ==1, NO,  ==2,3, NOx, ==4,5, NOy
#####################################################
##### put in calibrations ##### Units = ppbv ######
NOmode <- statNO > 0.5 & statNO < 1.5
NOxmode <- statNO > 1.5 & statNO < 3.5
NOymode <- statNO > 3.5

##### on June 9th (day 160), 2011, changed logger program to do 2 min. of NO, 2 min. of NOx, 1 min of NOy
if(month > 6)	{
	NOmode <- statNO > 0.5 & statNO < 2.5
	NOxmode <- statNO > 2.5 & statNO < 4.5
	NOymode <- statNO > 4.5
		}

#cat("\n test new statNO")
#browser()

## special case for June:
if(month == 6 & YEAR == 2011)	{
	## must make composite flag for NO status
	test1 <- dday > 160.45
	NOmode[test1] <- statNO[test] > 0.5 & statNO[test] < 2.5 
	NOxmode[test1] <- statNO[test] > 2.5 & statNO[test] < 4.5
	NOymode[test1] <- statNO[test] > 4.5
	}

if(month==8 & YEAR == 2010)	{
	## for Aug. data - read in a new time series of Z.no (which will also be Z.nox)
	## not so worried about Z.noy!
	d <- matrix(scan("c:/data/Manitou/Manitou010/ProfileData/Zno.aug010.txt", sep=" ", fill=T, list(dday=0,Zno=0)),byrow=T,ncol=2)
	NewZ.no <- as.ts(d[[2]])
	### before substituting - make sure it is the same length as our Z.no
	L1no <- length(Z.no)
	L2no <- length(NewZ.no)
		if(L1no == L2no)		{
			Z.no <- NewZ.no
			Z.nox <- Z.no
						}
		else				{
			cat("\n something wrong with Aug. 2010 Z.no - length not correct!")
			browser()
				}
		### also get rid of NOx measurements between 


		}		### end of if(month==8 and YEAR=2010)

### stop and plot the Z.no.
#par(mfrow=c(1,1))
#plot(dday,no,type="l",ylab="NO")
#points(dday,Z.no,lty=2,col=2,cex=0.3)
#cat("\n Plot Z.no again")
#############################################################
#browser()

######### For old data that does not have zero air offsets measured.
if(Z2 == 1)	{
	### put in "zero air" zeros for NOx and NOy  (will leave NO alone)
	### will "add" some signal to NOx and NOy zeros
	Z.nox <- Z.nox + 43.8			### Z.nox(prerxn) + [Z.nox(0air) - Z.no(0air)] = Z.nox(prerxn) + [59.9 - 15.1]
	Z.noy <- Z.noy + 27.4			### Z.noy(prerxn) + [Z.noy(0air) - Z.no(0air)] = Z.noy(prerxn) + [42.5 - 15.1]
		}

######## if we have offsets measured with zero air and want to use them:
if(UseNewZ == 1 & NewZ == 1)	{
	Z.no <- Z2.no
	Z.nox <- Z2.nox
	Z.noy <- Z2.noy
	############ if Pre-Rxn Z.no is OK, then....
	#Z.nox <- Z2.nox - (Z2.no - Z.no)		### Z(NOx)(ZAir) - [ Z(NO)(Zair) - Z(NO)(prerxn) ]  in case there is some small amt of NO in zero air. 
	#Z.noy <- Z2.noy	- (Z2.no - Z.no)		### Z.noy(prerxn) + [Z.noy(0air) - Z.no(0air)] = Z.noy(prerxn) + [59.9 - 15.1]
		}


no.ppb <- ifelse(NOmode==T,((no - Z.no)/G.no),NA)
nox.ppb <- ifelse(NOxmode==T,((no - Z.nox)/G.nox),NA) 
noy.ppb <- ifelse(NOymode==T,((no - Z.noy)/G.noy),NA)
### now have vectors of NO, NOx and NOy concentrations with other modes == -100 (to filter out later) 


### Now we need to put in the NOx and NOy flags:
nox.ppb <- ifelse(Fnox == 0,NA,nox.ppb)
noy.ppb <- ifelse(Fnoy == 0,NA,noy.ppb)

######### special Case filters
if(month == 7)	{
	## for July, 2011,  NOx conv.eff. tests, Low CO2
	test <- co2.ppm < 350
	so2.ppb <- ifelse(test==T,NA,so2.ppb)
	o3.ppb <- ifelse(test==T,NA,o3.ppb)
	no.ppb <- ifelse(test==T,NA,no.ppb)
	nox.ppb <- ifelse(test==T,NA,nox.ppb)
	noy.ppb <- ifelse(test==T,NA,noy.ppb)
	co2.ppm <- ifelse(test==T,NA,co2.ppm)
	h2o.mmol <- ifelse(test==T,NA,h2o.mmol)
		}

### WHERE'S PHOTOLYTIC CONVERSION IN HERE!!!!! - put in later before 5 min. averaging.
cat("\n Look at concentration data and make manual filters to data !!")
par(mfrow=c(2,1))
## Plot the raw voltages one at a time.
plot(date1[,6],so2.ppb,type="l",ylab="SO2")
#browser()
plot(date1[,6],o3.ppb,type="l",ylab="O3")
cat("\n Filter data manually:  o3.ppb and so2.ppb")
browser()

#browser()
plot(date1[,6],co2.ppm,type="l",ylab="CO2")
#browser()
plot(date1[,6],h2o.mmol,type="l",ylab="H2O")
cat("\n Filter data manually:  co2.ppm and h2o.mmol")
browser()
plot(date1[,6][NOmode],no.ppb[NOmode],type="l",ylab="NO", ylim=c(-.3,15))
lines(date1[,6][NOxmode],nox.ppb[NOxmode],col=2)
lines(date1[,6][NOymode],noy.ppb[NOymode],col=3)
cat("\n Look at filtered data!!")
##########################################
cat("\n Filter data manually:  no.ppb, nox.ppb and noy.ppb")
browser()

### then update the flags from manual filtering
p <- is.na(o3.ppb)
Fo3 <- ifelse(p==T,0,Fo3)
p <- is.na(so2.ppb)
Fso2 <- ifelse(p==T,0,Fso2)
p <- is.na(co2.ppm)
Fco2 <- ifelse(p==T,0,Fco2)
p <- is.na(h2o.mmol)
Fh2o <- ifelse(p==T,0,h2o.mmol)

#p <- is.na(no.ppb)
#Fno <- ifelse(p==T,0,Fno)
#p <- is.na(nox.ppb)
#Fnox <- ifelse(p==T,0,Fnox)
#p <- is.na(noy.ppb)
#Fnoy <- ifelse(p==T,0,Fnoy)

### Need to output 1 min. data 
### file for each trace gas (O3, SO2, CO2, H2O, etc.)
#######################################################

### now - need loop that goes through and makes 5 min. averages
### averages O3, So2, Co2 and H2O. 
### gets the proper mode for NO, NOx and NOy
### For SO2, O3, CO2 and H2O - averages minutes 2-5 (this assures complete line flushing).
### For NO - get statNO == 1,  For NOx - statNO == 3,  For NOy - statNO == 5
### After June 9th, 2011:  
### For NO - get statNO == 1-2,  For NOx - statNO == 3-4,  For NOy - statNO == 5
### ALSO - make sure Level stays constant!!!!!! 

### put everything into a matrix for averaging.
data <- cbind(index,date1,so2.ppb,o3.ppb,no.ppb,nox.ppb,noy.ppb,co2.ppm,h2o.mmol,Level,statNO)
flags <- cbind(index,date1,Fso2,Fo3,Fno,Fnox,Fnoy,Fco2,Fh2o)
### get output matrices to fill in.
NumRows <- (Lday/5) + 40
M <- matrix(NA,ncol=16,nrow=NumRows)
F <- matrix(0,ncol=14,nrow=NumRows)
### 
cat("\n REady for averaging loop")
cat("\n Check time tags for Sec != 0")
browser()
########## Note - this loop takes a long time to compute!!!! not sure why!!!?????
i <- 1
j <- i + 4 
k <- 1
looplength <- Lday - 5

while(i <= looplength)	{
	## Initially - must start on a 5 min. interval.
	test <- as.integer(date1[i,4]/5)
	test2 <- date1[i,4]/5
	if(test != test2)	{
		####increment i and move along until we hit a multiple of 5 min.
		i <- i + 1
		j <- j + 4
		cat("\n not on 5 min. interval",i,date1[i,4:5])
		if(i==6159)	browser()
			}
	else	{
		### get data and flags
		i <- i + 1		### need to start the next 5 min. period!! (note - time tag is at end of the minute!)
		j <- i + 4
		data.i <- data[i:j,]
		flags.i <- flags[i:j,]
			### first get the flags to make sure data is valid...
			## now get flags
			fso2.i <- sum(flags.i[2:5,8],na.rm=T)
			fo3.i <- sum(flags.i[2:5,9],na.rm=T)
			fco2.i <- sum(flags.i[2:5,13],na.rm=T)
			fh2o.i <- sum(flags.i[2:5,14],na.rm=T)
			fno.i <- flags.i[1,10]
			fnox.i <- sum(flags.i[2:3,11],na.rm=T)
			fnoy.i <- sum(flags.i[4:5,12],na.rm=T)
			
			### now put proper time tags, averages and flags into output matrix.
			## time tag and index number - beginning of period:
			date.i <- date1[(i-1),]		
			index.i <- data[(i-1),1]
			### get simple avg. concentrations (minutes 2-5)
				### (except for H2O - use minutes 3-5)
			if(month == 3)		ind.co2 <- 3
			else			ind.co2 <- 2

			if(fso2.i >= 1)		so2.i <- mean(data.i[2:5,8],na.rm=T)
			else			so2.i <- NA
			
			if(fo3.i >= 1)		o3.i <- mean(data.i[2:5,9],na.rm=T)
			else			o3.i <- NA

			if(fco2.i >= 1)		co2.i <- mean(data.i[ind.co2:5,13],na.rm=T)
			else			co2.i <- NA

			if(fh2o.i >= 1)		h2o.i <- mean(data.i[3:5,14],na.rm=T)
			else			h2o.i <- NA

			####### Note - this should also "miss' all the -100's in the time series here if sync'd properly.
			#if(k == 1347)	browser()
			#if(k == 1348)	browser()
			
			if(fno.i >= 0.9)	no.i <- data.i[1,10]
			else			no.i <- NA

			if(fnox.i >= 1)		nox.i <- data.i[3,11]
			else			nox.i <- NA

			if(fnoy.i >= 1)		noy.i <- data.i[5,12]
			else			noy.i <- NA


		if(month > 6)	{
			if(fno.i >= 0.9)	no.i <- mean(data.i[1:2,10],na.rm=T)
			else			no.i <- NA

			if(fnox.i >= 1)		nox.i <- mean(data.i[3:4,11],na.rm=T)
			else			nox.i <- NA

			if(fnoy.i >= 1)		noy.i <- data.i[5,12]
			else			noy.i <- NA
				}

			### get level (make sure all values are the same level!!!!)
			### see if 4th value of Level == 1st value of Level.
			
			

			Li <- data.i[,15]
			Ltest <- Li[4] - Li[1]
			level.i <- Li[1]
			Levflag.i <- 0

			if(Ltest == 0)	Levflag.i <- 1		
				
			
			#if(k == 1206)	browser()
			#if(k == 1207)	browser()
			

			data.avg <- c(date.i,index.i,so2.i,o3.i,co2.i,h2o.i,no.i,nox.i,noy.i,level.i, Levflag.i)
			M[k,] <- data.avg
			
			## create a 5 min. flag file.  (necessary???)
			fso2.i <- ifelse(fso2.i >= 1,1,0)
			fo3.i <- ifelse(fo3.i >= 1,1,0)
			fco2.i <- ifelse(fco2.i >= 1,1,0)
			fh2o.i <- ifelse(fh2o.i >= 1,1,0)
			fno.i <- ifelse(fno.i >= 1,1,0)
			fnox.i <- ifelse(fnox.i >= 1,1,0)
			fnoy.i <- ifelse(fnoy.i >= 1,1,0)
			flag.avg <- c(date.i,index.i,fso2.i,fo3.i,fno.i,fnox.i,fnoy.i,fco2.i,fh2o.i)
			F[k,] <- flag.avg
			
			#cat("\n index:",i,"mat.index:",k)

			### now increment i and move on
			i <- j			### or i = i + 4
			k <- k + 1		### next line in datafile.

			}		## end of else statement

	}	## end of while loop.


########## Now have 5 min. averages of concentrations #########################

cat("\n stop and test the 5min avg's for NOx!")
browser()

######### put in NOx conversion factor for BLC !!!!!!!!!!!!!!!  #######################################
#######################################################################################################
### put in photolytic conversion factor
PC <- 0.684			### 68.4 % conversion of NO2, 9/27/07 exp.
PC <- 0.604			### 60.4% conversion of NO2, 1/13/10 exp.
PC <- 0.423			### 42.3% conversion of NO2, 7/5/11 exp. - new glass photolysis cell.
M[,13] <- (M[,13] - M[,12])/PC		
### NO2 = (NOx meas - NO)/PC, subtract off NO conc. - then apply photolysis conversion.
########################################################################################################
### filter NO, NO2 and NOy for negative values (maybe from zero's)
p1 <- M[,12] < -.2
p2 <- M[,13] < -.1
p3 <- M[,14] < -.1
M[,12] <- ifelse(p1 == T,NA,M[,12])
M[,13] <- ifelse(p2 == T,NA,M[,13])
M[,14] <- ifelse(p3 == T,NA,M[,14])
## update flag file
Fno <- ifelse(p1 == T,0,Fno)
Fnox <- ifelse(p2 == T,0,Fnox)
Fnoy <- ifelse(p3 == T,0,Fnoy)
###################################

### Other stuff ###################
NOx <- M[,13] + M[,12]		### NOx = NO2 - NO
NOz <- M[,14] - NOx		### NOz = NOy - NOx
NO.NO2 <- M[,12]/M[,13]		### [NO]/[NO2] ratio


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

######################################################################################
##### have final concentrations now // should aleady be filtered.
######################################################################################

####
### need to filter data for flags again....
##### Make some plots of the data !!!!!!!!!!!!!!
### Only important in Profile calculations!!!!!!
### do it in that MEFprofile.R !!!!!!!!!!!!!
### but need to output statNO !!!!!!!!
#########################################################################################################
cat("\n Ready to output concentration data")
browser()
#### Now ready to output the data:

D <- cbind(M,NOx,NOz,NO.NO2)
Dflags <- F

par(mfrow=c(1,1))
p <- is.na(D[,8])
if(sum(p,na.rm=T) < length(p))	plot(D[,6],D[,8],type="l",ylab="SO2,ppbv",xlab="DOY")
browser()
p <- is.na(D[,9])
if(sum(p,na.rm=T) < length(p))	plot(D[,6],D[,9],type="l",ylab="O3,ppbv",xlab="DOY")
browser()
p <- is.na(D[,10])
if(sum(p,na.rm=T) < length(p))	plot(D[,6],D[,10],type="l",ylab="CO2,ppbv",xlab="DOY")
browser()
p <- is.na(D[,11])
if(sum(p,na.rm=T) < length(p))	plot(D[,6],D[,11],type="l",ylab="H2O,ppbv",xlab="DOY")
browser()
NOymax <- max(D[,14],na.rm=T)
p <- is.na(D[,14])
if(sum(p,na.rm=T) < length(p))	plot(D[,6],D[,14],type="l",ylab="NO cpds,ppbv",xlab="DOY", ylim=c(-1,NOymax))
lines(D[,6],D[,13],type="l", col=2)
lines(D[,6],D[,12],type="l",col=3)
browser()
lines(D[,6],D[,17],type="l", col=4)
lines(D[,6],D[,18],type="l",col=5)
browser()
if(sum(p,na.rm=T) < length(p))	plot(D[,6],D[,19],type="l",ylab="NO/NO2",xlab="DOY",ylim=c(-0.2,5))
browser()
cat("\n end of plots")


browser()

fun.asci(D,outfile)
fun.asci(Dflags,outflags)

###### Output file is 

#################################3
	}		### end of if(Lday == Lfday)
else	{
	cat("\n Length Flag file NOT equal to Length of data file")
	cat("\n Something is Wrong - System going down!!!")
	}





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


########################################################################################
##############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)

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