;look at ch4/co2 mixing ratios and flux relationships PRO wlef_nacp_ch4 ;Load WLEF flask CH4 and use it to gap-fill? ;non-linear fit? ;To run this ;also need to have myfit.pro and average_arr.pro in the path ;Both of these are in /data/code/util/idl ;Start IDL ;idl ;IN IDL TYPE: ;.compile wlef_nacp.pro ;wlef_nacp_ch4 ;which will do the following: ;READ the data file dta = (read_ascii(mydocs()+'cheas/tccon/ch4_co2.txt',data_start=1,header=h)).(0) ;Replace missing data (-999.0) with Not a Number dta[where(dta EQ -999)]=!values.f_nan ;Extract the four variables from 2005-2009 (discard 2004, I already ;removed 2010 from the file) tm = reform(dta[2,26:*]) ; ch4 = reform(dta[3,26:*]) ; co2 = reform(dta[4,26:*]) nee = reform(dta[5,26:*]) ;read in the new data restore,mydocs()+'cheas/tccon/newch4.xdf' xtm = 2005+(findgen(1825)/365.) ch4 = average_arr(trop_ch4/1000.,14,/nan) co2 = average_arr(xco2,14,/nan) ;This tells us which values of column CO2 are good gval = where(finite(co2)) gval2 = where(finite(ch4)) ;Linfit fits a line between x and y, and returns intercept and slope ;Here we fit a line between Time and CO2, and also Time and CH4 fos_co2 = linfit(tm[gval],co2[gval]) fos_ch4 = linfit(tm[gval2],ch4[gval2]) ;Using the slope and intercept, remove this linear trend in CO2 and ;CH4 ;The linear trend is the fossil fuel signal co2f = co2 - (fos_co2[0] + fos_co2[1]*tm) ch4f = ch4 - (fos_ch4[0] + fos_ch4[1]*tm) ;Gap-fill the CH4 and CO2 time series with a spline fit co2_fill = interpol(co2f[gval],tm[gval],tm,/spline) ch4_fill = interpol(ch4f[gval2],tm[gval2],tm,/spline) ;Remove the 2005 winter monthly NEE that is likely bad nee[18:25]=!values.f_nan ;Gap fill those values with the interpolate function gnee = where(finite(nee)) nee_f = interpol(nee[gnee],tm[gnee],tm,/spline) ;Compare the time derivative for CO2 and CH4, and then use a 5 point ;smoother on it, and then shift the time series forward by 1 to match ;up with the NEE time series dco2 = shift(ts_smooth(deriv(co2_fill),5),1) dch4 = shift(ts_smooth(deriv(ch4_fill),5),1) ;Estimate the slope of the fit between NEE and dCO2/dt ;Myfit is like linfit, but deals with missing data and gives you a sigma ; gvals_c = where(finite(nee) AND finite(co2)) cfit = myfit(dco2,nee_f,sigma=cfit_s,corr=cfit_r) ;*NOTE: PERHAPS WE SHOULD USE A CURVED FIT HERE? cfit2 = myfit(dco2,nee_f,sigma=cfit_s_2,corr=cfit_r_2,degree=2) cfit2_s = cfit2[1]+(2*dco2*cfit2[2]) ;Compute time series of CO2 NEE (from toweR) and CH4 NEE (from method) ;Convert units to useful values nee_co2 = nee_f * 1.0368 ;g C / day / m2 nee_ch4 = dch4*cfit[1] * 1036.8 ;mg C / day / m2 nee_ch4_2 = dch4*cfit2_s*1036.8 ;Compare the 5 year ensemble average monthly NEE for CO2 and CH4 ;This is a bit of IDL magic nee_ens_co2 = reform(nee_co2,26,5) nee_ens_co2 = total(nee_ens_co2,2,/nan)/total(finite(nee_ens_co2),2) nee_ens_ch4 = reform(nee_ch4,26,5) nee_ens_ch4 = total(nee_ens_ch4,2,/nan)/total(finite(nee_ens_ch4),2) nee_ens_ch4_2 = reform(nee_ch4_2,26,5) nee_ens_ch4_2 = total(nee_ens_ch4_2,2,/nan)/total(finite(nee_ens_ch4_2),2) nee_ens_tm = (findgen(26)*14)+7 ;read chamber flux CH4 data for the three wetlands sites = ['lc_means','wf_means','sf_means'] ; sites = ['lc_med','wf_med','sf_med'] FOR i = 0,2 DO BEGIN data = (read_ascii(mydocs()+'methane/chamber/'+sites[i])).(0) IF i GT 0 THEN chamb=[chamb,reform(data[3,*])] ELSE chamb = reform(data[3,*]) newctm = data[0,*] + (data[1,*]-1.0)/365. IF i GT 0 THEN ctm = [ctm,reform(newctm)] ELSE ctm = reform(newctm) ENDFOR ;Compute annual average flux and convert to units of gC/m2/yr or ;mgC/m2/yr depending on product ann_co2 = average_arr(nee_co2,26,/tot)*(365./26.) ann_ch4 = average_arr(nee_ch4,26,/tot)*(365./26.) ann_ch4_2 = average_arr(nee_ch4_2,26,/tot)*(365./26.) ann_tm = findgen(5)+2005 ;Enter debug mode stop ;gap-fill co2/ch4 first? ;correlation between nee and smoothed dc/dt lag 1 is high, better with ;smoother ;figures for poster ;Fig 1 ;TCCON CO2 and CH4 14-day avg (filled?) strat and trop ;Show actual daily CO2/CH4 and filled 14-day avg, also full ch4 ;EDDY TOWER NEE !p.multi = [0,1,2] plot,[0,0],[0,0],xrange=[2005,2010],yrange=[1.725,1.875],ytitle='ppm',xtitle='Year',xticks=5 oplot,xtm,xch4/1000.,psym=1,color=180 oplot,xtm,trop_ch4/1000.,psym=1 oplot,tm,interpol(ch4[gval2],tm[gval2],tm,/spline),thick=2,color=200 oplot,tm,fos_ch4[0] + fos_ch4[1]*tm,color=100 plot,[0,0],[0,0],xrange=[2005,2010],yrange=[370,395],ytitle='ppm',xtitle='Year',xticks=5 oplot,xtm,xco2,psym=1 oplot,tm,interpol(co2[gval],tm[gval],tm,/spline),thick=2,color=200 oplot,tm,fos_co2[0] + fos_co2[1]*tm,color=100 stop ;Fig 2 ;dco2 and dch4 !p.multi = [0,1,2] plot,[0,0],[0,0],xrange=[2005,2010],yrange=[-3.5,2],ytitle='umol m-2 s-1 / ppm 14_day-1 ',xtitle='Year',xticks=5 oplot,tm,nee,thick=2 oplot,tm,dco2,color=200,thick=2 plot,[0,0],[0,0],xrange=[2005,2010],yrange=[-4.5,4],ytitle='ppb 14_day-1 ',xtitle='Year',xticks=5 oplot,tm,dch4*1000 stop ;Fig 3 ;scatterplot of dCO2 to NEE (with correlation, use gvals only) !p.multi = [0,2,2] plot,dco2,nee,psym=1,xrange=[-3.5,2],yrange=[-3.5,2],xtitle='dc/dt',ytitle='NEE' oplot,findgen(8)-4,cfit[0]+cfit[1]*(findgen(8)-4),thick=2 xyouts,-3,1.5,'NEE = -0.26 + 0.89 dC/dt (r2=0.6)' stop ;Fig 4 ;time series of CH4 NEE (with chamber data pasted on top) ;comparison of annual CH4 to annnual CO2 NEE (in g or mg) !p.multi = [0,1,2] plot,[0,0],[0,0],xrange=[2005,2010],yrange=[-4,4],ytitle='CH4 NEE mgC m-2 day-1 CO2 NEE gC m-2 day-1',xtitle='Year',xticks=5 oplot,tm,nee_co2,color=180,thick=2 oplot,tm,nee_ch4,thick=2 oplot,ctm,chamb*1000,psym=1,color=100 plot,[0,0],[0,0],xrange=[2004.5,2009.5],yrange=[-200,75],ytitle='CH4 NEE mgC m-2 yr-1 CO2 NEE gC m-2 yr-1',xtitle='Year',xticks=4,xtickv=ann_tm oplot,ann_tm,ann_Co2,psym=1,color=180,symsize=2 oplot,ann_tm,ann_Ch4,psym=2,symsize=2 oplot,[2004,2010],[0,0] stop ;Fig 5 ;comparison of ensemble avg CO2 NEE and CH4 NEE !p.multi = [0,1,2] plot,[0,0],[0,0],xrange=[1,365],yrange=[-2.25,1.75],ytitle='CH4 NEE mgC m-2 day-1 CO2 NEE gC m-2 day-1',xtitle='Day of Year',xticks=25,xtickv=nee_ens_tm oplot,nee_ens_tm,nee_ens_co2,color=180,thick=2 oplot,nee_ens_tm,nee_ens_ch4,thick=2 oplot,[0,400],[0,0] stop END PRO wlef_nacp_readstrat ; hf = (read_ascii(mydocs()+'cheas/tccon/hf.txt',data_start=1,header=h)).(0) hf = read_xdf(mydocs()+'cheas/tccon/hf.xdf',header=h) ;2005-2009 gval = where(hf[0,*] GE 2005 AND hf[0,*] LE 2009 AND hf[2,*] GE 16.0 AND hf[2,*] LT 20.0) hfg = hf[*,gval] xch4 = make_Array(365.*5.,/float,value=nan()) xco2 = xch4 xhf = xch4 pout = xch4 hfg[16,*] = despike(hfg[16,*],nsig=1) hfg[18,*] = despike(hfg[18,*],nsig=1) hfg[28,*] = despike(hfg[28,*],nsig=1) FOR y = 2005,2009 DO BEGIN FOR d = 1,365 DO BEGIN gg = where(hfg[0,*] EQ y AND hfg[1,*] EQ d,ngg) IF ngg GT 0 THEN BEGIN xch4[((y-2005)*365)+(d-1)] = mean(hfg[16,gg],/nan)*1000. xco2[((y-2005)*365)+(d-1)] = mean(hfg[18,gg],/nan) xhf[((y-2005)*365)+(d-1)] = mean(hfg[28,gg],/nan) pout[((y-2005)*365)+(d-1)] = mean(hfg[8,gg],/nan) ENDIF ENDFOR ENDFOR xhf = despike(xhf,nsig=2) xch4 = despike(xch4,nsig=2) xco2 = despike(xco2,nsig=2) kAvogadro = 6.0221415e23 hf_column = xhf*pout*100/9.8/0.029*kAvogadro/10e4*10e-12 calc_slop = -8.343e-14*hf_column-577.4886 trop_ch4 = xch4 - (calc_slop * xhf)/1e3 stop ;save trop_ch4 save,trop_ch4,xch4,xco2,xhf,pout,filename=mydocs()+'cheas/tccon/newch4.xdf' stop ;HF_column=xhf_ppt*pout_hPa*100/9.8/0.029*kAvogadro/10e4*10e-12 [molec cm(-2)] ;calculated_slope=-8.343e-14*HF_column-577.4886 ;tropospheric_CH4_VMR=xch4_ppb - (calculated_slope*xhf_ppt)/10^3 [ppb] END FUNCTION wlef_ch4paper_Efit,x,a ;Y = A exp(Bx) term0 = exp(a[1]*reform(x)) term1 = a[0] * reform(x) * term0 efit = a[0] * term0 return,[[efit],[term0],[term1]] END PRO wlef_ch4paper ;Cumulative CH4 and CO2 flux, seasonal, spike based, spectral (EMD) ; p11 = read_xdf('/data/incoming/wlef/flux/2011/prefnee_2011.xdf',header=ph) ; p12 = read_xdf('/data/incoming/wlef/flux/2012/prefnee_2012.xdf',header=ph) ;Read CLEAN - get soil moisture, precip, temp, ch4 flux, co2 flux from here restore,'/data/incoming/wlef/alldata19952013-processed.sav' ch4dflux = average_arr(ch4flux_122[*,16:17],24,/nan)*1.0368 ch4d = average_arr(nee_ch4_122[*,16:17],24,/nan)*1.0368 td = average_arr(tair_30_filled[*,16:17],24,/nan)*1.0368 pard = average_arr(par_pref[*,16:17],24,/nan)*86400l/1e6/2.2 co2d = average_arr(fillnee[*,16:17],24,/nan)*1.0368 recod = average_arr(reco[*,16:17],24,/nan)*1.0368 gppd = average_arr((reco[*,16:17]-fillnee[*,16:17]),24,/nan)*1.0368 smd = average_arr(qsoil_surf[*,16:17],24,/nan)*100 precipd = average_arr(precip[*,16:17],24,/nan,/tot) preciphour = reform(precip[*,16:17],17568) ch4hourflux = reform(ch4flux_122[*,16:17],17568) ch4hour = reform(nee_ch4_122[*,16:17],17568) thour = reform(tair_30_filled[*,16:17],17568) smhour = reform(qsoil_surf[*,16:17],17568)*100 parhour = reform(par_pref[*,16:17],17568) neehour = reform(fillnee[*,16:17],17568) gpphour = reform(reco[*,16:17]-fillnee[*,16:17],17568) rehour = reform(reco[*,16:17],17568) ;read error (old version) e11 = read_xdf('/data/incoming/wlef/flux/2011/fluxch4err_final_2011.xdf',header=eh) e12 = read_xdf('/data/incoming/wlef/flux/2012/fluxch4err_final_2012.xdf',header=eh) ch4err12 = shift(e12[7,*],-6) ch4err11 = shift(e11[7,*],-6) ch4err = [reform(ch4err11),replicate(nan(),24),reform(Ch4err12)] ch4errd = (sqrt(average_Arr(ch4err^2,24,/nan,/tot))*(1.0368/24)) ch4errd2 = (average_arr(ch4err,24,/nan)*1.0368) ;new err test ch4errf = ch4err ch4errf[where(~finite(ch4err) OR (ch4err LE 0))]=median(ch4err) ch4errac = a_correlate(ch4errf,findgen(24)+1) ch4errf2 = ch4errf ch4errf3 = ch4errf ch4errf = [ch4errf,ch4errf[0:24]] FOR i = 0,17543 DO ch4errf2[i] = sqrt(total(replicate(ch4errf[i]^2,24) + (2 * ch4errac * replicate(ch4errf[i],24) * (ch4errf[i+1:i+24])))) FOR i = 0,17543 DO ch4errf3[i] = sqrt(ch4errf[i]^2+total(2 * ch4errac * replicate(ch4errf[i],24) * (ch4errf[i+1:i+24]))) ch4errd3 = (sqrt(average_arr(ch4errf3^2,24,/tot))*(1.0368/24)) gval = where(finite(ch4d) AND finite(td) AND finite(ch4errd3),ngv,comp=bval) ft = poly_fit(td[gval],ch4d[gval],2,measure_errors=ch4errd3[gval],yfit=ch4fit) ch4m = ft[0] + ft[1]*td + ft[2]*td*td ch4df = merge_array(ch4d,ch4m) Errmod = linfit(ch4d[gval],ch4fit-ch4d[gval]) ch4errmod = abs(errmod[0] + errmod[1] * ch4m) ch4errd[bval] = ch4errmod[bval] ch4errd3[bval] = ch4errmod[bval] ch4errdac = a_correlate(ch4errd3,findgen(24)+1) ch4errd4 = ch4errd3^2 ch4errd5 = [ch4errd3,ch4errd3[0:23]] FOR i = 0,729 DO ch4errd4[i] += total(2 * ch4errdac * replicate(ch4errd5[i],24) * (ch4errd5[i+1:i+24])) ch4errd4 = sqrt(ch4errd4) ; test = regress([transpose(gppd[gval]),transpose(recod[gval]),transpose(co2d[gval]),transpose(pard[gval]),transpose(td[gval])],ch4d[gval],const=ct,mcorrel=mc,correl=cr,yfit=yf) ; test = regress([transpose(average_arr(gppd,7,/nan)),transpose(average_arr(recod,7,/nan)),transpose(average_arr(co2d,7,/nan)),transpose(average_arr(pard,7,/nan)),transpose(average_arr(td,7,/nan))],average_arr(ch4df,7,/nan),const=ct,mcorrel=mc,correl=cr,yfit=yf) gppw = transpose(average_arr(gppd,7,/nan)) rew = transpose(average_arr(recod,7,/nan)) co2w = transpose(average_arr(co2d,7,/nan)) ch4w = (average_arr(ch4df,7,/nan)) ch4werr = sqrt(average_arr(ch4errd3^2,7,/tot,/nan)/7.0) ; ch4werr2 = sqrt((average_arr(ch4errd4^2,7,/tot,/nan))/7.0) tw = transpose(average_arr(td,7,/nan)) parw = transpose(average_arr(pard,7,/nan)) rmod = regress([tw,gppw],ch4w,mcorrel=mc,yfit=yf,measure_e=ch4werr,sigma=s,const=ct) seed = systime(/sec) ests = fltarr(104) FOR p = 0,103 DO ests[p] = stddev((((randomn(seed,1000)*2.0*s[0])+rmod[0])*tw[p])+(((randomn(seed,1000)*2.0*s[1])+rmod[1])*gppw[p])+ct) ;test = regress([transpose(gpphour[gg]),transpose(rehour[gg]),transpose(neehour[gg]),transpose(parhour[gg]),transpose(thour[gg])],ch4hour[gg],const=ct,mcorrel=mc,correl=cr,yfit=yf) ;errorbar doy = findgen(732)+1 ;GOTO,skip ;FIG 1 FLUXES stop ch4dfm = ch4df ch4dfm[where(finite(ch4d))]=nan() !p.multi = [0,1,4] avg_sm=7 red = fsc_color('red',202) blue=fsc_color('blue',201) plot,doy,ch4d,psym=1,ytitle='mgC m-2 day-1',yrange=[-10,20],/nodata,xticks=24,xtickn=['J','F','M','A','M','J','J','A','S','O','N','D','J','F','M','A','M','J','J','A','S','O','N','D',' '],charsize=2 FOR i = 0,729 DO oplot,[i+1,i+1],[ch4df[i]-ch4errd3[i],ch4df[i]+ch4errd3[i]],color=180,thick=3 oplot,doy,ch4d,psym=1,symsize=0.5 oplot,doy,ch4dfm,psym=1,symsize=0.5,color=red oplot,[0,1000],[0,0] oplot,[365,365],[-1000,1000] oplot,doy,ts_smooth(ch4df,avg_sm,order=10),thick=2,color=blue plot,doy,co2d,psym=1,ytitle='gC m-2 day-1',yrange=[-7,3],xticks=24,xtickn=['J','F','M','A','M','J','J','A','S','O','N','D','J','F','M','A','M','J','J','A','S','O','N','D',' '],charsize=2 oplot,[0,1000],[0,0] oplot,[365,365],[-1000,1000] oplot,doy,ts_smooth(co2d,avg_sm,order=10),thick=2,color=blue plot,doy,gppd,psym=1,ytitle='gC m-2 day-1',yrange=[0,20],xticks=24,xtickn=['J','F','M','A','M','J','J','A','S','O','N','D','J','F','M','A','M','J','J','A','S','O','N','D',' '],charsize=2 oplot,[365,365],[-1000,1000] oplot,doy,ts_smooth(gppd,avg_sm,order=10),thick=2,color=blue plot,doy,recod,psym=1,ytitle='gC m-2 day-1',yrange=[0,15],xticks=24,xtickn=['J','F','M','A','M','J','J','A','S','O','N','D','J','F','M','A','M','J','J','A','S','O','N','D',' '],charsize=2 oplot,[365,365],[-1000,1000] oplot,doy,ts_smooth(recod,avg_sm,order=10),thick=2,color=blue stop ;FIG 1b Diurnal ;Compare 2011 and 2012 Summer diurnal cycle for CO2 and CH4, also ;temp? dstart11 = (152-1)*24l dend11 = (243*24l)-1l dstart12 = dstart11+8784l+24l dend12 = (243*24l)-1l+8784l+24l ; dstart11 = (100-1)*24l ; dend11 = (276*24l)-1l ; dstart12 = dstart11+8784l+24l ; dend12 = (276*24l)-1l+8784l+24l ch4d11 = reform(ch4hour[dstart11:dend11],24,92) ch4d11_75 = average_arr(reform(transpose(ch4d11),24*92l),92,/nan,percent=75) ch4d11_25 = average_arr(reform(transpose(ch4d11),24*92l),92,/nan,percent=25) ch4d11 = total(ch4d11,2,/nan)/total(finite(ch4d11),2,/nan) ch4d12 = reform(ch4hour[dstart12:dend12],24,92) ch4d12_75 = average_arr(reform(transpose(ch4d12),24*92l),92,/nan,percent=75) ch4d12_25 = average_arr(reform(transpose(ch4d12),24*92l),92,/nan,percent=25) ch4d12 = total(ch4d12,2,/nan)/total(finite(ch4d12),2,/nan) need11 = reform(neehour[dstart11:dend11],24,92) need11_75 = average_arr(reform(transpose(need11),24*92l),92,/nan,percent=75) need11_25 = average_arr(reform(transpose(need11),24*92l),92,/nan,percent=25) need11 = total(need11,2,/nan)/total(finite(need11),2,/nan) need12 = reform(neehour[dstart12:dend12],24,92) need12_75 = average_arr(reform(transpose(need12),24*92l),92,/nan,percent=75) need12_25 = average_arr(reform(transpose(need12),24*92l),92,/nan,percent=25) need12 = total(need12,2,/nan)/total(finite(need12),2,/nan) gppd11 = reform(gpphour[dstart11:dend11],24,92) gppd11_75 = average_arr(reform(transpose(gppd11),24*92l),92,/nan,percent=75) gppd11_25 = average_arr(reform(transpose(gppd11),24*92l),92,/nan,percent=25) gppd11 = total(gppd11,2,/nan)/total(finite(gppd11),2,/nan) gppd12 = reform(gpphour[dstart12:dend12],24,92) gppd12_75 = average_arr(reform(transpose(gppd12),24*92l),92,/nan,percent=75) gppd12_25 = average_arr(reform(transpose(gppd12),24*92l),92,/nan,percent=25) gppd12 = total(gppd12,2,/nan)/total(finite(gppd12),2,/nan) red11 = reform(rehour[dstart11:dend11],24,92) red11_75 = average_arr(reform(transpose(red11),24*92l),92,/nan,percent=75) red11_25 = average_arr(reform(transpose(red11),24*92l),92,/nan,percent=25) red11 = total(red11,2,/nan)/total(finite(red11),2,/nan) red12 = reform(rehour[dstart12:dend12],24,92) red12_75 = average_arr(reform(transpose(red12),24*92l),92,/nan,percent=75) red12_25 = average_arr(reform(transpose(red12),24*92l),92,/nan,percent=25) red12 = total(red12,2,/nan)/total(finite(red12),2,/nan) td11 = reform(thour[dstart11:dend11],24,92) td11 = total(td11,2,/nan)/total(finite(td11),2,/nan) td12 = reform(thour[dstart12:dend12],24,92) td12 = total(td12,2,/nan)/total(finite(td12),2,/nan) !p.multi = [0,2,2] hrs = findgen(24) red = fsc_color('red',202) blue=fsc_color('blue',201) pink = fsc_color('pink',203) cyan = fsc_color('cyan',204) plot,[0,0],[0,0],/nodata,xrange=[0,23],yrange=[-10,30],ytitle='nmol CH4 m-2 s-1',charsize=1.5 errorbar,hrs,ch4d11,ch4d11_25,ch4d11_75,/over,/noline,fillcolor=cyan,/noadd errorbar,hrs,ch4d12,ch4d12_25,ch4d12_75,/over,/noline,fillcolor=pink,/noadd oplot,[0,25],[0,0] oplot,hrs,ch4d11,color=blue,thick=2 oplot,hrs,ch4d12,color=red,thick=2 plot,[0,0],[0,0],/nodata,xrange=[0,23],yrange=[-20,10],ytitle='umol CO2 m-2 s-1',charsize=1.5 errorbar,hrs,need12,need12_25,need12_75,/over,/noline,fillcolor=pink,/noadd errorbar,hrs,need11,need11_25,need11_75,/over,/noline,fillcolor=cyan,/noadd oplot,[0,25],[0,0] oplot,hrs,need11,color=blue,thick=2 oplot,hrs,need12,color=red,thick=2 plot,[0,0],[0,0],/nodata,xrange=[0,23],yrange=[-1,30],ytitle='umol CO2 m-2 s-1',charsize=1.5 errorbar,hrs,gppd12,gppd12_25,gppd12_75,/over,/noline,fillcolor=pink,/noadd errorbar,hrs,gppd11,gppd11_25,gppd11_75,/over,/noline,fillcolor=cyan,/noadd oplot,[0,25],[0,0] oplot,hrs,gppd11,color=blue,thick=2 oplot,hrs,gppd12,color=red,thick=2 plot,[0,0],[0,0],/nodata,xrange=[0,23],yrange=[0,15],ytitle='umol CO2 m-2 s-1',charsize=1.5 errorbar,hrs,red12,red12_25,red12_75,/over,/noline,fillcolor=pink,/noadd errorbar,hrs,red11,red11_25,red11_75,/over,/noline,fillcolor=cyan,/noadd oplot,[0,25],[0,0] oplot,hrs,red11,color=blue,thick=2 oplot,hrs,red12,color=red,thick=2 ;Fig 1c Flux error stop ;hourly and daily !p.multi = [0,2,2] red = fsc_color('red',200) blue = fsc_color('blue',201) errfit = myfit(abs(ch4hourflux),abs(ch4err),corr=errcorr,/lad) title = 'Err = '+string(errfit[0],format='(f0.2)')+' + '+string(errfit[1],format='(f0.2)')+' * Flux (r2='+string(errcorr^2,format='(f0.2)')+')' plot,[0,0],[0,0],xrange=[0,50],yrange=[0,20],xtitle='Abs(Flux) nmol CH4 m-2 s-1',ytitle='Error CH4 nmol m-2 s-1',title=title,charsize=1.25 oplot,abs(ch4hourflux),abs(ch4err),psym=1,color=80 avgerr = bin_avg(abs(ch4hourflux),abs(ch4err),min=0,max=50,binsize=10,median=mederr) oplot,findgen(5)*10+5,mederr,psym=-2,symsize=2,color=blue,thick=2 oplot,errfit[0]+errfit[1]*findgen(100),color=red,thick=2 print,errfit,errcorr errfitd = myfit(abs(ch4dflux),abs(ch4errd3),corr=errcorrd,/lad) title = 'Err = '+string(errfitd[0],format='(f0.2)')+' + '+string(errfitd[1],format='(f0.2)')+' * Flux (r2='+string(errcorrd^2,format='(f0.2)')+')' plot,[0,0],[0,0],xrange=[0,15],yrange=[0,5],xtitle='Abs(Flux) mgC m-2 day-1',ytitle='Error mgC m-2 day-1',title=title,charsize=1.25 oplot,abs(ch4d),abs(ch4errd3),psym=1,color=80 avgerr = bin_avg(abs(ch4dflux),abs(ch4errd3),min=0,max=20,binsize=5,median=mederr) oplot,findgen(4)*5+2.5,mederr,psym=-2,symsize=2,color=blue,thick=2 oplot,errfitd[0]+errfitd[1]*findgen(100),color=red,thick=2 print,errfitd,errcorrd ;Fig 1d fit of GPP (linear) and temperature (exp) to CH4, hourly and ;daily ;DO THIS stop hour = [fjday[*,16],fjday[*,17]] gv_h = where(finite(gpphour) AND finite(ch4hour) AND finite(thour) AND finite(ch4err) AND (gpphour GT 0) AND (ch4hour GT -50) AND (ch4hour LT 50) AND (thour GT 0) AND (gpphour LT 30)) gv_d = where(finite(gppd) AND finite(ch4d) AND finite(td) AND finite(ch4errd3) AND (gppd GT 0) AND (td GT 0)) gppfit_h = myfit(gpphour[gv_h],ch4hour[gv_h],degree=1,corr=gppfit_h_c,measure_error=abs(ch4err[gv_h]),sigma=gppfit_h_s) gppfit_h = myfit(loc,bin_avg(gpphour[gv_h],ch4hour[gv_h],binsize=2.0,loc=loc),degree=1,corr=gppfit_h_c2,sigma=gppfit_h_s,measure_error=sqrt(zapbadval(bin_avg(gpphour[gv_h],ch4err[gv_h]^2,binsize=2.0)))) gppfit_d = myfit(gppd[gv_d],ch4d[gv_d],degree=1,corr=gppfit_d_c,measure_error=abs(ch4errd3[gv_d]),sigma=gppfit_d_s) tfit_h = [0.36,0.2] tfit_d = [0.36,0.1] tstats_h = lmfit(thour[gv_h],ch4hour[gv_h],tfit_h,measure_error=abs(ch4err[gv_h]),sigma=tfit_h_s,FUNCTION='wlef_ch4paper_efit',/double,itmax=100) tstats_h2 = lmfit(loc,bin_avg(thour[gv_h],ch4hour[gv_h],binsize=1.0,loc=loc),tfit_h,measure_error=sqrt(zapbadval(bin_avg(thour[gv_h],ch4err[gv_h]^2,binsize=1.0))),sigma=tfit_h_s,FUNCTION='wlef_ch4paper_efit',/double) tstats_d = lmfit(td[gv_d],ch4d[gv_d],tfit_d,measure_error=abs(ch4errd3[gv_d]),sigma=tfit_d_s,FUNCTION='wlef_ch4paper_efit',/double,itmax=100) tfit_h_c = correl(Tstats_h,ch4hour[gv_h]) tfit_d_c = correl(tstats_d,ch4d[gv_d]) ;do same for t_d with exp fit gpprange_h = findgen(50) gppline_h = gppfit_h[0] + gppfit_h[1] * gpprange_h gpperr_h = fltarr(50) seed = systime(/sec) FOR p = 0,49 DO gpperr_h[p] = stddev((((randomn(seed,1000)*2.0*gppfit_h_s[0])+gppfit_h[0]))+(((randomn(seed,1000)*2.0*gppfit_h_s[1])+gppfit_h[1])*gpprange_h[p])) gpprange_d = findgen(50) gppline_d = gppfit_d[0] + gppfit_d[1] * gpprange_d gpperr_d = fltarr(50) seed = systime(/sec) FOR p = 0,49 DO gpperr_d[p] = stddev((((randomn(seed,1000)*2.0*gppfit_d_s[0])+gppfit_d[0]))+(((randomn(seed,1000)*2.0*gppfit_d_s[1])+gppfit_d[1])*gpprange_d[p])) trange_h = findgen(70)-30 tline_h = (wlef_ch4paper_efit(trange_h,tfit_h))[*,0] terr_h = fltarr(70) seed = systime(/sec) FOR p = 0,69 DO terr_h[p] = stddev(((randomn(seed,1000)*2.0*tfit_h_s[0])+tfit_h[0])*exp(((randomn(seed,1000)*2.0*tfit_h_s[1])+tfit_h[1]) * trange_h[p])) trange_d = findgen(70)-30 tline_d = (wlef_ch4paper_efit(trange_d,tfit_d))[*,0] terr_d = fltarr(70) seed = systime(/sec) FOR p = 0,69 DO terr_d[p] = stddev(((randomn(seed,1000)*2.0*tfit_d_s[0])+tfit_d[0])*exp(((randomn(seed,1000)*2.0*tfit_d_s[1])+tfit_d[1]) * trange_d[p])) !p.multi = [0,2,2] red=fsc_color('red',202) blue=fsc_color('blue',201) pink = fsc_color('pink',203) cyan = fsc_color('cyan',204) title = 'CH4 = '+string(gppfit_h[0],format='(f0.2)') + ' + ' + string(gppfit_h[1],format='(f0.2)') + ' GPP (r2=' + string(gppfit_h_c^2,format='(f0.2)') + ')' plot,[0,0],[0,0],xrange=[0,30],yrange=[-40,40],xtitle='umol CO2 m-2 s-1',ytitle='nmol CH4 m-2 s-1',title=title,/nodata,charsize=1.25 FOR i = 0,n_elements(gv_h)-1 DO oplot,[(gpphour[gv_h])[i],(gpphour[gv_h])[i]],[(ch4hour[gv_h])[i]-(ch4err[gv_h])[i],(ch4hour[gv_h])[i]+(ch4err[gv_h])[i]],color=180 errorbar,gpprange_h,gppline_h,gpperr_h,gpperr_h,/over,/noline,fillcolor=pink oplot,gpphour[gv_h],ch4hour[gv_h],psym=3,color=0 oplot,gpprange_h,gppline_h,color=red,thick=2 title = 'CH4 = '+string(gppfit_d[0],format='(f0.2)') + ' + ' + string(gppfit_d[1],format='(f0.2)') + ' GPP (r2=' + string(gppfit_d_c^2,format='(f0.2)') + ')' plot,[0,0],[0,0],xrange=[0,15],yrange=[-10,20],xtitle='gC m-2 day-1',ytitle='mgC m-2 day-1',title=title,/nodata,charsize=1.25 FOR i = 0,n_elements(gv_d)-1 DO oplot,[(gppd[gv_d])[i],(gppd[gv_d])[i]],[(ch4d[gv_d])[i]-(ch4err[gv_d])[i],(ch4d[gv_d])[i]+(ch4err[gv_d])[i]],color=180 errorbar,gpprange_d,gppline_d,gpperr_d,gpperr_d,/over,/noline,fillcolor=pink oplot,gppd[gv_d],ch4d[gv_d],psym=3,color=0 oplot,gpprange_d,gppline_d,color=red,thick=2 title = 'CH4 = '+string(tfit_h[0],format='(f0.2)') + 'exp(' + string(tfit_h[1],format='(f0.2)') + ' T) (r2=' + string(tfit_h_c^2,format='(f0.2)') + ')' plot,[0,0],[0,0],xrange=[0,30],yrange=[-40,40],xtitle='degrees C',ytitle='nmol m-2 s-1',title=title,/nodata,charsize=1.25 errorbar,trange_h,tline_h,terr_h,terr_h,/over,/noline,fillcolor=pink FOR i = 0,n_elements(gv_h)-1 DO oplot,[(thour[gv_h])[i],(thour[gv_h])[i]],[(ch4hour[gv_h])[i]-(ch4err[gv_h])[i],(ch4hour[gv_h])[i]+(ch4err[gv_h])[i]],color=180 oplot,thour[gv_h],ch4hour[gv_h],psym=3,color=0 oplot,trange_h,tline_h,color=red,thick=2 title = 'CH4 = '+string(tfit_d[0],format='(f0.2)') + 'exp(' + string(tfit_d[1],format='(f0.2)') + ' T) (r2=' + string(tfit_d_c^2,format='(f0.2)') + ')' plot,[0,0],[0,0],xrange=[0,30],yrange=[-10,20],xtitle='degrees C',ytitle='mgC m-2 day-1',title=title,/nodata,charsize=1.25 FOR i = 0,n_elements(gv_d)-1 DO oplot,[(td[gv_d])[i],(td[gv_d])[i]],[(ch4d[gv_d])[i]-(ch4err[gv_d])[i],(ch4d[gv_d])[i]+(ch4err[gv_d])[i]],color=180 errorbar,trange_d,tline_d,terr_d,terr_d,/over,/noline,fillcolor=pink oplot,td[gv_d],ch4d[gv_d],psym=3,color=0 oplot,trange_d,tline_d,color=red,thick=2 ;see wlef_efit lmfit ;q10 is daily is 2.28, base 1.45, for hourly is 2.61, base 1.29 ;q10 for Reco is 3.03 , base is 1.22 for daily, 2.79 for hourly, 1.25 base ;Fig 1e Recreate boxplot - include werner, read chambers ;boxplot function in coyote - ch4hour, werner from april-october, lf, ; sf, wf ;werner = [3.23,-1.29,6.36,10.16,15.55,17.04,18.56,14.84,13.58,11.26,1.04,3.94]/days_in_month(findgen(12)+1) ;ADD IN REST OF DATA (see peter's spreadsheet) stop ; sites = ['lc_means','wf_means','sf_means'] ; lcdata = ((read_ascii(mydocs()+'methane/chamber/'+sites[0])).(0))[3,*]*1000.0 ; wfdata = ((read_ascii(mydocs()+'methane/chamber/'+sites[1])).(0))[3,*]*1000.0 ; sfdata = ((read_ascii(mydocs()+'methane/chamber/'+sites[2])).(0))[3,*]*1000.0 chamberdata = read_csv(mydocs()+'methane/chamber/allchambers.csv') csite = strupcase(chamberdata.field1) cch4 = chamberdata.field5 * 1e9/86400. ctemp = chamberdata.field6 lcdata = cch4[where(csite EQ 'LC')] sfdata = cch4[where(csite EQ 'SF')] wfdata = cch4[where(csite EQ 'WF')] csdata = cch4[where(csite EQ 'CS')] rcdata = cch4[where(csite EQ 'RC')] tcdata = cch4[where(csite EQ 'TC')] wcdata = cch4[where(csite EQ 'WC')] wetdata = [lcdata,sfdata,wfdata,csdata] updata = [rcdata,tcdata,wcdata] ; upscale = (randomn(systime(/sec),10000)*stddev(wetdata)+mean(wetdata))*0.28*1.0368*179.0 ;*0.28*1.0368*179.0 ;179 days hour = [fjday[*,16],fjday[*,17]] gv_hbox = where(finite(ch4hour) and (ch4hour GT -50) AND (ch4hour LT 50) AND (hour GE 100) AND (hour LT 277)) towerd = ch4hour[gv_hbox] boxy = ptrarr(8) boxy[0] = ptr_new(lcdata) boxy[1] = ptr_new(wfdata) boxy[2] = ptr_new(sfdata) boxy[3] = ptr_new(csdata) boxy[4] = ptr_new(rcdata) boxy[5] = ptr_new(tcdata) boxy[6] = ptr_new(wcdata) boxy[7] = ptr_new(towerd) !p.multi = [0,1,2] boxplot,boxy,/fill,label=['LC','WF','SF','CS','RC','TC','WC','Tower'],ytitle='nmol m-2 s-1',yrange=[-50,50] oplot,[0,10],[0,0] oplot,[3.5,3.5],[-100,100] oplot,[7.5,7.5],[-100,100] xyouts,0.5,-20,'Wetland',charsize=2 xyouts,4,-20,'Upland',charsize=2 werndata = ([3.23,-1.29,6.36,10.16,15.55,17.04,18.56,14.84,13.58,11.26,1.04,3.94])*days_in_month(findgen(12)+1) ;/1.0368 ;APril-Oct ;split into separate figure mmo = reform(mo[0:*:24,16:17],732) ch4month = bin_Avg(mmo,ch4df)*days_in_month(findgen(12)+1) allmo = fltarr(3,12) allmo[0,*] = ch4month allmo[1,*] = werndata red = fsc_color('red',200) plot,[0,0],[0,0],/nodata,xrange=[0,36],yrange=[-40,600],xticks=1,xtickn=[' ',' '],charsize=1.25,ytitle='mg C m-2 mon-1' my_bar_plot,allmo,color=replicate_arr([0,200,255],12),background=255,barnames=['J',' ',' ','F',' ',' ','M',' ',' ','A',' ',' ','M',' ',' ','J',' ',' ','J',' ',' ','A',' ',' ','S',' ',' ','O',' ',' ','N',' ',' ','D',' ',' '],/over oplot,[0,1000],[0,0] ;Fig 1f from Ke - WPL ch4 vs in line CH4 ;FIG 2 FORCING stop mmo = mo[0:*:24,16] precipmo = [bin_avg(mmo,precipd[0:365]),bin_avg(mmo,precipd[366:*])] precipmo *= [days_in_month(findgen(12)+1),days_in_month(findgen(12)+1)] !p.multi = [0,1,4] red=fsc_color('red',202) blue=fsc_color('blue',201) pink = fsc_color('pink',203) cyan = fsc_color('cyan',204) plot,doy,td,psym=1,ytitle='degrees C',yrange=[-25,30],xticks=24,xtickn=['J','F','M','A','M','J','J','A','S','O','N','D','J','F','M','A','M','J','J','A','S','O','N','D',' '],charsize=2 oplot,[365,365],[-1000,1000] oplot,doy,ts_smooth(td,avg_sm,order=10),thick=2,color=blue plot,doy,pard,psym=1,ytitle='MJ day-1',yrange=[0,35],xticks=24,xtickn=['J','F','M','A','M','J','J','A','S','O','N','D','J','F','M','A','M','J','J','A','S','O','N','D',' '],charsize=2 oplot,[365,365],[-1000,1000] oplot,doy,ts_smooth(pard,avg_Sm,order=10),thick=2,color=blue ; oplot,[365,365],[-1000,1000] ; oplot,doy,ts_smooth(precipd,avg_Sm,order=10),thick=2,color=blue plot,doy,precipd,psym=1,ytitle='mm',yrange=[0,160],/nodata,xticks=24,xtickn=['J','F','M','A','M','J','J','A','S','O','N','D','J','F','M','A','M','J','J','A','S','O','N','D',' '],charsize=2 my_bar_plot,precipmo,/over,colors=[replicate(cyan,12),replicate(pink,12)] oplot,[365,365],[-1000,1000] ; oplot,doy[0:365],total(precipd[0:365],/cum)/10.,thick=2,color=blue ; oplot,doy[366:*],total(precipd[366:*],/cum)/10.,thick=2,color=red plot,doy,smd,psym=1,ytitle='% VWC',yrange=[0,50],xticks=24,xtickn=['J','F','M','A','M','J','J','A','S','O','N','D','J','F','M','A','M','J','J','A','S','O','N','D',' '],charsize=2 oplot,[365,365],[-1000,1000] oplot,doy,ts_smooth(smd,avg_Sm,order=10),thick=2,color=blue ;FIG 3 cumulative CH4 flux with quadrature error ;(sqrt(total(err^2,/cum,/nan))) compare stop !p.multi = [0,1,2] red=fsc_color('red',202) blue=fsc_color('blue',201) pink = fsc_color('pink',203) cyan = fsc_color('cyan',204) plot,[0,0],[0,0],xrange=[1,365],yrange=[0,1200],xtitle='Month',ytitle='mgC m-2',xticks=12,xtickn=['J','F','M','A','M','J','J','A','S','O','N','D',' '],charsize=1.25 ; errorbar,doy[0:364],total(ch4df[0:364],/cum,/nan),total(ch4errd[0:364],/cum,/nan),total(ch4errd[0:364],/cum,/nan),/over,/noline,fillcolor=cyan ; errorbar,doy[0:364],total(ch4df[365:*],/cum,/nan),total(ch4errd[365:*],/cum,/nan),total(ch4errd[365:*],/cum,/nan),/over,/noline,fillcolor=pink errorbar,doy[0:364],total(ch4df[0:364],/cum,/nan),sqrt(total(ch4errd4[0:364]^2,/cum,/nan)),sqrT(total(ch4errd4[0:364]^2,/cum,/nan)),/over,/noline,fillcolor=cyan errorbar,doy[0:364],total(ch4df[365:*],/cum,/nan),sqrt(total(ch4errd4[365:*]^2,/cum,/nan)),sqrt(total(ch4errd4[365:*]^2,/cum,/nan)),/over,/noline,fillcolor=pink oplot,doy[0:364],total(ch4df[0:364],/cum,/nan),color=blue,thick=2 oplot,doy[0:364],total(ch4df[365:*],/cum,/nan),color=red,thick=2 ;FIG 4 spectrum(co2d,1,freq=fq) , ;spectral smoothing with hanning window ;is short-term temporal variation in CH4 > CO2? ;or DEMD/HILBERT of normalized daily flux - marginal power spec emdch4 = demd(normalize(ch4d),/nodel) emdco2 = demd(normalize(co2d),/nodel) hch4 = hilbert_spec(emdch4,1,marginal=mch4,freq=fch4) hco2 = hilbert_spec(emdco2,1,marginal=mco2,freq=fco2) f_hch4 = filter_nd(hch4,11,'hanning') f_hco2 = filter_nd(hco2,11,'hanning') fhch4 = total(f_hch4,2,/nan) fhco2 = total(f_hco2,2,/nan) ; red = fsc_color('red') plot,1./fco2,fhco2*fco2,/ylog,/xlog,xrange=[2.,250],yrange=[0.03,4],thick=2,xtitle='Period (days)',ytitle='Normalized Power*Freq',charsize=1.25 oplot,1./fch4,fhch4*fch4,thick=2,color=red stop ;PART B ;FIG 5 Correlation analysis ;skip correctional analysis and regression analysis (replace with ;driver plots) GOTO,skip tcorr = lagcorr(thour,ch4hour,pval=tpval,avg=[1,24,24*7,24*30],lags=[0]) parcorr = lagcorr(parhour,ch4hour,pval=parpval,avg=[1,24,24*7,24*30],lags=[0]) precipcorr = lagcorr(preciphour,ch4hour,pval=precippval,avg=[1,24,24*7,24*30],lags=[0]) smcorr = lagcorr(smhour,ch4hour,pval=smpval,avg=[1,24,24*7,24*30],lags=[0]) neecorr = lagcorr(neehour,ch4hour,pval=neepval,avg=[1,24,24*7,24*30],lags=[0]) gppcorr = lagcorr(gpphour,ch4hour,pval=gpppval,avg=[1,24,24*7,24*30],lags=[0]) recorr = lagcorr(rehour,ch4hour,pval=repval,avg=[1,24,24*7,24*30],lags=[0]) allcorr = [tcorr,parcorr,smcorr,neecorr,gppcorr,recorr] allpval = [tpval,parpval,smpval,neepval,gpppval,repval] goodcorr = allcorr goodcorr[where(allpval GT 0.1)]=0 !p.multi = [0,2,2] loadct,33 tvlct,rd,gr,bl,/get rd[127] = 255 gr[127] = 255 bl[127] = 255 tvlct,rd,gr,bl !p.background = 127 tvplot,transpose(goodcorr),zrange=[-0.9,0.9],xticks=4,xtickv=numgen(0.0,3.0,0.75)+0.375,xtickn=['Hr','Dy','Wk','Mon'],xtitle='Averaging Period',yticks=6,ytickv=numgen(0.0,5.0,0.833)+0.416,ytickn=['T','PAR','SoilQ','NEE','GPP','Reco'],ntext=6,/order,charsize=2 stop loadct,0 !p.background=255 ; title = 'CH4 = '+string(ct,format='(f0.2)')+' + '+string(rmod[0],format='(f0.2)')+' * T + '+string(rmod[1],format='(f0.2)')+' * GPP'+' (r2='+string(mc^2,format='(f0.2)')+')' ; plot,ch4w,yf,/nodata,yrange=[-2,12],xrange=[-2,12],xtitle='Observed gC m-2 day-1',ytitle='Modeled gC m-2 day-1',title=title ; for i = 0,103 do oplot,[ch4w[i]-ch4werr[i],ch4w[i]+ch4werr[i]],[yf[i],yf[i]] ; for i = 0,103 do oplot,[ch4w[i],ch4w[i]],[yf[i]+ests[i],yf[i]-ests[i]] ; oplot,[-2,12],[-2,12] stop skip: ;other paper plots: ;ch4 diurnal cycle vs co2/t/par diurnal cycle ;compare separately for each year ;UPDATE THIS SECTION - only do DLEM comparison at two scales ;FIG 6 MODELS ;2011-2012 compared to models (inverse and forward and diagnostic) ;Read FUNG (interpolate to daily - spline) wlef_readfung,wetland_fung=wetland_fung,soil_fung=soil_fung ;kgC/m2/mo (maybe day?) wetland_fung(where(~finite(wetland_fung)))=0 fungch4 = ((wetland_fung*1e6)+(replicate(soil_fung*1e6/30.4167,12)))*(12./16.) ;Read WERNER (scan) werner = [3.23,-1.29,6.36,10.16,15.55,17.04,18.56,14.84,13.58,11.26,1.04,3.94]*days_in_month(findgen(12)+1) ;Read CT (interpolate?) ;save,filename='/data/projects/carbontracker/ct_ch4_wlef.sav',ch4ct_bio,ch4ct_ff,ch4ct_fire,ch4ct_ocn,ch4ct_ag restore,mydocs()+'carbontracker/ct_ch4_wlef.sav' ; ct_bio = (total(total(total(ch4ct_bio,2)/11.,2)/11.,2)/11.)*days_in_month(findgen(12)+1) ct_bio = (total(total(total(ch4ct_bio[*,*,3:7,3:7],2)/11.,2)/5.,2)/5.)*days_in_month(findgen(12)+1) ct_bios=stddev(total(total(ch4ct_bio[*,*,3:7,3:7],3)/5.,3)/5.,dim=2)*days_in_month(findgen(12)+1) ct_fire = (total(total(total(ch4ct_fire[*,*,3:7,3:7],2)/11.,2)/5.,2)/5.)*days_in_month(findgen(12)+1) ct_ff = (total(total(total(ch4ct_ff[*,*,3:7,3:7],2)/11.,2)/5.,2)/5.)*days_in_month(findgen(12)+1) ct_ocn = (total(total(total(ch4ct_ocn[*,*,3:7,3:7],2)/11.,2)/5.,2)/5.)*days_in_month(findgen(12)+1) ct_ag = (total(total(total(ch4ct_ag[*,*,3:7,3:7],2)/11.,2)/5.,2)/5.)*days_in_month(findgen(12)+1) stop ch4df = ch4df[0:729] ch4mo = mo_avg(total(reform(ch4df,365,2),2,/nan)/2.)*days_in_month(findgen(12)+1) ch4moerr = sqrt(mo_avg(total(reform(ch4errd4[0:729]^2,365,2),2)/2.)*days_in_month(findgen(12)+1)) mloc = (numgen(0.,335,30.4167))+15.2084 ;Read DLEM dlemsite = ((read_ascii(mydocs()+'methane/dlem/DLEM-site.csv')).(0))[0:729] dlemreg = ((read_ascii(mydocs()+'methane/dlem/DLEM-regional.csv')).(0)) dlemmo = mo_avg(total(reform(dlemsite,365,2),2,/nan)/2.)*days_in_month(findgen(12)+1) dlemr_av = total(reform(dlemreg,365,11),2)/11. dlemr_sd = stddev(reform(dlemreg,365,11),dim=2) dlemr_mo_all = reform(mo_avg(dlemreg),12,11) dlemr_mo = total(dlemr_mo_all,2)/11.*1e3*days_in_month(findgen(12)+1) dlemr_mo_sd = stddev(dlemr_mo_all,dim=2)*1e3*days_in_month(findgen(12)+1) stop !p.multi = [0,2,2] blue = fsc_color('blue',201) red = fsc_color('red',202) green = fsc_color('green',203) ; pur = fsc_color('purple',204) cyan = fsc_color('cyan',205) pink = fsc_color('pink',206) ; salmon = fsc_color('salmon',207) mos = findgen(12)+1 ; plot,[0,0],[0,0],xrange=[1,12],yrange=[-50,600],xtitle='Month',ytitle='mgC m-2 mo-1' ; errorbar,mos,ch4mo,ch4moerr,ch4moerr,/over,/noline,fillcolor=199 ; errorbar,mos,ct_bio,ct_bios,ct_bios,/over,/noline,fillcolor=cyan ; errorbar,mos,dlemr_mo,dlemr_mo_sd,dlemr_mo_sd,/over,/noline,fillcolor=pink ; oplot,mos,ch4mo,thick=3 ; oplot,mos,dlemr_mo,thick=3,color=red ; oplot,mos,ct_bio,thick=3,color=blue ; oplot,mos,werner,thick=3,color=green ; errorbar,mos,ch4mo,ch4moerr,ch4moerr,yrange=[-50,600],charsize=1,thick=2,ytitle='mgC m!U2!N mo!U-1!N',xtitle='Month',xrange=[0.5,12.5] ; oplot,mos,ct_bio,color=blue,thick=3 ; oplot,mos,dlemmo,color=red,thick=3 ; oplot,mos,werner,color=green,thick=3 !p.multi = [0,1,2] blue = fsc_color('blue',201) red = fsc_color('red',202) pink = fsc_color('pink',203) plot,doy,ch4d,psym=1,ytitle='mgC m-2 day-1',yrange=[-10,20],/nodata,xticks=24,xtickn=['J','F','M','A','M','J','J','A','S','O','N','D','J','F','M','A','M','J','J','A','S','O','N','D',' '],charsize=1.25 FOR i = 0,729 DO oplot,[i+1,i+1],[ch4df[i]-ch4errd3[i],ch4df[i]+ch4errd3[i]],color=180,thick=3 oplot,doy,ch4d,psym=1,symsize=0.5 oplot,doy,ch4dfm,psym=1,symsize=0.5,color=pink oplot,doy,dlemsite,color=red,thick=3 oplot,doy,dlemr_av*1e3,color=blue,thick=3 oplot,[0,800],[0,0] ; qq = total(reform(ch4df,365,2),2)/2.0 ; qq2 = (total(reform(ch4errd3,365,2),2)/2.0) ; ct_bio2 = (total(total(total(ch4ct_bio[*,*,3:7,3:7],2)/11.,2)/5.,2)/5.) ; ct_biod = interpol([ct_bio2[11],ct_bio2,ct_bio2[0]],[-16.0,mloc,380.0],findgen(365)+1,/spline) ; errorbar,findgen(365)+1,smooth(qq,5,/nan,/edge),smooth(qq2,5,/nan,/edge),smooth(qq2,5,/nan,/edge),charsize=2,thick=2,yrange=[-1,15],ytitle='mgC m!U2!N day!U-1!N',xtitle='Day of Year',title='Carbontracker interpolated=blue, weekly eddy flux=black' ; oplot,findgen(365)+1,ct_biod,color=blue,thick=3 ; oplot,findgen(365)+1,dlemsite,color=red,thick=3 stop END PRO wlefch4_paper2 restore,'/data/incoming/wlef/alldata19952013-processed.sav' gg = where(mo GE 6 AND mo LE 8 AND yr GE 2011 AND yr LE 2011,ngg) ndays = ngg / 24 test_c122 = reform(nee_co2_122_n[gg],24,ndays) test_p122 = reform(c_pref[gg],24,ndays) test_f122 = reform(co2flux_122[gg],24,ndays) test_s122 = reform(co2stor_122[gg],24,ndays) test_m122 = reform(nee_ch4_122[gg],24,ndays) test_ms122 = reform(ch4stor_122[gg],24,ndays) test_mf122 = reform(ch4flux_122[gg],24,ndays) test_mp122 = reform(ch4_prefflag[gg],24,ndays) test_pp122 = reform(c_prefflag[gg],24,ndays) test_ch4122 = reform(ch4_122[gg],24,ndays) test_ch430 = reform(ch4_30[gg],24,ndays) test_ch4396 = reform(ch4_396[gg],24,ndays) test_co2122 = reform(co2_122[gg],24,ndays) test_co230 = reform(co2_30[gg],24,ndays) test_co2396 = reform(co2_396[gg],24,ndays) badc = where(~finite(test_m122) OR ~finite(test_ms122) OR ~finite(test_mf122) OR test_mp122 LE 0) test_m122[badc] = nan() test_ms122[badc] = nan() test_mf122[badc] = nan() badcc = where(~finite(test_c122) OR ~finite(test_s122) OR ~finite(test_f122) OR ~finite(test_p122)) test_c122[badcc] = nan() test_s122[badcc] = nan() test_f122[badcc] = nan() test_p122[badcc] = nan() ;profile plot !p.multi=[0,2,2] blue = fsc_color('blue',200) red = fsc_color('red',201) plot,total(test_ch430,2,/nan)/total(finite(test_ch430),2,/nan),yrange=[1.9,2.0],charsize=1.5,xtitle='Local standard time',ytitle='ppm CH4',title='CH4 concentration black=30 blue=122 red=396',thick=3 oplot,total(test_ch4122,2,/nan)/total(finite(test_ch4122),2,/nan),color=blue,thick=3 oplot,total(test_ch4396,2,/nan)/total(finite(test_ch4396),2,/nan),color=red,thick=3 plot,total(test_co230,2,/nan)/total(finite(test_co230),2,/nan),yrange=[375,420],charsize=1.5,xtitle='Local standard time',ytitle='ppm CO2',title='CO2 concentration',thick=3 oplot,total(test_co2122,2,/nan)/total(finite(test_co2122),2,/nan),color=blue,thick=3 oplot,total(test_co2396,2,/nan)/total(finite(test_co2396),2,/nan),color=red,thick=3 plot,total(test_m122,2,/nan)/total(finite(test_m122),2,/nan),yrange=[-15,15],charsize=1.5,xtitle='Local standard time',ytitle='nmol CH4 m-2 s-1',title='CH4 flux black=NEE red=flux blue=storage',thick=3 oplot,total(test_ms122,2,/nan)/total(finite(test_ms122),2,/nan),color=blue,thick=3 oplot,total(test_mf122,2,/nan)/total(finite(test_mf122),2,/nan),color=red,thick=3 oplot,[0,24],[0,0] plot,total(test_c122,2,/nan)/total(finite(test_c122),2,/nan),yrange=[-15,10],charsize=1.5,xtitle='Local standard time',ytitle='umol CO2 m-2 s-1',title='CO2 flux',thick=3 oplot,total(test_s122,2,/nan)/total(finite(test_s122),2,/nan),color=blue,thick=3 oplot,total(test_f122,2,/nan)/total(finite(test_f122),2,/nan),color=red,thick=3 oplot,[0,24],[0,0] stop !p.multi = 0 ;pref vs c122 plot,total(test_c122,2,/nan)/total(finite(test_c122),2,/nan),yrange=[-15,10],charsize=1.5,xtitle='Local standard time',ytitle='umol CO2 m-2 s-1',title='CO2 flux',thick=3 oplot,total(test_p122,2,/nan)/total(finite(test_p122),2,/nan),color=blue,thick=3 oplot,[0,24],[0,0] END