pro plotc_lc,startdate,days,exclude,screen,plot,plotstart,plotend

;plotc_lc.pro
;IDL program opens 3 @ YYDDDHHc files from Lost Creek (day of interest and previous/following days),
;performs LI-6251 calibration, outputs 3 minute averaged ascii file for all levels, and plots data
;and regressions equations.
;
;WARNING: This program attempts to spline all gaps!  Be careful if large chunks are missing!
;You can exclude a single period of time from any or all levels if desired (-999 will appear in the *.ppm file)
;
;NOTE: If IDL returns NaN for a spline proceedure, check the raw data for clock resets that would
;cause the normally increasing time array to go backwards! (this is not typical, because the datalogger
;clock is synchronized every hour and reset if more than 1 second different)
;
;written by Bruce Cook
;modified from plotc_wc 20 September 2000
;modified for unix 8 April 2001

IF n_elements(startdate) EQ 0 THEN begin
   print,'Enter start date in yyddd format: '
   startdate=''
   read,startdate

   print,'Enter number of days to process (do not over overlap Jan 1): '
   days=1S
   read,days

   print,'Do you want to exclude a specific time period for any level?: '
   exclude=''
   read,exclude

   if (exclude eq 'y') then begin
      print,'Which level (enter "1" through "5", or "all")?: '
      level=''
      read,level
      level=str_sep(level,',')
      print,'Enter start time (hhmm): '
      starttime=0S
      read,starttime
      print,'Enter end time (hhmm): '
      endtime=0S
      read,endtime
   endif

   print,'Do you want to view screen plots?: '
   screen=''
   read,screen
   if (screen eq 'y') then plotstart=0S else plotstart=1S

   print,'Do you want to make postscript plots and output files?: '
   plot=''
   read,plot
   if (plot eq 'y') then plotend=1S else plotend=0S
endif

;extract file name information
yy=strtrim(strmid(startdate,0,2),2)
IF (fix(yy) GE 98) THEN yyyy='19'+yy ELSE yyyy='20'+yy
startday=fix(strmid(startdate,2,3))
mm=strmid(strtrim(jdatetodate(startday,yyyy),2),4,2)
dd=strmid(strtrim(jdatetodate(startday,yyyy),2),6,2)
jdn_start=julday(mm,dd,yyyy); determine Julian Day Number for startday

;start loop for multiple days
for x=1,days do begin
jdn=jdn_start+(x-1); advance to next day

skip=0; reset error flag (when there is not enough data for calculations)

;determine month, day, and year of files to open
caldat,jdn-1,m1,d1,y1
caldat,jdn,m2,d2,y2
caldat,jdn+1,m3,d3,y3

;determine yy and ddd for files to open
yy1=strmid(strtrim(string(y1),2),2,2)
yy2=strmid(strtrim(string(y2),2),2,2)
yy3=strmid(strtrim(string(y3),2),2,2)

ddd1=jdate2(y1,m1,d1)
ddd=jdate2(y2,m2,d2)
ddd3=jdate2(y3,m3,d3)

;add zeroes to Julian dates if needed
if (ddd1 lt 10) then ddd1='00'+strtrim(string(ddd1),1) else $
	if (ddd1 lt 100) then ddd1='0'+strtrim(string(ddd1),1) else ddd1=strtrim(string(ddd1),1)
if (ddd lt 10) then ddd2='00'+strtrim(string(ddd),1) else $
	if (ddd lt 100) then ddd2='0'+strtrim(string(ddd),1) else ddd2=strtrim(string(ddd),1)
if (ddd3 lt 10) then ddd3='00'+strtrim(string(ddd3),1) else $
	if (ddd3 lt 100) then ddd3='0'+strtrim(string(ddd3),1) else ddd3=strtrim(string(ddd3),1)

file1=yy+ddd1+'00c'
file2=yy+ddd2+'00c'
file3=yy+ddd3+'00c'

;setup array for data to be plotted
nmax=1600; maximum number of rows in input file
datafinal=dblarr(8,nmax); 8 columns x nmax rows
datafinal(*,*)=-999.
line=''

;begin loop to open three files (day of interest, and previous/following days)
i=0; reset line loop counter
for j=1,3 do begin; begin loop for each of 3 files
if j eq 1 then filename=file1 else $
if j eq 2 then filename=file2 else filename=file3

openr,lun1,filepath(filename+'.gz',root_dir='/eddy/s4/lcreek/rawc'),/get_lun,error=err,/compress; open input data file
if (err ne 0) then begin; error if no 'c.gz' file
   filename=strtrim(strmid(filename,0,7),2)+'C.GZ'
   openr,lun1,filepath(filename,root_dir='/eddy/s4/lcreek/rawc/'),/get_lun,/compress,error=err; open input data file
ENDIF; end if no 'c.gz' file

