# R program for processing up to one day of high rate (0.5 - 2 Hz) AIRCOA output data to produce ~2.5-minute data and diagnostics 
#
# Copyright 2006 University Corporation for Atmospheric Research
#
# Britton Stephens, National Center for Atmospheric Research, 1850 Table Mesa Dr., Boulder, CO 80305, USA, (303)497-1018, stephens@ucar.edu 
#

# define function name, options, and default values:
acrawday<-function(loc='NWR',unit='A',day='050901',outflag='',njog=150,skip=50,ico=0.0,eco=1.0,trh='humitter',png=F,datadir='/scr/raf/stephens/RCDAT',preserve=F) {
# loc is three letter location code
# unit is one letter AIRCOA unit name
# day is YYMMDD day of output to process
# outflag is output flag used when processing multiple parts of a day separately or the same part in different ways
# njog is the total number of data points for each line or cylinder analysis period (jog)
# skip is the number of data points to skip at the start of each jog
# ntrd is the number of points at start and end of nave period to average and then difference to get slopes (trends) for each jog - must be less than nave = njog - nskip
# ico and eco are initial and ending fractions of the day to allow selection of periods less than one day
# png selects png or bitmap processing for the output graphs (png processing is not available in batch mode, with bitmap processing output is still png files but of lower quality)
month<-substr(day,1,4) # YYMM month of data

print(paste(system('date',T),paste('RAWDAY starting',loc,day)))

# start output and error logging:
log<-file(paste(loc,"_",day,outflag,".log",sep=''),open="wt")
sink(log)
sink(log,type="message")
print(paste('aircoa_raw_day starting',loc,day))
print(system('date',T))

# if current month directory doesn't exist, then make it
#if(length(system(paste('ls -d ',month,sep=''),T))==0){
if(length(list.dirs(month))==0){
        system(paste('mkdir ',month,sep=''))
}

# if loc directory for current month doesn't exist, then make it
#if(length(system(paste('ls -d ',month,'/',loc,sep=''),T))==0){
if(length(list.dirs(paste(month,'/',loc,sep='')))==0){
        system(paste('mkdir ',month,'/',loc,sep=''))
}

# define input file prefix:
file<-paste(loc,'_',day,sep='')

# define number of data points to average at end of jog:
nave<-njog-skip

# define number of data points at start and end of nave period to average and then difference to get slopes (trends) for each jog, <= half of nave = njog - nskip
ntrd<-floor(nave/2)

# define leak check parameters:
nlpp<-njog/5 # 30 for 1 sec data, 15 for 2 sec data
nhpp<-njog/5 # 30 for 1 sec data, 15 for 2 sec data
#nlpp is the number of data points during the initial component of the low pressure leak check
#nhpp is the number of data points during the initial component of the high pressure leak check

#### set system parameters to access binary input files
#### must set RAWDATADIR somewhere (e.g. in autorun.r), example commands:
####  Sys.putenv("RAWDATADIR"="/usr/local/isff/aster/projects/LABTEST/raw_data/all")
####  Sys.putenv("RAWDATADIR"="/scr/isff/projects/CME04/raw_data/all")

# convert letters to numbers and capitals to lowercase to match input file naming convention for early units:
unum<-tolower(chartr('ABCDEFG','1234567',unit))

# unzip input files if necessary:
####zipfiles<-system(paste('ls $RAWDATADIR/a',unum,'_',day,'*.dat.gz',sep=''),T)
#zipfiles<-system(paste('ls ',datadir,'/a',unum,'_',day,'*.dat.gz',sep=''),T)
zipfiles<-list.files(datadir,glob2rx(paste('a',unum,'_',day,'*.dat.gz',sep='')))
if(length(zipfiles)>0){
	for(i in c(1:length(zipfiles))){
		print(paste('Unzipping data from file ',zipfiles[i]))
		#system(paste('gunzip -f',zipfiles[i]))
		system(paste('gunzip -f ',datadir,'/',zipfiles[i],sep=''))
	}
}
####files<-system(paste('ls $RAWDATADIR/a',unum,'_',day,'*.dat',sep=''),T)
####print(paste('ls $RAWDATADIR/a',unum,'_',day,'*.dat',sep=''))
#files<-system(paste('ls ',datadir,'/a',unum,'_',day,'*.dat',sep=''),T)
files<-list.files(datadir,glob2rx(paste('a',unum,'_',day,'*.dat',sep='')))
files=paste(datadir,'/',files,sep='')

#print(paste('ls ',datadir,'/a',unum,'_',day,'*.dat',sep=''))
print(files)

if(length(files)>0){ # do not run if no data

#set up data list
data<-list(time=NULL,co2=NULL,prs=NULL,ltmp=NULL,vlvhx=NULL,btmp=NULL,rh=NULL,sfl=NULL,fl1=NULL,fl2=NULL,fl3=NULL,fl4=NULL,fl5=NULL,fl6=NULL,sec=NULL)

print('Reading in data')

for (i in c(1:length(files))){

# decide if files are old-style with binary date, etc. or new-style all ascii
asctest<-!is.na(as.numeric(readChar(files[i],n=7)))
if(length(asctest)>0){ # otherwise crashes with an empty input file
if(asctest){ # time is in ascii format
        if(loc=='RBA'){ # RBA now has extra met data at end of row - temporarily flush
        input<-scan(files[i],sep="",what=list(time=0,co2=0,prs=0,ltmp=0,vlvhx="",btmp=0,rh=0,sfl=0,fl1=0,fl2=0,fl3=0,fl4=0,fl5=0,fl6=0),multi.line=F,fill=T,flush=T)
        input$fl5=rep(NA,length(input$fl5))
        input$fl6=rep(NA,length(input$fl6))
        } else {
        input<-scan(files[i],sep="",what=list(time=0,co2=0,prs=0,ltmp=0,vlvhx="",btmp=0,rh=0,sfl=0,fl1=0,fl2=0,fl3=0,fl4=0,fl5=0,fl6=0),multi.line=F,fill=T)
        }
	# time is in SSSSSS.S format
	input$sec<-input$time
} else { # time is in binary format
	# check if file already converted to ascii
	if(substr(readChar(files[i],n=7),3,3)!=':'|substr(readChar(files[i],n=7),6,6)!=':'){
	system('ignore previous error (from checking binary time format)')
	# convert binary data to ascii:
	print(paste('Dumping data from file ',files[i]))
	system(paste('../../data_dump -c 201 -f',files[i],'-A > temp.txt'))
	input<-scan('temp.txt',sep="",what=list(time="",port=0,char=0,co2=0,prs=0,ltmp=0,vlvhx="",btmp=0,rh=0,sfl=0,fl1=0,fl2=0,fl3=0,fl4=0,fl5=0,fl6=0),multi.line=F,fill=T)
	# time is in HH:MM:SS.S format
	# port is left over from ndaq system - not used
	# char is number of characters in line - not used
	system('rm temp.txt') 
	} else { # already converted
       	input<-scan(files[i],sep="",what=list(time="",co2=0,prs=0,ltmp=0,vlvhx="",btmp=0,rh=0,sfl=0,fl1=0,fl2=0,fl3=0,fl4=0,fl5=0,fl6=0),multi.line=F,fill=T)
	}
	# calculate seconds since midnight:
	input$sec<-(as.integer(substr(input$time,1,2))*60+as.integer(substr(input$time,4,5)))*60+as.real(substr(input$time,7,12))
} 
} else { # file is empty
	input<-scan(files[i],sep="",what=list(time=0,co2=0,prs=0,ltmp=0,vlvhx="",btmp=0,rh=0,sfl=0,fl1=0,fl2=0,fl3=0,fl4=0,fl5=0,fl6=0),multi.line=F,fill=T) 
}
# co2 is Li-820 CO2 in ppm
# prs is Li-820 cell pressure
# ltmp is Li-820 internal temperature
# vlvhx is hexadecimal valve status
# btmp is box temperature from Humitter or from custom system translated onto Humitter scale
# rh is sample relative humidity from Humitter or from custom system translated onto Humitter scale
# sfl is sample flow
# fl* is flow on lines 1 through 6

data$time<-c(data$time,input$time)
data$co2<-c(data$co2,input$co2)
data$prs<-c(data$prs,input$prs)
data$ltmp<-c(data$ltmp,input$ltmp)
data$vlvhx<-c(data$vlvhx,input$vlvhx)
data$btmp<-c(data$btmp,input$btmp)
data$rh<-c(data$rh,input$rh)
data$sfl<-c(data$sfl,input$sfl)
data$fl1<-c(data$fl1,input$fl1)
data$fl2<-c(data$fl2,input$fl2)
data$fl3<-c(data$fl3,input$fl3)
data$fl4<-c(data$fl4,input$fl4)
data$fl5<-c(data$fl5,input$fl5)
data$fl6<-c(data$fl6,input$fl6)
data$sec<-c(data$sec,input$sec)

} # looping through and cat-ing input files

# rezip input files if necessary
## if there was a mix of zipped and unzipped files, zip them all now
if(length(zipfiles)>0){
        #for(i in c(1:length(zipfiles))){
        #        print(paste('Zipping data to file ',zipfiles[i]))
        #        system(paste('gzip',substr(zipfiles[i],1,nchar(zipfiles[i])-3)))
	## if there was a mix of zipped and unzipped files, zip them all now
        for(i in c(1:length(files))){
                print(paste('Zipping data to file ',files[i]))
                system(paste('gzip',files[i]))
        }
}

print(paste('Analyzing ',file))

# determine data interval for labels
reso=round(median(data$sec[2:length(data$sec)]-data$sec[1:(length(data$sec)-1)],na.rm=T),1)

# define inlet line for LT surveillance gas:
ltlin<-6
# define vlvhx values corresponding to start of low pressure and high pressure leak checks
hpvh<-1060
lpvh<-2801
# if unit only has 3 sample lines, zero 2 flows and reset LT line and high pressure leak check vlvhx:
if(all(is.na(c(data$fl5,data$fl6)))){
	data$fl5<-rep(0,length(data$time))
	data$fl6<-rep(0,length(data$time))
	ltlin<-4
	hpvh<-1048
}

# convert text hexadecimal to number
data$vlvhx<-as.numeric(data$vlvhx)

print(paste('Approximately ',round(length(data$co2)/njog,0),' jogs'))

# handle any problems in serial communication to Li-820:
li820naok<-T  # accept occasional (less than 50% of nave period) NAs in Li820 data stream
if(li820naok){
data$prs[is.na(data$co2)]<-(-999.99)
data$ltmp[is.na(data$co2)]<-(-999.99)
data$co2[is.na(data$co2)]<-(-999.99)
data$co2[data$co2<100]<-(-999.99) # when system was having serial problems, it occasionallly put decimal in wrong place on CO2 value
data$co2[data$co2>3000]<-(-999.99) # when system was having serial problems, it occasionallly put decimal in wrong place on CO2 value
} else {
data$co2[data$co2<100]<-NA
data$co2[data$co2>3000]<-NA
}

# define line removal reporting function:
remrep<-function(){
	if(length(rem)>1){
		if((rem[2]-rem[1])!=1){
			print(paste('at:',rem[1]))
		}else{
			print(paste('from:',rem[1]))
		}
		if(length(rem)>2){
			for(i in c(2:(length(rem)-1))){
				if((rem[i]-rem[i-1])!=1&(rem[i+1]-rem[i])!=1){
					print(paste('at:',rem[i]))
				}else if((rem[i]-rem[i-1])!=1){
					print(paste('from:',rem[i]))
				}else if((rem[i+1]-rem[i])!=1){
					print(paste('to:',rem[i]))
				}
			}
			if((rem[length(rem)]-rem[length(rem)-1])!=1){
				print(paste('at:',rem[length(rem)]))
			}else{
				print(paste('to:',rem[length(rem)]))
			}
		} else {
			if((rem[2]-rem[1])!=1){
				print(paste('at:',rem[2]))
			}else{
				print(paste('to:',rem[2]))
			}
		}
	} else {
		print(paste('at:',rem[1]))
	}
}

# remove any lines with NAs in them (currently no way to accept co2 data when diagnostic data missing):
data<-data.frame(data)
data<-na.omit(data)
if(length(attributes(data)$na.action)>=1){
	rem<-unique(attributes(data)$na.action)
	print(paste('Removed ',length(rem),' records with NA in them'))
	remrep()
}
attributes(data)$na.action<-NULL

if(length(data$vlvhx)>0){ # do not run if all data screened

# fix problem of having last value in file from 1st second of next day:
if(data$sec[length(data$sec)]==0){
	data$sec[length(data$sec)]<-NA
	print('Removed last value at 00:00:00')
}
if(data$sec[1]==1){
	data<-rbind(data[1,],data)	
	print('Added dummy first value at 00:00:00')
}

# remove lines before ico or after eco:
data$co2[data$sec/86400<ico]<-NA
data$co2[data$sec/86400>eco]<-NA
data<-na.omit(data)
if(length(attributes(data)$na.action)>=1){
	rem<-unique(attributes(data)$na.action)
	print(paste('Removed ',length(rem),' record(s) > eco or < ico'))
	remrep()
}
attributes(data)$na.action<-NULL

} # do not run if all data screened

if(length(data$vlvhx)>0){ # do not run if all data screened

# remove any lines from manual control periods:
data$co2[data$vlvhx>=80000]<-NA
data<-na.omit(data)
if(length(attributes(data)$na.action)>=1){
	rem<-unique(attributes(data)$na.action)
	print(paste('Removed ',length(rem),' records from manual control period(s)'))
	remrep()
}
attributes(data)$na.action<-NULL

} # do not run if all data screened

if(length(data$vlvhx)>(njog-1)){ # do not run if all (or almost all) data screened

# convert flows from volts to sccm
# for AWM3300, use a 3rd-order approximation over region 0-500 sccm
# for AWM43600, use empirical fit from NCAR/VUV2/AO2/awm43600cal.xls for 1-3 V (0 - ~ 1000 sccm) and linear interp of factory cals above that
print('Converting flows')
# first read in mass-flowmeter information
# type 1 = AWM3300, type 2 = AWM43600 (not using serial number now, but will implement flow-meter specific calibrations in future)
mfms<-scan('../../aircoa_mfms.txt',what=list(loc='',unit='',date=0,time='',SLTYPE=0,SLSERNUM=0,L1TYPE=0,L1SERNUM=0,L2TYPE=0,L2SERNUM=0,L3TYPE=0,L3SERNUM=0,L4TYPE=0,L4SERNUM=0,L5TYPE=0,L5SERNUM=0,L6TYPE=0,L6SERNUM=0),skip=1)
mfms$dayf<-as.numeric(substr(mfms$time,1,2))/24+as.numeric(substr(mfms$time,3,4))/60/24
#print(mfms)
hires<-c(0:100)/50+1
calvolts<-c(0,hires,3.1,3.8,4.4,4.7,4.89,5.0)
calsccm<-c(0,53.298*hires^4-422.09*hires^3+1343.9*hires^2-1565.4*hires+590.66,1000,2000,3000,4000,5000,6000)
#print(cbind(calvolts,calsccm))
# find most recent previous values
for (i in c(1:length(mfms$loc))){
	if(mfms$loc[i]==loc&(mfms$date[i]+mfms$dayf[i])<=(as.numeric(day)+ico)){
		sltype<-mfms$SLTYPE[i]	
		slsernum<-mfms$SLSERNUM[i]	
		l1type<-mfms$L1TYPE[i]	
		l1sernum<-mfms$L1SERNUM[i]	
		l2type<-mfms$L2TYPE[i]	
		l2sernum<-mfms$L2SERNUM[i]	
		l3type<-mfms$L3TYPE[i]	
		l3sernum<-mfms$L3SERNUM[i]	
		l4type<-mfms$L4TYPE[i]	
		l4sernum<-mfms$L4SERNUM[i]	
		l5type<-mfms$L5TYPE[i]	
		l5sernum<-mfms$L5SERNUM[i]	
		l6type<-mfms$L6TYPE[i]	
		l6sernum<-mfms$L6SERNUM[i]	
	}
}
# apply calibrations
# read in zero and span for every calibrated flow meter and apply here
mfmcals<-scan('../../aircoa_mfmcals.txt',what=list(date='',time='',type=0,sernum=0,zero=0,span=0),skip=1)
mfmcalsdatetime=strptime(paste(mfmcals$date,mfmcals$time),format='%y%m%d %H%M',tz="GMT")
mfmcalsyrdays=rep(365,length(mfmcalsdatetime$year))
mfmcalsyrdays[(mfmcalsdatetime$year+1900)/4-trunc((mfmcalsdatetime$year+1900)/4)==0]=366
mfmcals$yrfrac=round(mfmcalsdatetime$year+1900+as.numeric(difftime(mfmcalsdatetime,strptime(paste(mfmcalsdatetime$year+1900,'-01-01 00:00:00',sep=''),format='%Y-%m-%d %H:%M:%S',tz="GMT"),units='days'))/mfmcalsyrdays,7)
datadatetime=round.POSIXt(as.POSIXlt(strptime(day,format='%y%m%d',tz="GMT")+data$sec),units="secs")
if((datadatetime$year[1]+1900)/4-trunc((datadatetime$year[1]+1900)/4)==0){ datayrdays=366 } else { datayrdays=365 }
datayrfrac=round(datadatetime$year+1900+as.numeric(difftime(datadatetime,strptime(paste(datadatetime$year+1900,'-01-01 00:00:00',sep=''),format='%Y-%m-%d %H:%M:%S',tz="GMT"),units='days'))/datayrdays,8) # 0.3 sec resolution - needed?
if(sltype==1){
	data$sfl<-8.9902*(data$sfl-1)^3+4.7735*(data$sfl-1)^2+58.728*(data$sfl-1)
}else if(sltype==2){
	data$sfl<-approx(calvolts,calsccm,data$sfl)$y
}
if(any(is.element(slsernum,mfmcals$sernum))){
	if(sum(is.element(slsernum,mfmcals$sernum))>1){
		zeros=approx(mfmcals$yrfrac[mfmcals$sernum==slsernum],mfmcals$zero[mfmcals$sernum==slsernum],datayrfrac,rule=2)$y
		spans=approx(mfmcals$yrfrac[mfmcals$sernum==slsernum],mfmcals$span[mfmcals$sernum==slsernum],datayrfrac,rule=2)$y
	} else {
		zeros=rep(mfmcals$zero[mfmcals$sernum==slsernum],length(datayrfrac))
		spans=rep(mfmcals$span[mfmcals$sernum==slsernum],length(datayrfrac))
	}
	print(paste('Applying sfl mfm calibration: span =',median(spans,na.rm=T),', zero =',median(zeros,na.rm=T)))
	data$sfl=data$sfl*spans+zeros
}
if(l1type==1){
	data$fl1<-8.9902*(data$fl1-1)^3+4.7735*(data$fl1-1)^2+58.728*(data$fl1-1)
}else if(l1type==2){
	data$fl1<-approx(calvolts,calsccm,data$fl1)$y
}
if(any(is.element(l1sernum,mfmcals$sernum))){
	if(sum(is.element(l1sernum,mfmcals$sernum))>1){
		zeros=approx(mfmcals$yrfrac[mfmcals$sernum==l1sernum],mfmcals$zero[mfmcals$sernum==l1sernum],datayrfrac,rule=2)$y
		spans=approx(mfmcals$yrfrac[mfmcals$sernum==l1sernum],mfmcals$span[mfmcals$sernum==l1sernum],datayrfrac,rule=2)$y
	} else {
		zeros=rep(mfmcals$zero[mfmcals$sernum==l1sernum],length(datayrfrac))
		spans=rep(mfmcals$span[mfmcals$sernum==l1sernum],length(datayrfrac))
	}
	print(paste('Applying fl1 mfm calibration: span =',median(spans,na.rm=T),', zero =',median(zeros,na.rm=T)))
	data$fl1=data$fl1*spans+zeros
}
if(l2type==1){
	data$fl2<-8.9902*(data$fl2-1)^3+4.7735*(data$fl2-1)^2+58.728*(data$fl2-1)
}else if(l2type==2){
	data$fl2<-approx(calvolts,calsccm,data$fl2)$y
}
if(any(is.element(l2sernum,mfmcals$sernum))){
	if(sum(is.element(l2sernum,mfmcals$sernum))>1){
		zeros=approx(mfmcals$yrfrac[mfmcals$sernum==l2sernum],mfmcals$zero[mfmcals$sernum==l2sernum],datayrfrac,rule=2)$y
		spans=approx(mfmcals$yrfrac[mfmcals$sernum==l2sernum],mfmcals$span[mfmcals$sernum==l2sernum],datayrfrac,rule=2)$y
	} else {
		zeros=rep(mfmcals$zero[mfmcals$sernum==l2sernum],length(datayrfrac))
		spans=rep(mfmcals$span[mfmcals$sernum==l2sernum],length(datayrfrac))
	}
	print(paste('Applying fl2 mfm calibration: span =',median(spans,na.rm=T),', zero =',median(zeros,na.rm=T)))
	data$fl2=data$fl2*spans+zeros
}
if(l3type==1){
	data$fl3<-8.9902*(data$fl3-1)^3+4.7735*(data$fl3-1)^2+58.728*(data$fl3-1)
}else if(l3type==2){
	data$fl3<-approx(calvolts,calsccm,data$fl3)$y
}
if(any(is.element(l3sernum,mfmcals$sernum))){
	if(sum(is.element(l3sernum,mfmcals$sernum))>1){
		zeros=approx(mfmcals$yrfrac[mfmcals$sernum==l3sernum],mfmcals$zero[mfmcals$sernum==l3sernum],datayrfrac,rule=2)$y
		spans=approx(mfmcals$yrfrac[mfmcals$sernum==l3sernum],mfmcals$span[mfmcals$sernum==l3sernum],datayrfrac,rule=2)$y
	} else {
		zeros=rep(mfmcals$zero[mfmcals$sernum==l3sernum],length(datayrfrac))
		spans=rep(mfmcals$span[mfmcals$sernum==l3sernum],length(datayrfrac))
	}
	print(paste('Applying fl3 mfm calibration: span =',median(spans,na.rm=T),', zero =',median(zeros,na.rm=T)))
	data$fl3=data$fl3*spans+zeros
}
if(l4type==1){
	data$fl4<-8.9902*(data$fl4-1)^3+4.7735*(data$fl4-1)^2+58.728*(data$fl4-1)
}else if(l4type==2){
	data$fl4<-approx(calvolts,calsccm,data$fl4)$y
}
if(any(is.element(l4sernum,mfmcals$sernum))){
	if(sum(is.element(l4sernum,mfmcals$sernum))>1){
		zeros=approx(mfmcals$yrfrac[mfmcals$sernum==l4sernum],mfmcals$zero[mfmcals$sernum==l4sernum],datayrfrac,rule=2)$y
		spans=approx(mfmcals$yrfrac[mfmcals$sernum==l4sernum],mfmcals$span[mfmcals$sernum==l4sernum],datayrfrac,rule=2)$y
	} else {
		zeros=rep(mfmcals$zero[mfmcals$sernum==l4sernum],length(datayrfrac))
		spans=rep(mfmcals$span[mfmcals$sernum==l4sernum],length(datayrfrac))
	}
	print(paste('Applying fl4 mfm calibration: span =',median(spans,na.rm=T),', zero =',median(zeros,na.rm=T)))
	data$fl4=data$fl4*spans+zeros
}
if(l5type==1){
	data$fl5<-8.9902*(data$fl5-1)^3+4.7735*(data$fl5-1)^2+58.728*(data$fl5-1)
}else if(l5type==2){
	data$fl5<-approx(calvolts,calsccm,data$fl5)$y
}
if(any(is.element(l5sernum,mfmcals$sernum))){
	if(sum(is.element(l5sernum,mfmcals$sernum))>1){
		zeros=approx(mfmcals$yrfrac[mfmcals$sernum==l5sernum],mfmcals$zero[mfmcals$sernum==l5sernum],datayrfrac,rule=2)$y
		spans=approx(mfmcals$yrfrac[mfmcals$sernum==l5sernum],mfmcals$span[mfmcals$sernum==l5sernum],datayrfrac,rule=2)$y
	} else {
		zeros=rep(mfmcals$zero[mfmcals$sernum==l5sernum],length(datayrfrac))
		spans=rep(mfmcals$span[mfmcals$sernum==l5sernum],length(datayrfrac))
	}
	print(paste('Applying fl5 mfm calibration: span =',median(spans,na.rm=T),', zero =',median(zeros,na.rm=T)))
	data$fl5=data$fl5*spans+zeros
}
if(l6type==1){
	data$fl6<-8.9902*(data$fl6-1)^3+4.7735*(data$fl6-1)^2+58.728*(data$fl6-1)
}else if(l6type==2){
	data$fl6<-approx(calvolts,calsccm,data$fl6)$y
}
if(any(is.element(l6sernum,mfmcals$sernum))){
	if(sum(is.element(l6sernum,mfmcals$sernum))>1){
		zeros=approx(mfmcals$yrfrac[mfmcals$sernum==l6sernum],mfmcals$zero[mfmcals$sernum==l6sernum],datayrfrac,rule=2)$y
		spans=approx(mfmcals$yrfrac[mfmcals$sernum==l6sernum],mfmcals$span[mfmcals$sernum==l6sernum],datayrfrac,rule=2)$y
	} else {
		zeros=rep(mfmcals$zero[mfmcals$sernum==l6sernum],length(datayrfrac))
		spans=rep(mfmcals$span[mfmcals$sernum==l6sernum],length(datayrfrac))
	}
	print(paste('Applying fl6 mfm calibration: span =',median(spans,na.rm=T),', zero =',median(zeros,na.rm=T)))
	data$fl6=data$fl6*spans+zeros
}

# if using new custom T/RH (Honeywell HIH-3610 humidity sensor + ?) 
#	convert new (0.8 - 4.07 V = 0 - 100 % RH) to old Humitter scale (0 - 1 V = 0 - 100% RH)
#	convert new (0 - 1 V = 0 - 100 C) to old Humitter scale (0 - 1 V = -40 - 60 C)
if(trh=='custom') {
	data$rh<-(data$rh-0.8)/(4.07-0.8)
	data$btmp<-data$btmp+0.4
}

# convert hexadecimal valve status to nlin (line selected), ngas (cal or sample gas selected), lchk (0 for open, 1 for closed), pump (0 for off, 1 for on) 
print('Converting hexadecimal')
nlin<-rep(0,length(data$co2))
ngas<-rep(0,length(data$co2))
lchk<-rep(0,length(data$co2))
pump<-rep(0,length(data$co2))
hex<-unique(data$vlvhx)[!is.na(unique(data$vlvhx))]
for(i in c(1:length(hex))){
	temp<-hex[i]
	if(temp>=4000){ # bit 14
		pump[data$vlvhx==hex[i]]<-1 # micro pump on
		temp<-temp-4000
	}
	if(temp>=2000){ # bit 13
		pump[data$vlvhx==hex[i]]<-1 # micro pump on
		temp<-temp-2000
	}
	if(temp>=1000){ # bit 12 
		lchk[data$vlvhx==hex[i]]<-1 # leak check valve closed
		temp<-temp-1000
	}
	if(temp>=800){ # bit 11
		ngas[data$vlvhx==hex[i]]<-6 # guest cal valve selected
		temp<-temp-800
	}
	if(temp>=400){ # bit 10
		ngas[data$vlvhx==hex[i]]<-5 # hs2 cal gas selected 
		temp<-temp-400
	}
	if(temp>=200){ # bit 9
		ngas[data$vlvhx==hex[i]]<-4 # hs1 cal gas selected 
		temp<-temp-200
	}
	if(temp>=100){ # bit 8
		ngas[data$vlvhx==hex[i]]<-3 # ls1 cal gas selected 
		temp<-temp-100
	}
	if(temp>=80){ # bit 7
		ngas[data$vlvhx==hex[i]]<-2 # ls2 cal gas selected 
		temp<-temp-80
	}
	if(temp>=40){ # bit 6
		ngas[data$vlvhx==hex[i]]<-1 # sample gas selected 
		temp<-temp-40
	}
	if(temp>=20){ # bit 5
		nlin[data$vlvhx==hex[i]]<-6 # inlet 6 selected 
		temp<-temp-20
	}
	if(temp>=10){ # bit 4
		nlin[data$vlvhx==hex[i]]<-5 # inlet 5 selected 
		temp<-temp-10
	}
	if(temp>=8){ # bit 3
		nlin[data$vlvhx==hex[i]]<-4 # inlet 4 selected 
		temp<-temp-8
	}
	if(temp>=4){ # bit 2
		nlin[data$vlvhx==hex[i]]<-3 # inlet 3 selected 
		temp<-temp-4
	}
	if(temp>=2){ # bit 1
		nlin[data$vlvhx==hex[i]]<-2 # inlet 2 selected 
		temp<-temp-2
	}
	if(temp>=1){ # bit 0
		nlin[data$vlvhx==hex[i]]<-1 # inlet 1 selected 
		temp<-temp-1
	}
}
data$nlin<-nlin
data$ngas<-ngas
data$lchk<-lchk
data$pump<-pump
rm(nlin,ngas,lchk,pump)

# if unit is only using lines 2 and 3 (e.g. HDP), reset low pressure leak check vlvhx:
if(any(data$nlin==2&data$ngas==6)){ # using line 2 instead of line 1 as nafion purge during low pressure leak check
	lpvh<-2802
}

# find start of first complete jog:
if(all(data$vlvhx[1:njog]==data$vlvhx[1])|(all(data$vlvhx[1:nlpp]==data$vlvhx[1])&all(data$lchk[(nlpp+1):njog]==1))|(all(data$vlvhx[1:nhpp]==data$vlvhx[1])&all(data$lchk[1:njog]==1))) { # starts on an even jog
	pos1<-1 
} else { # find first change in vlvhx
	pos1<-c(1:length(data$vlvhx))[cumsum((data$vlvhx[1:length(data$vlvhx)]!=data$vlvhx[1])*1)==1]
}
# check again:
if(length(data$vlvhx)>=(pos1+njog-1)){
	if(all(data$vlvhx[pos1:(pos1+njog-1)]==data$vlvhx[pos1])|(all(data$vlvhx[pos1:(pos1+nlpp-1)]==data$vlvhx[pos1])&all(data$lchk[(pos1+nlpp):(pos1+njog-1)]==1))|
		(all(data$vlvhx[pos1:(pos1+nhpp-1)]==data$vlvhx[pos1])&all(data$lchk[pos1:(pos1+njog-1)]==1))) { # pos1 is at start of an even jog
	} else { # find next change in vlvhx
		pos1<-c(pos1:length(data$vlvhx))[cumsum((data$vlvhx[pos1:length(data$vlvhx)]!=data$vlvhx[pos1])*1)==1]
	}
 	print(paste('Skipping initial ',pos1-1,' records.'))
} else {
	print('No complete records')
}

# check for any incomplete jogs:
post<-pos1
for(i in c(1:1000)){ # can handle 1000 jogs (only 576 2.5 min jogs in a day)
    if(post<=(length(data$vlvhx)-njog+1)){ # not at end of file, ok to proceed
	if(data$vlvhx[post]==hpvh) { # high pressure leak check
		if(all(data$vlvhx[post:(post+nhpp-1)]==data$vlvhx[post])&all(data$lchk[post:(post+njog-1)]==1)) { # good jog
			post<-post+njog
		} else {
			print(paste('bad high pressure leak check at ',post))
			# go to change in vlvhx and check again
			posn<-c(post:length(data$vlvhx))[cumsum((data$vlvhx[post:length(data$vlvhx)]!=data$vlvhx[post])*1)==1]
			if(length(posn)==0){ posn<-length(data$vlvhx)+1 } # hit end of file
			data$co2[post:(posn-1)]<-NA			
			post<-posn
		}
	} else if(data$vlvhx[post]==lpvh) { # low pressure leak check
		if(all(data$vlvhx[post:(post+nlpp-1)]==data$vlvhx[post])&all(data$lchk[(post+nlpp):(post+njog-1)]==1)) { # good jog
			post<-post+njog
		} else {
			print(paste('bad ambient pressure leak check at ',post))
			# go to change in vlvhx and check again
			posn<-c(post:length(data$vlvhx))[cumsum((data$vlvhx[post:length(data$vlvhx)]!=data$vlvhx[post])*1)==1]
			if(length(posn)==0){ posn<-length(data$vlvhx)+1 } # hit end of file
			data$co2[post:(posn-1)]<-NA			
			post<-posn
		}
	} else { # sample or cal jog
		if(all(data$vlvhx[post:(post+njog-1)]==data$vlvhx[post])) { # good jog
				post<-post+njog
		} else {
			print(paste('bad sample or cal. period at ',post))
			# go to change in vlvhx and check again
			posn<-c(post:length(data$vlvhx))[cumsum((data$vlvhx[post:length(data$vlvhx)]!=data$vlvhx[post])*1)==1]
			if(length(posn)==0){ posn<-length(data$vlvhx)+1 } # hit end of file
			data$co2[post:(posn-1)]<-NA			
			post<-posn
		}
	}
    }
}

# remove any incomplete jogs:
data<-na.omit(data)
if(length(attributes(data)$na.action)>=1){
	rem<-unique(attributes(data)$na.action)
	print(paste('Removed ',length(rem),' incomplete jog records'))
	remrep()
}
attributes(data)$na.action<-NULL

if(li820naok){ # accept occasional NAs in Li820 data stream

	# remove any jogs with more than 50% NAs for Li820 data in averaging period:
	post<-pos1
	for(i in c(1:1000)){ # can handle 1000 jogs (only 576 2.5 min jogs in a day)
    		if(post<=(length(data$vlvhx)-njog+1)){ # not at end of file
			if(sum(rep(1,nave)[data$co2[(post+skip):(post+njog-1)]!=-999.99])>(nave*0.5)){ # good jog
				post<-post+njog
			} else {
				print(paste('more than 50% Li820 NAs at ',post))
				data$co2[post:(post+njog-1)]<-NA
				post<-post+njog
			}
		}
	}

	data<-na.omit(data)
	if(length(attributes(data)$na.action)>=1){
		rem<-unique(attributes(data)$na.action)
		print(paste('Removed ',length(rem),' records from jogs with too many Li820 NAs'))
		remrep()
	}
	attributes(data)$na.action<-NULL

	data$prs[data$co2==-999.99]<-NA
	data$ltmp[data$co2==-999.99]<-NA
	data$co2[data$co2==-999.99]<-NA

}

# set last possible jog start point:
pos2<-length(data$co2)-njog+1 

# set up matrices and vectors for averaging variables over jog cycle: 
shapeslfl<-matrix(0,18,njog) # 18 = 6 lines + 2 leak checks + 4 cals + 2 cals * 2 diff preceding gases + 1 all sample + 1 all cal
shapessfl<-matrix(0,18,njog)
shapesprs<-matrix(0,18,njog)
shapesco2<-matrix(0,18,njog)
# shapeslfl is line flow, shapessfl is sample flow, shapesprs is sample pressure, and shapesco2 is co2
shapescnt<-rep(0,18) # counter
shapesngas<-rep(0,18) # to check that bins are all one gas
shapesnlin<-rep(0,18) # to check that bins are all one line

# shapes binned according to:
# BIN  STATUS
#  1   line 1
#  2   line 2
#  3   line 3
#  4   line 4
#  5   line 5
#  6   line 6
#  7   ls2
#  8   ls1
#  9   hs1
# 10   hs2
# 11   high pressure leak check
# 12   ambient pressure leak check
# 13   ls1 preceded by hs1
# 14   ls1 preceded by ls2
# 15   hs1 preceded by hs2
# 16   hs1 preceded by ls1
# 17   all sample lines averaged
# 18   all cal lines averaged

# set up output list, with 42 components corresponding to:       		   		component
# AVERAGES of sec, co2, prs, ltmp, btmp, rh, sfl, fl1, fl2,        				[ 1- 9]
#             fl3, fl4, fl5, fl6, nlin, ngas, lchk, pump 					[10-17] 
# STD DEVS of co2, prs, ltmp, btmp, rh, sfl, fl1, fl2, fl3, fl4, fl5, fl6       		[18-29]  # add 'sd' to start of var names
# SLOPES of co2, prs, ltmp, btmp, rh, sfl, fl1, fl2, fl3, fl4, fl5, fl6         		[30-41]  # add 'sl' to start of var names
# shapesbin                		                                         		[42]
output<-list(sec=NULL,co2=NULL,prs=NULL,ltmp=NULL,btmp=NULL,rh=NULL,sfl=NULL,fl1=NULL,fl2=NULL,
fl3=NULL,fl4=NULL,fl5=NULL,fl6=NULL,nlin=NULL,ngas=NULL,lchk=NULL,pump=NULL,
sdco2=NULL,sdprs=NULL,sdltmp=NULL,sdbtmp=NULL,sdrh=NULL,sdsfl=NULL,sdfl1=NULL,sdfl2=NULL,sdfl3=NULL,sdfl4=NULL,sdfl5=NULL,sdfl6=NULL,
slco2=NULL,slprs=NULL,slltmp=NULL,slbtmp=NULL,slrh=NULL,slsfl=NULL,slfl1=NULL,slfl2=NULL,slfl3=NULL,slfl4=NULL,slfl5=NULL,slfl6=NULL,
shpbn=NULL)


## loop through all jogs, calculating shapes, means, std devs, and slopes:

pos<-pos1
jog<-0
for (i in c(1:1000)) {  # can handle 1000 jogs (only 576 2.5 min jogs in a day) 

	if(pos<=pos2){ # not at end of file, ok to proceed

		jog<-jog+1

		print(paste('Processing jog',jog,'at position',pos))

		# determine bin for shape averaging:
		shapesbin<-0
		if(all(data$ngas[pos:(pos+njog-1)]==1)&all(data$lchk[pos:(pos+njog-1)]!=1)) { # line sampling
			shapesbin<-data$nlin[pos]
		} else if (all(data$ngas[pos:(pos+nlpp-1)]==6)&all(data$lchk[(pos+nlpp):(pos+njog-1)]==1)){ # low p leak check
			shapesbin<-12
		} else if (all(data$lchk[pos:(pos+njog-1)]==1)){ # high p leak check
			shapesbin<-11
		} else { # cal gas
			shapesbin<-data$ngas[pos]+5
			if(i>1){ # not first jog
				if((data$ngas[pos]==3)&(data$ngas[pos-1]==4)) shapesbin<-13 # ls1 preceded by hs1
				if((data$ngas[pos]==3)&(data$ngas[pos-1]==2)) shapesbin<-14 # ls1 preceded by ls2
				if((data$ngas[pos]==4)&(data$ngas[pos-1]==5)) shapesbin<-15 # hs1 preceded by hs2
				if((data$ngas[pos]==4)&(data$ngas[pos-1]==3)) shapesbin<-16 # hs1 preceded by ls1
			} 
		}
		output$shpbn<-c(output$shpbn,shapesbin)

		# calculate shapes:
		if (shapesbin!=0){
			co2window<-data$co2[pos:(pos+njog-1)]
			if(any(is.na(co2window))){
				co2window<-approx(c(1:njog)[!is.na(co2window)],co2window[!is.na(co2window)],c(1:njog),rule=2)$y # interpolate to replace any CO2 NAs
			}
			shapesco2[shapesbin,(1:njog)]<-shapesco2[shapesbin,(1:njog)]+co2window
			prswindow<-data$prs[pos:(pos+njog-1)]
			if(any(is.na(prswindow))){
				prswindow<-approx(c(1:njog)[!is.na(prswindow)],prswindow[!is.na(prswindow)],c(1:njog),rule=2)$y # interpolate to replace any pressure NAs
			}
			shapesprs[shapesbin,(1:njog)]<-shapesprs[shapesbin,(1:njog)]+prswindow
			shapessfl[shapesbin,(1:njog)]<-shapessfl[shapesbin,(1:njog)]+data$sfl[pos:(pos+njog-1)]
			if(data$nlin[pos]==1){ shapeslfl[shapesbin,(1:njog)]<-shapeslfl[shapesbin,(1:njog)]+data$fl1[pos:(pos+njog-1)]}
			if(data$nlin[pos]==2){ shapeslfl[shapesbin,(1:njog)]<-shapeslfl[shapesbin,(1:njog)]+data$fl2[pos:(pos+njog-1)]}
			if(data$nlin[pos]==3){ shapeslfl[shapesbin,(1:njog)]<-shapeslfl[shapesbin,(1:njog)]+data$fl3[pos:(pos+njog-1)]}
			if(data$nlin[pos]==4){ shapeslfl[shapesbin,(1:njog)]<-shapeslfl[shapesbin,(1:njog)]+data$fl4[pos:(pos+njog-1)]}
			if(data$nlin[pos]==5){ shapeslfl[shapesbin,(1:njog)]<-shapeslfl[shapesbin,(1:njog)]+data$fl5[pos:(pos+njog-1)]}
			if(data$nlin[pos]==6){ shapeslfl[shapesbin,(1:njog)]<-shapeslfl[shapesbin,(1:njog)]+data$fl6[pos:(pos+njog-1)]}
			shapesngas[shapesbin]<-shapesngas[shapesbin]+mean(data$ngas[pos:(pos+njog-1)])
			shapesnlin[shapesbin]<-shapesnlin[shapesbin]+mean(data$nlin[pos:(pos+njog-1)])
			shapescnt[shapesbin]<-shapescnt[shapesbin]+1
		} else { 
			print(paste('No shapebin at position ',pos))
		}

		# calculate jog averages, std devs, and slopes:
		output$sec<-c(output$sec,mean(data$sec[(pos+skip):(pos+skip+nave-1)]))
		output$co2<-c(output$co2,mean(data$co2[(pos+skip):(pos+skip+nave-1)],na.rm=T))
		output$prs<-c(output$prs,mean(data$prs[(pos+skip):(pos+skip+nave-1)],na.rm=T))
		output$ltmp<-c(output$ltmp,mean(data$ltmp[(pos+skip):(pos+skip+nave-1)],na.rm=T))
		output$btmp<-c(output$btmp,mean(data$btmp[(pos+skip):(pos+skip+nave-1)]))
		output$rh<-c(output$rh,mean(data$rh[(pos+skip):(pos+skip+nave-1)]))
		output$sfl<-c(output$sfl,mean(data$sfl[(pos+skip):(pos+skip+nave-1)]))
		output$fl1<-c(output$fl1,mean(data$fl1[(pos+skip):(pos+skip+nave-1)]))
		output$fl2<-c(output$fl2,mean(data$fl2[(pos+skip):(pos+skip+nave-1)]))
		output$fl3<-c(output$fl3,mean(data$fl3[(pos+skip):(pos+skip+nave-1)]))
		output$fl4<-c(output$fl4,mean(data$fl4[(pos+skip):(pos+skip+nave-1)]))
		output$fl5<-c(output$fl5,mean(data$fl5[(pos+skip):(pos+skip+nave-1)]))
		output$fl6<-c(output$fl6,mean(data$fl6[(pos+skip):(pos+skip+nave-1)]))
		output$nlin<-c(output$nlin,mean(data$nlin[(pos+skip):(pos+skip+nave-1)]))
		output$ngas<-c(output$ngas,mean(data$ngas[(pos+skip):(pos+skip+nave-1)]))
		output$lchk<-c(output$lchk,mean(data$lchk[(pos+skip):(pos+skip+nave-1)]))
		output$pump<-c(output$pump,mean(data$pump[(pos+skip):(pos+skip+nave-1)]))
		output$sdco2<-c(output$sdco2,sqrt(var(data$co2[(pos+skip):(pos+skip+nave-1)],na.rm=T)))
		output$sdprs<-c(output$sdprs,sqrt(var(data$prs[(pos+skip):(pos+skip+nave-1)],na.rm=T)))
		output$sdltmp<-c(output$sdltmp,sqrt(var(data$ltmp[(pos+skip):(pos+skip+nave-1)],na.rm=T)))
		output$sdbtmp<-c(output$sdbtmp,sqrt(var(data$btmp[(pos+skip):(pos+skip+nave-1)],na.rm=T)))
		output$sdrh<-c(output$sdrh,sqrt(var(data$rh[(pos+skip):(pos+skip+nave-1)],na.rm=T)))
		output$sdsfl<-c(output$sdsfl,sqrt(var(data$sfl[(pos+skip):(pos+skip+nave-1)],na.rm=T)))
		output$sdfl1<-c(output$sdfl1,sqrt(var(data$fl1[(pos+skip):(pos+skip+nave-1)],na.rm=T)))
		output$sdfl2<-c(output$sdfl2,sqrt(var(data$fl2[(pos+skip):(pos+skip+nave-1)],na.rm=T)))
		output$sdfl3<-c(output$sdfl3,sqrt(var(data$fl3[(pos+skip):(pos+skip+nave-1)],na.rm=T)))
		output$sdfl4<-c(output$sdfl4,sqrt(var(data$fl4[(pos+skip):(pos+skip+nave-1)],na.rm=T)))
		output$sdfl5<-c(output$sdfl5,sqrt(var(data$fl5[(pos+skip):(pos+skip+nave-1)],na.rm=T)))
		output$sdfl6<-c(output$sdfl6,sqrt(var(data$fl6[(pos+skip):(pos+skip+nave-1)],na.rm=T)))
		output$slco2<-c(output$slco2,mean(data$co2[(pos+skip+ntrd):(pos+skip+nave-1)],na.rm=T)-mean(data$co2[(pos+skip):(pos+skip+ntrd-1)],na.rm=T))
		output$slprs<-c(output$slprs,mean(data$prs[(pos+skip+ntrd):(pos+skip+nave-1)],na.rm=T)-mean(data$prs[(pos+skip):(pos+skip+ntrd-1)],na.rm=T))
		output$slltmp<-c(output$slltmp,mean(data$ltmp[(pos+skip+ntrd):(pos+skip+nave-1)],na.rm=T)-mean(data$ltmp[(pos+skip):(pos+skip+ntrd-1)],na.rm=T))
		output$slbtmp<-c(output$slbtmp,mean(data$btmp[(pos+skip+ntrd):(pos+skip+nave-1)])-mean(data$btmp[(pos+skip):(pos+skip+ntrd-1)]))
		output$slrh<-c(output$slrh,mean(data$rh[(pos+skip+ntrd):(pos+skip+nave-1)])-mean(data$rh[(pos+skip):(pos+skip+ntrd-1)]))
		output$slsfl<-c(output$slsfl,mean(data$sfl[(pos+skip+ntrd):(pos+skip+nave-1)])-mean(data$sfl[(pos+skip):(pos+skip+ntrd-1)]))
		output$slfl1<-c(output$slfl1,mean(data$fl1[(pos+skip+ntrd):(pos+skip+nave-1)])-mean(data$fl1[(pos+skip):(pos+skip+ntrd-1)]))
		output$slfl2<-c(output$slfl2,mean(data$fl2[(pos+skip+ntrd):(pos+skip+nave-1)])-mean(data$fl2[(pos+skip):(pos+skip+ntrd-1)]))
		output$slfl3<-c(output$slfl3,mean(data$fl3[(pos+skip+ntrd):(pos+skip+nave-1)])-mean(data$fl3[(pos+skip):(pos+skip+ntrd-1)]))
		output$slfl4<-c(output$slfl4,mean(data$fl4[(pos+skip+ntrd):(pos+skip+nave-1)])-mean(data$fl4[(pos+skip):(pos+skip+ntrd-1)]))
		output$slfl5<-c(output$slfl5,mean(data$fl5[(pos+skip+ntrd):(pos+skip+nave-1)])-mean(data$fl5[(pos+skip):(pos+skip+ntrd-1)]))
		output$slfl6<-c(output$slfl6,mean(data$fl6[(pos+skip+ntrd):(pos+skip+nave-1)])-mean(data$fl6[(pos+skip):(pos+skip+ntrd-1)]))
		 
		pos<-pos+njog
		
	} # not at end of file

} # looping through jogs

# estimate span and offset:
co2scl<-1
co2off<-0 # because Li-820 gives approximate ppm, use these values for diagnostics plotting

# finish processing shapes for composite bins:
shapesco2[8,]<-shapesco2[8,]+shapesco2[13,]+shapesco2[14,] # all LS1s
shapesprs[8,]<-shapesprs[8,]+shapesprs[13,]+shapesprs[14,]
shapessfl[8,]<-shapessfl[8,]+shapessfl[13,]+shapessfl[14,]
shapeslfl[8,]<-shapeslfl[8,]+shapeslfl[13,]+shapeslfl[14,]
shapesngas[8]<-shapesngas[8]+shapesngas[13]+shapesngas[14]
shapesnlin[8]<-shapesnlin[8]+shapesnlin[13]+shapesnlin[14]
shapescnt[8]<-shapescnt[8]+shapescnt[13]+shapescnt[14]

shapesco2[9,]<-shapesco2[9,]+shapesco2[15,]+shapesco2[16,] # all HS1s
shapesprs[9,]<-shapesprs[9,]+shapesprs[15,]+shapesprs[16,]
shapessfl[9,]<-shapessfl[9,]+shapessfl[15,]+shapessfl[16,]
shapeslfl[9,]<-shapeslfl[9,]+shapeslfl[15,]+shapeslfl[16,]
shapesngas[9]<-shapesngas[9]+shapesngas[15]+shapesngas[16]
shapesnlin[9]<-shapesnlin[9]+shapesnlin[15]+shapesnlin[16]
shapescnt[9]<-shapescnt[9]+shapescnt[15]+shapescnt[16]

shapesco2[17,]<-shapesco2[1,]+shapesco2[2,]+shapesco2[3,]+shapesco2[4,]+shapesco2[5,]+shapesco2[6,] # all sample lines (including LT)
shapesprs[17,]<-shapesprs[1,]+shapesprs[2,]+shapesprs[3,]+shapesprs[4,]+shapesprs[5,]+shapesprs[6,]
shapessfl[17,]<-shapessfl[1,]+shapessfl[2,]+shapessfl[3,]+shapessfl[4,]+shapessfl[5,]+shapessfl[6,]
shapeslfl[17,]<-shapeslfl[1,]+shapeslfl[2,]+shapeslfl[3,]+shapeslfl[4,]+shapeslfl[5,]+shapeslfl[6,]
shapesngas[17]<-shapesngas[1]+shapesngas[2]+shapesngas[3]+shapesngas[4]+shapesngas[5]+shapesngas[6]
shapesnlin[17]<-shapesnlin[1]+shapesnlin[2]+shapesnlin[3]+shapesnlin[4]+shapesnlin[5]+shapesnlin[6]
shapescnt[17]<-shapescnt[1]+shapescnt[2]+shapescnt[3]+shapescnt[4]+shapescnt[5]+shapescnt[6]

shapesco2[18,]<-shapesco2[7,]+shapesco2[8,]+shapesco2[9,]+shapesco2[10,] # all cals
shapesprs[18,]<-shapesprs[7,]+shapesprs[8,]+shapesprs[9,]+shapesprs[10,]
shapessfl[18,]<-shapessfl[7,]+shapessfl[8,]+shapessfl[9,]+shapessfl[10,]
shapeslfl[18,]<-shapeslfl[7,]+shapeslfl[8,]+shapeslfl[9,]+shapeslfl[10,]
shapesngas[18]<-shapesngas[7]+shapesngas[8]+shapesngas[9]+shapesngas[10]
shapesnlin[18]<-shapesnlin[7]+shapesnlin[8]+shapesnlin[9]+shapesnlin[10]
shapescnt[18]<-shapescnt[7]+shapescnt[8]+shapescnt[9]+shapescnt[10]

# convert shapes from sums to averages:
for(i in c(1:18)){
	if(shapescnt[i]!=0){
		shapesco2[i,]<-shapesco2[i,]/shapescnt[i]
		shapesprs[i,]<-shapesprs[i,]/shapescnt[i]
		shapessfl[i,]<-shapessfl[i,]/shapescnt[i]
		shapeslfl[i,]<-shapeslfl[i,]/shapescnt[i]
		shapesngas[i]<-shapesngas[i]/shapescnt[i]
		shapesnlin[i]<-shapesnlin[i]/shapescnt[i]
	}
}

if(!is.null(output$sec)){ # do not run if all jogs bad

# write out processed data:
outmat<-rbind(round(output$sec,1),round(output$co2,3),round(output$prs,3),round(output$ltmp,3),round(output$btmp,3),round(output$rh,3),
round(output$sfl,2),round(output$fl1,2),round(output$fl2,2),round(output$fl3,2),round(output$fl4,2),round(output$fl5,2),round(output$fl6,2),
output$nlin,output$ngas,output$lchk,output$pump,
round(output$sdco2,3),round(output$sdprs,3),round(output$sdltmp,3),round(output$sdbtmp,3),round(output$sdrh,3),
round(output$sdsfl,2),round(output$sdfl1,2),round(output$sdfl2,2),round(output$sdfl3,2),round(output$sdfl4,2),round(output$sdfl5,2),round(output$sdfl6,2),
round(output$slco2,3),round(output$slprs,3),round(output$slltmp,3),round(output$slbtmp,3),round(output$slrh,3),
round(output$slsfl,2),round(output$slfl1,2),round(output$slfl2,2),round(output$slfl3,2),round(output$slfl4,2),round(output$slfl5,2),round(output$slfl6,2),output$shpbn)
write(outmat,paste(month,'/',loc,'/',file,outflag,".out",sep=''),ncol=42)
system(paste('chmod 664 ',month,'/',loc,'/',file,outflag,".out",sep='')) # allows program to be run and output to be used by different users 

# convert seconds into day of month + day fraction and handle midnight:
domfrachr<-as.integer(substr(day,5,6))+data$sec/86400 # high rate (0.5 - 2 Hz) values
hodfrachr<-data$sec/60/60 # high rate (0.5 - 2 Hz) values
domfrac<-as.integer(substr(day,5,6))+output$sec/86400 # ~ 2.5-minute values
hodfrac<-output$sec/60/60 # ~ 2.5-minute values
daystep<-output$sec-output$sec[1] # positive except for if last value of day is recorded as 00:00:00 of the next day 
domfrac[daystep<0]<-domfrac[daystep<0]+1
hodfrac[daystep<0]<-hodfrac[daystep<0]+1


## make plots: 

if(png){png(paste(file,outflag,"_co2shapes.png",sep=''),width=950,height=550)} else {
bitmap(paste(file,outflag,"_co2shapes.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)}

## plot all shapes:

# set up title list:
tlist<-c('Line 1','Line 2','Line 3','Line 4','Line 5','Line 6','LS2','LS1','HS1','HS2','High P Check',
'Low P Check','LS1 after HS1','LS1 after LS2','HS1 after HS2','HS1 after LS1')
tlist[ltlin]<-paste(tlist[ltlin],' (LT)',sep='')

# set up vector of shape bins:
binorder<-NULL
for(i in c(1:6)){
	if(shapescnt[i]!=0&any(!is.na(shapesco2[i,]))){
		binorder<-c(binorder,i)
	}
}
binorder<-c(binorder,7:10,13:16,12:11)

# co2 (LOC_YYMMDDf_1.png):
par(mfrow=c(4,4))
par(oma=c(2.5,0,2,0))
par(mar=c(2,3,2,1)+.1)
par(mgp=c(2,1,0))
for(i in binorder){
	plot(c(1:njog),shapesco2[i,]*co2scl+co2off,main=paste(tlist[i],' gas=',round(shapesngas[i],2),
     ' lin=',round(shapesnlin[i],2),' n=',shapescnt[i],sep=''),xlab='',ylab='approx. ppm',type='n',
     ylim=c(-5,5)+mean(shapesco2[i,(skip+1):njog],na.rm=T)*co2scl+co2off)
	points(c(1:njog),shapesco2[i,]*co2scl+co2off,col=2)
	temp1<-round(mean(output$slco2[output$shpbn==i]/ntrd*60,na.rm=T),2)
	temp2<-round(mean(output$sdco2[output$shpbn==i],na.rm=T),2)
	if(i==8){ temp1<-round(mean(output$slco2[output$shpbn==8|output$shpbn==13|output$shpbn==14]/ntrd*60,na.rm=T),2)
	temp2<-round(mean(output$sdco2[output$shpbn==8|output$shpbn==13|output$shpbn==14],na.rm=T),2)}
	if(i==9){ temp1<-round(mean(output$slco2[output$shpbn==9|output$shpbn==15|output$shpbn==16]/ntrd*60,na.rm=T),2)
	temp2<-round(mean(output$sdco2[output$shpbn==9|output$shpbn==15|output$shpbn==16],na.rm=T),2)}
	text(njog*.67,mean(shapesco2[i,(skip+1):njog],na.rm=T)*co2scl+co2off-2.5,paste('Ave. Slope = ',temp1,' ppm/min'))
	text(njog*.67,mean(shapesco2[i,(skip+1):njog],na.rm=T)*co2scl+co2off-4,paste('Ave. Std. Dev. = ',temp2,' ppm'))
	abline(mean(shapesco2[i,(skip+1):njog],na.rm=T)*co2scl+co2off+1,0)
	abline(mean(shapesco2[i,(skip+1):njog],na.rm=T)*co2scl+co2off-1,0)
	abline(v=skip+.5)
}
mtext(paste(file,outflag,'  CO2 Shapes',sep=''),3,0.5,T)
mtext(paste('Jog Position (',reso,' second resolution)',sep=''),1,1,T)

dev.off()

if(png){png(paste(file,outflag,"_co2shapes_ave.png",sep=''),width=950,height=550)} else {
bitmap(paste(file,outflag,"_co2shapes_ave.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)}

# co2 averaged by sample or cal (LOC_YYMMDDf_2.png):
par(mfrow=c(2,2))
par(oma=c(2.5,0,2,0))
par(mar=c(2,3,2,1)+.1)
par(mgp=c(2,1,0))
if(shapescnt[17]!=0&any(!is.na(shapesco2[17,]))){
	plot(c(1:njog),shapesco2[17,]*co2scl+co2off,main=paste('All Sample Lines, gas=',round(shapesngas[17],2),
     ' lin=',round(shapesnlin[17],2),' n=',shapescnt[17],sep=''),xlab='',ylab='approx. ppm',type='n',
     ylim=c(-5,5)+mean(shapesco2[17,(skip+1):njog],na.rm=T)*co2scl+co2off)
	points(c(1:njog),shapesco2[17,]*co2scl+co2off,col=2)
	abline(v=skip+.5)
	abline(mean(shapesco2[17,(skip+1):njog],na.rm=T)*co2scl+co2off+1,0)
	abline(mean(shapesco2[17,(skip+1):njog],na.rm=T)*co2scl+co2off-1,0)

	plot(c(1:njog),shapesco2[17,]*co2scl+co2off,main=paste('All Sample Lines, gas=',round(shapesngas[17],2),
     ' lin=',round(shapesnlin[17],2),' n=',shapescnt[17],' (zoomed)',sep=''),xlab='',ylab='approx. ppm',type='n',
     ylim=c(-.5,.5)+mean(shapesco2[17,(skip+1):njog],na.rm=T)*co2scl+co2off)
	points(c(1:njog),shapesco2[17,]*co2scl+co2off,col=2)
	temp1<-round(mean(output$slco2[output$shpbn<7]/ntrd*60,na.rm=T),2)
	temp2<-round(mean(output$sdco2[output$shpbn<7],na.rm=T),2)
	text(njog*.67,mean(shapesco2[17,(skip+1):njog],na.rm=T)*co2scl+co2off-.25,paste('Ave. Slope = ',temp1,' ppm/min'))
	text(njog*.67,mean(shapesco2[17,(skip+1):njog],na.rm=T)*co2scl+co2off-.4,paste('Ave. Std. Dev. = ',temp2,' ppm'))
	abline(mean(shapesco2[17,(skip+1):njog],na.rm=T)*co2scl+co2off+.1,0)
	abline(mean(shapesco2[17,(skip+1):njog],na.rm=T)*co2scl+co2off-.1,0)
	abline(v=skip+.5)

}
if(shapescnt[18]!=0&any(!is.na(shapesco2[18,]))){
	plot(c(1:njog),shapesco2[18,]*co2scl+co2off,main=paste('All Cal. Lines, gas=',round(shapesngas[18],2),
     ' lin=',round(shapesnlin[18],2),' n=',shapescnt[18],sep=''),xlab='',ylab='approx. ppm',type='n',
     ylim=c(-5,5)+mean(shapesco2[18,(skip+1):njog],na.rm=T)*co2scl+co2off)
	points(c(1:njog),shapesco2[18,]*co2scl+co2off,col=2)
	abline(v=skip+.5)
	abline(mean(shapesco2[18,(skip+1):njog],na.rm=T)*co2scl+co2off+1,0)
	abline(mean(shapesco2[18,(skip+1):njog],na.rm=T)*co2scl+co2off-1,0)

	plot(c(1:njog),shapesco2[18,]*co2scl+co2off,main=paste('All Cal. Lines, gas=',round(shapesngas[18],2),
     ' lin=',round(shapesnlin[18],2),' n=',shapescnt[18],' (zoomed)',sep=''),xlab='',ylab='approx. ppm',type='n',
     ylim=c(-.5,.5)+mean(shapesco2[18,(skip+1):njog],na.rm=T)*co2scl+co2off)
	points(c(1:njog),shapesco2[18,]*co2scl+co2off,col=2)
	temp1<-round(mean(output$slco2[output$shpbn>6&output$shpbn<11]/ntrd*60,na.rm=T),2)
	temp2<-round(mean(output$sdco2[output$shpbn>6&output$shpbn<11],na.rm=T),2)
	text(njog*.67,mean(shapesco2[18,(skip+1):njog],na.rm=T)*co2scl+co2off-.25,paste('Ave. Slope = ',temp1,' ppm/min'))
	text(njog*.67,mean(shapesco2[18,(skip+1):njog],na.rm=T)*co2scl+co2off-.8,paste('Ave. Std. Dev. = ',temp2,' ppm'))
	abline(mean(shapesco2[18,(skip+1):njog],na.rm=T)*co2scl+co2off+.1,0)
	abline(mean(shapesco2[18,(skip+1):njog],na.rm=T)*co2scl+co2off-.1,0)
	abline(v=skip+.5)
}
mtext(paste(file,outflag,'  CO2 Shapes Averaged by Sample or Calibration',sep=''),3,0.5,T)
mtext(paste('Jog Position (',reso,' second resolution)',sep=''),1,1,T)

dev.off()

if(png){png(paste(file,outflag,"_prsshapes.png",sep=''),width=950,height=550)} else {
bitmap(paste(file,outflag,"_prsshapes.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)}

# reset vector of shape bins for pressure:
binorder<-NULL
for(i in c(1:6)){
	if(shapescnt[i]!=0&any(!is.na(shapesprs[i,]))){
		binorder<-c(binorder,i)
	}
}
binorder<-c(binorder,7:10,13:16,12:11)

# pressure (LOC_YYMMDDf_3.png):
par(mfrow=c(4,4))
par(oma=c(2.5,0,2,0))
par(mar=c(2,3,2,1)+.1)
par(mgp=c(2,1,0))
for(i in binorder){
	plot(c(1:njog),shapesprs[i,],main=paste(tlist[i],' gas=',round(shapesngas[i],2),
     ' lin=',round(shapesnlin[i],2),' n=',shapescnt[i],sep=''),xlab='',ylab='kPa',type='n',
     ylim=c(-2.5,5)+mean(shapesprs[i,(skip+1):njog],na.rm=T))
	points(c(1:njog),shapesprs[i,],col=3)
	temp<-round(mean(output$slprs[output$shpbn==i]/ntrd*60,na.rm=T),3)
	if(i==8){ temp<-round(mean(output$slprs[output$shpbn==8|output$shpbn==13|output$shpbn==14]/ntrd*60,na.rm=T),3)}
	if(i==9){ temp<-round(mean(output$slprs[output$shpbn==9|output$shpbn==15|output$shpbn==16]/ntrd*60,na.rm=T),3)}
	text(njog*.67,mean(shapesprs[i,(skip+1):njog],na.rm=T)+2.5,paste('Ave. Slope = ',temp,' kPa/min'))
	temp2<-round(mean(output$prs[output$shpbn==i],na.rm=T),2)
	if(i==8){ temp2<-round(mean(output$prs[output$shpbn==8|output$shpbn==13|output$shpbn==14],na.rm=T),2) }
	if(i==9){ temp2<-round(mean(output$prs[output$shpbn==9|output$shpbn==15|output$shpbn==16],na.rm=T),2) }
	text(njog*.67,mean(shapesprs[i,(skip+1):njog],na.rm=T)+3.5,paste('Ave. Press. = ',temp2,' kPa'))
	abline(mean(shapesprs[i,(skip+1):njog],na.rm=T)+1,0)
	abline(mean(shapesprs[i,(skip+1):njog],na.rm=T)-1,0)
	abline(v=skip+.5)
}
mtext(paste(file,outflag,'  Pressure Shapes',sep=''),3,0.5,T)
mtext(paste('Jog Position (',reso,' second resolution)',sep=''),1,1,T)

dev.off()

if(png){png(paste(file,outflag,"_sflshapes.png",sep=''),width=950,height=550)} else {
bitmap(paste(file,outflag,"_sflshapes.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)}

# reset vector of shape bins for sample flow:
binorder<-NULL
for(i in c(1:6)){
	if(shapescnt[i]!=0){
		binorder<-c(binorder,i)
	}
}
binorder<-c(binorder,7:10,13:16,12:11)

# sample flow (LOC_YYMMDDf_4.png):
par(mfrow=c(4,4))
par(oma=c(2.5,0,2,0))
par(mar=c(2,3,2,1)+.1)
par(mgp=c(2,1,0))
for(i in binorder){
	plot(c(1:njog),shapessfl[i,],main=paste(tlist[i],' gas=',round(shapesngas[i],2),
     ' lin=',round(shapesnlin[i],2),' n=',shapescnt[i],sep=''),xlab='',ylab='ml/min',type='n',ylim=c(-25,200))
	points(c(1:njog),shapessfl[i,],col=4)
	temp2<-round(mean(output$sfl[output$shpbn==i],na.rm=T),1)
	if(i==8){ temp2<-round(mean(output$sfl[output$shpbn==8|output$shpbn==13|output$shpbn==14],na.rm=T),1) }
	if(i==9){ temp2<-round(mean(output$sfl[output$shpbn==9|output$shpbn==15|output$shpbn==16],na.rm=T),1) }
	text(njog*.67,150,paste('Ave. Flow = ',temp2,' ml/min'))
	abline(mean(shapessfl[i,(skip+1):njog],na.rm=T)+10,0)
	abline(mean(shapessfl[i,(skip+1):njog],na.rm=T)-10,0)
	abline(v=skip+.5)
}
mtext(paste(file,outflag,'  Sample Flow Shapes',sep=''),3,0.5,T)
mtext(paste('Jog Position (',reso,' second resolution)',sep=''),1,1,T)

dev.off()

if(png){png(paste(file,outflag,"_lflshapes.png",sep=''),width=950,height=550)} else {
bitmap(paste(file,outflag,"_lflshapes.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)}

# reset vector of shape bins for line flow:
binorder<-NULL
for(i in c(1:6)){
	if(shapescnt[i]!=0){
		binorder<-c(binorder,i)
	}
}

# line flow (LOC_YYMMDDf_5.png):
par(mfrow=c(2,3))
par(oma=c(2.5,0,2,0))
par(mar=c(2,3,2,1)+.1)
par(mgp=c(2,1,0))
for(i in binorder){
	plot(c(1:njog),shapeslfl[i,],main=paste(tlist[i],' gas=',round(shapesngas[i],2),
     ' lin=',round(shapesnlin[i],2),' n=',shapescnt[i],sep=''),xlab='',ylab='ml/min',type='n',ylim=c(-25,200))
	points(c(1:njog),shapeslfl[i,],col=5)
	if(i==1){ text(njog*.67,150,paste('Ave. Flow = ',round(mean(output$fl1[output$shpbn==i],na.rm=T),1),' ml/min')) }
	if(i==2){ text(njog*.67,150,paste('Ave. Flow = ',round(mean(output$fl2[output$shpbn==i],na.rm=T),1),' ml/min')) }
	if(i==3){ text(njog*.67,150,paste('Ave. Flow = ',round(mean(output$fl3[output$shpbn==i],na.rm=T),1),' ml/min')) }
	if(i==4){ text(njog*.67,150,paste('Ave. Flow = ',round(mean(output$fl4[output$shpbn==i]),1),' ml/min')) }
	if(i==5){ text(njog*.67,150,paste('Ave. Flow = ',round(mean(output$fl5[output$shpbn==i]),1),' ml/min')) }
	if(i==6){ text(njog*.67,150,paste('Ave. Flow = ',round(mean(output$fl6[output$shpbn==i]),1),' ml/min')) }
	abline(mean(shapeslfl[i,(skip+1):njog],na.rm=T)+10,0)
	abline(mean(shapeslfl[i,(skip+1):njog],na.rm=T)-10,0)
	abline(v=skip+.5)
}
mtext(paste(file,outflag,'  Line Flow Shapes',sep=''),3,0.5,T)
mtext(paste('Jog Position (',reso,' second resolution)',sep=''),1,1,T)

dev.off()

if(png){png(paste(file,outflag,"_hrdiag.png",sep=''),width=950,height=550)} else {
bitmap(paste(file,outflag,"_hrdiag.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)}

## plot high rate diagnostics:

# plot full day (LOC_YYMMDDf_6.png):
par(mfrow=c(4,3))
par(oma=c(2.5,0,2,0))
par(mar=c(2,3,2,1)+.1)
par(mgp=c(2,1,0))

# setup title list for flows:
tlistf<-c('Line 1','Line 2','Line 3','Line 4','Line 5','Line 6')
tlistf[ltlin]<-paste(tlistf[ltlin],' (LT)',sep='')
tlistf<-paste(tlistf,' Flow',sep='')

plot(hodfrachr,data$co2,main='Li-820 CO2',ylab='Licor ppm',type='n',xlim=c(0,24),xaxp=c(0,24,6))
lines(hodfrachr,data$co2,col=2)
plot(hodfrachr,data$prs,main='Li-820 Pressure',ylab='kPa',type='n',xlim=c(0,24),xaxp=c(0,24,6))
lines(hodfrachr,data$prs,col=3)
plot(hodfrachr,data$ltmp,main='Li-820 Temperature',ylab='degrees C',type='n',xlim=c(0,24),xaxp=c(0,24,6))
lines(hodfrachr,data$ltmp,col=6)
plot(hodfrachr,data$btmp*100-40,main='Ambient Box Temperature',ylab='degrees C',type='n',xlim=c(0,24),xaxp=c(0,24,6))
lines(hodfrachr,data$btmp*100-40,col=6)
plot(hodfrachr,data$rh*100,main='Gas Stream RH',ylab='percent',type='n',xlim=c(0,24),xaxp=c(0,24,6))
lines(hodfrachr,data$rh*100,col='orange')
plot(hodfrachr,data$sfl,main='Sample Flow',ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6),ylim=c(-25,200))
lines(hodfrachr,data$sfl,col=4)
if(any(data$nlin==1&data$ngas==1)){
	plot(hodfrachr,data$fl1,main=tlistf[1],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6),ylim=c(-100,1000))
	lines(hodfrachr,data$fl1,col=5)
}
if(any(data$nlin==2&data$ngas==1)){
	plot(hodfrachr,data$fl2,main=tlistf[2],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6),ylim=c(-100,1000))
	lines(hodfrachr,data$fl2,col=5)
}
if(any(data$nlin==3&data$ngas==1)){
	plot(hodfrachr,data$fl3,main=tlistf[3],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6),ylim=c(-100,1000))
	lines(hodfrachr,data$fl3,col=5)
}
if(any(data$nlin==4&data$ngas==1)){
	plot(hodfrachr,data$fl4,main=tlistf[4],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6),ylim=c(-100,1000))
	lines(hodfrachr,data$fl4,col=5)
}
if(any(data$nlin==5&data$ngas==1)){
	plot(hodfrachr,data$fl5,main=tlistf[5],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6),ylim=c(-100,1000))
	lines(hodfrachr,data$fl5,col=5)
}
if(any(data$nlin==6&data$ngas==1)){
	plot(hodfrachr,data$fl6,main=tlistf[6],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6),ylim=c(-100,1000))
	lines(hodfrachr,data$fl6,col=5)
}
mtext(paste(file,outflag,'  ',1/reso,' Hz Diagnostics',sep=''),3,0.5,T)
mtext('Hour of Day (UTC)',1,1,T)

dev.off()

if(png){png(paste(file,outflag,"_hrdiag_zmam.png",sep=''),width=950,height=550)} else {
bitmap(paste(file,outflag,"_hrdiag_zmam.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)}

## plot zoomed looks at data:

# determine zoomed periods in AM and PM:
ammin<-(max(domfrachr)-min(domfrachr))*11.5/24+min(domfrachr) # 11:30 GMT (4:30 AM MDT) in full day record
ammax<-(max(domfrachr)-min(domfrachr))*12/24+min(domfrachr) # 12:00 GMT (5:00 AM MDT) in full day record
pmmin<-(max(domfrachr)-min(domfrachr))*23.5/24+min(domfrachr) # 23:30 GMT (4:30 PM MDT) in full day record
pmmax<-(max(domfrachr)-min(domfrachr))*24/24+min(domfrachr) # 24:00 GMT (5:00 PM MDT) in full day record
## 6/09 converted plots to hour of day, could convert AM PM filters too

# plot zoomed am period (LOC_YYMMDDf_7.png):
if(any(domfrachr>ammin&domfrachr<ammax)){ # data available from zoomed am period, ok to proceed

	par(mfrow=c(4,3))
	par(oma=c(2.5,0,2,0))
	par(mar=c(2,3,2,1)+.1)
	par(mgp=c(2,1,0))

	if(any(!is.na(data$co2[domfrachr>ammin&domfrachr<ammax]))){ # Li-820 data available from zoomed am period
		plot(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$co2[domfrachr>ammin&domfrachr<ammax],main='Li-820 CO2',ylab='Licor ppm',type='n')
		lines(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$co2[domfrachr>ammin&domfrachr<ammax],col=2)
		plot(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$prs[domfrachr>ammin&domfrachr<ammax],main='Li-820 Pressure',ylab='kPa',type='n')
		lines(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$prs[domfrachr>ammin&domfrachr<ammax],col=3)
		plot(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$ltmp[domfrachr>ammin&domfrachr<ammax],main='Li-820 Temperature',ylab='degrees C',type='n')
		lines(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$ltmp[domfrachr>ammin&domfrachr<ammax],col=6)
	}
	plot(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$btmp[domfrachr>ammin&domfrachr<ammax]*100-40,main='Ambient Box Temperature',ylab='degrees C',type='n')
	lines(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$btmp[domfrachr>ammin&domfrachr<ammax]*100-40,col=6)
	plot(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$rh[domfrachr>ammin&domfrachr<ammax]*100,main='Gas Stream RH',ylab='percent',type='n')
	lines(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$rh[domfrachr>ammin&domfrachr<ammax]*100,col='orange')
	plot(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$sfl[domfrachr>ammin&domfrachr<ammax],main='Sample Flow',ylab='ml/min',type='n',ylim=c(-25,200))
	lines(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$sfl[domfrachr>ammin&domfrachr<ammax],col=4)
	if(any(data$nlin==1&data$ngas==1)){
		plot(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$fl1[domfrachr>ammin&domfrachr<ammax],main=tlistf[1],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$fl1[domfrachr>ammin&domfrachr<ammax],col=5)
	}
	if(any(data$nlin==2&data$ngas==1)){
		plot(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$fl2[domfrachr>ammin&domfrachr<ammax],main=tlistf[2],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$fl2[domfrachr>ammin&domfrachr<ammax],col=5)
	}
	if(any(data$nlin==3&data$ngas==1)){
		plot(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$fl3[domfrachr>ammin&domfrachr<ammax],main=tlistf[3],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$fl3[domfrachr>ammin&domfrachr<ammax],col=5)
	}
	if(any(data$nlin==4&data$ngas==1)){
		plot(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$fl4[domfrachr>ammin&domfrachr<ammax],main=tlistf[4],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$fl4[domfrachr>ammin&domfrachr<ammax],col=5)
	}
	if(any(data$nlin==5&data$ngas==1)){
		plot(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$fl5[domfrachr>ammin&domfrachr<ammax],main=tlistf[5],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$fl5[domfrachr>ammin&domfrachr<ammax],col=5)
	}
	if(any(data$nlin==6&data$ngas==1)){
		plot(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$fl6[domfrachr>ammin&domfrachr<ammax],main=tlistf[6],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$fl6[domfrachr>ammin&domfrachr<ammax],col=5)
	}
} else {
	plot(c(1:3),c(1:3),type='n',bty='n',axes=F,xlab='',ylab='') # create a blank plot
}
#mtext(paste(file,outflag,'  ',1/reso,' Hz Diagnostics Zoomed ',round((ammin-trunc(ammin))*24,2),'-',round((ammax-trunc(ammax))*24,2),sep=''),3,0.5,T)
mtext(paste(file,outflag,'  ',1/reso,' Hz Diagnostics Zoomed AM',sep=''),3,0.5,T)
mtext('Hour of Day (UTC)',1,1,T)

dev.off()

if(png){png(paste(file,outflag,"_hrdiag_zmpm.png",sep=''),width=950,height=550)} else {
bitmap(paste(file,outflag,"_hrdiag_zmpm.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)}

# plot zoomed pm period (LOC_YYMMDDf_8.png):
if(any(domfrachr>pmmin&domfrachr<pmmax)){ # data available from zoomed pm period, ok to proceed

	par(mfrow=c(4,3))
	par(oma=c(2.5,0,2,0))
	par(mar=c(2,3,2,1)+.1)
	par(mgp=c(2,1,0))

	if(any(!is.na(data$co2[domfrachr>pmmin&domfrachr<pmmax]))){
		plot(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$co2[domfrachr>pmmin&domfrachr<pmmax],main='Li-820 CO2',ylab='Licor ppm',type='n')
		lines(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$co2[domfrachr>pmmin&domfrachr<pmmax],col=2)
		plot(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$prs[domfrachr>pmmin&domfrachr<pmmax],main='Li-820 Pressure',ylab='kPa',type='n')
		lines(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$prs[domfrachr>pmmin&domfrachr<pmmax],col=3)
		plot(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$ltmp[domfrachr>pmmin&domfrachr<pmmax],main='Li-820 Temperature',ylab='degrees C',type='n')
		lines(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$ltmp[domfrachr>pmmin&domfrachr<pmmax],col=6)
	}
	plot(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$btmp[domfrachr>pmmin&domfrachr<pmmax]*100-40,main='Ambient Box Temperature',ylab='degrees C',type='n')
	lines(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$btmp[domfrachr>pmmin&domfrachr<pmmax]*100-40,col=6)
	plot(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$rh[domfrachr>pmmin&domfrachr<pmmax]*100,main='Gas Stream RH',ylab='percent',type='n')
	lines(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$rh[domfrachr>pmmin&domfrachr<pmmax]*100,col='orange')
	plot(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$sfl[domfrachr>pmmin&domfrachr<pmmax],main='Sample Flow',ylab='ml/min',type='n',ylim=c(-25,200))
	lines(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$sfl[domfrachr>pmmin&domfrachr<pmmax],col=4)
	if(any(data$nlin==1&data$ngas==1)){
		plot(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$fl1[domfrachr>pmmin&domfrachr<pmmax],main=tlistf[1],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$fl1[domfrachr>pmmin&domfrachr<pmmax],col=5)
	}
	if(any(data$nlin==2&data$ngas==1)){
		plot(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$fl2[domfrachr>pmmin&domfrachr<pmmax],main=tlistf[2],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$fl2[domfrachr>pmmin&domfrachr<pmmax],col=5)
	}
	if(any(data$nlin==3&data$ngas==1)){
		plot(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$fl3[domfrachr>pmmin&domfrachr<pmmax],main=tlistf[3],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$fl3[domfrachr>pmmin&domfrachr<pmmax],col=5)
	}
	if(any(data$nlin==4&data$ngas==1)){
		plot(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$fl4[domfrachr>pmmin&domfrachr<pmmax],main=tlistf[4],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$fl4[domfrachr>pmmin&domfrachr<pmmax],col=5)
	}
	if(any(data$nlin==5&data$ngas==1)){
		plot(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$fl5[domfrachr>pmmin&domfrachr<pmmax],main=tlistf[5],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$fl5[domfrachr>pmmin&domfrachr<pmmax],col=5)
	}
	if(any(data$nlin==6&data$ngas==1)){
		plot(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$fl6[domfrachr>pmmin&domfrachr<pmmax],main=tlistf[6],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$fl6[domfrachr>pmmin&domfrachr<pmmax],col=5)
	}
} else {
	plot(c(1:3),c(1:3),type='n',bty='n',axes=F,xlab='',ylab='') # create a blank plot
}
#mtext(paste(file,outflag,'  ',1/reso,' Hz Diagnostics Zoomed ',round((pmmin-trunc(pmmin))*24,2),'-',round((pmmax-trunc(pmmax))*24,2),sep=''),3,0.5,T)
mtext(paste(file,outflag,'  ',1/reso,' Hz Diagnostics Zoomed PM',sep=''),3,0.5,T)
mtext('Hour of Day (UTC)',1,1,T)

dev.off()

if(png){png(paste(file,outflag,"_hrdiag_zmco2.png",sep=''),width=950,height=550)} else {
bitmap(paste(file,outflag,"_hrdiag_zmco2.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)}

# plot zoomed am and pm raw co2 data (LOC_YYMMDDf_9.png):
par(mfrow=c(2,1))
par(oma=c(2.5,0,2,0))
par(mar=c(2,3,2,1)+.1)
par(mgp=c(2,1,0))

if(any(domfrachr>ammin&domfrachr<ammax)&any(!is.na(data$co2[domfrachr>ammin&domfrachr<ammax]))){ # data available from zoomed am period, ok to proceed
#	plot(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$co2[domfrachr>ammin&domfrachr<ammax],main=paste(round((ammin-trunc(ammin))*24,2),'-',round((ammax-trunc(ammax))*24,2)),ylab='Licor ppm',type='n')
	plot(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$co2[domfrachr>ammin&domfrachr<ammax],main='AM',ylab='Licor ppm',type='n')
	points(hodfrachr[domfrachr>ammin&domfrachr<ammax],data$co2[domfrachr>ammin&domfrachr<ammax],col=2)
} else {
	plot(c(1:3),c(1:3),type='n',bty='n',axes=F,xlab='',ylab='') # create a blank plot
}
if(any(domfrachr>pmmin&domfrachr<pmmax)&any(!is.na(data$co2[domfrachr>pmmin&domfrachr<pmmax]))){ # data available from zoomed pm period, ok to proceed
#	plot(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$co2[domfrachr>pmmin&domfrachr<pmmax],main=paste(round((pmmin-trunc(pmmin))*24,2),'-',round((pmmax-trunc(pmmax))*24,2)),ylab='Licor ppm',type='n')
	plot(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$co2[domfrachr>pmmin&domfrachr<pmmax],main='PM',ylab='Licor ppm',type='n')
	points(hodfrachr[domfrachr>pmmin&domfrachr<pmmax],data$co2[domfrachr>pmmin&domfrachr<pmmax],col=2)
} else {
	plot(c(1:3),c(1:3),type='n',bty='n',axes=F,xlab='',ylab='') # create a blank plot
}
mtext(paste(file,outflag,'  ',1/reso,' Hz CO2 Zoomed',sep=''),3,0.5,T)
mtext('Hour of Day (UTC)',1,1,T)

dev.off()

if(png){png(paste(file,outflag,"_avediag.png",sep=''),width=950,height=550)} else {
bitmap(paste(file,outflag,"_avediag.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)}

## plot ~2.5-minute average diagnostics:

# plot averages (LOC_YYMMDDf_10.png):
par(mfrow=c(4,3))
par(oma=c(2.5,0,2,0))
par(mar=c(2,3,2,1)+.1)
par(mgp=c(2,1,0))

plot(hodfrac[output$lchk==0],output$co2[output$lchk==0],main='Li-820 CO2',ylab='Licor ppm',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1&output$lchk==0],output$co2[output$ngas==1&output$lchk==0],col=2)
points(hodfrac[output$ngas>1],output$co2[output$ngas>1],col=1)
plot(hodfrac[output$lchk==0],output$prs[output$lchk==0],main='Li-820 Pressure',ylab='kPa',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1&output$lchk==0],output$prs[output$ngas==1&output$lchk==0],col=3)
points(hodfrac[output$ngas>1],output$prs[output$ngas>1],col=1)
plot(hodfrac[output$lchk==0],output$ltmp[output$lchk==0],main='Li-820 Temperature',ylab='degrees C',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1&output$lchk==0],output$ltmp[output$ngas==1&output$lchk==0],col=6)
points(hodfrac[output$ngas>1],output$ltmp[output$ngas>1],col=1)
plot(hodfrac[output$lchk==0],output$btmp[output$lchk==0]*100-40,main='Ambient Box Temperature',ylab='degrees C',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1&output$lchk==0],output$btmp[output$ngas==1&output$lchk==0]*100-40,col=6)
points(hodfrac[output$ngas>1],output$btmp[output$ngas>1]*100-40,col=1)
plot(hodfrac[output$lchk==0],output$rh[output$lchk==0]*100,main='Gas Stream RH',ylab='percent',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1&output$lchk==0],output$rh[output$ngas==1&output$lchk==0]*100,col='orange')
points(hodfrac[output$ngas>1],output$rh[output$ngas>1]*100,col=1)
plot(hodfrac[output$lchk==0],output$sfl[output$lchk==0],main='Sample Flow',ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6),ylim=c(-25,200))
points(hodfrac[output$ngas==1&output$lchk==0],output$sfl[output$ngas==1&output$lchk==0],col=4)
points(hodfrac[output$ngas>1],output$sfl[output$ngas>1],col=1)
if(any(data$nlin==1&data$ngas==1)){
	plot(hodfrac[output$lchk==0],output$fl1[output$lchk==0],main=tlistf[1],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6),ylim=c(-100,1000))
	points(hodfrac[output$ngas==1&output$lchk==0],output$fl1[output$ngas==1&output$lchk==0],col=5)
	points(hodfrac[output$ngas>1],output$fl1[output$ngas>1],col=1)
}
if(any(data$nlin==2&data$ngas==1)){
	plot(hodfrac[output$lchk==0],output$fl2[output$lchk==0],main=tlistf[2],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6),ylim=c(-100,1000))
	points(hodfrac[output$ngas==1&output$lchk==0],output$fl2[output$ngas==1&output$lchk==0],col=5)
	points(hodfrac[output$ngas>1],output$fl2[output$ngas>1],col=1)
}
if(any(data$nlin==3&data$ngas==1)){
	plot(hodfrac[output$lchk==0],output$fl3[output$lchk==0],main=tlistf[3],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6),ylim=c(-100,1000))
	points(hodfrac[output$ngas==1&output$lchk==0],output$fl3[output$ngas==1&output$lchk==0],col=5)
	points(hodfrac[output$ngas>1],output$fl3[output$ngas>1],col=1)
}
if(any(data$nlin==4&data$ngas==1)){
	plot(hodfrac[output$lchk==0],output$fl4[output$lchk==0],main=tlistf[4],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6),ylim=c(-100,1000))
	points(hodfrac[output$ngas==1&output$lchk==0],output$fl4[output$ngas==1&output$lchk==0],col=5)
	points(hodfrac[output$ngas>1],output$fl4[output$ngas>1],col=1)
}
if(any(data$nlin==5&data$ngas==1)){
	plot(hodfrac[output$lchk==0],output$fl5[output$lchk==0],main=tlistf[5],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6),ylim=c(-100,1000))
	points(hodfrac[output$ngas==1&output$lchk==0],output$fl5[output$ngas==1&output$lchk==0],col=5)
	points(hodfrac[output$ngas>1],output$fl5[output$ngas>1],col=1)
}
if(any(data$nlin==6&data$ngas==1)){
	plot(hodfrac[output$lchk==0],output$fl6[output$lchk==0],main=tlistf[6],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6),ylim=c(-100,1000))
	points(hodfrac[output$ngas==1&output$lchk==0],output$fl6[output$ngas==1&output$lchk==0],col=5)
	points(hodfrac[output$ngas>1],output$fl6[output$ngas>1],col=1)
}
mtext(paste(file,outflag,'  Average Diagnostics',sep=''),3,0.5,T)
mtext('Hour of Day (UTC)',1,1,T)

dev.off()

if(png){png(paste(file,outflag,"_sddiag.png",sep=''),width=950,height=550)} else {
bitmap(paste(file,outflag,"_sddiag.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)}

# plot standard deviations (LOC_YYMMDDf_11.png):
par(mfrow=c(4,3))
par(oma=c(2.5,0,2,0))
par(mar=c(2,3,2,1)+.1)
par(mgp=c(2,1,0))

# set up title list for standard deviations:
tlistsd<-paste('Stan. Dev. ',tlistf,sep='')

if(any(!is.na(output$sdco2[output$ngas!=0]))){
plot(hodfrac[output$ngas!=0],output$sdco2[output$ngas!=0],main='Stan. Dev. Li-820 CO2',ylab='Licor ppm',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1],output$sdco2[output$ngas==1],col=2)
points(hodfrac[output$ngas>1],output$sdco2[output$ngas>1],col=1)
}
if(any(!is.na(output$sdprs[output$ngas!=0]))){
plot(hodfrac[output$ngas!=0],output$sdprs[output$ngas!=0],main='Stan. Dev. Li-820 Pressure',ylab='kPa',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1],output$sdprs[output$ngas==1],col=3)
points(hodfrac[output$ngas>1],output$sdprs[output$ngas>1],col=1)
}
if(any(!is.na(output$sdltmp[output$ngas!=0]))){
plot(hodfrac[output$ngas!=0],output$sdltmp[output$ngas!=0],main='Stan. Dev. Li-820 Temp.',ylab='degrees C',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1],output$sdltmp[output$ngas==1],col=6)
points(hodfrac[output$ngas>1],output$sdltmp[output$ngas>1],col=1)
}
if(any(!is.na(output$sdbtmp[output$ngas!=0]))){
plot(hodfrac[output$ngas!=0],output$sdbtmp[output$ngas!=0]*100,main='Stan. Dev. Amb. Box Temp.',ylab='degrees C',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1],output$sdbtmp[output$ngas==1]*100,col=6)
points(hodfrac[output$ngas>1],output$sdbtmp[output$ngas>1]*100,col=1)
}
if(any(!is.na(output$sdrh[output$ngas!=0]))){
plot(hodfrac[output$ngas!=0],output$sdrh[output$ngas!=0]*100,main='Stan. Dev. Gas Stream RH',ylab='percent',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1],output$sdrh[output$ngas==1]*100,col='orange')
points(hodfrac[output$ngas>1],output$sdrh[output$ngas>1]*100,col=1)
}
if(any(!is.na(output$sdsfl[output$ngas!=0]))){
plot(hodfrac[output$ngas!=0],output$sdsfl[output$ngas!=0],main='Stan. Dev. Sample Flow',ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1],output$sdsfl[output$ngas==1],col=4)
points(hodfrac[output$ngas>1],output$sdsfl[output$ngas>1],col=1)
}
if(any(data$nlin==1&data$ngas==1)){
	plot(hodfrac[output$ngas!=0],output$sdfl1[output$ngas!=0],main=tlistsd[1],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6))
	points(hodfrac[output$ngas==1],output$sdfl1[output$ngas==1],col=5)
	points(hodfrac[output$ngas>1],output$sdfl1[output$ngas>1],col=1)
}
if(any(data$nlin==2&data$ngas==1)){
	plot(hodfrac[output$ngas!=0],output$sdfl2[output$ngas!=0],main=tlistsd[2],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6))
	points(hodfrac[output$ngas==1],output$sdfl2[output$ngas==1],col=5)
	points(hodfrac[output$ngas>1],output$sdfl2[output$ngas>1],col=1)
}
if(any(data$nlin==3&data$ngas==1)){
	plot(hodfrac[output$ngas!=0],output$sdfl3[output$ngas!=0],main=tlistsd[3],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6))
	points(hodfrac[output$ngas==1],output$sdfl3[output$ngas==1],col=5)
	points(hodfrac[output$ngas>1],output$sdfl3[output$ngas>1],col=1)
}
if(any(data$nlin==4&data$ngas==1)){
	plot(hodfrac[output$ngas!=0],output$sdfl4[output$ngas!=0],main=tlistsd[4],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6))
	points(hodfrac[output$ngas==1],output$sdfl4[output$ngas==1],col=5)
	points(hodfrac[output$ngas>1],output$sdfl4[output$ngas>1],col=1)
}
if(any(data$nlin==5&data$ngas==1)){
	plot(hodfrac[output$ngas!=0],output$sdfl5[output$ngas!=0],main=tlistsd[5],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6))
	points(hodfrac[output$ngas==1],output$sdfl5[output$ngas==1],col=5)
	points(hodfrac[output$ngas>1],output$sdfl5[output$ngas>1],col=1)
}
if(any(data$nlin==6&data$ngas==1)){
	plot(hodfrac[output$ngas!=0],output$sdfl6[output$ngas!=0],main=tlistsd[6],ylab='ml/min',type='n',xlim=c(0,24),xaxp=c(0,24,6))
	points(hodfrac[output$ngas==1],output$sdfl6[output$ngas==1],col=5)
	points(hodfrac[output$ngas>1],output$sdfl6[output$ngas>1],col=1)
}
mtext(paste(file,outflag,'  Std. Dev. Diagnostics',sep=''),3,0.5,T)
mtext('Hour of Day (UTC)',1,1,T)

dev.off()

if(png){png(paste(file,outflag,"_sldiag.png",sep=''),width=950,height=550)} else {
bitmap(paste(file,outflag,"_sldiag.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)}

# plot slopes (LOC_YYMMDDf_12.png):
par(mfrow=c(4,3))
par(oma=c(2.5,0,2,0))
par(mar=c(2,3,2,1)+.1)
par(mgp=c(2,1,0))

# set up title list for slopes:
tlistsl<-paste('Slope ',tlistf,sep='')

plot(hodfrac[output$ngas!=0],output$slco2[output$ngas!=0]/ntrd*60,main='Slope Li-820 CO2',ylab='Licor ppm per min.',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1],output$slco2[output$ngas==1]/ntrd*60,col=2)
points(hodfrac[output$ngas>1],output$slco2[output$ngas>1]/ntrd*60,col=1)
plot(hodfrac[output$ngas!=0],output$slprs[output$ngas!=0]/ntrd*60,main='Slope Li-820 Pressure',ylab='kPa per min.',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1],output$slprs[output$ngas==1]/ntrd*60,col=3)
points(hodfrac[output$ngas>1],output$slprs[output$ngas>1]/ntrd*60,col=1)
plot(hodfrac[output$ngas!=0],output$slltmp[output$ngas!=0]/ntrd*60,main='Slope Li-820 Temp.',ylab='degrees C per min.',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1],output$slltmp[output$ngas==1]/ntrd*60,col=6)
points(hodfrac[output$ngas>1],output$slltmp[output$ngas>1]/ntrd*60,col=1)
plot(hodfrac[output$ngas!=0],output$slbtmp[output$ngas!=0]*100/ntrd*60,main='Slope Amb. Box Temp.',ylab='degrees C per min.',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1],output$slbtmp[output$ngas==1]*100/ntrd*60,col=6)
points(hodfrac[output$ngas>1],output$slbtmp[output$ngas>1]*100/ntrd*60,col=1)
plot(hodfrac[output$ngas!=0],output$slrh[output$ngas!=0]*100/ntrd*60,main='Slope Gas Stream RH',ylab='percent per min.',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1],output$slrh[output$ngas==1]*100/ntrd*60,col='orange')
points(hodfrac[output$ngas>1],output$slrh[output$ngas>1]*100/ntrd*60,col=1)
plot(hodfrac[output$ngas!=0],output$slsfl[output$ngas!=0]/ntrd*60,main='Slope Sample Flow',ylab='ml/min per min.',type='n',xlim=c(0,24),xaxp=c(0,24,6))
points(hodfrac[output$ngas==1],output$slsfl[output$ngas==1]/ntrd*60,col=4)
points(hodfrac[output$ngas>1],output$slsfl[output$ngas>1]/ntrd*60,col=1)
if(any(data$nlin==1&data$ngas==1)){
	plot(hodfrac[output$ngas!=0],output$slfl1[output$ngas!=0]/ntrd*60,main=tlistsl[1],ylab='ml/min per min.',type='n',xlim=c(0,24),xaxp=c(0,24,6))
	points(hodfrac[output$ngas==1],output$slfl1[output$ngas==1]/ntrd*60,col=5)
	points(hodfrac[output$ngas>1],output$slfl1[output$ngas>1]/ntrd*60,col=1)
}
if(any(data$nlin==2&data$ngas==1)){
	plot(hodfrac[output$ngas!=0],output$slfl2[output$ngas!=0]/ntrd*60,main=tlistsl[2],ylab='ml/min per min.',type='n',xlim=c(0,24),xaxp=c(0,24,6))
	points(hodfrac[output$ngas==1],output$slfl2[output$ngas==1]/ntrd*60,col=5)
	points(hodfrac[output$ngas>1],output$slfl2[output$ngas>1]/ntrd*60,col=1)
}
if(any(data$nlin==3&data$ngas==1)){
	plot(hodfrac[output$ngas!=0],output$slfl3[output$ngas!=0]/ntrd*60,main=tlistsl[3],ylab='ml/min per min.',type='n',xlim=c(0,24),xaxp=c(0,24,6))
	points(hodfrac[output$ngas==1],output$slfl3[output$ngas==1]/ntrd*60,col=5)
	points(hodfrac[output$ngas>1],output$slfl3[output$ngas>1]/ntrd*60,col=1)
}
if(any(data$nlin==4&data$ngas==1)){
	plot(hodfrac[output$ngas!=0],output$slfl4[output$ngas!=0]/ntrd*60,main=tlistsl[4],ylab='ml/min per min.',type='n',xlim=c(0,24),xaxp=c(0,24,6))
	points(hodfrac[output$ngas==1],output$slfl4[output$ngas==1]/ntrd*60,col=5)
	points(hodfrac[output$ngas>1],output$slfl4[output$ngas>1]/ntrd*60,col=1)
}
if(any(data$nlin==5&data$ngas==1)){
	plot(hodfrac[output$ngas!=0],output$slfl5[output$ngas!=0]/ntrd*60,main=tlistsl[5],ylab='ml/min per min.',type='n',xlim=c(0,24),xaxp=c(0,24,6))
	points(hodfrac[output$ngas==1],output$slfl5[output$ngas==1]/ntrd*60,col=5)
	points(hodfrac[output$ngas>1],output$slfl5[output$ngas>1]/ntrd*60,col=1)
}
if(any(data$nlin==6&data$ngas==1)){
	plot(hodfrac[output$ngas!=0],output$slfl6[output$ngas!=0]/ntrd*60,main=tlistsl[6],ylab='ml/min per min.',type='n',xlim=c(0,24),xaxp=c(0,24,6))
	points(hodfrac[output$ngas==1],output$slfl6[output$ngas==1]/ntrd*60,col=5)
	points(hodfrac[output$ngas>1],output$slfl6[output$ngas>1]/ntrd*60,col=1)
}
mtext(paste(file,outflag,'  Slope Diagnostics',sep=''),3,0.5,T)
mtext('Hour of Day (UTC)',1,1,T)

dev.off()

if(png){png(paste(file,outflag,"_co2ts.png",sep=''),width=950,height=550)} else {
bitmap(paste(file,outflag,"_co2ts.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)}

# plot co2 timeseries (LOC_YYMMDDf_13.png): 
par(mfrow=c(2,1))
par(oma=c(2.5,0,2,0))
par(mar=c(2,3,2,1)+.1)
par(mgp=c(2,1,0))

plot(hodfrachr,data$co2*co2scl+co2off,type='n',ylab='appox. ppm',xlab='day of month',xlim=c(0,24),xaxp=c(0,24,6),ylim=c(min(output$co2,na.rm=T)*co2scl+co2off,max(output$co2,na.rm=T)*co2scl+co2off))
lines(hodfrachr,data$co2*co2scl+co2off,col=8)
points(hodfrac[output$ngas>1],output$co2[output$ngas>1]*co2scl+co2off,col=1)
if(any(output$nlin==1&output$ngas==1)){
	points(hodfrac[output$nlin==1&output$ngas==1],output$co2[output$nlin==1&output$ngas==1]*co2scl+co2off,col=2)
}
if(any(output$nlin==2&output$ngas==1)){
	points(hodfrac[output$nlin==2&output$ngas==1],output$co2[output$nlin==2&output$ngas==1]*co2scl+co2off,col=3)
}
if(any(output$nlin==3&output$ngas==1)){
	points(hodfrac[output$nlin==3&output$ngas==1],output$co2[output$nlin==3&output$ngas==1]*co2scl+co2off,col=4)
}
if(any(output$nlin==4&output$ngas==1)){
	points(hodfrac[output$nlin==4&output$ngas==1],output$co2[output$nlin==4&output$ngas==1]*co2scl+co2off,col=5)
}
if(any(output$nlin==5&output$ngas==1)){
	points(hodfrac[output$nlin==5&output$ngas==1],output$co2[output$nlin==5&output$ngas==1]*co2scl+co2off,col=6)
}
if(any(output$nlin==6&output$ngas==1)){
	points(hodfrac[output$nlin==6&output$ngas==1],output$co2[output$nlin==6&output$ngas==1]*co2scl+co2off,col=7)
}

if(any(!is.na(output$co2[output$nlin>0&output$nlin<4&output$ngas==1]))){ # sample jogs exist, ok to proceed
	plot(hodfrac[output$nlin==1],output$co2[output$nlin==1]*co2scl+co2off,xlab='day of month',ylab='approx. ppm',type='n',ylim=
	c(min(output$co2[output$nlin>0&output$nlin<4&output$ngas==1],na.rm=T)*co2scl+co2off,max(output$co2[output$nlin>0&output$nlin<4&output$ngas==1],na.rm=T)*co2scl+co2off),
#	xlim=c(min(hodfrac,na.rm=T),max(hodfrac,na.rm=T)))
	xlim=c(0,24),xaxp=c(0,24,6))
	if(any(output$nlin==1&output$ngas==1&ltlin!=1)){
		points(hodfrac[output$nlin==1&output$ngas==1],output$co2[output$nlin==1&output$ngas==1]*co2scl+co2off,type='b',pch='1',col=2)
	}
	if(any(output$nlin==2&output$ngas==1&ltlin!=2)){
		points(hodfrac[output$nlin==2&output$ngas==1],output$co2[output$nlin==2&output$ngas==1]*co2scl+co2off,type='b',pch='2',col=3)
	}
	if(any(output$nlin==3&output$ngas==1&ltlin!=3)){
		points(hodfrac[output$nlin==3&output$ngas==1],output$co2[output$nlin==3&output$ngas==1]*co2scl+co2off,type='b',pch='3',col=4)
	}
	if(any(output$nlin==4&output$ngas==1&ltlin!=4)){
		points(hodfrac[output$nlin==4&output$ngas==1],output$co2[output$nlin==4&output$ngas==1]*co2scl+co2off,type='b',pch='4',col=5)
	}
	if(any(output$nlin==5&output$ngas==1&ltlin!=5)){
		points(hodfrac[output$nlin==5&output$ngas==1],output$co2[output$nlin==5&output$ngas==1]*co2scl+co2off,type='b',pch='5',col=6)
	}
	if(any(output$nlin==6&output$ngas==1&ltlin!=6)){
		points(hodfrac[output$nlin==6&output$ngas==1],output$co2[output$nlin==6&output$ngas==1]*co2scl+co2off,type='b',pch='6',col=7)
	}
}

mtext(paste(file,outflag,'  CO2 Timeseries',sep=''),3,0.5,T)
mtext('Hour of Day (UTC)',1,1,T)

dev.off()

if(png){png(paste(file,outflag,"_cal_lc.png",sep=''),width=950,height=550)} else {
bitmap(paste(file,outflag,"_cal_lc.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)}

# plot calibration and LT values and leak check slopes (LOC_YYMMDDf_14.png):
par(mfrow=c(2,3))
par(oma=c(2.5,0,2,0))
par(mar=c(2,3,2,1)+.1)
par(mgp=c(2,1,0))

if(any(!is.na(output$co2[output$ngas==5]))){
	plot(hodfrac[output$ngas==5],output$co2[output$ngas==5]*co2scl+co2off,main='HS2',xlim=c(0,24),xaxp=c(0,24,6),ylab='approx. ppm',type='n', 
	xlab='sec',ylim=c(mean(output$co2[output$ngas==5],na.rm=T)*co2scl+co2off-2.5,mean(output$co2[output$ngas==5],na.rm=T)*co2scl+co2off+2.5))
	points(hodfrac[output$ngas==5],output$co2[output$ngas==5]*co2scl+co2off,col=2)
	abline(mean(output$co2[output$ngas==5],na.rm=T)*co2scl+co2off+1,0)
	abline(mean(output$co2[output$ngas==5],na.rm=T)*co2scl+co2off-1,0)
}
if(any(!is.na(output$co2[output$ngas==4]))){
	plot(hodfrac[output$ngas==4],output$co2[output$ngas==4]*co2scl+co2off,main='HS1',xlim=c(0,24),xaxp=c(0,24,6),ylab='approx. ppm',type='n', 
	xlab='sec',ylim=c(mean(output$co2[output$ngas==4],na.rm=T)*co2scl+co2off-2.5,mean(output$co2[output$ngas==4],na.rm=T)*co2scl+co2off+2.5))
	points(hodfrac[output$ngas==4],output$co2[output$ngas==4]*co2scl+co2off,col=2)
	abline(mean(output$co2[output$ngas==4],na.rm=T)*co2scl+co2off+1,0)
	abline(mean(output$co2[output$ngas==4],na.rm=T)*co2scl+co2off-1,0)
}
if(any(!is.na(output$co2[output$ngas==3]))){
	plot(hodfrac[output$ngas==3],output$co2[output$ngas==3]*co2scl+co2off,main='LS1',xlim=c(0,24),xaxp=c(0,24,6),ylab='approx. ppm',type='n', 
	xlab='sec',ylim=c(mean(output$co2[output$ngas==3],na.rm=T)*co2scl+co2off-2.5,mean(output$co2[output$ngas==3],na.rm=T)*co2scl+co2off+2.5))
	points(hodfrac[output$ngas==3],output$co2[output$ngas==3]*co2scl+co2off,col=2)
	abline(mean(output$co2[output$ngas==3],na.rm=T)*co2scl+co2off+1,0)
	abline(mean(output$co2[output$ngas==3],na.rm=T)*co2scl+co2off-1,0)
}
if(any(!is.na(output$co2[output$ngas==2]))){
	plot(hodfrac[output$ngas==2],output$co2[output$ngas==2]*co2scl+co2off,main='LS2',xlim=c(0,24),xaxp=c(0,24,6),ylab='approx. ppm',type='n', 
	xlab='sec',ylim=c(mean(output$co2[output$ngas==2],na.rm=T)*co2scl+co2off-2.5,mean(output$co2[output$ngas==2],na.rm=T)*co2scl+co2off+2.5))
	points(hodfrac[output$ngas==2],output$co2[output$ngas==2]*co2scl+co2off,col=2)
	abline(mean(output$co2[output$ngas==2],na.rm=T)*co2scl+co2off+1,0)
	abline(mean(output$co2[output$ngas==2],na.rm=T)*co2scl+co2off-1,0)
}
if(any(!is.na(output$co2[output$nlin==ltlin&output$ngas==1]))){
	plot(hodfrac[output$nlin==ltlin&output$ngas==1],output$co2[output$nlin==ltlin&output$ngas==1]*co2scl+co2off,main='LT',xlim=c(0,24),xaxp=c(0,24,6),ylab='approx. ppm',type='n', 
	xlab='sec',ylim=c(mean(output$co2[output$nlin==ltlin&output$ngas==1],na.rm=T)*co2scl+co2off-2.5,mean(output$co2[output$nlin==ltlin&output$ngas==1],na.rm=T)*co2scl+co2off+2.5))
	points(hodfrac[output$nlin==ltlin&output$ngas==1],output$co2[output$nlin==ltlin&output$ngas==1]*co2scl+co2off,col=2)
	abline(mean(output$co2[output$nlin==ltlin&output$ngas==1],na.rm=T)*co2scl+co2off+1,0)
	abline(mean(output$co2[output$nlin==ltlin&output$ngas==1],na.rm=T)*co2scl+co2off-1,0)
}
if(any(!is.na(output$slprs[output$lchk==1]))){
	plot(hodfrac[output$lchk==1],output$slprs[output$lchk==1]/ntrd*60,main='Leak Check Pressure Slopes',xlim=c(0,24),xaxp=c(0,24,6),ylab='kPa per min.',xlab='sec',type='n',ylim=c(min(c(-0.075,output$slprs[output$lchk==1]/ntrd*60),na.rm=T),max(c(0.075,output$slprs[output$lchk==1]/ntrd*60),na.rm=T)))
	abline(.05,0)
	abline(-.05,0)
	points(hodfrac[output$shpbn==11],output$slprs[output$shpbn==11]/ntrd*60,col=3)
	points(hodfrac[output$shpbn==12],output$slprs[output$shpbn==12]/ntrd*60,col=4)
} else {
	plot(hodfrac,output$slprs/ntrd*60,main='Leak Check Pressure Slopes',xlim=c(0,24),xaxp=c(0,24,6),ylab='kPa per min.',xlab='sec',type='n') # make a blank plot
}
legend(par('usr')[1],-0.05,c('High P','Low P'),pch=0,col=c(3,4))

mtext(paste(file,outflag,'  Calibration and LT Values and Leak Check Slopes',sep=''),3,0.5,T)
mtext('Hour of Day (UTC)',1,1,T)

dev.off()

} else { # do not run if all jogs bad

print('All jogs bad')
write(" ",paste(month,'/',loc,'/',file,outflag,".out",sep='')) # must output empty *.out file to prevent automated scripts from rerunning 
system(paste('chmod 664 ',month,'/',loc,'/',file,outflag,".out",sep='')) 

if(length(data$sec)>0){ # still plot high rate diagnostics if some raw data present

	#convert seconds into day of month + day fraction and handle midnight
	domfrachr<-as.integer(substr(day,5,6))+data$sec/86400
	hodfrachr<-data$sec/60/60
	# bitmap command necessary to run in batch mode (png gets better resolution)
	if(png){png(paste(file,outflag,"_hrdiag.png",sep=''),width=950,height=550)} else {
	bitmap(paste(file,outflag,"_hrdiag.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)}

	# setup title list for flows:
	tlistf<-c('Line 1','Line 2','Line 3','Line 4','Line 5','Line 6')
	tlistf[ltlin]<-paste(tlistf[ltlin],' (LT)',sep='')
	tlistf<-paste(tlistf,' Flow',sep='')

	# plot full day (LOC_YYMMDDf_6.png):
	par(mfrow=c(4,3))
	par(oma=c(2.5,0,2,0))
	par(mar=c(2,3,2,1)+.1)
	par(mgp=c(2,1,0))
	if(any(!is.na(data$co2))){
		plot(hodfrachr,data$co2,main='Li-820 CO2',ylab='Licor ppm',type='n')
		lines(hodfrachr,data$co2,col=2)
	}
	if(any(!is.na(data$prs))){
		plot(hodfrachr,data$prs,main='Li-820 Pressure',ylab='kPa',type='n')
		lines(hodfrachr,data$prs,col=3)
	}
	if(any(!is.na(data$ltmp))){
		plot(hodfrachr,data$ltmp,main='Li-820 Temperature',ylab='degrees C',type='n')
		lines(hodfrachr,data$ltmp,col=6)
	}
	if(any(!is.na(data$btmp))){
		plot(hodfrachr,data$btmp*100-40,main='Ambient Box Temperature',ylab='degrees C',type='n')
		lines(hodfrachr,data$btmp*100-40,col=6)
	}
	if(any(!is.na(data$rh))){
		plot(hodfrachr,data$rh*100,main='Gas Stream RH',ylab='percent',type='n')
		lines(hodfrachr,data$rh*100,col='orange')
	}
	if(any(!is.na(data$sfl))){
		plot(hodfrachr,data$sfl,main='Sample Flow',ylab='ml/min',type='n',ylim=c(-25,200))
		lines(hodfrachr,data$sfl,col=4)
	}
	if(any(!is.na(data$fl1)&data$nlin==1&data$ngas==1)){
		plot(hodfrachr,data$fl1,main=tlistf[1],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr,data$fl1,col=5)
	}
	if(any(!is.na(data$fl2)&data$nlin==2&data$ngas==1)){
		plot(hodfrachr,data$fl2,main=tlistf[2],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr,data$fl2,col=5)
	}
	if(any(!is.na(data$fl3)&data$nlin==3&data$ngas==1)){
		plot(hodfrachr,data$fl3,main=tlistf[3],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr,data$fl3,col=5)
	}
	if(any(!is.na(data$fl4)&data$nlin==4&data$ngas==1)){
		plot(hodfrachr,data$fl4,main=tlistf[4],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr,data$fl4,col=5)
	}
	if(any(!is.na(data$fl5)&data$nlin==5&data$ngas==1)){
		plot(hodfrachr,data$fl5,main=tlistf[5],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr,data$fl5,col=5)
	}
	if(any(!is.na(data$fl6)&data$nlin==6&data$ngas==1)){
		plot(hodfrachr,data$fl6,main=tlistf[6],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr,data$fl6,col=5)
	}
	mtext(paste(file,outflag,'  ',1/reso,' Hz Diagnostics',sep=''),3,0.5,T)
	mtext('Hour of Day (UTC)',1,1,T)
	dev.off()

} # do if some raw data present

} # do not run if all jogs bad

} else { # do not run if all data screened

	write(" ",paste(month,'/',loc,'/',file,outflag,".out",sep='')) # must output empty *.out file to prevent automated scripts from rerunning 
	system(paste('chmod 664 ',month,'/',loc,'/',file,outflag,".out",sep='')) 

if(length(data$vlvhx)>0){ # still plot high rate diagnostics if between 0 and njog values

	print('Almost all data screened')

	#convert seconds into day of month + day fraction and handle midnight
	domfrachr<-as.integer(substr(day,5,6))+data$sec/86400
	hodfrachr<-data$sec/60/60
	# bitmap command necessary to run in batch mode (png gets better resolution)
	if(png){png(paste(file,outflag,"_hrdiag.png",sep=''),width=950,height=550)} else {
	bitmap(paste(file,outflag,"_hrdiag.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)}

	# setup title list for flows:
	tlistf<-c('Line 1','Line 2','Line 3','Line 4','Line 5','Line 6')
	tlistf[ltlin]<-paste(tlistf[ltlin],' (LT)',sep='')
	tlistf<-paste(tlistf,' Flow',sep='')

	# plot diagnostics 
	par(mfrow=c(4,3))
	par(oma=c(2.5,0,2,0))
	par(mar=c(2,3,2,1)+.1)
	par(mgp=c(2,1,0))
	if(any(!is.na(data$co2))){
		plot(hodfrachr,data$co2,main='Li-820 CO2',ylab='Licor ppm',type='n')
		lines(hodfrachr,data$co2,col=2)
	}
	if(any(!is.na(data$prs))){
		plot(hodfrachr,data$prs,main='Li-820 Pressure',ylab='kPa',type='n')
		lines(hodfrachr,data$prs,col=3)
	}
	if(any(!is.na(data$ltmp))){
		plot(hodfrachr,data$ltmp,main='Li-820 Temperature',ylab='degrees C',type='n')
		lines(hodfrachr,data$ltmp,col=6)
	}
	if(any(!is.na(data$btmp))){
		plot(hodfrachr,data$btmp*100-40,main='Ambient Box Temperature',ylab='degrees C',type='n')
		lines(hodfrachr,data$btmp*100-40,col=6)
	}
	if(any(!is.na(data$rh))){
		plot(hodfrachr,data$rh*100,main='Gas Stream RH',ylab='percent',type='n')
		lines(hodfrachr,data$rh*100,col='orange')
	}
	if(any(!is.na(data$sfl))){
		plot(hodfrachr,data$sfl,main='Sample Flow',ylab='ml/min',type='n',ylim=c(-25,200))
		lines(hodfrachr,data$sfl,col=4)
	}
	if(any(!is.na(data$fl1)&data$nlin==1&data$ngas==1)){
		plot(hodfrachr,data$fl1,main=tlistf[1],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr,data$fl1,col=5)
	}
	if(any(!is.na(data$fl2)&data$nlin==2&data$ngas==1)){
		plot(hodfrachr,data$fl2,main=tlistf[2],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr,data$fl2,col=5)
	}
	if(any(!is.na(data$fl3)&data$nlin==3&data$ngas==1)){
		plot(hodfrachr,data$fl3,main=tlistf[3],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr,data$fl3,col=5)
	}
	if(any(!is.na(data$fl4)&data$nlin==4&data$ngas==1)){
		plot(hodfrachr,data$fl4,main=tlistf[4],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr,data$fl4,col=5)
	}
	if(any(!is.na(data$fl5)&data$nlin==5&data$ngas==1)){
		plot(hodfrachr,data$fl5,main=tlistf[5],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr,data$fl5,col=5)
	}
	if(any(!is.na(data$fl6)&data$nlin==6&data$ngas==1)){
		plot(hodfrachr,data$fl6,main=tlistf[6],ylab='ml/min',type='n',ylim=c(-100,1000))
		lines(hodfrachr,data$fl6,col=5)
	}
	mtext(paste(file,outflag,'  ',1/reso,' Hz Diagnostics',sep=''),3,0.5,T)
	mtext('Hour of Day (UTC)',1,1,T)
	dev.off()

} else {
	print('All data screened')
} # do if between 0 and njog values

} # do not run if all data screened

} else { # do not run if no data
	print('No raw data found')
}

# finish output and error logging:
print('aircoa_raw_day ending')
print(system('date',T))
sink(type="message")
sink()
close(log)

print(paste(system('date',T),paste('RAWDAY ending',loc,day)))

} # end of function