load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;************************************************
; Conduct EOF analysis and make plots by YM Cheng
;************************************************

; Descriptiosn:
; Calculate EOFs of filtered data with period and season of interest
; Generate the associated PCs for the entire length of the dataset
; Output the EOFs and PCs in a nc file
; Plot the EOFs

begin
;;=============================================
;;=========== beginning of settings ===========
;;=============================================
neval = 8; number of EOFs to retain
latRange=(/-10,35/); domain to conduct EOF analysis
lonRange=(/-40,40/); domain to conduct EOF analysis
yrs=1984; year start, set the time period and season of interest for EOF
yre=2013; year end
mns=7; month start
mne=9; month end

diri="./filtered"; input directory
fin="2-6d.filtered.vwnd.700.nc"; input file
varname="vwnd" ; variable name in the filtered field

diro="."; output directory
fout="EOF.nc"; output EOF file name
EOFcor="False"; Not using correlation matrix

wksType="x11"; "png"
wksName=diro+"/"+"eof"; output fig name
;;=============================================
;;============== end of settings ==============
;;========== can leave the rest unchanged======
;;=============================================


print("varname       = "+varname)
print("filtered file = "+diri+"/"+fin)
print("EOF file      = "+diro+"/"+fout)
print("reading data...")

inf=addfile(diri+"/"+fin,"r")
vartmp=inf->$varname$
; test out if lonFlip is needed, i.e., applied when lon=[0,360)
if(any(getvardims(vartmp) .eq. "lon" .or. getvardims(vartmp) .eq. "LON"))
        if(vartmp&lon(dimsizes(vartmp&lon)-1) .gt. 180.)
                vartmp = lonFlip(vartmp)
        end if
end if
varTotal=vartmp(:,{latRange(0):latRange(1)},{lonRange(0):lonRange(1)})
; read in all the time for PC calculation of the full length of the dataset
; only read in the domain of interest

time=varTotal&time
cd      = cd_calendar(time,0)
indx    = ind(cd(:,0) .ge. yrs .and. cd(:,0) .le. yre .and. \
          cd(:,1) .ge. mns .and. cd(:,1) .le. mne)
          ; get the index for the period and season of interst for EOF
var     = vartmp(indx,{latRange(0):latRange(1)},{lonRange(0):lonRange(1)})
; read in the targeted domain, seasons and period for EOFs
delete(vartmp)
; printVarSummary(varTotal)
; printVarSummary(var)
print("finishing reading data")
print("=================================")
print("=========Calculating EOF=========")
print("=================================")
;*******************************************
; weight the data by area
;*******************************************

rad    = 4.*atan(1.)/180.
clat   = var&lat ; get the lat vector
clat   = sqrt( cos(rad*clat) )
varWeighted	= var ; copy meta data from var
varWeighted	= var*conform(var, clat, 1)
varWeighted := varWeighted(lat|:,lon|:,time|:)
; reorder the dimension to make time the rightmost for EOF computation

varTotalWeighted	= varTotal ; copy meta data from varTotal
varTotalWeighted	= varTotal*conform(varTotal, clat, 1)
varTotalWeighted  := varTotalWeighted(lat|:,lon|:,time|:)
; reorder the dimension to make time the rightmost for PC calculation


optEOF      = True
optETS      = True
if(EOFcor .eq. "True")
  optEOF@jopt = 1
  ; Use correlation matrix to compute EOF, default is convariance matrix
  optETS@jopt = 1
  ; Use the standardized data matrix to compute the time series.
  ; The default is to use data and evec.
  print("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")
  print("!!!!!!!!!!!!!!!! Using correlation matrix !!!!!!!!!!!!!!!!")
  print("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")
else
  optEOF@jopt=0
  optETS@jopt=0
end if
;*******************************************
; compute EOFs
;*******************************************
eof = eofunc_Wrap(varWeighted,neval,optEOF)
; apply EOF on varWeighted
eof_ts	= eofunc_ts_Wrap(varTotalWeighted,eof,optETS)
; get the associated PC for the lenght of the dataset
eof@domain="("+latRange(0)+","+latRange(1)+")/("+lonRange(0)+","+lonRange(1)+")"

; printVarSummary(eof)
; printVarSummary(eof_ts)
;*******************************************
; North significance test
;*******************************************
ndim   = dimsizes(varWeighted)
ntim   = ndim(2); take in the time dimension

prinfo = True; Print the calculated 'delta lambda', the eigenvalues, \
;the lower and upper test bounds and the 'separated significance' (sig).
print("==== North et al. 1982 test ====")
print("====Testing if eigenvalues are significantly separated====")
sig_ev  = eofunc_north(eof@eval, ntim, prinfo)
;*******************************************
; writing out EOF file
;*******************************************
system("/bin/rm -f " + diro +"/"+ fout)
outf=addfile(diro+"/"+fout,"c")
outf->eof=eof
outf->eof_ts=eof_ts

;*******************************************
; ============ plotting EOF ============
;*******************************************

plot = new((/neval,2/),graphic)
wks	=	gsn_open_wks(wksType,wksName)
; resources for maps
mapres = True
mapres@gsnDraw = False
mapres@gsnFrame = False
mapres@tmXBLabelFontHeightF	=	0.014
mapres@tmYLLabelFontHeightF	=	0.014
mapres@mpMaxLatF = latRange(1)
mapres@mpMinLatF = latRange(0)
mapres@mpMaxLonF = lonRange(1)
mapres@mpMinLonF = lonRange(0)
mapres@mpOutlineOn = True
mapres@mpGridAndLimbOn      = True ; default is every 15 deg
mapres@mpGeophysicalLineThicknessF = 5
mapres@mpOutlineDrawOrder = "PostDraw"
mapres@mpGridLineDashPattern  = 2  ; lat/lon lines dashed
mapres@mpGridLatSpacingF      = 10
mapres@mpGridLonSpacingF      = 15

; resources for shadings
xyres=True
xyres@gsnFrame =       False
xyres@gsnDraw  =       False
xyres@cnLineLabelPlacementMode =       "Constant"
xyres@cnFillOn = True
xyres@cnFillPalette = "BlWhRe"
xyres@cnLinesOn     = False ; turn off contour line
xyres@cnInfoLabelOn = False
xyres@lbLabelBarOn  = False ; turn off individual lb's
symMinMaxPlt(eof, 16, False, xyres); set symmetric plot min/max


; resources for panelling
resP                     = True         ; modify the panel plot
resP@gsnMaximize         = True         ; large format
resP@gsnPanelLabelBar    = True         ; add common colorbar
resP@lbLabelAutoStride   = True         ; auto stride on labels
resP@lbLabelFontHeightF  = 0.012

;; plotting the maps!
do n=0,neval-1
	xyres@gsnLeftString  = "EOF "+(n+1)
	xyres@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%"
	plot(n,0)	=	gsn_csm_map(wks,mapres)
	plot(n,1) = gsn_csm_contour(wks,eof(n,:,:),xyres)
	overlay(plot(n,0),plot(n,1))
end do
gsn_panel(wks,plot(:,0),(/4,2/),resP)     ; draw all eofs as one plot
end