;Create the 3 figures for the paper ;1. Temperature ;2. Winds ;3. Relationship FUNCTION weibull_dist,x,param,_EXTRA=ex IF n_elements(param) EQ 0 THEN BEGIN beta = 1.0 gamma = 0.0 eta = 1.0 ENDIF ELSE BEGIN beta = param[0] ;1.0 gamma = param[1] ;0.0 eta = param[2] ;1.0 ENDELSE ;Pavia and O'brien 1986 JClim arg = (x-float(gamma))/eta R = exp(-arg^beta) lambda = (float(beta)/eta)*arg^(beta-1.) return,lambda*R END FUNCTION weibull_fit,data x = transpose(findgen(21)) hdata = transpose(histogram(data,min=0,max=20,binsize=1)/float(n_elements(data))) param = {name: ['beta','gamma','eta'], value: [1.0,0.0,1.0], max: [10.0,0.0,5.0], min : [0.0,0.0,0.0], knob: [0.0,0.0,0.0], changeable: [1,0,1]} mcmc,x,hdata,param,/fast,outputll=outputll,outputvalue=outputvalue,outputy=y,model='weibull_dist',/quiet bestparam = [mean(outputvalue[0,*]),mean(outputvalue[2,*])] range = [stddev(outputvalue[0,*]),stddev(outputvalue[2,*])] return,[bestparam,range] END PRO superiorplots ;Read temp and wind data from superior.pro and superiornarr.pro restore,mydocs()+'superior/Uopt2.sav' restore,mydocs()+'superior/qscat_annvars.sav' restore,mydocs()+'superior/buoy_annvars.sav' restore,mydocs()+'superior/narr_annvars.sav' restore,mydocs()+'superior/narr_annvars_850.sav' ;REVISION replace 250 with 850 sum_current = [0.0006,0.0004,-0.0013,0.0007,-0.0001,0.0012,0.0037,-0.0015,0.0008,0.001,-0.0018,-0.0011,0.0022,0.0001,0,0.0003,0.0007,-0.0018,-0.0016,0.0015,0.0028,-0.0004,0.0012,0.0024,0.0017,0.0063,0.007,0.0069]+0.0237 sum_current = [0.0366,0.0370,0.0297,0.0355,0.0367,0.0375,0.0398,0.0340,0.0381,0.0363,0.0297,0.0308,0.0381,0.0363,0.0337,0.0350,0.0342,0.0321,0.0290,0.0373,0.0400,0.0345,0.0365,0.0419,0.0384,0.0426,0.0486,0.0453] sum_current*=1000.0 sum_mld = [-15.0935,-15.8429,-13.9336,-15.5239,-13.1158,-14.4860,-16.3001,-16.0349,-14.6434,-14.7419,-12.7891,-14.1636,-15.3152,-16.5052,-14.7215,-14.1968,-14.6595,-13.4612,-13.6099,-14.3437,-14.3160,-14.6325,-14.5022,-14.0894,-13.5045,-14.8319,-14.3576,-15.1848] sum_mld*=(-1.0) sum_strat = [157.5428,160.0544,168.4317,155.2975,149.6383,154.0055,151.0985,145.9324,140.3108,163.0140,141.9425,140.6872,147.7279,173.2544,158.9280,168.9368,169.5509,153.1214,150.2041,155.2105,162.4019,166.2102] yrs_n = numgen(1979,2008) yrs_b = numgen(1980,2008) yrs_q = numgen(1999,2008) yrs_s = numgen(1979,2006) yrs_st = numgen(1985,2006) ; yrstart = 1981 yrstart = 1985 nyrs = 2009-yrstart bb = yrstart-1980 nn = yrstart-1979 qsstart = 2000 qsnyrs = 2009-qsstart qq = qsstart-1999 bb_qq = qsstart-1980 nn_qq = qsstart-1979 ;temp trends trendline,yrs_b[bb:*],airt_sum_av[bb:*],xplot=ybx,yplot=at_trend,fit=at_fit,r=at_r,sigma=at_sig trendline,yrs_b[bb:*],laket_sum_av[bb:*],xplot=ybx,yplot=LT_trend,fit=LT_fit,r=LT_r,sigma=LT_sig trendline,yrs_n[nn:*],t2s_lnd[nn:*]-273.15,xplot=ynx,yplot=t2_trend,fit=t2_fit,r=t2_r,sigma=t2_sig trendline,yrs_n[nn:*],tsfcs_lnd[nn:*]-273.15,xplot=ynx,yplot=ts_trend,fit=ts_fit,r=ts_r,sigma=ts_sig trendline,yrs_b[bb:*],deltat_sum_av[bb:*],xplot=ybx,yplot=dt_trend,fit=dt_fit,r=dt_r,sigma=dt_sig trendline,yrs_n[nn:*],tsfcs_lnd[nn:*]-t2s_lnd[nn:*],xplot=ynx,yplot=dtl_trend,fit=dtl_fit,r=dtl_r,sigma=dtl_sig ;wind trends trendline,yrs_st,sum_strat,xplot=ystartx,yplot=strat_trend,fit=strat_fit,r=start_r,sigma=strat_sig trendline,yrs_s[nn:*],sum_current[nn:*],xplot=ynxx,yplot=cur_trend,fit=cur_fit,r=cur_r,sigma=cur_sig trendline,yrs_s[nn:*],sum_mld[nn:*],xplot=ynxx,yplot=mld_trend,fit=mld_fit,r=mld_r,sigma=mld_sig trendline,yrs_b[bb:*],wspd_sum_av[bb:*],xplot=ybx,yplot=wlk_trend,fit=wlk_fit,r=wlk_r,sigma=wlk_sig trendline,yrs_n[nn:*],wspd_lnd[nn:*],xplot=ynx,yplot=wlnd_trend,fit=wlnd_fit,r=wlnd_r,sigma=wlnd_sig trendline,yrs_n[nn:*],wspd_lk[nn:*],xplot=ynx,yplot=wlkn_trend,fit=wlkn_fit,r=wlkn_r,sigma=wlkn_sig trendline,yrs_n[nn:*],wspd_250_av[nn:*],xplot=ynx,yplot=w250_trend,fit=w250_fit,r=w250_r,sigma=w250_sig trendline,yrs_n[0:*],wspd_850_av[0:*],xplot=ynx2,yplot=w850_trend,fit=w850_fit,r=w850_r,sigma=w850_sig ;bootstrapping trends at_bs = make_array(24,24,/float,value=nan()) FOR i = 5,23 DO FOR j = bb,28-i DO at_bs[i,j] = (linfit(yrs_b[j:(j+i)],airt_sum_av[j:(j+i)]))[1] lt_bs = make_array(24,24,/float,value=nan()) FOR i = 5,23 DO FOR j = bb,28-i DO lt_bs[i,j] = (linfit(yrs_b[j:(j+i)],laket_sum_av[j:(j+i)]))[1] ws_bs = make_array(24,24,/float,value=nan()) FOR i = 5,23 DO FOR j = bb,28-i DO ws_bs[i,j] = (linfit(yrs_b[j:(j+i)],wspd_sum_av[j:(j+i)]))[1] wlnd_bs = make_array(24,24,/float,value=nan()) FOR i = 5,23 DO FOR j = nn,29-i DO wlnd_bs[i,j-1] = (linfit(yrs_n[j:(j+i)],wspd_lnd[j:(j+i)]))[1] w850_bs = make_array(24,24,/float,value=nan()) FOR i = 5,23 DO FOR j = nn,29-i DO w850_bs[i,j-1] = (linfit(yrs_n[j:(j+i)],wspd_850_av[j:(j+i)]))[1] dt_bs = make_array(24,24,/float,value=nan()) FOR i = 5,23 DO FOR j = bb,28-i DO dt_bs[i,j] = (linfit(yrs_b[j:(j+i)],deltat_sum_av[j:(j+i)]))[1] ;REVISION WE NEED DT_SIG and WLK_SIG ;wind-temp correlations ; trendline,deltat2[bb:*],modU[bb:*]-mean(modU[bb:*]),xplot=dtplot,yplot=modUplot,fit=dtModU_fit,r=dtmodu_r ; dtModU_fit = [0.294195,0.294195] ;5.5 m/s at 8 m ; dtmodU_fit = [0.394028,0.315223] ;6.0 m/s at 20m ; dtmodU_fit = [5.32423,0.283105] ;6.0 m/s at 20 m dtModU_fit = [4.93545,0.410301]-[mean(wspd_lnd[nn:*]),0.0] dtmodU_sigma = [0.0272827,0.0184465] ; trendline,deltat2[bb:*],wspd_sum_av[bb:*]-mean(wspd_sum_av[bb:*]),xplot=dtplot,yplot=wsdtplot,fit=wsdt_fit,r=wsdt_r ; trendline,deltat2[bb_qq:*],ls_sum_all[qq:*]-mean(ls_sum_all[qq:*]),xplot=dtplot2,yplot=qsdtplot,fit=qsdt_fit,r=qsdt_r trendline,deltat2[bb:*],wspd_sum_av[bb:*]-wspd_lnd[nn:*],xplot=dtplot,yplot=wsdtplot,fit=wsdt_fit,r=wsdt_r,sigma=wsdt_sig,pval=wsdt_p trendline,deltat2[bb_qq:*],ls_sum_all[qq:*]-wspd_lnd[nn_qq:*],xplot=dtplot2,yplot=qsdtplot,fit=qsdt_fit,r=qsdt_r,sigma=qsdt_sig,pval=qsdt_p dtplot = (-1.0)*findgen(31)/10.0 qsdtplot = qsdt_fit[0] + qsdt_fit[1]*dtplot wsdtplot = wsdt_fit[0] + wsdt_fit[1]*dtplot modUplot = dtmodU_fit[0] + dtmodU_fit[1]*dtplot modu = interpol(Uopt2[*,49],dtrange,deltat2) ;Variables of interest ;Quickscat winds ls_sum_all ;Buoy data laket_sum_av, airt_sum_Av, deltat_sum_Av, deltat2, wspd_sum_av ;NARR winds wspd_lnd,wdir_lnd,t2s_lnd,tsfcs_lnd,wspd_250_av,wdir_250_av ;Model winds modU,Uopt,zrange,dtrange stop ; lowess,x,y,smooth,yout ;Fig. 1a. Lake and air temperature and trendline !p.multi = [0,2,2] plot,[0,0],[0,0],xrange=[yrstart,2008],yrange=[5,20],ytitle='Temperature (degrees C)',/nodata oplot,yrs_b,airt_sum_av,psym=-1 oplot,yrs_b,laket_sum_av,psym=-2 oplot,ybx,at_trend oplot,ybx,LT_trend ;Fig 1b. dT plot,[0,0],[0,0],xrange=[yrstart,2008],yrange=[-3,0],ytitle='Delta T (degrees C)',/nodata oplot,yrs_b,deltat_sum_av,psym=-1 oplot,ybx,dt_trend stop ;Fig 2. Wind speeds - Lake land upper ;REVISION USE 850MB ;trend lines !p.multi = [0,2,3] ;,0,1] plot,[0,0],[0,0],/nodata,xrange=[yrstart,2008],yrange=[3,6],ytitle='Surface Wind Speed (m s-1)' oplot,yrs_b,wspd_sum_av ; oplot,yrs_b,wspd_lk oplot,yrs_q[1:*],ls_sum_all[1:*],color=200,thick=3 oplot,ybx,wlk_trend ; plot,[0,0],[0,0],/nodata,xrange=[yrstart,2008],yrange=[3,6],title='Wind Speed - land',ytitle='m s-1' oplot,yrs_n,wspd_lnd,linestyle=2 oplot,ynx,wlnd_trend ; oplot,yrs_n,wspd_850_av,linestyle=1 ; oplot,ynx,w850_trend plot,[0,0],[0,0],/nodata,xrange=[yrstart,2008],yrange=[2,10],ytitle='850mb Wind Speed (m s-1)' oplot,yrs_n,wspd_850_av oplot,ynx,w850_trend ; plot,[0,0],[0,0],/nodata,xrange=[yrstart,2008],yrange=[20,40],ytitle='250mb Wind Speed (m s-1)' ; oplot,yrs_n,wspd_250_av ; oplot,ynx,w250_trend stop ;Fig. 3 DT to wind speed correlation !p.multi = [0,2,2] plot,[0,0],[0,0],/nodata,xrange=[-2.5,0.0],yrange=[0,1.5],xtitle='Temperature gradient (degrees C)',ytitle='Wind speed gradient (m s-1)' oplot,deltat2[bb:*],wspd_sum_Av[bb:*]-wspd_lnd[nn:*],psym=1 oplot,deltat2[bb_qq:*],ls_sum_all[qq:*]-wspd_lnd[nn_qq:*],psym=2 oplot,dtplot,wsdtplot,linestyle=1 oplot,dtplot,qsdtplot,linestyle=2 ; oplot,dtplot,modUplot ; oplot,deltat2[bb:*],wspd_sum_av[bb:*]-mean(wspd_sum_av[bb:*]),psym=1 ; oplot,deltat2[bb_qq:*],ls_sum_all[qq:*]-mean(ls_sum_all[qq:*]),psym=2 ; xyouts,deltat2[bb:*],wspd_sum_av[bb:*]-mean(Wspd_sum_av[bb:*]),yrs_b[bb:*] stop ;Fig. 4 Model ;Perhaps a separate figure (Fig 4) !p.multi = [0,2,2] plot,[0,0],[0,0],/nodata,xrange=[yrstart,2006],yrange=[10,50],ytitle='Current Speed (mm s-1) Mixed layer depth (m)',ystyle=8 oplot,yrs_s,sum_current oplot,ynxx,cur_trend oplot,yrs_s,sum_mld,linestyle=1 oplot,ynxx,mld_trend axis,/yaxis,/save,yrange=[125,225],ytitle='Stratified length (days)' oplot,yrs_st,sum_strat,linestyle=2 oplot,ystartx,strat_trend stop ;Trends and numbers here print,'Temperature' print,'Air temp (degrees/decade) ',at_fit[1]*10.0,' r2 =',at_r^2,' p =',corr_ttest(corr=at_r,n=nyrs),' sigma =',at_sig[1]*10.0 print,'Lake temp (degrees/decade) ',lt_fit[1]*10.0,' r2 =',lt_r^2,' p =',corr_ttest(corr=lt_r,n=nyrs),' sigma =',LT_sig[1]*10.0 print,'Land temp (degrees/decade) ',t2_fit[1]*10.0,' r2 =',t2_r^2,' p =',corr_ttest(corr=t2_r,n=nyrs), 'sigma =',t2_sig[1]*10.0 print,'Lake deltaT (degrees/decade) ',dt_fit[1]*10.0,' r2 =',dt_r^2,' p =',corr_ttest(corr=dt_r,n=nyrs),' sigma =',dt_sig[1]*10.0 print,'Land deltat (degrees/decade) ',dtl_fit[1]*10.0,' r2 =',dtl_r^2,' p =',corr_ttest(corr=dtl_r,n=nyrs),' sigma =',dtl_sig[1]*10.0 print,'Wind Speed' print,'Lake wind (m/s/decade) ',wlk_fit[1]*10.0,' r2 =',wlk_r^2,' p =',corr_ttest(corr=wlk_r,n=nyrs),' sigma =',wlk_sig[1]*10.0 print,'Land wind (m/s/decade) ',wlnd_fit[1]*10.0,' r2 =',wlnd_r^2,' p =',corr_ttest(corr=wlnd_r,n=nyrs),' sigma =',wlnd_sig[1]*10.0 print,'250mb wind (m/s/decade) ',w250_fit[1]*10.0,' r2 =',w250_r^2,' p =',corr_ttest(corr=w250_r,n=nyrs), 'sigma =',w250_sig[1]*10.0 print,'850mb wind (m/s/decade) ',w850_fit[1]*10.0,' r2 =',w850_r^2,' p =',corr_ttest(corr=w850_r,n=nyrs), 'sigma =',w850_sig[1]*10.0 print,'Trends' print,'Model (m/s/degree) ',dtmodu_fit[1] print,'Buoy (m/s/degree) ',wsdt_fit[1],' r2 =',wsdt_r^2,' p =',corr_ttest(corr=wsdt_r,n=nyrs),' sigma =',wsdt_sig[1] print,'Satellite (m/s/degree) ',qsdt_fit[1],' r2 =',qsdt_r^2,' p =',corr_ttest(corr=qsdt_r,n=qsnyrs),' sigma =',qsdt_sig[1] stop ;Supplemental Figures ;Supp Fig 1. Model profiles of wind speed as f of dT and z !p.multi = [0,2,2] ustar = 0.2 z0 = 6.5e-5 plot,Uopt2[0,*]/ustar,zrange/z0,xrange=[15,30],yrange=[1500,300000],ytitle='z/z0',xtitle='U/u*',thick=2 oplot,Uopt2[4,*]/ustar,zrange/z0,thick=2,color=100 oplot,Uopt2[8,*]/ustar,zrange/z0,thick=2,color=200 xyouts,22,120000,'0' xyouts,22,100000,'-1' xyouts,22,60000,'-2' stop ;Supp Fig 2. Land temperature and deltaT !p.multi = [0,2,2] plot,[0,0],[0,0],xrange=[yrstart,2008],yrange=[14,20],title='Temperature',ytitle='Degrees C',/nodata oplot,yrs_n,t2s_lnd-273.15,psym=-4 oplot,yrs_n,tsfcs_lnd-273.15,psym=-5 oplot,ynx,t2_trend oplot,ynx,ts_trend plot,[0,0],[0,0],xrange=[yrstart,2008],yrange=[0,3],title='Delta_T',ytitle='Degrees C',/nodata oplot,ynx,dtl_trend oplot,yrs_n,tsfcs_lnd-t2s_lnd,psym=-4 stop ;Supp Fig 3. Histograms ;compare 1982-1986 2004-2008 ;compare 1985-1989 2004-2008 ;compare 1985-1994 1999-2008 restore,mydocs()+'superior/buoy_data.sav' ws_1 = reform(wspd_sum[*,5:14],2232.*10.) ws_2 = reform(wspd_sum[*,19:28],2232.*10.) ws_1 = ws_1[where(finite(Ws_1),nws1)] ws_2 = ws_2[where(finite(ws_2),nws2)] wd_1 = reform(wdir_sum[*,5:14],2232.*10.) wd_2 = reform(wdir_sum[*,19:28],2232.*10.) wd_1 = wd_1[where(finite(Wd_1),nwd1)] wd_2 = wd_2[where(finite(wd_2),nwd2)] gs_1 = reform(gst_sum[*,5:14],2232.*10.) gs_2 = reform(gst_sum[*,19:28],2232.*10.) gs_1 = gs_1[where(finite(Gs_1),ngs1)] gs_2 = gs_2[where(finite(gs_2),ngs2)] dt_1 = reform(deltat_sum[*,5:14],2232.*10.) dt_2 = reform(deltat_sum[*,19:28],2232.*10.) dt_1 = dt_1[where(finite(Dt_1),ndt1)] dt_2 = dt_2[where(finite(dt_2),ndt2)] !p.multi = [0,2,2] plot,[0,0],[0,0],/nodata,yrange=[0,0.4],xrange=[0,15],xtitle='Wind speed (m s-1)',ytitle='Relative frequency' oplot,histogram(ws_1,min=0,max=30)/float(nws1),linestyle=2 oplot,histogram(ws_2,min=0,max=30)/float(nws2) plot,[0,0],[0,0],/nodata,yrange=[0,0.2],xrange=[0,360],xtitle='Wind Direction (degree)',ytitle='Relative frequency' oplot,findgen(18)*20,histogram(wd_1,min=0,max=360,bin=20)/float(nwd1),linestyle=2 oplot,findgen(18)*20,histogram(wd_2,min=0,max=360,bin=20)/float(nwd2) plot,[0,0],[0,0],/nodata,yrange=[0,0.4],xrange=[0,15],xtitle='Wind Gust (m s-1)',ytitle='Relative frequency' oplot,histogram(gs_1,min=0,max=30)/float(ngs1),linestyle=2 oplot,histogram(gs_2,min=0,max=30)/float(ngs2) plot,[0,0],[0,0],/nodata,yrange=[0,0.3],xrange=[-10,10],xtitle='Delta T (degree C)',ytitle='Relative frequency' oplot,findgen(20)-10,histogram(dt_1,min=-10,max=10)/float(ngs1),linestyle=2 oplot,findgen(20)-10,histogram(dt_2,min=-10,max=10)/float(ngs2) ;Weibull games ;first for ws_1 and ws_2 ; print,'Shape/Mean 1985-1994: ',weibull_fit(ws_1) ; print,'Shape/Mean 1999-2008: ',weibull_fit(ws_2) ;for each year wspd_sum plot A and C stop ;Supp Fig 4. CyCLeS Model ;if any stop ;significance ;corr_ttest(corr=r,n=nvals END