;This script gives a generic example of how to use ncl to compute simple statistics on RegCM3 output. In this example
;the input array is dimensioned time x lat x lon, with statistics calculated in lat x lon space. The array that
;is used consisists of 12 monthly averages from RegCM3, on the lat x lon grid,
;but array sizes are dynamically assigned so ARRAYS OF ANY DIMENSIONALITY CAN BE USED (for instance, daily
;or 6-hourly values on the lat x lon grid,a series of july averages from each year of a simulation, etc.
;Some of the statistics included below will be more or less usefull depending on the array supplied).
;
;Note that this script produces a netcdf output file containing the results of the statistical analysis. These
;new variables can be plotted using the plotting routines on the website.
;
;********************************************
load "/usr/local/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "/usr/local/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "/usr/local/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   
load "/usr/local/lib/ncarg/nclscripts/csm/contributed.ncl"   
;********************************************

begin

;user-specied buffer on each side of model domain
;four variables are used to provide maximum flexibility, though all values may be the same
;depending on user intention.
east_buffer = 8
west_buffer = 8
north_buffer = 8
south_buffer = 8

diri = "/vol1/msnyder/RegCM3/oct6/RegCM/Cal_Ncep_opt/"
fili = "ATM1960AVG.nc"

diri1 = "/vol1/msnyder/RegCM3/oct6/RegCM/Cal_Ncep_opt/"
fili1 = "HEAD_OUT.nc"

dirout = "/vol1/suresh/"
filout = "ATM1960AVG_stats.nc"

f = addfile(diri+fili,"r")

f1 = addfile(diri1+fili1,"r")

fout = addfile(dirout+filout,"c")

; Read in variables
ts_in     = f->TK(:,0,:,:)
lat2d_in  = f1->XLAT(0,:,:)
lon2d_in  = f1->XLON(0,:,:)

; Get dimensions from input variable
dimvar_in = dimsizes(ts_in)
ntim_in   = dimvar_in(0)
jlat_in   = dimvar_in(1)
ilon_in   = dimvar_in(2)

;trim buffer (buffer size is set by the user above)
ts = ts_in(:,south_buffer:jlat_in-north_buffer-1,west_buffer:ilon_in-north_buffer-1)
lat2d = lat2d_in(south_buffer:jlat_in-north_buffer-1,west_buffer:ilon_in-north_buffer-1)
lon2d = lon2d_in(south_buffer:jlat_in-north_buffer-1,west_buffer:ilon_in-north_buffer-1)

;get dimensions of trimmed variable
dimvar = dimsizes(ts)
ntim = dimvar(0)
jlat = dimvar(1)
ilon = dimvar(2)

;use ncl functions to compute simple statistics on the data
;whole domain for all time planes
d_avg = avg(ts)						;domain average over all time planes
d_max = max(ts)						;domain max over all time planes
d_min = min(ts)						;domain min over all time planes
d_stddev = stddev(ts)					;domain standard deviation over all time planes
d_var = variance(ts)					;domain variance over all time planes

;whole domain at each time plane
;allocate memory for variables
t_avg = new((/ntim/),typeof(ts))
t_max = new((/ntim/),typeof(ts))
t_min = new((/ntim/),typeof(ts))
t_stddev = new((/ntim/),typeof(ts))
t_var = new((/ntim/),typeof(ts))

do i=0,ntim-1
	t_avg(i) = avg(ts(i,:,:))				;domain average for each time plane
	t_max(i) = max(ts(i,:,:))				;domain max for each time plane
	t_min(i) = min(ts(i,:,:))				;domain min for each time plane
	t_stddev(i) = stddev(ts(i,:,:))			;domain standard deviation for each time plane
	t_var(i) = variance(ts(i,:,:))			;domain variance for each time plane
end do

;lat x lon for all time planes
time_avg = dim_avg(ts(lat|:,lon|:,time|:))		;lat x lon average at each grid point
time_med = dim_median(ts(lat|:,lon|:,time|:))		;lat x lon median at each grid point
time_min = dim_min(ts(lat|:,lon|:,time|:))		;lat x lon min at each grid point 
time_max = dim_max(ts(lat|:,lon|:,time|:))		;lat x lon max at each grid point
time_stddev = dim_stddev(ts(lat|:,lon|:,time|:))	;lat x lon standard deviation at each grid point
time_var = dim_variance(ts(lat|:,lon|:,time|:))		;lat x lon variance at each grid point

;assign some meta-data so it's easier to tell what the variables are in the netcdf file
d_avg@long_name = "domain average over all time planes"
d_max@long_name = "domain max over all time planes"
d_min@long_name	= "domain min over all time planes"
d_stddev@long_name = "domain standard deviation over all time planes"
d_var@long_name = "domain variance over all time planes"

t_avg@long_name = "domain average for each time plane"
t_max@long_name = "domain max for each time plane"
t_min@long_name = "domain min for each time plane"
t_stddev@long_name = "domain standard deviation for each time plane"
t_var@long_name = "domain variance for each time plane"

time_avg@long_name = "lat x lon average at each grid point"
time_med@long_name = "lat x lon median at each grid point"
time_min@long_name = "lat x lon min at each grid point"
time_max@long_name = "lat x lon max at each grid point"
time_stddev@long_name = "lat x lon standard deviation at each grid point"
time_var@long_name = "lat x lon variance at each grid point"

;assign coordinate variables
;scalar
d_avg!0 = "avg"
d_max!0 = "avg"
d_min!0 = "avg"
d_stddev!0 = "avg"
d_var!0 = "avg"

;time
t_avg!0 = "time"
t_max!0 = "time"
t_min!0 = "time"
t_stddev!0 = "time"
t_var!0 = "time"

;lat x lon
time_avg!0 = "lat"
time_avg!1 = "lon"

time_med!0 = "lat"
time_med!1 = "lon"

time_min!0 = "lat"
time_min!1 = "lon"

time_max!0 = "lat"
time_max!1 = "lon"

time_stddev!0 = "lat"
time_stddev!1 = "lon"

time_var!0 = "lat"
time_var!1 = "lon"

;output variables to file
fout->d_avg = d_avg
fout->d_max = d_max
fout->d_min = d_min
fout->d_stddev = d_stddev
fout->d_var = d_var

fout->t_avg = t_avg
fout->t_max = t_max
fout->t_min = t_min
fout->t_stddev = t_stddev
fout->t_var = t_var

fout->time_avg = time_avg
fout->time_med = time_med
fout->time_max = time_max
fout->time_min = time_min
fout->time_stddev = time_stddev
fout->time_var = time_var

fout->xlat = lat2d
fout->xlon = lon2d

end