fun.MEFprofile011 <- function(PLOT=F, FILTER=T,NO2ratio=T,YEAR=2011)
	{

## Function designed to take in concentration data from MEF and :
## (1) deal with timing (throw out the first minute on some sensors) day
## (2) Group by profile level
## (3) Make some plots of the profiles   
## (4) output files by level.

## originally written for CHATs
## modified for AMAZE 2008
## modified for Beachon-SRM08

## PLOT - shows data plots
## GRAD == T -- we are computing profiles
## RAWPLOT == T, plots raw data, == F, plots filtered data
## FILTER == T, filter for "extreme" events in NO and SO2.  (i.e, vehicle traffic, etc.)
## NO2ratio == T,  from Aug. 2010 on - added extra column to input data (no/no2 ratio).

### Input text file with names of input/output files in it.
### will create monthly files (5 min. averages)
if(YEAR==2009)	datafile <- scan("c:/data/Manitou/Manitou09/Datafiles09.txt",what="character")
if(YEAR==2011)  datafile <- scan("c:/data/Manitou/Manitou011/Datafiles011.txt",what="character")

## Data flags// is there data??
P.o3 <- T
P.so2 <- T
P.no <- T
P.co2 <- T
P.h2o <- T
P.nox <- T
P.noy <- T
## Parse "datafile" into the input and output filenames
infile <- datafile[2]
### flag file
inflags <- datafile[17]
outflags <- datafile[17]
######################## Output files ###############################
### output file for each level
out1 <- datafile[10]
out2 <- datafile[11]
out3 <- datafile[12]
out4 <- datafile[13]
out5 <- datafile[14]
out6 <- datafile[15]
MONTH <- datafile[30]		## to use for specific filtering.
cat("\n Year:",YEAR,"Month:",MONTH)
#######################
### output profile data - species
so2prof <- datafile[4]
o3prof <- datafile[5]
co2prof <- datafile[6]
h2oprof <- datafile[7]
NOprof <- datafile[8]
NOxprof <- datafile[9]
NOyprof <- datafile[22]
NO2prof <- datafile[23]

so2diel <- datafile[38]
o3diel <- datafile[39]
co2diel <- datafile[40]
h2odiel <- datafile[41]
NOdiel <- datafile[42]
NOxdiel <- datafile[43]
NOydiel <- datafile[44]
NO2diel <- datafile[45]

######## 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 == "sep")	month <- 9
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 #######
if(NO2ratio==T)	{
d <- matrix(scan(infile,fill=T,list(yr=0,day=0,hr=0,min=0,sec=0,doy=0,index=0,so2=0,o3=0,co2=0,h2o=0,NO=0,NOx=0,NOy=0,Level=0,Lflag=0,NO2=0,NOz=0,NO.NO2=0)),byrow=T,ncol=19)
		}
else		{
d <- matrix(scan(infile,fill=T,list(yr=0,day=0,hr=0,min=0,sec=0,doy=0,index=0,so2=0,o3=0,co2=0,h2o=0,NO=0,NOx=0,NOy=0,Level=0,Lflag=0,NO2=0,NOz=0)),byrow=T,ncol=18)
		}

#### Read in Flags ######
Flags <- matrix(scan(inflags,fill=T,list(yr=0, day=0, hr=0,min=0,sec=0,doy=0,index=0, Fso2=0, Fo3=0, Fno=0,Fnox=0,Fnoy=0, Fco2=0,Fh2o=0)),byrow=T,ncol=14)

#### Note - we have switched NOy and NOyZ in the sampling - so we switch here as well.
doy <- as.ts(d[[6]])
day <- as.ts(d[[2]])
year <- as.ts(d[[1]])
min <- as.ts(d[[4]])
index <- as.ts(d[[7]])
hour <- as.ts(d[[3]])
############################
so2 <- as.ts(d[[8]])
o3 <- as.ts(d[[9]])
co2 <- as.ts(d[[10]])
h2o <- as.ts(d[[11]])
NO <- as.ts(d[[12]])
NOx <- as.ts(d[[13]])
NOy <- as.ts(d[[14]])
NOz <- as.ts(d[[18]])
Level <- as.ts(d[[15]])
Level <- round(Level,d=0)		### make sure we have integers!
NO2 <- as.ts(d[[17]])
Lflag <- as.ts(d[[16]])

Fdoy <- as.ts(Flags[[6]])
#Fttag <- as.ts(Flags[[2]])
#Flevel <- as.ts(Flags[[3]])
Fso2 <- as.ts(Flags[[8]])
Fo3 <- as.ts(Flags[[9]])
Fco2 <- as.ts(Flags[[13]])
Fh2o <- as.ts(Flags[[14]])
FNO <- as.ts(Flags[[10]])
FNOx <- as.ts(Flags[[11]])
FNOy <- as.ts(Flags[[12]])
########## Flag Codes:
## 1 == Okay, sampling
## 0 == bad data (for whatever reason// including calibration)

### make an hrmin vector
hrmin <- hour/100 + min/60

cat("\n check on input data")
browser()
### make sure that Fday and day are the same size
Lday <- length(index)
Lfday <- length(Fdoy)

if(MONTH == "feb")	{
	o3.test <- o3 < 15
	o3 <- ifelse(o3.test == T,NA,o3)
			}


### for 2009 data - should all be 5 min. avg. data!!!!
##############################################################################
#########33 now should have 5 min. time series - avg values for each level!!!!
######### combine into single matrix
M <-  cbind(year,day,hour,min,doy,so2,o3,co2,h2o,NO,NOx,NOy,NO2,NOz,Level, hrmin)
##############################################################################
##############################################################################

## now sort by levels ##############
### get set of indices for each level
Level <- M[,15]
L1 <- Level == 1
L2 <- Level == 2
L3 <- Level == 3
L4 <- Level == 4
L5 <- Level == 5
L6 <- Level == 6

### test - make sure we have complete profiles
cat("\n Level 1:",sum(L1,na.rm=T))
cat("\n Level 2:",sum(L2,na.rm=T))
cat("\n Level 3:",sum(L3,na.rm=T))
cat("\n Level 4:",sum(L4,na.rm=T))
cat("\n Level 5:",sum(L5,na.rm=T))
cat("\n Level 6:",sum(L6,na.rm=T))
cat("Number of samples in each level should be the same!")


########## Not sure why we need complete profiles here ##############
########## May need them eventually for Inv. Langragian Method ######

#######################################################################
###M <-  cbind(year,day,hour,min,doy,so2,o3,co2,h2o,NO,NOx,NOy,NO2,NOz,Level,hrmin)
##############################################################################
### time tags first:
### just use the day of year here. 
t.1 <- M[L1,5]
t.2 <- M[L2,5]
t.3 <- M[L3,5]
t.4 <- M[L4,5]
t.5 <- M[L5,5]
t.6 <- M[L6,5]
t.L <- cbind(t.1,t.2,t.3,t.4,t.5,t.6)

### ozone 
o3.1 <- M[L1,7]
o3.2 <- M[L2,7]
o3.3 <- M[L3,7]
o3.4 <- M[L4,7]
o3.5 <- M[L5,7]
o3.6 <- M[L6,7]
### NO
NO.1 <- M[L1,10]
NO.2 <- M[L2,10]
NO.3 <- M[L3,10]
NO.4 <- M[L4,10]
NO.5 <- M[L5,10]
NO.6 <- M[L6,10]
### SO2:
so2.1 <- M[L1,6]
so2.2 <- M[L2,6]
so2.3 <- M[L3,6]
so2.4 <- M[L4,6]
so2.5 <- M[L5,6]
so2.6 <- M[L6,6]
### NOx
NOx.1 <- M[L1,11]
NOx.2 <- M[L2,11]
NOx.3 <- M[L3,11]
NOx.4 <- M[L4,11]
NOx.5 <- M[L5,11]
NOx.6 <- M[L6,11]
### NOy
NOy.1 <- M[L1,12]
NOy.2 <- M[L2,12]
NOy.3 <- M[L3,12]
NOy.4 <- M[L4,12]
NOy.5 <- M[L5,12]
NOy.6 <- M[L6,12]
### NOz
NOz.1 <- M[L1,14]
NOz.2 <- M[L2,14]
NOz.3 <- M[L3,14]
NOz.4 <- M[L4,14]
NOz.5 <- M[L5,14]
NOz.6 <- M[L6,14]
### CO2
co2.1 <- M[L1,8]
co2.2 <- M[L2,8]
co2.3 <- M[L3,8]
co2.4 <- M[L4,8]
co2.5 <- M[L5,8]
co2.6 <- M[L6,8]
### h2o
h2o.1 <- M[L1,9]
h2o.2 <- M[L2,9]
h2o.3 <- M[L3,9]
h2o.4 <- M[L4,9]
h2o.5 <- M[L5,9]
h2o.6 <- M[L6,9]
### NO2
NO2.1 <- M[L1,13]
NO2.2 <- M[L2,13]
NO2.3 <- M[L3,13]
NO2.4 <- M[L4,13]
NO2.5 <- M[L5,13]
NO2.6 <- M[L6,13]

cat("\n check out consistency of levels")
browser()
### now have divided by level

################################################################
########### Data by Level (Level files) ########################
### now put into 6 matrices = each level
### Data should already be filtered (from MEFtracegas.R routine)
M1 <- cbind(t.1,so2.1,o3.1,co2.1,h2o.1,NO.1,NOx.1,NOy.1,NOz.1,NO2.1)
M2 <- cbind(t.2,so2.2,o3.2,co2.2,h2o.2,NO.2,NOx.2,NOy.2,NOz.2,NO2.2)
M3 <- cbind(t.3, so2.3,o3.3,co2.3,h2o.3, NO.3,NOx.3,NOy.3,NOz.3,NO2.3)
M4 <- cbind(t.4,so2.4,o3.4,co2.4,h2o.4, NO.4,NOx.4,NOy.4,NOz.4,NO2.4)
M5 <- cbind(t.5,so2.5,o3.5,co2.5,h2o.5,NO.5,NOx.5,NOy.5,NOz.5,NO2.5)
M6 <- cbind(t.6,so2.6,o3.6,co2.6,h2o.6,NO.6,NOx.6,NOy.6,NOz.6,NO2.6)

############## Quick plots #########################################
###### plot levels 1, 3 and 6 ######################################
par(mfrow=c(2,1))
plot(M1[,1],M1[,2],type="l",ylab="SO2")
lines(M3[,1], M3[,2],col=2)
lines(M6[,1], M6[,2],col=3)
plot(M1[,1],M1[,3],type="l",ylab="O3")
lines(M3[,1], M3[,3],col=2)
lines(M6[,1], M6[,3],col=3)
browser()
plot(M1[,1],M1[,4],type="l",ylab="CO2")
lines(M3[,1], M3[,4],col=2)
lines(M6[,1], M6[,4],col=3)
plot(M1[,1],M1[,5],type="l",ylab="H2O")
lines(M3[,1], M3[,5],col=2)
lines(M6[,1], M6[,5],col=3)
browser()
#### skip NOx for now - they are confusing anyway. ##############
####### Output these files ######################################
fun.asci(M1, out1)
fun.asci(M2, out2)
fun.asci(M3, out3)
fun.asci(M4, out4)
fun.asci(M5, out5)
fun.asci(M6, out6)
####################################################################

### Now get profiles by species ####################################
######### here is where we need to use only complete profiles for each species (== length vectors) ###
######################################################################################################
### now - need to get complete profiles.....1 every 30 min.
#######################################################################
###M <-  cbind(year,day,hour,min,doy,so2,o3,co2,h2o,NO,NOx,NOy,NO2,NOz,Level,hrmin)
cat("\n check out consistency of levels")
browser()
############First - get rid of NAs (sometimes have rows of NAs at end!!)
##### get rid of all NAs in "Level" (don't know where you are in profile!)
test <- is.na(M[,15])
Mnew <- M[test==F,]
M <- Mnew
Lm <- length(M[,1])
## get an index (to figure out where things go bad
ind <- seq(from=1,to=Lm,by=1)
i <- 1
k <- 0
p <- 1
#browser()
while(i < Lm)	{
	### find first (next) where Level = 1
	j <- i + 5
	#cat("\n i=",i,"j=",j)
	if(j < Lm)	{
		 Mi <- M[i:j,]
		if(Mi[1,15] == 1 & Mi[6,15] == 6)	{
			#cat("\n Index:",ind[i],"Good")

			if(k == 0)	Ni <- Mi
			else		Ni <- rbind(Ni,Mi)

			k <- k + 1
			i <- i + 6
				}		## end of if(Mi[i,3] == 1)
		else	{
			i <- i + 1
			cat("\n Incomplete profile - #",p,"Index=",ind[i])
			#browser()
			p <- p + 1
			}
	}	## end of if j < Lm
	else	{
		i <- i + 1
		cat("\n Incomplete profile - #",p,"Index=",ind[i])
		p <- p + 1
		}
	
	}	## end of while loop

cat("Have only complete profiles now!")
#browser()
###Now have filtered for incomplete profiles
### now "M" has only complete profiles
M <- Ni

if(FILTER == T)	{
	## Filter for "extreme" events in SO2 and NO (vehicles passing etc.)
	####### use these two lines to set high cut-offs !!!!
	NOcut <- 15
	SO2cut <- 15
	p1 <- M[,6] > SO2cut
	p2 <- M[,10] > NOcut
	p1 <- sum(p1,na.rm=T)
	p2 <- sum(p2,na.rm=T)
	if(p1 > 0)		cat("\nSee ",p1,"SO2 events")
	if(p2 > 0)		cat("\nSee ",p2,"NO events")
	M[,6] <- ifelse(M[,6] > SO2cut,NA,M[,6])
	M[,10] <- ifelse(M[,10] > SO2cut,NA,M[,10])
			}

##### parse by Level again.
Level <- M[,15]
L1 <- Level == 1
L2 <- Level == 2
L3 <- Level == 3
L4 <- Level == 4
L5 <- Level == 5
L6 <- Level == 6
### time tags first:
### just use the day of year here. 
t.1 <- M[L1,5]
t.2 <- M[L2,5]
t.3 <- M[L3,5]
t.4 <- M[L4,5]
t.5 <- M[L5,5]
t.6 <- M[L6,5]
t.L <- cbind(t.1,t.2,t.3,t.4,t.5,t.6)
### ozone 
o3.1 <- M[L1,7]
o3.2 <- M[L2,7]
o3.3 <- M[L3,7]
o3.4 <- M[L4,7]
o3.5 <- M[L5,7]
o3.6 <- M[L6,7]
### NO
NO.1 <- M[L1,10]
NO.2 <- M[L2,10]
NO.3 <- M[L3,10]
NO.4 <- M[L4,10]
NO.5 <- M[L5,10]
NO.6 <- M[L6,10]
### SO2:
so2.1 <- M[L1,6]
so2.2 <- M[L2,6]
so2.3 <- M[L3,6]
so2.4 <- M[L4,6]
so2.5 <- M[L5,6]
so2.6 <- M[L6,6]
### NOx
NOx.1 <- M[L1,11]
NOx.2 <- M[L2,11]
NOx.3 <- M[L3,11]
NOx.4 <- M[L4,11]
NOx.5 <- M[L5,11]
NOx.6 <- M[L6,11]
### NOy
NOy.1 <- M[L1,12]
NOy.2 <- M[L2,12]
NOy.3 <- M[L3,12]
NOy.4 <- M[L4,12]
NOy.5 <- M[L5,12]
NOy.6 <- M[L6,12]
### NOz
NOz.1 <- M[L1,14]
NOz.2 <- M[L2,14]
NOz.3 <- M[L3,14]
NOz.4 <- M[L4,14]
NOz.5 <- M[L5,14]
NOz.6 <- M[L6,14]
### CO2
co2.1 <- M[L1,8]
co2.2 <- M[L2,8]
co2.3 <- M[L3,8]
co2.4 <- M[L4,8]
co2.5 <- M[L5,8]
co2.6 <- M[L6,8]
### h2o
h2o.1 <- M[L1,9]
h2o.2 <- M[L2,9]
h2o.3 <- M[L3,9]
h2o.4 <- M[L4,9]
h2o.5 <- M[L5,9]
h2o.6 <- M[L6,9]
### NO2
NO2.1 <- M[L1,13]
NO2.2 <- M[L2,13]
NO2.3 <- M[L3,13]
NO2.4 <- M[L4,13]
NO2.5 <- M[L5,13]
NO2.6 <- M[L6,13]
### hrmin
hm.1 <- M[L1,16]
hm.2 <- M[L2,16]
hm.3 <- M[L3,16]
hm.4 <- M[L4,16]
hm.5 <- M[L5,16]
hm.6 <- M[L6,16]
hm.1 <- round(hm.1,d=1)
hm.1 <- ifelse(hm.1 > 23.75,0.0,hm.1)

### test - make sure we have complete profiles
cat("\n Level 1:",sum(L1,na.rm=T))
cat("\n Level 2:",sum(L2,na.rm=T))
cat("\n Level 3:",sum(L3,na.rm=T))
cat("\n Level 4:",sum(L4,na.rm=T))
cat("\n Level 5:",sum(L5,na.rm=T))
cat("\n Level 6:",sum(L6,na.rm=T))
cat("Number of samples in each level should be the same!")


cat("\n check out consistency of levels")
browser()
### now divide by species
############## Data by species #################################
##### use the time tag at the end of the 30 min. period. #######
S.o3 <- cbind(t.1,hm.1,o3.1,o3.2,o3.3,o3.4,o3.5,o3.6)
S.so2 <- cbind(t.1,hm.1,so2.1,so2.2,so2.3,so2.4,so2.5,so2.6)
S.co2 <- cbind(t.1,hm.1,co2.1,co2.2,co2.3,co2.4,co2.5,co2.6)
S.h2o <- cbind(t.1,hm.1,h2o.1,h2o.2,h2o.3,h2o.4,h2o.5,h2o.6)
S.NO <- cbind(t.1,hm.1,NO.1,NO.2,NO.3,NO.4,NO.5,NO.6)
S.NOx <- cbind(t.1,hm.1,NOx.1,NOx.2,NOx.3,NOx.4,NOx.5,NOx.6)
S.NOy <- cbind(t.1,hm.1,NOy.1,NOy.2,NOy.3,NOy.4,NOy.5,NOy.6)
S.NOz <- cbind(t.1,hm.1,NOz.1,NOz.2,NOz.3,NOz.4,NOz.5,NOz.6)
S.NO2 <- cbind(t.1,hm.1,NO2.1,NO2.2,NO2.3,NO2.4,NO2.5,NO2.6)

####### make sure there is data for each species (use level 6)
p.o3 <- sum(is.na(S.o3[,8]),na.rm=T)
if(p.o3 == length(S.o3[,8]))		P.o3 <- F

p.so2 <- sum(is.na(S.so2[,8]),na.rm=T)
if(p.so2 == length(S.so2[,8]))		P.so2 <- F

p.no <- sum(is.na(S.NO[,8]),na.rm=T)
if(p.no == length(S.NO[,8]))		P.no <- F

p.co2 <- sum(is.na(S.co2[,8]),na.rm=T)
if(p.co2 == length(S.co2[,8]))		P.co2 <- F

p.h2o <- sum(is.na(S.h2o[,8]),na.rm=T)
if(p.h2o == length(S.h2o[,8]))		P.h2o <- F

p.nox <- sum(is.na(S.NOx[,8]),na.rm=T)
if(p.nox == length(S.NOx[,8]))		P.nox <- FALSE	

p.noy <- sum(is.na(S.NOy[,8]),na.rm=T)
if(p.noy == length(S.NOy[,8]))		P.noy <- F

Dflags <- c(P.o3,P.so2,P.co2,P.h2o,P.no,P.nox,P.noy)
cat("\n Data Flags:",Dflags)
browser()
par(mfrow=c(1,1))
if(P.o3 == T)	{
o3.lim <- range(S.o3[,3:8],na.rm=T)
plot(S.o3[,8],S.o3[,7],pch=1,ylim=o3.lim,xlim=o3.lim,xlab="O3 at 25m")
points(S.o3[,8],S.o3[,6],col=2)
points(S.o3[,8],S.o3[,5],col=3)
points(S.o3[,8],S.o3[,4],col=4)
points(S.o3[,8],S.o3[,3],col=5)
cat("\n Colors: Black=L5,Red=L4,Green=L3,Blue=L2,Cyan=L1")
browser()
		}
if(P.so2 == T)	{
so2.lim <- range(S.so2[,3:8],na.rm=T)
plot(S.so2[,8],S.so2[,7],pch=1,ylim=so2.lim,xlim=so2.lim,xlab="so2 at 25m")
points(S.so2[,8],S.so2[,6],col=2)
points(S.so2[,8],S.so2[,5],col=3)
points(S.so2[,8],S.so2[,4],col=4)
points(S.so2[,8],S.so2[,3],col=5)
cat("\n Colors: Black=L5,Red=L4,Green=L3,Blue=L2,Cyan=L1")
browser()
		}
if(P.co2 == T)	{
co2.lim <- range(S.co2[,3:8],na.rm=T)
plot(S.co2[,8],S.co2[,7],pch=1,ylim=co2.lim,xlim=co2.lim,xlab="co2 at 25m")
points(S.co2[,8],S.co2[,6],col=2)
points(S.co2[,8],S.co2[,5],col=3)
points(S.co2[,8],S.co2[,4],col=4)
points(S.co2[,8],S.co2[,3],col=5)
cat("\n Colors: Black=L5,Red=L4,Green=L3,Blue=L2,Cyan=L1")
browser()
		}
if(P.h2o == T)	{
h2o.lim <- range(S.h2o[,3:8],na.rm=T)
plot(S.h2o[,8],S.h2o[,7],pch=1,ylim=h2o.lim,xlim=h2o.lim,xlab="h2o at 25m")
points(S.h2o[,8],S.h2o[,6],col=2)
points(S.h2o[,8],S.h2o[,5],col=3)
points(S.h2o[,8],S.h2o[,4],col=4)
points(S.h2o[,8],S.h2o[,3],col=5)
cat("\n Colors: Black=L5,Red=L4,Green=L3,Blue=L2,Cyan=L1")
browser()
			}
if(P.no == T)	{
NO.lim <- range(S.NO[,3:8],na.rm=T)
plot(S.NO[,8],S.NO[,7],pch=1,ylim=NO.lim,xlim=NO.lim,xlab="NO at 25m")
points(S.NO[,8],S.NO[,6],col=2)
points(S.NO[,8],S.NO[,5],col=3)
points(S.NO[,8],S.NO[,4],col=4)
points(S.NO[,8],S.NO[,3],col=5)
cat("\n Colors: Black=L5,Red=L4,Green=L3,Blue=L2,Cyan=L1")
browser()
			}
if(P.nox == T)	{
NOx.lim <- range(S.NOx[,3:8],na.rm=T)
plot(S.NOx[,8],S.NOx[,7],pch=1,ylim=NOx.lim,xlim=NOx.lim,xlab="NOx at 25m")
points(S.NOx[,8],S.NOx[,6],col=2)
points(S.NOx[,8],S.NOx[,5],col=3)
points(S.NOx[,8],S.NOx[,4],col=4)
points(S.NOx[,8],S.NOx[,3],col=5)
cat("\n Colors: Black=L5,Red=L4,Green=L3,Blue=L2,Cyan=L1")
browser()
			}
if(P.noy == T)	{
NOy.lim <- range(S.NOy[,3:8],na.rm=T)
plot(S.NOy[,8],S.NOy[,7],pch=1,ylim=NOy.lim,xlim=NOy.lim,xlab="NOy at 25m")
points(S.NOy[,8],S.NOy[,6],col=2)
points(S.NOy[,8],S.NOy[,5],col=3)
points(S.NOy[,8],S.NOy[,4],col=4)
points(S.NOy[,8],S.NOy[,3],col=5)
cat("\n Colors: Black=L5,Red=L4,Green=L3,Blue=L2,Cyan=L1")
browser()
			}
cat("\n end of plots to check on data ")



###############################################################################
## Make some plots:

if(PLOT==T) {
### set some early plotting limits:
		#x1 <- range(t.L)
		#x.co2 <- c(135.5,max(t.L,na.rm=T))
		##3 ylimits:
		y.o3 <- c(0,25)
		#y.NO <- c(0,10)
		#y.NOx <- c(0,22)
		#y.NOy <- c(0,22)
		#y.co2 <- c(360,480)
		#y.h2o <- c(7,13.5)

		#par(mfrow=c(2,1))
		#browser()

o3.L <- cbind(o3.1,o3.2,o3.3,o3.4,o3.5,o3.6)
so2.L <- cbind(so2.1,so2.2,so2.3,so2.4,so2.5,so2.6)
NO.L <- cbind(NO.1,NO.2,NO.2,NO.4,NO.5,NO.6)
NOx.L <- cbind(NOx.1,NOx.2,NOx.2,NOx.4,NOx.5,NOx.6)
NOy.L <- cbind(NOy.1,NOy.2,NOy.2,NOy.4,NOy.5,NOy.6)
co2.L <- cbind(co2.1,co2.2,co2.3,co2.4,co2.5,co2.6)
h2o.L <- cbind(h2o.1,h2o.2,h2o.3,h2o.4,h2o.5,h2o.6)

fun.fastplot(t.L,o3.L,Label="O3")
	browser()
fun.fastplot(t.L,so2.L,Label="SO2")
	browser()
fun.fastplot(t.L,co2.L,Label="CO2")
	browser()
fun.fastplot(t.L,h2o.L,Label="H2O")
	browser()
fun.fastplot(t.L,NO.L,Label="NO")
	browser()
#fun.fastplot(t.L,NOx.L,Label="NOx")
	browser()
fun.fastplot(t.L,NOy.L,Label="NOy")
browser()
#fun.fastplot(t.L,NO2.L,Label="NO2")
cat("\n Last plot - will then output files")
			} ### end of if(PLOT==T)

#####################################################################################
fun.asci(S.o3, o3prof)
fun.asci(S.so2, so2prof)
fun.asci(S.co2, co2prof)
fun.asci(S.h2o, h2oprof)
fun.asci(S.NO, NOprof)
fun.asci(S.NOx, NOxprof)
fun.asci(S.NOy, NOyprof)
fun.asci(S.NO2, NO2prof)

####################
########################## OK to here !!!!! #######################

cat("\n check on diurnal avg. time tags ")
######## get diurnal average for each level ##########################
#########################################################################
DA.o3 <- fun.diurnalavg(S.o3,month)
DA.so2 <- fun.diurnalavg(S.so2,month)
DA.co2 <- fun.diurnalavg(S.co2,month)
DA.h2o <- fun.diurnalavg(S.h2o,month)
DA.NO <- fun.diurnalavg(S.NO,month)
DA.NOx <- fun.diurnalavg(S.NOx,month)
DA.NOy <- fun.diurnalavg(S.NOy,month)
DA.NOz <- fun.diurnalavg(S.NOz,month)
DA.NO2 <- fun.diurnalavg(S.NO2,month)

par(mfrow=c(1,1))
### plot diurnal avgs.
if(P.o3 == T)	{
o3.lim <- range(DA.o3[,8:13],na.rm=T)
plot(DA.o3[,1],DA.o3[,8],type="l",ylab="Median O3",ylim=o3.lim)
lines(DA.o3[,1],DA.o3[,9],col=2)
lines(DA.o3[,1],DA.o3[,10],col=3)
lines(DA.o3[,1],DA.o3[,11],col=4)
lines(DA.o3[,1],DA.o3[,12],col=5)
lines(DA.o3[,1],DA.o3[,13],col=6)
cat("\n L1 = black,L2=red, L3=green, L4=blue, L5 = cyan, L6=Magenta")
browser()
			}
if(P.so2 == T)	{
so2.lim <- range(DA.so2[,8:13],na.rm=T)
plot(DA.so2[,1],DA.so2[,8],type="l",ylab="Median SO2",ylim=so2.lim)
lines(DA.so2[,1],DA.so2[,9],col=2)
lines(DA.so2[,1],DA.so2[,10],col=3)
lines(DA.so2[,1],DA.so2[,11],col=4)
lines(DA.so2[,1],DA.so2[,12],col=5)
lines(DA.so2[,1],DA.so2[,13],col=6)
browser()
			}
if(P.co2 == T)	{
co2.lim <- range(DA.co2[,8:13],na.rm=T)
plot(DA.co2[,1],DA.co2[,8],type="l",ylab="Median CO2",ylim=co2.lim)
lines(DA.co2[,1],DA.co2[,9],col=2)
lines(DA.co2[,1],DA.co2[,10],col=3)
lines(DA.co2[,1],DA.co2[,11],col=4)
lines(DA.co2[,1],DA.co2[,12],col=5)
lines(DA.co2[,1],DA.co2[,13],col=6)
browser()
			}
if(P.h2o == T)	{
h2o.lim <- range(DA.h2o[,8:13],na.rm=T)
plot(DA.h2o[,1],DA.h2o[,8],type="l",ylab="Median H2O",ylim=h2o.lim)
lines(DA.h2o[,1],DA.h2o[,9],col=2)
lines(DA.h2o[,1],DA.h2o[,10],col=3)
lines(DA.h2o[,1],DA.h2o[,11],col=4)
lines(DA.h2o[,1],DA.h2o[,12],col=5)
lines(DA.h2o[,1],DA.h2o[,13],col=6)
browser()
			}
if(P.no == T)	{
NO.lim <- range(DA.NO[,8:13],na.rm=T)
plot(DA.NO[,1],DA.NO[,8],type="l",ylab="Median NO",ylim=NO.lim)
lines(DA.NO[,1],DA.NO[,9],col=2)
lines(DA.NO[,1],DA.NO[,10],col=3)
lines(DA.NO[,1],DA.NO[,11],col=4)
lines(DA.NO[,1],DA.NO[,12],col=5)
lines(DA.NO[,1],DA.NO[,13],col=6)
browser()
			}
if(P.nox == T)	{
NOx.lim <- range(DA.NOx[,8:13],na.rm=T)
plot(DA.NOx[,1],DA.NOx[,8],type="l",ylab="Median NOx",ylim=NOx.lim)
lines(DA.NOx[,1],DA.NOx[,9],col=2)
lines(DA.NOx[,1],DA.NOx[,10],col=3)
lines(DA.NOx[,1],DA.NOx[,11],col=4)
lines(DA.NOx[,1],DA.NOx[,12],col=5)
lines(DA.NOx[,1],DA.NOx[,13],col=6)
browser()
			}
if(P.noy == T)	{
NOy.lim <- range(DA.NOy[,8:13],na.rm=T)
plot(DA.NOy[,1],DA.NOy[,8],type="l",ylab="Median NOy",ylim=NOy.lim)
lines(DA.NOy[,1],DA.NOy[,9],col=2)
lines(DA.NOy[,1],DA.NOy[,10],col=3)
lines(DA.NOy[,1],DA.NOy[,11],col=4)
lines(DA.NOy[,1],DA.NOy[,12],col=5)
lines(DA.NOy[,1],DA.NOy[,13],col=6)
			}
cat("\n End of Diel Plots")
browser()

### Output the mean profiles (DA.xx files)
fun.asci(DA.o3, o3diel)
fun.asci(DA.so2, so2diel)
fun.asci(DA.co2, co2diel)
fun.asci(DA.h2o, h2odiel)
fun.asci(DA.NO, NOdiel)
fun.asci(DA.NOx, NOxdiel)
fun.asci(DA.NOy, NOydiel)
fun.asci(DA.NO2, NO2diel)



#fun.asci(Profiles,"c:/data/AMAZE08/O3Profiles.txt")


####### get average profile at each 1/2 hr from the diurnal averages above #############
########################################################################################

### measurement heights
z.heights <- c(2.0,3.5,7,10,15.5,25.3)			### check this with tape measure!!! 
		### In my old book that i couldn't find -

### transpose the DA matrices
tDA.o3 <- t(DA.o3[,8:13])		## get the means only
tDA.so2 <- t(DA.so2[,8:13])
tDA.co2 <- t(DA.co2[,8:13])
tDA.h2o <- t(DA.h2o[,8:13])
tDA.NO <- t(DA.NO[,8:13])
tDA.NOx <- t(DA.NOx[,8:13])
tDA.NOy <- t(DA.NOy[,8:13])
#### combine with heights
tDA.o3 <- cbind(z.heights,tDA.o3)		## get the means only
tDA.so2 <- cbind(z.heights,tDA.so2)
tDA.co2 <- cbind(z.heights,tDA.co2)
tDA.h2o <- cbind(z.heights,tDA.h2o)
tDA.NO <- cbind(z.heights,tDA.NO)
tDA.NOx <- cbind(z.heights,tDA.NOx)
tDA.NOy <- cbind(z.heights,tDA.NOy)
### should be matrices with 6 rows (each level), col. of heights, each 1/2 hour of the day (ncol = 49)

#### what should we do with this data ????????????
##################################################################################################


### Can probably leave this going 
######################################################################

###################################### calculate and plot some average gradients:
### find day and night times
### make up index time series:
cat("\n starting to get avg. profiles")
#browser()
### look at specific species:
### However, all time tags are the same - can just use one to get day/nights flags.
tt.prof <- S.so2[,1]
### time tags first:
tt.day <- floor(tt.prof)
tt.index <- tt.prof - tt.day

morning <- tt.index > 0.2 & tt.index < 0.34			## denote: M
daytime <- tt.index > 0.36 & tt.index < 0.63			## denote: D
evening <- tt.index > 0.67 & tt.index < 0.84			## denote: E
night <- tt.index > 0.84 | tt.index < 0.17			## denote: N
cat("\n Test indexing for day/night")
#browser()

## (1) SO2
########### Daytime Gradient:
M.so2 <- fun.Pskewness(S.so2[morning,],NAME="SO2 morning")
D.so2 <- fun.Pskewness(S.so2[daytime,],NAME="SO2 daytime")
E.so2 <- fun.Pskewness(S.so2[evening,],NAME="SO2 evening")
N.so2 <- fun.Pskewness(S.so2[night,],NAME="SO2 night")
## (2) O3
########### Daytime Gradient:
M.o3 <- fun.Pskewness(S.o3[morning,],NAME="O3 morning")
D.o3 <- fun.Pskewness(S.o3[daytime,],NAME="O3 daytime")
E.o3 <- fun.Pskewness(S.o3[evening,],NAME="O3 evening")
N.o3 <- fun.Pskewness(S.o3[night,],NAME="O3 night")
## (3) CO2
########### Daytime Gradient:
M.co2 <- fun.Pskewness(S.co2[morning,],NAME="co2 morning")
D.co2 <- fun.Pskewness(S.co2[daytime,],NAME="co2 daytime")
E.co2 <- fun.Pskewness(S.co2[evening,],NAME="co2 evening")
N.co2 <- fun.Pskewness(S.co2[night,],NAME="co2 night")
## (4) H2O
########### Daytime Gradient:
M.h2o <- fun.Pskewness(S.h2o[morning,],NAME="h2o morning")
D.h2o <- fun.Pskewness(S.h2o[daytime,],NAME="h2o daytime")
E.h2o <- fun.Pskewness(S.h2o[evening,],NAME="h2o evening")
N.h2o <- fun.Pskewness(S.h2o[night,],NAME="h2o night")
## (5) NO
########### Daytime Gradient:
M.NO <- fun.Pskewness(S.NO[morning,],NAME="NO morning")
D.NO <- fun.Pskewness(S.NO[daytime,],NAME="NO daytime")
E.NO <- fun.Pskewness(S.NO[evening,],NAME="NO evening")
N.NO <- fun.Pskewness(S.NO[night,],NAME="NO night")
## (6) NOx
########### Daytime Gradient:
M.NOx <- fun.Pskewness(S.NOx[morning,],NAME="NOx morning")
D.NOx <- fun.Pskewness(S.NOx[daytime,],NAME="NOx daytime")
E.NOx <- fun.Pskewness(S.NOx[evening,],NAME="NOx evening")
N.NOx <- fun.Pskewness(S.NOx[night,],NAME="NOx night")
## (7) NOy
########### Daytime Gradient:
M.NOy <- fun.Pskewness(S.NOy[morning,],NAME="NOy morning")
D.NOy <- fun.Pskewness(S.NOy[daytime,],NAME="NOy daytime")
E.NOy <- fun.Pskewness(S.NOy[evening,],NAME="NOy evening")
N.NOy <- fun.Pskewness(S.NOy[night,],NAME="NOy night")


DayProf.mean <- rbind(D.so2$xbar,D.o3$xbar,D.co2$xbar,D.h2o$xbar,D.NO$xbar,D.NOx$xbar,D.NOy$xbar)
DayProf.med <- rbind(D.so2$xmed,D.o3$xmed,D.co2$xmed,D.h2o$xmed,D.NO$xmed,D.NOx$xmed,D.NOy$xmed)
DayProf.std <- rbind(D.so2$xstd,D.o3$xstd,D.co2$xstd,D.h2o$xstd,D.NO$xstd,D.NOx$xstd,D.NOy$xstd)
DayProf.N <- rbind(D.so2$N,D.o3$N,D.co2$N,D.h2o$N,D.NO$N,D.NOx$N,D.NOy$N)

MornProf.mean <- rbind(M.so2$xbar,M.o3$xbar,M.co2$xbar,M.h2o$xbar,M.NO$xbar,M.NOx$xbar,M.NOy$xbar)
MornProf.med <- rbind(M.so2$xmed,M.o3$xmed,M.co2$xmed,M.h2o$xmed,M.NO$xmed,M.NOx$xmed,M.NOy$xmed)
MornProf.std <- rbind(M.so2$xstd,M.o3$xstd,M.co2$xstd,M.h2o$xstd,M.NO$xstd,M.NOx$xstd,M.NOy$xstd)
MornProf.N <- rbind(M.so2$N,M.o3$N,M.co2$N,M.h2o$N,M.NO$N,M.NOx$N,M.NOy$N)

EveProf.mean <- rbind(E.so2$xbar,E.o3$xbar,E.co2$xbar,E.h2o$xbar,E.NO$xbar,E.NOx$xbar,E.NOy$xbar)
EveProf.med <- rbind(E.so2$xmed,E.o3$xmed,E.co2$xmed,E.h2o$xmed,E.NO$xmed,E.NOx$xmed,E.NOy$xmed)
EveProf.std <- rbind(E.so2$xstd,E.o3$xstd,E.co2$xstd,E.h2o$xstd,E.NO$xstd,E.NOx$xstd,E.NOy$xstd)
EveProf.N <- rbind(E.so2$N,E.o3$N,E.co2$N,E.h2o$N,E.NO$N,E.NOx$N,E.NOy$N)

NightProf.mean <- rbind(N.so2$xbar,N.o3$xbar,N.co2$xbar,N.h2o$xbar,N.NO$xbar,N.NOx$xbar,N.NOy$xbar)
NightProf.med <- rbind(N.so2$xmed,N.o3$xmed,N.co2$xmed,N.h2o$xmed,N.NO$xmed,N.NOx$xmed,N.NOy$xmed)
NightProf.std <- rbind(N.so2$xstd,N.o3$xstd,N.co2$xstd,N.h2o$xstd,N.NO$xstd,N.NOx$xstd,N.NOy$xstd)
NightProf.N <- rbind(N.so2$N,N.o3$N,N.co2$N,N.h2o$N,N.NO$N,N.NOx$N,N.NOy$N)

cat("\n final avg profile plots")
browser()
par(mfrow=c(2,2))
plot(NightProf.mean[2,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Night O3, ppbv")
points(NightProf.med[2,2:7],z.heights,type="o",col=2,pch=5)
plot(MornProf.mean[2,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Morning O3, ppbv")
points(MornProf.med[2,2:7],z.heights,type="o",pch=5,col=2)
plot(DayProf.mean[2,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Midday O3, ppbv")
points(DayProf.med[2,2:7],z.heights,type="o",pch=5,col=2)
plot(EveProf.mean[2,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Evening O3, ppbv")
points(EveProf.med[2,2:7],z.heights,type="o",pch=5,col=2)
browser()

plot(NightProf.mean[1,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Night SO2, ppbv")
points(NightProf.med[1,2:7],z.heights,type="o",col=2,pch=5)
plot(MornProf.mean[1,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Morning SO2, ppbv")
points(MornProf.med[1,2:7],z.heights,type="o",pch=5,col=2)
plot(DayProf.mean[1,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Midday SO2, ppbv")
points(DayProf.med[1,2:7],z.heights,type="o",pch=5,col=2)
plot(EveProf.mean[1,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Evening SO2, ppbv")
points(EveProf.med[1,2:7],z.heights,type="o",pch=5,col=2)
browser()

plot(NightProf.mean[3,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Night CO2, ppbv")
points(NightProf.med[3,2:7],z.heights,type="o",col=2,pch=5)
plot(MornProf.mean[3,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Morning CO2, ppbv")
points(MornProf.med[3,2:7],z.heights,type="o",pch=5,col=2)
plot(DayProf.mean[3,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Midday CO2, ppbv")
points(DayProf.med[3,2:7],z.heights,type="o",pch=5,col=2)
plot(EveProf.mean[3,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Evening CO2, ppbv")
points(EveProf.med[3,2:7],z.heights,type="o",pch=5,col=2)
browser()

plot(NightProf.mean[4,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Night H2O, ppbv")
points(NightProf.med[4,2:7],z.heights,type="o",col=2,pch=5)
plot(MornProf.mean[4,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Morning H2O, ppbv")
points(MornProf.med[4,2:7],z.heights,type="o",pch=5,col=2)
plot(DayProf.mean[4,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Midday H2O, ppbv")
points(DayProf.med[4,2:7],z.heights,type="o",pch=5,col=2)
plot(EveProf.mean[4,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Evening H2O, ppbv")
points(EveProf.med[4,2:7],z.heights,type="o",pch=5,col=2)
browser()

plot(NightProf.mean[5,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Night NO, ppbv")
points(NightProf.med[5,2:7],z.heights,type="o",col=2,pch=5)
plot(MornProf.mean[5,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Morning NO, ppbv")
points(MornProf.med[5,2:7],z.heights,type="o",pch=5,col=2)
plot(DayProf.mean[5,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Midday NO, ppbv")
points(DayProf.med[5,2:7],z.heights,type="o",pch=5,col=2)
plot(EveProf.mean[5,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Evening NO, ppbv")
points(EveProf.med[5,2:7],z.heights,type="o",pch=5,col=2)
browser()

#plot(NightProf.mean[6,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Night NOx, ppbv")
#points(NightProf.med[6,2:7],z.heights,type="o",col=2,pch=5)
#plot(MornProf.mean[6,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Morning NOx, ppbv")
#points(MornProf.med[6,2:7],z.heights,type="o",pch=5,col=2)
#plot(DayProf.mean[6,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Midday NOx, ppbv")
#points(DayProf.med[6,2:7],z.heights,type="o",pch=5,col=2)
#plot(EveProf.mean[6,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Evening NOx, ppbv")
#points(EveProf.med[6,2:7],z.heights,type="o",pch=5,col=2)
#browser()

plot(NightProf.mean[7,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Night NOy, ppbv")
points(NightProf.med[7,2:7],z.heights,type="o",col=2,pch=5)
plot(MornProf.mean[7,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Morning NOy, ppbv")
points(MornProf.med[7,2:7],z.heights,type="o",pch=5,col=2)
plot(DayProf.mean[7,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Midday NOy, ppbv")
points(DayProf.med[7,2:7],z.heights,type="o",pch=5,col=2)
plot(EveProf.mean[7,2:7],z.heights,type="o",pch=1,ylab="z,m",xlab="Evening NOy, ppbv")
points(EveProf.med[7,2:7],z.heights,type="o",pch=5,col=2)

cat("\n End of function")
browser()



}

##############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.Pskewness <- function(x,NAME)
	{
## function designed to take a matrix x and 
## calculate the mean, median, std. dev., skewness and
## kurtosis for each column.
## x = matrix:  time tag, Conc. at L1, L2,L3,L4,L5,L6
# number of columns (should be 7, ttag + 6 levels):

Lc <- length(x[1,])

### find  NA's
p <- is.na(x[,2:7])
p1 <- sum(p[,1],na.rm=T)
p2 <- sum(p[,2],na.rm=T)
p3 <- sum(p[,3],na.rm=T)
p4 <- sum(p[,4],na.rm=T)
p5 <- sum(p[,5],na.rm=T)
p6 <- sum(p[,6],na.rm=T)
N1 <- length(x[,1]) - p1
N2 <- length(x[,1]) - p2
N3 <- length(x[,1]) - p3
N4 <- length(x[,1]) - p4
N5 <- length(x[,1]) - p5
N6 <- length(x[,1]) - p6

#if(NAME=="SO2 morning")	browser()
N <-  cbind(N1,N2,N3,N4,N5,N6)
xbar <- apply(x,2,mean,na.rm=T)
xmed <- apply(x,2,median,na.rm=T)
x.min <- apply(x,2,min,na.rm=T)
x.max <- apply(x,2,max,na.rm=T)

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

Residuals <- x - xbar		### time series (matrix??)

Skx <- (apply(Residuals, 2, mean, na.rm=T)^3) / (sdx^3)

Krx <- (apply(Residuals, 2, mean, na.rm=T)^4) / (sdx^4)

	cat("\n Species = ",NAME)
	cat("\n N = ",N)
	cat("\n xbar = ",signif(xbar,d=3))
	cat("\n xmed = ",signif(xmed,d=3))
	cat("\n sdx = ",signif(sdx,d=3))

	
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.diurnalavg <- function(X,month)
	{
## Subroutine to get diurnal avg from time series of profiles for a single species
#  X = time series with columns:  doy, hrmin, Levels 1 - 6  (8 columns)

i <- 0.0
j <- 1
Y <- matrix(NA,ncol=31,nrow=48)
hrmin <- X[,2]
#browser()

while(i <= 23.5)	{
	#if(i == 7.0)	browser()
	test <- hrmin == i
	### should index all the points at 0.0 hours
	Xtest <- X[test,]

	## for May, 2009 - calibrations are off by 1 hour for a while - get rid of 
	if(month == 5)	{
		if(i == 0)	Xtest[,5] <- NA
		if(i == 4)	Xtest[,5] <- NA
		if(i == 8)	Xtest[,5] <- NA
		if(i == 12)	Xtest[,5] <- NA
		if(i == 16)	Xtest[,5] <- NA
		if(i == 20)	Xtest[,5] <- NA
				}
	### get means, medians and std. dev.and # of points (minus the NAs)
	P <- is.na(Xtest[,3:8])
	P[,1:6] <- ifelse(P[,1:6] == 1,0,1)		### switch around na matrix.
	N <- apply(P,2,sum,na.rm=T)		### number of pts.
	xbar <- apply(Xtest[,3:8], 2,mean,na.rm=T)
	xmed <- apply(Xtest[,3:8],2,median,na.rm=T)
	xvar <- apply(Xtest[,3:8],2,var,na.rm=T)
	xstd <- (xvar)^0.5
	xsterr <- xstd/(N^.5)
	############################
	### put into a row of output matrix 
	Y[j,] <- c(i,N,xbar,xmed,xstd,xsterr)
	###########################
	## now move to the next 1/2 hour
	i <- i + 0.5
	j <- j + 1
	}		## end of the while loop

## output Y
Y

}

# end function
##########################################################################