PRO plot_m,data,days,screen=screen ;need averages - maybe? yr = fix(data[0,0]) MOD 100 doys = fix(data[1,*]) plotfilename = co2_const('file_head')+string(yr,format='(i2.2)')+'.m' ;days give the start and end point of the data IF NOT keyword_set(days) THEN tdays = [0,days_in_year(yr)] ELSE tdays = days time = data[2,*] / 24.0 time = time + data[1,*] ;day of year + fraction of day labels = numgen(float(tdays[0]),float(tdays[1])+1,numvals=7) labels = strtrim(string(labels,format='(f10.2)'),2) IF NOT keyword_set(Screen) THEN BEGIN plotdir = co2_basedir('plots')+'mfiles/' create_dir,plotdir print,'Outputting postscript plots to ',plotdir+plotfilename+'.ps' printon,plotdir+plotfilename,/nodir ENDIF oldmulti = !p.multi !p.multi = [0,1,5] ;plot graphs (Temp,Ta,Ts) (Rh,Dp,Q) (Gpar1,Gpar2,Throughfall) ;(Wspd_20, Wdir_20) ;air temp at 20,10,2 = 12,13,14 range = [min(data[12:14,*],/nan),max(data[12:14,*],/nan)] IF NOT isallnan(data[12,*]) THEN plot,time,data[12,*],title='Air temperature at 20,10,2 m',ytitle='!uo!nC',xtitle='Day of Year',charthick=2,xtickname=labels,xticks=6,xminor=4,yrange=range IF NOT isallnan(data[13,*]) THEN oplot,time,data[13,*] IF NOT isallnan(data[14,*]) THEN oplot,time,data[14,*] ;air temp at 100,50,25,5,0 = 15,16,17,18,19 range = [min(data[15:19,*],/nan),max(data[15:19,*],/nan)] plot,time,data[15,*],title='Air temperature above soil 100,50,25,5,0 cm',ytitle='!uo!nC',xtitle='Day of Year',charthick=2,xtickname=labels,xticks=6,xminor=4,yrange=range oplot,time,data[16,*] oplot,time,data[17,*] oplot,time,data[18,*] oplot,time,data[19,*] ;soil temp at 5,10,25,50,100 = 20,21,22,23,24 range = [min(data[20:24,*],/nan),max(data[20:24,*],/nan)] IF NOT isallnan(data[20,*]) THEN plot,time,data[20,*],title='Soil temperature 5,10,25,50,100 cm',ytitle='!uo!nC',xtitle='Day of Year',charthick=2,xtickname=labels,xticks=6,xminor=4,yrange=range IF NOT isallnan(data[21,*]) THEN oplot,time,data[21,*] IF NOT isallnan(data[22,*]) THEN oplot,time,data[22,*] IF NOT isallnan(data[23,*]) THEN oplot,time,data[23,*] IF NOT isallnan(data[24,*]) THEN oplot,time,data[24,*] ;RH at 20,10,2 = 3,4,5 range = [min(data[3:5,*],/nan),max(data[3:5,*],/nan)] IF NOT isallnan(data[3,*]) THEN plot,time,data[3,*],title='RH at 20,10,2 m',ytitle='%',xtitle='Day of Year',charthick=2,xtickname=labels,xticks=6,xminor=4,yrange=range IF NOT isallnan(data[4,*]) THEN oplot,time,data[4,*] IF NOT isallnan(data[5,*]) THEN oplot,time,data[5,*] ;DP at 20,10,2 = 6,7,8 range = [min(data[6:8,*],/nan),max(data[6:8,*],/nan)] IF NOT isallnan(data[6,*]) THEN plot,time,data[6,*],title='Dewpoint at 20,10,2 m',ytitle='!uo!nC',xtitle='Day of Year',charthick=2,xtickname=labels,xticks=6,xminor=4,yrange=range IF NOT isallnan(data[7,*]) THEN oplot,time,data[7,*] IF NOT isallnan(data[8,*]) THEN oplot,time,data[8,*] ;Q at 20,10,2 = 9,10,11 range = [min(data[9:11,*],/nan),max(data[9:11,*],/nan)] IF NOT isallnan(data[9,*]) THEN plot,time,data[9,*],title='H2O mixing ratio at 20,10,2 m',ytitle='g/kg',xtitle='Day of Year',charthick=2,xtickname=labels,xticks=6,xminor=4,yrange=range IF NOT isallnan(data[10,*]) THEN oplot,time,data[10,*] IF NOT isallnan(data[11,*]) THEN oplot,time,data[11,*] ;Ground PAR = 25,26 range = [min(data[25:26,*],/nan),max(data[9:11,*],/nan)] IF NOT isallnan(data[25,*]) THEN plot,time,data[25,*],title='PAR at 1 m',ytitle='umol/m2s',xtitle='Day of Year',charthick=2,xtickname=labels,xticks=6,xminor=4,yrange=range IF NOT isallnan(data[26,*]) THEN oplot,time,data[26,*] ;Windspeed = 27 IF NOT isallnan(data[27,*]) THEN plot,time,data[27,*],title='Subcanopy wind speed 20m',ytitle='m/s',xtitle='Day of Year',charthick=2,xtickname=labels,xticks=6,xminor=4 ;Winddir = 28 IF NOT isallnan(data[28,*]) THEN plot,time,data[28,*],title='Subcanopy wind direction 20m',ytitle='deg',xtitle='Day of Year',charthick=2,xtickname=labels,xticks=6,xminor=4,yrange=[0,360] ;Throughfall = 29 IF NOT isallnan(data[29,*]) THEN plot,time,data[29,*],title='1 hour surface throughfall',ytitle='mm',xtitle='Day of Year',charthick=2,xtickname=labels,xticks=6,xminor=4 IF NOT keyword_set(Screen) THEN printoff !p.multi = oldmulti END PRO avg_m,yr,nocalc=nocalc,screen=screen ;clean micromet file ;calc dewpoint and mixing ratio IF NOT keyword_set(yr) THEN yr = co2_const('cur_year') year = string(yr MOD 100,format='(i2.2)') nxtyear = string(yr+1 MOD 100,format='(i2.2)') outdir = co2_basedir('micromet',yr=yr) halfname = co2_const('file_head')+year+'.m' fullname = outdir + halfname + '.xdf' poutdir = co2_basedir('qfiles',yr=yr) poutdir2 = co2_basedir('qfiles',yr=yr+1) halfpressname = co2_const('file_head')+year+'.q' halfpressname2 = co2_const('file_head')+nxtyear+'.q' pressname = poutdir + halfpressname + '.xdf' pressname2 = poutdir2 + halfpressname2 + '.xdf' IF file_test(fullname,/read) THEN data = read_xdf(fullname,header=header) ELSE BEGIN print,'Cannot read ',fullname GOTO,leave ENDELSE ;clean data i = where((data EQ -999) OR (data EQ -9999) OR (data EQ -6999) OR (data EQ 6999) OR (data EQ -99999) OR (data EQ 9999) OR (data EQ 99999),num) IF num GT 0 THEN data[i] = nan() ;calculate dewpoint and mixing ratio IF NOT keyword_set(nocalc) THEN BEGIN temp20 = data[12,*] temp10 = data[13,*] temp02 = data[14,*] rh20 = data[3,*] rh10 = data[4,*] rh02 = data[5,*] des20 = 0.6112*exp((17.67*temp20)/(243.5+temp20)) des10 = 0.6112*exp((17.67*temp10)/(243.5+temp10)) des02 = 0.6112*exp((17.67*temp02)/(243.5+temp02)) vp20 = (des20 * rh20 * 0.01) vp10 = (des10 * rh10 * 0.01) vp02 = (des02 * rh02 * 0.01) wr20 = alog(1.6373 * vp20) wr10 = alog(1.6373 * vp10) wr02 = alog(1.6373 * vp02) dew20 = (241.88 * wr20) / (17.558 - wr20) dew10 = (241.88 * wr10) / (17.558 - wr10) dew02 = (241.88 * wr02) / (17.558 - wr02) data[6,*] = dew20 data[7,*] = dew10 data[8,*] = dew02 dea20 = 0.6112*exp((17.67*dew20)/(243.5+dew20)) dea10 = 0.6112*exp((17.67*dew10)/(243.5+dew10)) dea02 = 0.6112*exp((17.67*dew02)/(243.5+dew02)) IF file_test(pressname,/read) THEN BEGIN pdata = read_xdf(pressname,header=pheader) press36 = shift(pdata[2,*],-6) ;shift UTC data to LST np36 = n_elements(press36) press36[np36-6:np36-1] = nan() IF file_test(pressname2,/read) THEN BEGIN pdata2 = read_xdf(pressname2,header=pheader2) press36_2 = pdata2[2,0:5] press36[np36-6:np36-1] = press36_2 ENDIF badp36 = where((press36 GT 100.0) OR (press36 LT 90.0),nbp36) IF nbp36 GT 0 THEN press36[badp36] = nan() press36 = zapbadval(average_arr(press36,2,/nan)) badp36 = where((press36 GT 100.0) OR (press36 LT 90.0),nbp36) IF nbp36 GT 0 THEN press36[badp36] = mean(press36,/nan) ;dp = (101.325 * (1 -((1-(dz/44307.69231))^5.25328))) press20 = press36 + 0.192054 press10 = press36 + 0.311961 press02 = press36 + 0.407783 qmix20 = (dea20*0.622)/(press20-dea20)*10^3 qmix10 = (dea10*0.622)/(press10-dea10)*10^3 qmix02 = (dea02*0.622)/(press02-dea02)*10^3 data[9,*] = qmix20 data[10,*] = qmix10 data[11,*] = qmix02 ENDIF ENDIF ;write new datafile print,'Writing cleaned ',fullname write_xdf,fullname,data,header=header ;get subdata IF NOT keyword_set(days) THEN BEGIN ;clip to good data gooddat = where(isnotnan(data[3,*]),numw) IF numw GT 1 THEN BEGIN tdays = [data[1,gooddat[0]],data[1,gooddat[numw-1]]] ENDIF ELSE tdays = [1,days_in_year(yr)] ENDIF ELSE BEGIN IF n_elements(days) EQ 1 THEN tdays = [days,days] ELSE tdays = days ENDELSE startday = long(tdays[0]) endday = long(tdays[1]) startloc = (startday - 1) * 24 endloc = ((endday - 1) * 24)+23 data = data[*,startloc:endloc] plot_m,data,tdays,screen=screen leave: print,'avg_m complete' END