;NCL test script to estimate statistical significance and optionally, t-values
;of wet season 2 meter temperature for 1985-2010 for modern and potential natural vegetation land covers 
;rewritten Fall 2013

;Load necessary ncl code

load "/opt/local/ncl/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "/opt/local/ncl/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "/opt/local/ncl/lib/ncarg/nclscripts/csm/contributed.ncl"
load "/opt/local/ncl/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin

;User defined variables

	sFileName	= "Ttest_t2mwet85_10"  ;set the output/plot file name

;	sFileName2	= "Ttest_t2mwet85_10sig0.05"

;set the directories when the RegCM data are stored
	
	sDataDir1	= "/home/vol5/hstauffe/SEARegCM4.1.1/postprocMOD/" 
				;MOD postproc directory

	sDataDir2	= "/home/vol5/hstauffe/SEAR4.1.1NAT/postprocNAT/"
				;NAT postproc directory

;set the data files; seasonal averages for 25 years, 25 records
	
	sDataFile1	= "SRF.wet8510.nc" ;MOD seasonal average file

	sDataFile2	= "SRF.wet8510N.nc" ;NAT seasonal average file

	;Set the Mercator projection parameters
	;(defined in /.../PreProc/Terrain/domain.param)
	fCenterLon		= 120.	
	fTrueLatL		= 25.
	fTrueLatH		= 50.

;variable reading code

;read in coordinates and temperature variable for specific year
;from MOD file

	;Open SRF MOD files, which hold header info including xlan and xlon
	;same xlat, xlon for all files

	fSRF	= addfile(sDataDir1 + sDataFile1, "r") ; MOD file
	fSRF2	= addfile(sDataDir2 + sDataFile2, "r") ; NAT file

;read in lat/lon for the map's data projection with buffer
lat2d	= fSRF->xlat(8:122,8:166)
lon2d	= fSRF->xlon(8:122,8:166)

;use dimsizes() to set the number of grid points in the 
;horizontal (nlon) and vertical (nlat) direction
iDimSizes1	=dimsizes(lat2d)
nlat		=iDimSizes1(0)
nlon		=iDimSizes1(1)

;read in variable from MOD file

var1	=fSRF->t2m(:,:,8:122,8:166) ;read in variable  with buffer
					;25 records per grid cell

;read in variable from NAT file

var2	=fSRF2->t2m(:,:,8:122,8:166) ;read in variable with buffer
					;25 records per grid cell

;reordering the dimensions temporarily so that they are (lat,lon,time) instead of
;(time, m2, iy, jx)

var1tmp = var1(jx|:,iy|:,m2|:,time|:)

var2tmp = var2(jx|:,iy|:,m2|:,time|:)

;number of values for each variable, xN, yN should be 25

;the means for each grid point for the entire time frame (mean of 25 records)
;using dim_avg

v1ave	= dim_avg(var1tmp)
v2ave	= dim_avg(var2tmp)

;need variance; averaged over time for each grid cell 
v1var = dim_variance(var1tmp)
v2var = dim_variance(var2tmp) 

;set the critical significance level
sigr = 0.05

;two samples have different population variances
iflag=True

;return only probabilities=False
;return value will be same dimensions as var1ave

tval_opt=False

prob=ttest(v1ave,v1var,25,v2ave,v2var,25,iflag,False)

;using where function to plot 
;using less than/greater than function
;for prob values less than 0.05, using where function

;xprob = prob
;xprob = where(prob.gt.0.05,0,prob)

;plot the results
;plot settings

res	= True

res@cnFillOn = True ; color fill
res@cnLinesOn = False ; no contour lines
res@cnLineLabelsOn = False ; no contour labels
res@gsnSpreadColors = True ; use total colormap

res@cnLevelSelectionMode	="ManualLevels"
res@cnMinLevelValF		=0.05
res@cnMaxLevelValF		=0.5
res@cnLevelSpacingF		=0.05

res@pmTickMarkDisplayMode  = "Always"           ; turn on tickmarks
        
res@tiMainFontHeightF      = 0.015             ; Set a small title
res@gsnAddCyclic           = False             ; Data is not cyclic; it is regional data
res@tfDoNDCOverlay         = True              ; Data is already on a map projection
							; Don't transform it ( =True)
res@mpGridAndLimbOn        = True               ; lat/lon lines on
res@mpGridSpacingF         = 10.0               ; lat/lon lines every 10 deg
res@mpGridLineDashPattern  = 2                  ; lat/lon lines as dashed
res@mpOutlineBoundarySets  = "GeophysicalAndUSStates" ; state boundaries

res@mpProjection           = "mercator"		;Set the map projection

res@mpCenterLatF	= 0.0
res@mpCenterLonF	= 120.0

;Set the boundaries of the map based on the corners of lat2d and lon2d
res@mpLimitMode            = "Corners"
res@mpLeftCornerLatF    =       min(lat2d)
res@mpLeftCornerLonF    =       min(lon2d)
res@mpRightCornerLatF   =       max(lat2d);(nlat-1,nlon-1)
res@mpRightCornerLonF   =       max(lon2d);(nlat-1,nlon-1)

;set small font

res@lbAutoManage	= False

res@lbLabelFontHeightF	= 0.01

;Plotting code

wks= gsn_open_wks("eps",sFileName)

plot= new(1,graphic)

gsn_define_colormap(wks,"gui_default")

;set up first plot
;set plot title

res@gsnLeftString = "Ttest results for Wet Season T2m raw prob values"

;set up plot

plot(0) =gsn_csm_contour_map(wks,prob(:,:,0),res)

end