if (err ne 0) then begin; error if file not found
	print,filename+ ' not found'
        IF (filename EQ file2) THEN BEGIN
           print,'WARNING!! NO DATA FOUND, SKIPPING CALCULATIONS...'
           skip=1
           GOTO, skipcalcs
        endif 
endif else begin; begin loading if file found
print,'Loading '+filename

;read data and create file for ploting
while (not eof(lun1)) do begin; begin loop
	readf,lun1,line; read first line of comma delimited data
	data=str_sep(line,','); separate line into data array

	;change selected data from string to floating point
 	jday=double(data(2)); Julian day
	time=double(data(3)); HHMM (0 to 2459)

	;create fractional julian date
	if time le 59 then begin; IF 2 digit value (i.e., <1 hr)
		fjday=jday+(time/1440)
	endif else begin
		if time le 959 then begin; IF 3 digit value
			fjday=jday+(double(strmid(data(3),0,1))/24)+(double(strmid(data(3),1,2))/1440)
		endif else begin; IF 4 digit value
			fjday=jday+(double(strmid(data(3),0,2))/24)+(double(strmid(data(3),2,2))/1440)
		endelse
	endelse

	datafinal(*,i)=[fjday,double(data(1:7))]; write selected data to file

	i=i+1; advance counter for next row

endwhile

;close files
close,lun1
free_lun,lun1

endelse; end if no file found

endfor; end open 3 files loop

