2009-12-06, 15:49
#1
Jag tycker vi p FB ska underska den verkomna programmeringskod frn CRU. CRU r ttt sammankopplade med FN:s klimatpanel IPCC. CRU:s privata email och programkod ligger nu fritt p Internet dr autenticiteten i deras forskning verkligen kan ifrgasttas.
P olika bloggar kommenteras fljande kod frn CRU. Ls frst
http://wattsupwiththat.com/2009/12/0...-smoking-code/
http://wattsupwiththat.com/2009/12/0...g-code-part-2/
och ge sen er syn p koden. Kan vi ana fult spel frn CRU?
Koden:
001 1. ;
002 2. ; PLOTS 'ALL' REGION MXD timeseries from age banded and from hugershoff
003 3. ; standardised datasets.
004 4. ; Reads Harry's regional timeseries and outputs the 1600-1992 portion
005 5. ; with missing values set appropriately. Uses mxd, and just the
006 6. ; "all band" timeseries
007 7. ;****** APPLIES A VERY ARTIFICIAL CORRECTION FOR DECLINE*********
008 8. ;
009 9. yrloc=[1400,findgen(19)*5.+1904]
010 10. valadj=[0.,0.,0.,0.,0.,-0.1,-0.25,-0.3,0.,-0.1,0.3,0.8,1.2,1.7,2.5,2.6,2.6,$
011 11. 2.6,2.6,2.6]*0.75 ; fudge factor
012 12. if n_elements(yrloc) ne n_elements(valadj) then message,'Oooops!'
013 13. ;
014 14. loadct,39
015 15. def_1color,20,color='red'
016 16. plot,[0,1]
017 17. multi_plot,nrow=4,layout='large'
018 18. if !d.name eq 'X' then begin
019 19. window, ysize=800
020 20. !p.font=-1
021 21. endif else begin
022 22. !p.font=0
023 23. device,/helvetica,/bold,font_size=18
024 24. endelse
025 25. ;
026 26. ; Get regional tree lists and rbar
027 27. ;
028 28. restore,filename='reglists.idlsave'
029 29. harryfn=['nwcan','wnam','cecan','nweur','sweur','nsib','csi b','tib',$
030 30. 'esib','allsites']
031 31. ;
032 32. rawdat=fltarr(4,2000)
033 33. for i = nreg-1 , nreg-1 do begin
034 34. fn='mxd.'+harryfn(i)+'.pa.mean.dat'
035 35. print,fn
036 36. openr,1,fn
037 37. readf,1,rawdat
038 38. close,1
039 39. ;
040 40. densadj=reform(rawdat(2:3,*))
041 41. ml=where(densadj eq -99.999,nmiss)
042 42. densadj(ml)=!values.f_nan
043 43. ;
044 44. x=reform(rawdat(0,*))
045 45. kl=where((x ge 1400) and (x le 1992))
046 46. x=x(kl)
047 47. densall=densadj(1,kl) ; all bands
048 48. densadj=densadj(0,kl) ; 2-6 bands
049 49. ;
050 50. ; Now normalise w.r.t. 1881-1960
051 51. ;
052 52. mknormal,densadj,x,refperiod=[1881,1960],refmean=refmean,refsd=refsd
053 53. mknormal,densall,x,refperiod=[1881,1960],refmean=refmean,refsd=refsd
054 54. ;
055 55. ; APPLY ARTIFICIAL CORRECTION
056 56. ;
057 57. yearlyadj=interpol(valadj,yrloc,x)
058 58. densall=densall+yearlyadj
059 59. ;
060 60. ; Now plot them
061 61. ;
062 62. filter_cru,20,tsin=densall,tslow=tslow,/nan
063 63. cpl_barts,x,densall,title='Age-banded MXD from all sites',$
064 64. xrange=[1399.5,1994.5],xtitle='Year',/xstyle,$
065 65. zeroline=tslow,yrange=[-7,3]
066 66. oplot,x,tslow,thick=3
067 67. oplot,!x.crange,[0.,0.],linestyle=1
068 68. ;
069 69. endfor
070 70. ;
071 71. ; Restore the Hugershoff NHD1 (see Nature paper 2)
072 72. ;
073 73. xband=x
074 74. restore,filename='../tree5/densadj_MEAN.idlsave'
075 75. ; gets: x,densadj,n,neff
076 76. ;
077 77. ; Extract the post 1600 part
078 78. ;
079 79. kl=where(x ge 1400)
080 80. x=x(kl)
081 81. densadj=densadj(kl)
082 82. ;
083 83. ; APPLY ARTIFICIAL CORRECTION
084 84. ;
085 85. yearlyadj=interpol(valadj,yrloc,x)
086 86. densadj=densadj+yearlyadj
087 87. ;
088 88. ; Now plot it too
089 89. ;
090 90. filter_cru,20,tsin=densadj,tslow=tshug,/nan
091 91. cpl_barts,x,densadj,title='Hugershoff-standardised MXD from all sites',$
092 92. xrange=[1399.5,1994.5],xtitle='Year',/xstyle,$
093 93. zeroline=tshug,yrange=[-7,3],bar_color=20
094 94. oplot,x,tshug,thick=3,color=20
095 95. oplot,!x.crange,[0.,0.],linestyle=1
096 96. ;
097 97. ; Now overplot their bidecadal components
098 98. ;
099 99. plot,xband,tslow,$
100 100. xrange=[1399.5,1994.5],xtitle='Year',/xstyle,$
101 101. yrange=[-6,2],thick=3,title='Low-pass (20-yr) filtered comparison'
102 102. oplot,x,tshug,thick=3,color=20
103 103. oplot,!x.crange,[0.,0.],linestyle=1
104 104. ;
105 105. ; Now overplot their 50-yr components
106 106. ;
107 107. filter_cru,50,tsin=densadj,tslow=tshug,/nan
108 108. filter_cru,50,tsin=densall,tslow=tslow,/nan
109 109. plot,xband,tslow,$
110 110. xrange=[1399.5,1994.5],xtitle='Year',/xstyle,$
111 111. yrange=[-6,2],thick=3,title='Low-pass (50-yr) filtered comparison'
112 112. oplot,x,tshug,thick=3,color=20
113 113. oplot,!x.crange,[0.,0.],linestyle=1
114 114. ;
115 115. ; Now compute the full, high and low pass correlations between the two
116 116. ; series
117 117. ;
118 118. perst=1400.
119 119. peren=1992.
120 120. ;
121 121. openw,1,'corr_age2hug.out'
122 122. thalf=[10.,30.,50.,100.]
123 123. ntry=n_elements(thalf)
124 124. printf,1,'Correlations between timeseries'
125 125. printf,1,'Age-banded vs. Hugershoff-standardised'
126 126. printf,1,' Region Full <10 >10 >30 >50 >100'
127 127. ;
128 128. kla=where((xband ge perst) and (xband le peren))
129 129. klh=where((x ge perst) and (x le peren))
130 130. ts1=densadj(klh)
131 131. ts2=densall(kla)
132 132. ;
133 133. r1=correlate(ts1,ts2)
134 134. rall=fltarr(ntry)
135 135. for i = 0 , ntry-1 do begin
136 136. filter_cru,thalf(i),tsin=ts1,tslow=tslow1,tshigh=t shi1,/nan
137 137. filter_cru,thalf(i),tsin=ts2,tslow=tslow2,tshigh=t shi2,/nan
138 138. if i eq 0 then r2=correlate(tshi1,tshi2)
139 139. rall(i)=correlate(tslow1,tslow2)
140 140. endfor
141 141. ;
142 142. printf,1,'ALL SITES',r1,r2,rall,$
143 143. format='(A11,2X,6F6.2)'
144 144. ;
145 145. printf,1,' '
146 146. printf,1,'Correlations carried out over the period ',perst,peren
147 147. ;
148 148. close,1
149 149. ;
150 150. end
P olika bloggar kommenteras fljande kod frn CRU. Ls frst
http://wattsupwiththat.com/2009/12/0...-smoking-code/
http://wattsupwiththat.com/2009/12/0...g-code-part-2/
och ge sen er syn p koden. Kan vi ana fult spel frn CRU?
Koden:
001 1. ;
002 2. ; PLOTS 'ALL' REGION MXD timeseries from age banded and from hugershoff
003 3. ; standardised datasets.
004 4. ; Reads Harry's regional timeseries and outputs the 1600-1992 portion
005 5. ; with missing values set appropriately. Uses mxd, and just the
006 6. ; "all band" timeseries
007 7. ;****** APPLIES A VERY ARTIFICIAL CORRECTION FOR DECLINE*********
008 8. ;
009 9. yrloc=[1400,findgen(19)*5.+1904]
010 10. valadj=[0.,0.,0.,0.,0.,-0.1,-0.25,-0.3,0.,-0.1,0.3,0.8,1.2,1.7,2.5,2.6,2.6,$
011 11. 2.6,2.6,2.6]*0.75 ; fudge factor
012 12. if n_elements(yrloc) ne n_elements(valadj) then message,'Oooops!'
013 13. ;
014 14. loadct,39
015 15. def_1color,20,color='red'
016 16. plot,[0,1]
017 17. multi_plot,nrow=4,layout='large'
018 18. if !d.name eq 'X' then begin
019 19. window, ysize=800
020 20. !p.font=-1
021 21. endif else begin
022 22. !p.font=0
023 23. device,/helvetica,/bold,font_size=18
024 24. endelse
025 25. ;
026 26. ; Get regional tree lists and rbar
027 27. ;
028 28. restore,filename='reglists.idlsave'
029 29. harryfn=['nwcan','wnam','cecan','nweur','sweur','nsib','csi b','tib',$
030 30. 'esib','allsites']
031 31. ;
032 32. rawdat=fltarr(4,2000)
033 33. for i = nreg-1 , nreg-1 do begin
034 34. fn='mxd.'+harryfn(i)+'.pa.mean.dat'
035 35. print,fn
036 36. openr,1,fn
037 37. readf,1,rawdat
038 38. close,1
039 39. ;
040 40. densadj=reform(rawdat(2:3,*))
041 41. ml=where(densadj eq -99.999,nmiss)
042 42. densadj(ml)=!values.f_nan
043 43. ;
044 44. x=reform(rawdat(0,*))
045 45. kl=where((x ge 1400) and (x le 1992))
046 46. x=x(kl)
047 47. densall=densadj(1,kl) ; all bands
048 48. densadj=densadj(0,kl) ; 2-6 bands
049 49. ;
050 50. ; Now normalise w.r.t. 1881-1960
051 51. ;
052 52. mknormal,densadj,x,refperiod=[1881,1960],refmean=refmean,refsd=refsd
053 53. mknormal,densall,x,refperiod=[1881,1960],refmean=refmean,refsd=refsd
054 54. ;
055 55. ; APPLY ARTIFICIAL CORRECTION
056 56. ;
057 57. yearlyadj=interpol(valadj,yrloc,x)
058 58. densall=densall+yearlyadj
059 59. ;
060 60. ; Now plot them
061 61. ;
062 62. filter_cru,20,tsin=densall,tslow=tslow,/nan
063 63. cpl_barts,x,densall,title='Age-banded MXD from all sites',$
064 64. xrange=[1399.5,1994.5],xtitle='Year',/xstyle,$
065 65. zeroline=tslow,yrange=[-7,3]
066 66. oplot,x,tslow,thick=3
067 67. oplot,!x.crange,[0.,0.],linestyle=1
068 68. ;
069 69. endfor
070 70. ;
071 71. ; Restore the Hugershoff NHD1 (see Nature paper 2)
072 72. ;
073 73. xband=x
074 74. restore,filename='../tree5/densadj_MEAN.idlsave'
075 75. ; gets: x,densadj,n,neff
076 76. ;
077 77. ; Extract the post 1600 part
078 78. ;
079 79. kl=where(x ge 1400)
080 80. x=x(kl)
081 81. densadj=densadj(kl)
082 82. ;
083 83. ; APPLY ARTIFICIAL CORRECTION
084 84. ;
085 85. yearlyadj=interpol(valadj,yrloc,x)
086 86. densadj=densadj+yearlyadj
087 87. ;
088 88. ; Now plot it too
089 89. ;
090 90. filter_cru,20,tsin=densadj,tslow=tshug,/nan
091 91. cpl_barts,x,densadj,title='Hugershoff-standardised MXD from all sites',$
092 92. xrange=[1399.5,1994.5],xtitle='Year',/xstyle,$
093 93. zeroline=tshug,yrange=[-7,3],bar_color=20
094 94. oplot,x,tshug,thick=3,color=20
095 95. oplot,!x.crange,[0.,0.],linestyle=1
096 96. ;
097 97. ; Now overplot their bidecadal components
098 98. ;
099 99. plot,xband,tslow,$
100 100. xrange=[1399.5,1994.5],xtitle='Year',/xstyle,$
101 101. yrange=[-6,2],thick=3,title='Low-pass (20-yr) filtered comparison'
102 102. oplot,x,tshug,thick=3,color=20
103 103. oplot,!x.crange,[0.,0.],linestyle=1
104 104. ;
105 105. ; Now overplot their 50-yr components
106 106. ;
107 107. filter_cru,50,tsin=densadj,tslow=tshug,/nan
108 108. filter_cru,50,tsin=densall,tslow=tslow,/nan
109 109. plot,xband,tslow,$
110 110. xrange=[1399.5,1994.5],xtitle='Year',/xstyle,$
111 111. yrange=[-6,2],thick=3,title='Low-pass (50-yr) filtered comparison'
112 112. oplot,x,tshug,thick=3,color=20
113 113. oplot,!x.crange,[0.,0.],linestyle=1
114 114. ;
115 115. ; Now compute the full, high and low pass correlations between the two
116 116. ; series
117 117. ;
118 118. perst=1400.
119 119. peren=1992.
120 120. ;
121 121. openw,1,'corr_age2hug.out'
122 122. thalf=[10.,30.,50.,100.]
123 123. ntry=n_elements(thalf)
124 124. printf,1,'Correlations between timeseries'
125 125. printf,1,'Age-banded vs. Hugershoff-standardised'
126 126. printf,1,' Region Full <10 >10 >30 >50 >100'
127 127. ;
128 128. kla=where((xband ge perst) and (xband le peren))
129 129. klh=where((x ge perst) and (x le peren))
130 130. ts1=densadj(klh)
131 131. ts2=densall(kla)
132 132. ;
133 133. r1=correlate(ts1,ts2)
134 134. rall=fltarr(ntry)
135 135. for i = 0 , ntry-1 do begin
136 136. filter_cru,thalf(i),tsin=ts1,tslow=tslow1,tshigh=t shi1,/nan
137 137. filter_cru,thalf(i),tsin=ts2,tslow=tslow2,tshigh=t shi2,/nan
138 138. if i eq 0 then r2=correlate(tshi1,tshi2)
139 139. rall(i)=correlate(tslow1,tslow2)
140 140. endfor
141 141. ;
142 142. printf,1,'ALL SITES',r1,r2,rall,$
143 143. format='(A11,2X,6F6.2)'
144 144. ;
145 145. printf,1,' '
146 146. printf,1,'Correlations carried out over the period ',perst,peren
147 147. ;
148 148. close,1
149 149. ;
150 150. end
