;pro example1 ;Program: Gives an example of linear regression trend analysis ; of a time series using the IDL procedure "ztrend.pro". ; The trend model includes seasonal cycle, linear trend, ; QBO, and SOLAR proxy terms. (please read documentation ; in procedure "ZTREND.PRO") ; This example program also plots out the various ; quantities generated by ztrend.pro (as a postscript file). ;The data: ; Nimbus 7 TOMS version 7 1979-1992 monthly zonal mean ; total column ozone time series averaged between 40N and 45N. ; (Data Source: see http://toms.gsfc.nasa.gov/) ;To run this example program: ; 1) Make sure the following files are in your current ; directory: ; ztrend.pro ; data.txt ; example1.pro (this program) ; 2) Get into IDL and just type ; .r example1 ; 3) The output graphics postscript file will be: idl.ps ;Define necessary constant parameters for procedure ztrend.pro: IERROR=1 ; <== Do multivariate statistical analysis for coefficients N=168 ; <== Total number of months in time series (1979-1992) NUM=12 ; <== Number of samples per year in time series TMIN=50 ; <== Ozone values less than 50 Dobson units are flagged missing TMAX=600 ; <== Ozone values greater than 600 Dobson units are flagged missing M1=9 ; <== 9 constants for seasonal fit coefficient M2=7 ; <== 7 " " linear trend coefficient M3=7 ; <== 7 " " QBO coefficient M4=7 ; <== 7 " " SOLAR coefficient M5=0 ; <== 0 " " EXTRA1 coefficient (EXTRA1 not used) M6=0 ; <== 0 " " EXTRA2 coefficient (EXTRA2 not used) M=M1+M2+M3+M4+M5+M6 ;Define arrays for procedure ztrend.pro: T=fltarr(N) & QBO=fltarr(N) & SOLAR=fltarr(N) EXTRA1=fltarr(N) & EXTRA2=fltarr(N) ALPHA=FLTARR(NUM) & BETA=FLTARR(NUM) & GAMMA=FLTARR(NUM) DELTA=FLTARR(NUM) & EPS1=FLTARR(NUM) & EPS2=FLTARR(NUM) VAR_ALPHA=FLTARR(NUM) & VAR_BETA=FLTARR(NUM) VAR_GAMMA=FLTARR(NUM) & VAR_DELTA=FLTARR(NUM) VAR_EPS1=FLTARR(NUM) & VAR_EPS2=FLTARR(NUM) ALPHA_FIT=FLTARR(N) & BETA_FIT=FLTARR(N) & GAMMA_FIT=FLTARR(N) DELTA_FIT=FLTARR(N) & EPS1_FIT=FLTARR(N) & EPS2_FIT=FLTARR(N) RESIDUAL_FIT=FLTARR(N) VAR_FACTOR=FLTARR(NUM) ;Now read in Ozone time series to be analyzed, and QBO and SOLAR proxies: stuff=strarr(1) openr,2,'data.txt' readf,2,stuff & readf,2,t readf,2,stuff & readf,2,qbo readf,2,stuff & readf,2,solar close,2 ;Now do the trend analysis: ZTREND, $ IERROR,N,M,NUM,T,TMIN,TMAX, $ M1,M2,M3,M4,M5,M6, $ QBO,SOLAR,EXTRA1,EXTRA2, $ ALPHA,BETA,GAMMA,DELTA,EPS1,EPS2, $ VAR_ALPHA,VAR_BETA,VAR_GAMMA,VAR_DELTA,VAR_EPS1,VAR_EPS2, $ ALPHA_FIT,BETA_FIT,GAMMA_FIT,DELTA_FIT,EPS1_FIT,EPS2_FIT, $ RESIDUAL_FIT,VAR_FACTOR, $ F,WK,A,AINV,S,C,COV,CS ;Now print out the derived coefficients with their 2-sigma values: print,'Seasonal Cycle Coefficient (Dobson Units):' print,' Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec' print,alpha,format='(12I5)' print,'Seasonal Cycle Coefficient 2-SIGMA (Dobson Units):' print,' Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec' print,2*sqrt(var_alpha),format='(12f5.1)' print,' ' print,'Trend Coefficient (% per decade):' print,' Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec' print,beta*12000/alpha,format='(12f5.1)' print,'Trend Coefficient 2-SIGMA (% per decade):' print,' Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec' print,2*sqrt(var_beta)*12000/alpha,format='(12f5.1)' print,' ' print,'QBO Coefficient (DU per 10m/s):' print,' Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec' print,gamma*10,format='(12f5.1)' print,'QBO Coefficient 2-SIGMA (DU per 10m/s):' print,' Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec' print,2*sqrt(var_gamma)*10,format='(12f5.1)' print,' ' print,'SOLAR Coefficient (DU per 100F10.7):' print,' Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec' print,delta*100,format='(12f5.1)' print,'SOLAR Coefficient 2-SIGMA (DU per 100F10.7):' print,' Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec' print,2*sqrt(var_delta)*100,format='(12f5.1)' ;Plotting setup: !p.font=1 !p.title='Nimbus 7 TOMS 1979-1992 40N to 45N Avg Zonal Mean Data' !p.multi=[0,1,2] set_plot,'ps' device,/landscape xstf=['J','F','M','A','M','J','J','A','S','O','N','D'] xstf2=['79','80','81','82','83','84','85','86','87','88',$ '89','90','91','92','93'] ;Plot original time series and regression model fit: plot,t,/xs,charsize=1.1,xcharsize=1.1,ycharsize=1.1,$ xthick=3,ythick=3,thick=3,psym=-6,symsize=0.7,$ xtitle='Year',ytitle='Dobson Units',$ xticks=14,xtickname=xstf2,/ys,yrange=[250,450],min_value=1 q=alpha_fit+beta_fit+gamma_fit+delta_fit+eps1_fit+eps2_fit oplot,q,linestyle=1,thick=6 oplot,[20,50],[430,430],linestyle=0,thick=3,psym=-6,symsize=0.7 xyouts,55,430,'Original Time Series',font=1,charsize=1.3 oplot,[20,50],[410,410],linestyle=1,thick=6 xyouts,55,410,'Regression Fit',font=1,charsize=1.3 ;Plot residual error time series: plot,residual_fit,/xs,charsize=1.1,xcharsize=1.1,ycharsize=1.1,$ xthick=3,ythick=3,thick=3,$ xtitle='Year',ytitle='Model Residual (Dobson Units)',$ xticks=14,xtickname=xstf2,/ys,yrange=[-30,30] oplot,[0,N-1],[0,0],thick=3 !p.multi=[0,2,2] ;Plot seasonal cycle in Dobson units: t=alpha plot,t,/xs,charsize=0.8,xcharsize=1.3,ycharsize=1.3,$ xthick=3,ythick=3,thick=3,$ xtitle='Month',ytitle='Seasonal Cycle (Dobson Units)',$ xticks=11,xtickname=xstf,/ys,yrange=[280,420] dt=2*sqrt(var_alpha) errplot,t-dt,t+dt ;Plot trend in percent per decade: t=beta*12000/alpha plot,t,/xs,charsize=0.8,xcharsize=1.3,ycharsize=1.3,$ xthick=3,ythick=3,thick=3,$ xtitle='Month',ytitle='Trend (% Per Decade)',$ xticks=11,xtickname=xstf,/ys,yrange=[-10,2] dt=2*sqrt(var_beta)*12000/alpha errplot,t-dt,t+dt oplot,[0,11],[0,0],thick=3 ;Plot QBO coefficient in Dobson units per 10 m/s change in QBO wind: t=gamma*10 plot,t,/xs,charsize=0.8,xcharsize=1.3,ycharsize=1.3,$ xthick=3,ythick=3,thick=3,$ xtitle='Month',ytitle='QBO Coeff (DU Per 10m/s)',$ xticks=11,xtickname=xstf,/ys,yrange=[-6,2] dt=2*sqrt(var_gamma)*10 errplot,t-dt,t+dt oplot,[0,11],[0,0],thick=3 ;Plot SOLAR coefficient in Dobson units per 100 units change in F10.7: t=delta*100 plot,t,/xs,charsize=0.8,xcharsize=1.3,ycharsize=1.3,$ xthick=3,ythick=3,thick=3,$ xtitle='Month',ytitle='SOLAR Coeff (DU Per 100F10.7)',$ xticks=11,xtickname=xstf,/ys,yrange=[-10,15] dt=2*sqrt(var_delta)*100 errplot,t-dt,t+dt oplot,[0,11],[0,0],thick=3 ;Plot seasonal cycle fit in Dobson units: t=alpha_fit plot,t,/xs,charsize=0.8,xcharsize=1.3,ycharsize=1.3,$ xthick=3,ythick=3,thick=3,$ xtitle='Year',ytitle='Seasonal Fit (DU)',$ xticks=14,xtickname=xstf2,/ys,yrange=[280,420] ;Plot trend fit in Dobson units: t=beta_fit plot,t,/xs,charsize=0.8,xcharsize=1.3,ycharsize=1.3,$ xthick=3,ythick=3,thick=3,$ xtitle='Year',ytitle='Trend Fit (DU)',$ xticks=14,xtickname=xstf2,/ys,yrange=[-40,10] oplot,[0,N-1],[0,0],thick=3 ;Plot QBO fit in Dobson units: t=gamma_fit plot,t,/xs,charsize=0.8,xcharsize=1.3,ycharsize=1.3,$ xthick=3,ythick=3,thick=3,$ xtitle='Year',ytitle='QBO Fit (DU)',$ xticks=14,xtickname=xstf2,/ys,yrange=[-10,10] oplot,[0,N-1],[0,0],thick=3 ;Plot SOLAR fit in Dobson units: t=delta_fit plot,t,/xs,charsize=0.8,xcharsize=1.3,ycharsize=1.3,$ xthick=3,ythick=3,thick=3,$ xtitle='Year',ytitle='Solar Fit (DU)',$ xticks=14,xtickname=xstf2,/ys,yrange=[-10,10] oplot,[0,N-1],[0,0],thick=3 device,/close end