;standard gas concentrations
;ORIGINAL STANDARDS (low, Toll #T389171; ref, Toll #T420298/22719724; high, Toll #HP-T032547) 
;   - Calibrated 10 August 1999 (low, high)
;   - Recalibrated 26 May 2000 (low, ref, high)
stds=[(345.68+344.87)/2,442.46,(551.39+551.51)/2]

;REFERENCE GAS REPLACED (Toll 31777) - 21 Oct 2001
;   - Initially calibrated ?? (439.66 ppm)
;   - Recalibrated 9 July 2002 (436.24 ppm)
IF ((yy eq 01) AND (ddd GE 294)) OR ((yy EQ 02) AND (ddd LT 314)) THEN $
   stds=[(345.68+344.87)/2.,(439.66+436.24)/2.,(551.39+551.51)/2.]

;REFERENCE GAS REPLACED (AGA C148663) - 10 Nov 2002
;   - Initially calibrated 11 Mar 2003 (426.44 ppm) using old method
;        stds=[(345.68+344.87)/2.,426.44,(551.39+551.51)/2.]
;   - All standards recalibrated 31 Aug 2003 w/4 @ NOAA stds
IF ((yy EQ 02) AND (ddd GE 314)) OR (yy eq 03) OR (yy EQ 04 AND ddd LT 127) THEN $
   stds=[345.28,426.44,551.45]

;REFERENCE GAS REPLACED (AGA C44429) - 6 May 2004
;   - calibrate 19 June 2004
IF (yy EQ 04 AND (ddd GE 127 AND ddd LT 171)) THEN $
   stds=[345.28,396.66,551.45]

;HIGH STANDARD REPLACED (AGA 503116) - 19 June 2004
;   - All standards calibrated 19 June 2004 w/4 @ NOAA stds
IF (yy EQ 04 AND (ddd GE 171 AND ddd LT 349)) THEN $
   stds=[344.23,396.66,528.92]

;LOW STANDARD REPLACED (AGA T95033644) - 14 December 2004
;   - Initially uncalibrated; 345 ppm on label
IF ((yy EQ 04 AND ddd GE 349) OR (yy EQ 05 AND ddd LT 175)) THEN $
   stds=[342.86,396.66,528.92]

;ALL STANDARDS RECALIBRATED - 27 June 2005
;   - calibrated using 4 @ NOAA stds
IF (yy EQ 05 AND (ddd GE 175 AND ddd LE 2005)) then $
   stds=[342.86,404.73,529.14]

;REFERENCE GAS REPLACED (C44429) - 27 July 2005
;   - ref gas calibrated 27 June 2005 using 4 @ NOAA stds
IF ((yy EQ 05 AND ddd GT 205) OR (yy EQ 06 AND ddd LE 161)) THEN $
   stds=[342.86,444.00,529.14]

;REFERENCE GAS REPLACED (AGA002473) - 10 June 2006
;   - calibrated using 4 @ NOAA standards 16 Aug 2006
IF (yy EQ 06 AND ddd GT 161 AND ddd LE 229) THEN $
   stds=[342.86,401.89,529.14]

;HIGH STANDARD REPLACED (AGA C148649) - 16 Aug 2006
;   - calibrated using 4 @ NOAA standards 16 Aug 2006
IF (yy EQ 06 AND ddd Gt 229) THEN $
   stds=[341.41,401.89,534.41]

lowstd=stds(0)
refgas=stds(1)
highstd=stds(2)

;set jday in previous year to 0, and jday in following year to 366 or 367
IF (ddd EQ 1) THEN BEGIN
   index=where(datafinal(0,*) GT ddd1)
   IF (index(0) NE -1) THEN datafinal(0,index)=datafinal(0,index)-ddd1
ENDIF

IF (fix(ddd3) EQ 1) THEN BEGIN
   index=where(datafinal(0,*) GT ddd1 AND datafinal(0,*) NE -999.)
   IF (index(0) NE -1) THEN datafinal(0,index)=datafinal(0,index)+ddd
endif

;locate most recent standards on previous and following days and trim dataset
id=datafinal(4,*)
jday=datafinal(2,*)
stds1=where(id eq 8 and jday eq double(strmid(file1,2,3)),count1)
stds2=where(id eq 8 and jday eq double(strmid(file2,2,3)),count2)
stds3=where(id eq 8 and jday eq double(strmid(file3,2,3)),count3)

IF count2 LT 1 THEN BEGIN; do not perform calcs if <1 std set for selected day
   print,'ERROR!! NO STANDARDS FOR SELECTED DAY, SKIPPING CALCULATIONS...'
   skip=1
   goto, skipcalcs
endif

if stds1(0) ne -1 then first_std=stds1(count1-1)-2 else first_std=stds2(0)-2; locate first standard set
if stds3(0) ne -1 then last_std=stds3(0) else last_std=stds2(count2-1)-2; located last standard set
datafinal=datafinal(*,first_std:last_std); keep only first/last standards and data inbetween
allstds=where(datafinal(4,*) eq 8,allcount); total number of standards

if allcount lt 3 then begin; do not perform calcs if <3 std sets for all days
   print,'ERROR!! NOT ENOUGH STANDARDS (n='+strtrim(string(allcount),2)+'), SKIPPING CALCULATIONS...'
   skip=1
   goto, skipcalcs
endif

;identify variables
jday=datafinal(2,*)
id=datafinal(4,*)
li_t=datafinal(5,*)
li_p=datafinal(6,*)
IF (yy EQ '01') AND ((ddd GE 193) AND (ddd LE 233)) THEN li_p=li_p+100
co2=datafinal(7,*)

;screen for extreme values (including -99999) and replace with -999
ili_t=where((li_t lt -50) or (li_t gt 75))
if (ili_t(0) ne -1) then li_t(ili_t)=-999

ili_p=where((li_p lt 94) or (li_p gt 120))
if (ili_p(0) ne -1) then li_p(ili_p)=-999

ico2=where((co2 lt -2500) or (co2 gt 2500))
if (ico2(0) ne -1) then co2(ico2)=-999

;locate bad values and change id to -999 to prevent use in calibrations and calculations
badval=where(li_t eq -999 or li_p eq -999 or co2 eq -999,cbadval)
if (cbadval ne 0) then id(badval)=-999

;locate samples, zero gas, standards, beginning and end of day
index1=where(id eq 1,count1);  4 feet above the soil surface
index2=where(id eq 2,count2);  6 feet
index3=where(id eq 3,count3); 10 feet
index4=where(id eq 4,count4); 18 feet
index5=where(id eq 5,count5); 33 feet, 8 inches (9.93 m)
index8=where(id eq 6,count8); reference gas
index9=where(id eq 7,count9); low standard
index10=where(id eq 8,count10); high standard
index8b=index10-2; reference gas when followed by standards
begind=where(jday eq double(strmid(file1,2,3)),begcount)
endind=where(jday eq double(strmid(file1,2,3)),endcount)

;skip calculations if no good standards
IF (count8 LT 3) OR (count9 LT 3) OR (count10 LT 3) THEN BEGIN
   print,'ERROR!! NOT ENOUGH GOOD STANDARDS, SKIPPING CALCULATIONS...'
   skip=1
   GOTO,skipcalcs
endif   

;create individual vectors for fractional Julian date (d#), LI-6251 temperature (t#), pressure (p#), and CO2 (c#)
if (count1 ge 3) then begin
	d1=datafinal(0,index1)
	t1=datafinal(5,index1)
	p1=datafinal(6,index1)
	c1=datafinal(7,index1)
endif else begin
	d1=[min(jday),median(jday),max(jday)]
	t1=[-999.,-999.,-999.]
	p1=[-999.,-999.,-999.]
	c1=[-999.,-999.,-999.]
endelse

if (count2 ge 3) then begin
	d2=datafinal(0,index2)
	t2=datafinal(5,index2)
	p2=datafinal(6,index2)
	c2=datafinal(7,index2)
endif else begin
	d2=[min(jday),median(jday),max(jday)]
	t2=[-999.,-999.,-999.]
	p2=[-999.,-999.,-999.]
	c2=[-999.,-999.,-999.]
endelse

if (count3 ge 3) then begin
	d3=datafinal(0,index3)
	t3=datafinal(5,index3)
	p3=datafinal(6,index3)
	c3=datafinal(7,index3)
endif else begin
	d3=[min(jday),median(jday),max(jday)]
	t3=[-999.,-999.,-999.]
	p3=[-999.,-999.,-999]
	c3=[-999.,-999.,-999]
endelse

if (count4 ge 3) then begin
	d4=datafinal(0,index4)
	t4=datafinal(5,index4)
	p4=datafinal(6,index4)
	c4=datafinal(7,index4)
endif else begin
	d4=[min(jday),median(jday),max(jday)]
	t4=[-999.,-999.,-999.]
	p4=[-999.,-999.,-999]
	c4=[-999.,-999.,-999]
endelse

if (count5 ge 3) then begin
	d5=datafinal(0,index5)
	t5=datafinal(5,index5)
	p5=datafinal(6,index5)
	c5=datafinal(7,index5)
endif else begin
	d5=[min(jday),median(jday),max(jday)]
	t5=[-999.,-999.,-999.]
	p5=[-999.,-999.,-999]
	c5=[-999.,-999.,-999]
endelse

d8=datafinal(0,index8)
t8=datafinal(5,index8)
p8=datafinal(6,index8)
c8=datafinal(7,index8)

d9=datafinal(0,index9)
t9=datafinal(5,index9)
p9=datafinal(6,index9)
c9=datafinal(7,index9)

d10=datafinal(0,index10)
t10=datafinal(5,index10)
p10=datafinal(6,index10)
c10=datafinal(7,index10)

d8b=datafinal(0,index8b)
t8b=datafinal(5,index8b)
p8b=datafinal(6,index8b)
c8b=datafinal(7,index8b)

;set time values for interpolation (every 3 minutes)
rows=n_elements(id)
n=(datafinal(0,rows-1)-datafinal(0,0))*24*60/3
t=datafinal(0,0)+(findgen(n+1)/n)*(datafinal(0,rows-1)-datafinal(0,0))

;calculate jday and hhmm for interpolated values
jday=fix(t)
hh=fix((t-jday)*24)
mm=(((t-jday)*24)-hh)*60
for z=0,n do begin
	if (float(mm(z)) eq 60) then begin
	mm(z)=0
	hh(z)=hh(z)+1
	endif
endfor
hhmm=hh*100+mm

;fit data using spline function, and interpolate
t1spline=spline(d1,t1,t)
p1spline=spline(d1,p1,t)
c1spline=spline(d1,c1,t)

t2spline=spline(d2,t2,t)
p2spline=spline(d2,p2,t)
c2spline=spline(d2,c2,t)

t3spline=spline(d3,t3,t)
p3spline=spline(d3,p3,t)
c3spline=spline(d3,c3,t)

t4spline=spline(d4,t4,t)
p4spline=spline(d4,p4,t)
c4spline=spline(d4,c4,t)

t5spline=spline(d5,t5,t)
p5spline=spline(d5,p5,t)
c5spline=spline(d5,c5,t)

t8spline=spline(d8,t8,t)
p8spline=spline(d8,p8,t)
c8spline=spline(d8,c8,t)

;apply pressure correction (P0/P) to mV readings
if (count1 ge 3) then c1spline=c1spline*100./p1spline
if (count2 ge 3) then c2spline=c2spline*100./p2spline
if (count3 ge 3) then c3spline=c3spline*100./p3spline
if (count4 ge 3) then c4spline=c4spline*100./p4spline
if (count5 ge 3) then c5spline=c5spline*100./p5spline
c8spline=c8spline*100./p8spline
c8=c8*100./p8
c9=c9*100./p9
c10=c10*100./p10
c8b=c8b*100./p8b

;perform 2nd order polynomial calibrations
polyfunc=dblarr(3,count10); setup array for variables (3 columns x n polynomial fits)
stdconc=dblarr(3,count10); setup array for temperature corrected standard concentrations (standard minus ref concentration * (T0/T))
for k=0,count10-1 do begin; do for as many standard were run during this time period
stdconc(*,k)=[(lowstd-refgas)*(273.15)/(t9(k)+273.15),0,(highstd-refgas)*(273.15)/(t10(k)+273.15)]; standard minus ref conc. (ppm) * (T0/T)
mv=[c9(k),c8b(k),c10(k)]; =mV*(P0/P) measurements
polyfunc(*,k)=poly_fit(mv,stdconc(*,k),2)
endfor

;manual input of coefficients for period when calibration proceedure was not working
IF (yy EQ 01) AND ((ddd GE 102) AND (ddd LE 126)) THEN begin
   polyfunc(0,*)=1.947     ; C0 = avg from 69-72 (4 d before) and 127-130 (4 d after)
   polyfunc(1,*)=0.2184     ; C1 = avg from 69-72 (4 d before) and 127-130 (4 d after)
   polyfunc(2,*)=3.133E-05 ; C2 = avg from 69-72 (4 d before) and 127-130 (4 d after)
endif

;calculate average values for plotting
avg1=moment(polyfunc(0,*))
avg2=moment(polyfunc(1,*))
avg3=moment(polyfunc(2,*))
polyavg1='Averages (n='+strtrim(string(count10),2)+')'
polyavg2='C0='+strtrim(string(avg1(0)),1)+'   C1='+strtrim(string(avg2(0)),1)+'   C2='+strtrim(string(avg3(0)),1)
polyavg3='Low standard='+strtrim(string(lowstd),1)+' ppm'
polyavg4='Reference gas='+strtrim(string(refgas),1)+' ppm'
polyavg5='High standard='+strtrim(string(highstd),1)+' ppm'

;create array for plotting polynomial fits
polyplot=dblarr(count10+1,100); setup array for plotting as many columns as fits plus one for co2 concs
polyplot(0,*)=min(c9)+(findgen(100)/(100-1))*(max(c10)-min(c9)); create 100 data points between min and max CO2 values
for l=1,count10 do begin
polyplot(l,*)=polyfunc(0,l-1)+polyfunc(1,l-1)*polyplot(0,*)+polyfunc(2,l-1)*polyplot(0,*)^2
endfor

;fit polynomial data using spline function and interpolate
polyspline=dblarr(3,n+1)
polyspline(0,*)=spline(d9,polyfunc(0,*),t)
polyspline(1,*)=spline(d9,polyfunc(1,*),t)
polyspline(2,*)=spline(d9,polyfunc(2,*),t)

;calculate ppm for individual and interpolated reference gas values
ppm8=(polyspline(0,*)+polyspline(1,*)*c8spline+polyspline(2,*)*c8spline^2)*(t8spline+273.15)/273.15
refpts=ppm8(*,index8)

;calculate ppm for interpolated air sampling heights
;(calculate ppm using polynomial, multiply by T/TO, subtract offset, add to known concentration of ref gas)
if (count1 ge 3) then ppm1=(polyspline(0,*)+polyspline(1,*)*c1spline+polyspline(2,*)*c1spline^2)*((t1spline+273.15)/273.15)+refgas-ppm8 $
	else ppm1=c1spline
if (count2 ge 3) then ppm2=(polyspline(0,*)+polyspline(1,*)*c2spline+polyspline(2,*)*c2spline^2)*((t2spline+273.15)/273.15)+refgas-ppm8 $
	else ppm2=c2spline
if (count3 ge 3) then ppm3=(polyspline(0,*)+polyspline(1,*)*c3spline+polyspline(2,*)*c3spline^2)*((t3spline+273.15)/273.15)+refgas-ppm8 $
	else ppm3=c3spline
if (count4 ge 3) then ppm4=(polyspline(0,*)+polyspline(1,*)*c4spline+polyspline(2,*)*c4spline^2)*((t4spline+273.15)/273.15)+refgas-ppm8 $
	else ppm4=c4spline
if (count5 ge 3) then ppm5=(polyspline(0,*)+polyspline(1,*)*c5spline+polyspline(2,*)*c5spline^2)*((t5spline+273.15)/273.15)+refgas-ppm8 $
	else ppm5=c5spline

;exclude range of values (specified by user at start of program)
if exclude eq 'y' then BEGIN
        index=where(level(*) EQ '1' OR level(*) EQ 'all')
	if index(0) NE -1 then begin
		badrange=where((hhmm ge starttime) and (hhmm le endtime))
		ppm1(badrange)=-999.
		t1spline(badrange)=-999.
		p1spline(badrange)=-999.
		c1spline(badrange)=-999.
	endif

        index=where(level(*) EQ '2' OR level(*) EQ 'all')
	if index(0) NE -1 then begin
		badrange=where((hhmm ge starttime) and (hhmm le endtime))
		ppm2(badrange)=-999.
		t2spline(badrange)=-999.
		p2spline(badrange)=-999.
		c2spline(badrange)=-999.
	endif

        index=where(level(*) EQ '3' OR level(*) EQ 'all')
	if index(0) NE -1 then begin
		badrange=where((hhmm ge starttime) and (hhmm le endtime))
		ppm3(badrange)=-999.
		t3spline(badrange)=-999.
		p3spline(badrange)=-999.
		c3spline(badrange)=-999.
	endif

        index=where(level(*) EQ '4' OR level(*) EQ 'all')
	if index(0) NE -1 then begin
		badrange=where((hhmm ge starttime) and (hhmm le endtime))
		ppm4(badrange)=-999.
		t4spline(badrange)=-999.
		p4spline(badrange)=-999.
		c4spline(badrange)=-999.
	endif

        index=where(level(*) EQ '5' OR level(*) EQ 'all')
	if index(0) NE -1 then begin
		badrange=where((hhmm ge starttime) and (hhmm le endtime))
		ppm5(badrange)=-999.
		t5spline(badrange)=-999.
		p5spline(badrange)=-999.
		c5spline(badrange)=-999.
	endif
endif

;combine vectors and trim file to 24 hours
ppmall=dblarr(7,n+1)
ppmall(*,*)=-999.
ppmall(0,*)=jday
ppmall(1,*)=hhmm
ppmall(2,*)=ppm1(0,*)
ppmall(3,*)=ppm2(0,*)
ppmall(4,*)=ppm3(0,*)
ppmall(5,*)=ppm4(0,*)
ppmall(6,*)=ppm5(0,*)
dddkeep=where(jday eq ddd)
ppmall=ppmall(*,dddkeep)
help,ppmall; print array information

;compute stats for graphing
allc=[c1spline(dddkeep),c2spline(dddkeep),c3spline(dddkeep),c4spline(dddkeep), $
	c5spline(dddkeep)]
avgc=moment(allc(where(allc ne -999)))
allt=[t1spline(dddkeep),t2spline(dddkeep),t3spline(dddkeep),t4spline(dddkeep), $
	t5spline(dddkeep)]
avgt=moment(allt(where(allt ne -999)))
allp=[p1spline(dddkeep),p2spline(dddkeep),p3spline(dddkeep),p4spline(dddkeep), $
	p5spline(dddkeep)]
avgp=moment(allp(where(allp ne -999)))

;make sure there are enough good values to compute statistics (after removing levels)
newindex1=where(ppmall(2,*) NE -999,newcount1);  2 feet
newindex2=where(ppmall(3,*) NE -999,newcount2);  5 feet
newindex3=where(ppmall(4,*) NE -999,newcount3); 10 feet
newindex4=where(ppmall(5,*) NE -999,newcount4); 25 feet
newindex5=where(ppmall(6,*) NE -999,newcount5); 45 feet

if (newcount1 ge 3) then avgppm1=moment(ppmall(2,where(ppmall(2,newindex1) ne -999))) else avgppm1=[-999.,-999.,-999.,-999.]
if (newcount2 ge 3) then avgppm2=moment(ppmall(3,where(ppmall(3,newindex2) ne -999))) else avgppm2=[-999.,-999.,-999.,-999.]
if (newcount3 ge 3) then avgppm3=moment(ppmall(4,where(ppmall(4,newindex3) ne -999))) else avgppm3=[-999.,-999.,-999.,-999.]
if (newcount4 ge 3) then avgppm4=moment(ppmall(5,where(ppmall(5,newindex4) ne -999))) else avgppm4=[-999.,-999.,-999.,-999.]
if (newcount5 ge 3) then avgppm5=moment(ppmall(6,where(ppmall(6,newindex5) ne -999))) else avgppm5=[-999.,-999.,-999.,-999.]
avgppm8=moment(ppm8(dddkeep))

;begin plotting loop

for plot_loop=plotstart,plotend do BEGIN

fname='lc'+strmid(file2,0,5)+'.ppm'

if (plot_loop eq 0) then begin
	window,/free,xsize=382.5,ysize=495,title=fname,xpos=5,ypos=200
endif

if (plot eq 'y') and (plot_loop eq 1) then begin; set plotting to postscript second time through loop
	set_plot,'ps'
	device,file='/eddy/s4/lcreek/ppm/ps/'+file2+'.ps',/inches,xoffset=.5,yoffset=.5,xsize=7.5,ysize=10
endif

if (plot eq 'y') or (screen eq 'y') then begin; plot only if screen or ps files requested
	;plot graphs (note:0 empty fields,2 columns, 3 rows, 0 z-dim plots, starts in top-left then down)
	!p.multi=[0,2,3,0,1]

	;make x-axis labels
	labels=['0','4','8','12','16','20','24']

	;plot fractional Julian date vs. mV*p0/p for 7 levels
	maxy=max(allc(where(allc ne -999)))
	miny=min(allc(where(allc ne -999)))
	plot,t(dddkeep),c1spline(dddkeep),title='Raw CO!d2!n voltages',ytitle='mV*(P0/P)',xtitle='Time (UT)',xrange=[ddd,ddd+1],charthick=2,yrange=[miny,maxy],xticks=6,xtickname=labels,xminor=4,xstyle=1,min_value=-998
	oplot,t(dddkeep),c2spline(dddkeep)
	oplot,t(dddkeep),c3spline(dddkeep)
	oplot,t(dddkeep),c4spline(dddkeep)
	oplot,t(dddkeep),c5spline(dddkeep)
endif

if (plot eq 'y') and (plot_loop eq 1) then begin; other titles for ps plot only
	xyouts,.48,1.02,fname,/normal,charthick=2; main title=ouput file name

	xyouts,.58,.97,polyavg1,size=0.5,/normal,charthick=2; polynominal stats
	xyouts,.58,.96,polyavg2,size=0.5,/normal,charthick=2
	xyouts,.58,.94,polyavg3,size=0.5,/normal,charthick=2
	xyouts,.58,.93,polyavg4,size=0.5,/normal,charthick=2
	xyouts,.58,.92,polyavg5,size=0.5,/normal,charthick=2

	xyouts,.08,.97,'Level',size=0.5,/normal,charthick=2; mean/stdev mV*(P0/P) (7 levels)
	xyouts,.08,.96,'All',size=0.5,/normal,charthick=2
	xyouts,.14,.97,'Mean',size=0.5,/normal,charthick=2
	xyouts,.1,.96,avgc(0),size=0.5,/normal,charthick=2
	xyouts,.19,.97,'Std dev',size=0.5,/normal,charthick=2
	xyouts,.16,.96,avgc(1)^.5,size=0.5,/normal,charthick=2

	xyouts,.08,.635,'Level',size=0.5,/normal,charthick=2; mean/stdev LI-6251 temperature (7 levels)
	xyouts,.08,.625,'All',size=0.5,/normal,charthick=2
	xyouts,.14,.635,'Mean',size=0.5,/normal,charthick=2
	xyouts,.1,.625,avgt(0),size=0.5,/normal,charthick=2
	xyouts,.19,.635,'Std dev',size=0.5,/normal,charthick=2
	xyouts,.16,.625,avgt(1)^.5,size=0.5,/normal,charthick=2

	xyouts,.08,.30,'Level',size=0.5,/normal,charthick=2; mean/stdev LI-6251 pressure (7 levels)
	xyouts,.08,.29,'All',size=0.5,/normal,charthick=2
	xyouts,.14,.30,'Mean',size=0.5,/normal,charthick=2
	xyouts,.1,.29,avgp(0),size=0.5,/normal,charthick=2
	xyouts,.19,.30,'Std dev',size=0.5,/normal,charthick=2
	xyouts,.16,.29,avgp(1)^.5,size=0.5,/normal,charthick=2

	xyouts,.58,.635,'Mean='+strtrim(string(avgppm8(0)),2),size=0.5,/normal,charthick=2; reference gas (ppm)
	xyouts,.58,.625,'Std dev='+strtrim(string(avgppm8(0)),2),size=0.5,/normal,charthick=2

	xyouts,.58,.30,'Level',size=0.5,/normal,charthick=2; ppm Mixing ratios  (7 levels)
	xyouts,.58,.29,' 2 ft',size=0.5,/normal,charthick=2
	xyouts,.58,.28,' 5 ft',size=0.5,/normal,charthick=2
	xyouts,.58,.27,'10 ft',size=0.5,/normal,charthick=2
	xyouts,.58,.26,'25 ft',size=0.5,/normal,charthick=2
	xyouts,.58,.25,'45 ft',size=0.5,/normal,charthick=2
	xyouts,.64,.3,'Mean',size=0.5,/normal,charthick=2
	xyouts,.6,.29,avgppm1(0),size=0.5,/normal,charthick=2
	xyouts,.6,.28,avgppm2(0),size=0.5,/normal,charthick=2
	xyouts,.6,.27,avgppm3(0),size=0.5,/normal,charthick=2
	xyouts,.6,.26,avgppm4(0),size=0.5,/normal,charthick=2
	xyouts,.6,.25,avgppm5(0),size=0.5,/normal,charthick=2
	xyouts,.69,.30,'Std dev',size=0.5,/normal,charthick=2
	xyouts,.66,.29,avgppm1(1)^.5,size=0.5,/normal,charthick=2
	xyouts,.66,.28,avgppm2(1)^.5,size=0.5,/normal,charthick=2
	xyouts,.66,.27,avgppm3(1)^.5,size=0.5,/normal,charthick=2
	xyouts,.66,.26,avgppm4(1)^.5,size=0.5,/normal,charthick=2
	xyouts,.66,.25,avgppm5(1)^.5,size=0.5,/normal,charthick=2
endif

if (plot eq 'y') or (screen eq 'y') then begin; plot only if screen or ps files requested
	;plot fractional Julian date vs. LI-6251 temperature for 7 levels
	maxy=max(allt(where(allt ne -999)))
	miny=min(allt(where(allt ne -999)))
	plot,t(dddkeep),t1spline(dddkeep),title='LI-6262 Temperature',ytitle='!uo!nC',yrange=[miny,maxy],xtitle='Time (UT)',xrange=[ddd,ddd+1],charthick=2,xticks=6,xtickname=labels,xminor=4,xstyle=1,min_value=-998
	oplot,t(dddkeep),t2spline(dddkeep)
	oplot,t(dddkeep),t3spline(dddkeep)
	oplot,t(dddkeep),t4spline(dddkeep)
	oplot,t(dddkeep),t5spline(dddkeep)

	;plot fractional Julian date vs. LI-6251 pressure for 7 levels
	maxy=max(allp(where(allp ne -999)))
	miny=min(allp(where(allp ne -999)))
	plot,t(dddkeep),p1spline(dddkeep),title='LI-6262 Pressure',ytitle='kPa',xtitle='Time (UT)',yrange=[miny,maxy],xrange=[ddd,ddd+1],charthick=2,xticks=6,xtickname=labels,xminor=4,xstyle=1,min_value=-998
	oplot,t(dddkeep),p2spline(dddkeep)
	oplot,t(dddkeep),p3spline(dddkeep)
	oplot,t(dddkeep),p4spline(dddkeep)
	oplot,t(dddkeep),p5spline(dddkeep)

	;plot polynomial calibrations (mV*P0/P vs. ppm (Cs-Cr)*T0/T)
	plot,polyplot(0,*),polyplot(1,*),title='Calibration Curves',ytitle='(Cs-Cr)*T0/T',xtitle='mV*P0/P',charthick=2
	oplot,c9,stdconc(0,*),psym=1,symsize=0.5
	oplot,c8b,stdconc(1,*),psym=1,symsize=0.5
	oplot,c10,stdconc(2,*),psym=1,symsize=0.5
	for m=2,count10 do begin
		oplot,polyplot(0,*),polyplot(m,*)
	endfor

	;plot fractional Julian date vs. reference gas (ppm CO2)
	plot,t(dddkeep),ppm8(dddkeep),title='Reference gas',ytitle='ppm CO!d2',xtitle='Time (UT)',xrange=[ddd,ddd+1],charthick=2,xticks=6,xtickname=labels,xminor=4,xstyle=1
	;oplot,d8,refpts,psym=1,symsize=0.5

	;plot fractional Julian date vs. ppm CO2 at seven sampling heights
	ppm=[ppm1(dddkeep),ppm2(dddkeep),ppm3(dddkeep),ppm4(dddkeep),ppm5(dddkeep)]
	maxy=max(ppm(where(ppm ne -999)))
	miny=min(ppm(where(ppm ne -999)))
	plot,t(dddkeep),ppm1(dddkeep),title='CO!d2!n Mixing ratios',ytitle='ppm CO!d2',xtitle='Time (UT)',yrange=[miny,maxy],xrange=[ddd,ddd+1],charthick=2,xticks=6,xtickname=labels,xminor=4,xstyle=1,min_value=-998
	oplot,t(dddkeep),ppm2(dddkeep)
	oplot,t(dddkeep),ppm3(dddkeep)
	oplot,t(dddkeep),ppm4(dddkeep)
	oplot,t(dddkeep),ppm5(dddkeep)
	endif

if (plot eq 'y') and (plot_loop eq 1) then begin; close ps file, reset plotting to 1 per page, return plotting to unix device
	device,/close
	!p.multi=0
	set_plot,'X'
endif

endfor; end plot loop

if (plot eq 'y') then begin
	cmd="spawn,'gzip -f /eddy/s4/lcreek/ppm/ps/'+file2+'.ps'"
	r=execute(cmd)
endif

skipcalcs:

;set up 24 hour array to include all possible times
ppmday=dblarr(7,480); 3 minute intervals
ppmday(0,*)=ddd; fill day column

;fill hhmm column
i=0
for j=0,23 do begin
	for k=0,19 do begin
		hhmm=(j*100)+(k*3)
		ppmday(1,i)=hhmm
		i=i+1
	endfor
endfor
ppmday(2:6,*)=-999

IF skip NE 1 THEN begin; fill data where available
   for i=0,479 do begin
      index=where(float(ppmall(1,*)) eq float(ppmday(1,i)))
      if (index(0) ne -1) then ppmday(2:6,i)=ppmall(2:6,index)
   ENDFOR
endif

if (plot eq 'y') then begin
	;save 24 hour array to ascii file
        fname='lc'+strmid(file2,0,5)+'.ppm'
	print,'Saving '+fname
	openw,lun1,filepath(fname,root_dir='/eddy/s4/lcreek/ppm'),/get_lun
	thisformat='(i3,TR3,i4,TR1,7(d9.3,TR1))'
	for i=0,479 do begin
		printf,lun1,ppmday(*,i),format=thisformat
	endfor
	close,lun1
	free_lun,lun1
endif

endfor; end multiple day loop

print,'Completed!'

return
end