* writes out 95th percentile of 30yr anomalies for present day (1979-2005) relative to historic (1850-1879)
* mask uses Student t-test distribution
 'reinit'

 'xdfopen /large/obs2/cmip5/mo/bcc/hist/pr_Amon_bcc-csm1-1_historical_r1i1p1_185001-201212.ctl'

 'set lon 0 360'
 'set lat -60 60'

* define clim (1850-79)
 'set time jan1850 dec1879'
 'define prc=ave(pr,t+0,time=dec1879,1yr)'
 'modify prc seasonal'

* define anom for full time series (independent mean)
 'set time jan1850 dec2005'
 'define pa=(pr-prc)*86400'

* define std for independent mean. Place avg in each of the first 12 time steps.
 'set t 1 13'
 'define pastd=ave(pa*pa,time=jan1850,time=dec1879)'
* define Student t-test distribution parameters. \
* Right tail percentage points for n=30yrs and p (t >= t_alpha) = 0.025. 
* factor of 2.0 might mean 2 sample test and 95th percentile?
* see Wilke for table parameters
 'define np=30'
 'define tfac=2.042'
 'define sigma=tfac*sqrt(2.0*pastd/np)'

* define 30yr avg for end of 20th cent. modify for seasonal avg
 'set time jan1976 dec2005'
 'define p30a=ave(pa,t+0,time=dec2005,1yr)'
 'modify p30a seasonal'

* Color scale (blue(50) to red(68))
 'set rgb 43 0 0 0'
 'set rgb 49 0 13 129'
 'set rgb 50 0 30 164'
 'set rgb 52 21 79 161'
 'set rgb 54 24 150 153'
 'set rgb 56 67 167 41'
 'set rgb 57 156 210 41'
 'set rgb 58 187 240 0'
 'set rgb 59 219 255 68'
 'set rgb 60 255 242 187'
 'set rgb 61 255 242 106'
 'set rgb 62 255 182 0'
 'set rgb 64 247 107 0'
 'set rgb 66 232 78 0'
 'set rgb 68 204 61 0'
 'set rgb 69 179 0 0'
 'set rgb 70 109 5 0'

 'set t 1'
 'set vpage 0 8.5 5.5 11'
 'set grads off'
 'set gxout grfill'
 'set clevs -3 -2.5 -2 -1.5 -1 -0.5 -0.25 0 0.25 0.5 1 1.5 2 2.5 3'
 'set ccols 70 69 68 66 64 62 61 60 59 58 57 56 54 52 50 49'
 'define p30DJFa=ave(p30a,t=12,t=14)'
 'd maskout(p30DJFa,abs(p30DJFa)-sigma)'
 'draw title BCC DJF Pra(76-05) wrt 1850-79 (95%)'
 'run /cbar'


 'set t 7'
 'set vpage 0 8.5 0 5.5'
 'set grads off'
 'set gxout grfill'
 'set clevs -3 -2.5 -2 -1.5 -1 -0.5 -0.25 0 0.25 0.5 1 1.5 2 2.5 3'
 'set ccols 70 69 68 66 64 62 61 60 59 58 57 56 54 52 50 49'
 'define p30JJAa=ave(p30a,t=6,t=8)'
 'd maskout(p30JJAa,abs(p30JJAa)-sigma)'
 'draw title BCC JJA Pra(76-05) wrt 1850-79 (95%)'
 'run /cbar'

 'print'
 '!mv grads.eps bcc_pra30yr_djfjja_79-05_95sig.eps'

return
 'd ave(sigma,t-1,t+1)'
 'draw title CanESM2 JJA Preca sigma(t-test) (1850-1879)'
 'run /cbar'
 'print'

 '!mv grads.eps CanESM2_pra_jja_sigma_1850-79.eps'

**********