; ==============================================================
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
; ==============================================================
; User defined parameters that specify region of globe and
; ==============================================================

; Latitudes: min=-47.9383   max=49.2609
; Longitude: min=-28.7721   max=64.7721

  latS   = -12.
  latN   =  23. 
  lonL   =  20.
  lonR   =  60.

  yrStrt = 1989
  yrLast = 2008

  season = "MAM"   

  neof   = 3       
  optEOF = True       
  optEOF@jopt = 0  
  optETS = False

; ==============================================================
; Open the file: Read only the user specified period
; ==============================================================
 ;;dir   = "/scratch/DATA/CORDEX/UC-WRF/"
   f  = addfile (dir+"AFRICA_UC-WRF311_CTL_ERAINT_MM_50km_1989-2008_pr.nc", "r")

   TIME = f->time
   YYYY = cd_calendar(TIME, -1)/100
   iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)

   pr = f->pr             ; (:,:,:)
   pr = (/ pr*86400 /)
   pr@units = "Kg/m2-day"
   printVarSummary(pr)     ; variable overview

   pr&lat   = f->lat(:,0)
   pr&lon   = f->lon(0,:)

   x = pr(:,{latS:latN},{lonL:lonR})
   printVarSummary(x)
   printMinMax(x&lat, 0)
   printMinMax(x&lon, 0)

   delete(pr)

; ==============================================================
;compute desired monthly climatology and monthly anomalies
;===============================================================

  xClm = clmMonTLL(x)
  printVarSummary(xClm)

  xAnom = calcMonAnomTLL(x,xClm)
  printVarSummary(xAnom)

; ================================================================
; Reorder (lat,lon,time) the *weighted* input data
; Access the area of interest via coordinate subscripting
; ================================================================
  xr      = xAnom(lat|:,lon|:,time|:)
  
  eof    = eofunc_Wrap(xr, neof, optEOF) 
; eof    =-1*eof
  eof_ts = eofunc_ts_Wrap (xr, eof, optETS)

  printVarSummary( eof )                         
  printVarSummary( eof_ts )

;============================================================
; PLOTS
;============================================================
  plot = new(neof,graphic)                ; create graphic array
  plot2 = plot

  xwks = gsn_open_wks("x11","eof")
  pswks = gsn_open_wks("eps","Africa-WRF-Cordex_eof-mam1989-2008") 
  gsn_define_colormap(pswks,"BlWhRe")
                                          ; only needed if paneling
; EOF patterns

  res                      = True         
  res@gsnDraw              = False        ; don't draw yet
  res@gsnFrame             = False        ; don't advance frame yet
  res@gsnSpreadColors      = True         ; spread out color table
  res@gsnAddCyclic         = False        ; plotted dataa are not cyclic
 
  res@mpFillOn                    = False        ; turn off map fill
  res@mpOutlineBoundarySets       = "National" 
  res@mpGeophysicalLineColor      = "Navy"
;  res@mpGeophysicalLineThickness = 1.5
  
  res@mpMinLatF            = latS         ; zoom in on map
  res@mpMaxLatF            = latN
  res@mpMinLonF            = lonL
  res@mpMaxLonF            = lonR

  res@cnFillOn             = True         ; turn on color fill
  res@cnLinesOn            = False        ; True is default
;  res@cnLineLabelsOn       = False        ; True is default
  res@lbLabelBarOn         = False        ; turn off individual lb's

                                          ; set symmetric plot min/max
  symMinMaxPlt(eof, 16, False, res)       ; contributed.ncl

; panel plot only resources
  resP                     = True         ; modify the panel plot
  resP@gsnMaximize         = True         ; large format
  resP@gsnPanelLabelBar    = True         ; add common colorbar
  resP@lbLabelAutoStride   = True         ; auto stride on 

;*******************************************
; first plot
;*******************************************
  do n=0,neof-1
     res@gsnLeftString  = "EOF "+(n+1)
     res@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%"
     plot(n)=gsn_csm_contour_map_ce(xwks,eof(n,:,:),res)
	    plot2(n)=gsn_csm_contour_map_ce(pswks,eof(n,:,:),res)
  end do
  gsn_panel(xwks,plot,(/neof,1/),resP)     ; now draw as one plot
  gsn_panel(pswks,plot2,(/neof,1/),resP)

;*******************************************
; second plot
;*******************************************
;; EOF time series  [bar form]

   time = eof_ts&time
  do n=0,neof-1
    ascii_file1 = "WRF-cordex_eof1-mam1989-2008.txt" 
	ascii_file2 = "WRF-cordex_eof2-mam1989-2008.txt"
	ascii_file3 = "WRF-cordex_eof3-mam1989-2008.txt"

	asciiwrite(ascii_file1,eof_ts(0,:))
	asciiwrite(ascii_file2,eof_ts(1,:))
	asciiwrite(ascii_file3,eof_ts(2,:))
  end do
end