; ==============================================================
; eof_1.ncl
;
; Concepts illustrated:
;   - Calculating EOFs
;   - Using coordinate subscripting to read a specified geographical region
;   - Rearranging longitude data to span -180 to 180
;   - Calculating symmetric contour intervals
;   - Drawing filled bars above and below a given reference line
;   - Drawing subtitles at the top of a plot
;
; ==============================================================
; Calculate EOFs of the Sea Level Pressure over the North Atlantic.
; ==============================================================
; The slp.mon.mean file can be downloaded from:
; http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.surface.html
; ==============================================================

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
; ==============================================================
  latS   =  -14.
  latN   =  17. 
  lonW   =  22.
  lonE   =  54.

  yrStrt = 1950
  yrLast = 2008

  season = "MAM"    ; choose Mar-Apr-May seasonal mean

  neof   = 1       ; number of EOFs
  optEOF = True       
  optEOF@jopt = 0   ; This is the default; most commonly used; no need to specify.
;; optEOF@jopt = 1   ; **only** if the correlation EOF is desired

  optETS = False

; ==============================================================
; Open the file: Read only the user specified period
; ==============================================================
  f      = addfile ("cru.nc", "w")
 ;print(f)

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

  cru    = f->pre(iYYYY,{latS:latN},{lonW:lonE})
  printVarSummary(cru)                              ; variable overview


; ==============================================================
; compute desired global seasonal mean: month_to_season (contributed.ncl) 
; ==============================================================
  SLP    = month_to_season (cru, season)
  nyrs   = dimsizes(SLP&time)
  printVarSummary(SLP)

; =================================================================
; create weights:  sqrt(cos(lat))   [or sqrt(gw) ]
; =================================================================
  rad    = 4.*atan(1.)/180.
  clat   = SLP&lat           
  clat   = sqrt( cos(rad*clat) )                 ; gw for gaussian grid

; =================================================================
; weight all observations 
; =================================================================
  wSLP   = SLP                                   ; copy meta data
  wSLP   = SLP*conform(SLP, clat, 1)
  wSLP@long_name = "Wgt: "+wSLP@long_name

; =================================================================
; Reorder (lat,lon,time) the *weighted* input data
; Access the area of interest via coordinate subscripting
; =================================================================
  x      = wSLP({lat|latS:latN},{lon|lonW:lonE},time|:)

  eof    = eofunc_Wrap(x, neof, optEOF)      
  eof_ts = eofunc_ts_Wrap (x, eof, optETS)
  ; print(eof_ts)
  printVarSummary( eof )                         ; examine EOF variables
  printVarSummary( eof_ts )

;;prinfo = True
;;sig_ev = eofunc_north(eof@eval, nyrs, prinfo)  ; 6.3.0 onward

; =================================================================
; Normalize time series: Sum spatial weights over the area of used
; =================================================================
  dimx   = dimsizes( x )
  mln    = dimx(1)
  sumWgt = mln*sum( clat )
  eof_ts = eof_ts/sumWgt
  print(eof_ts)
 
; ==============================================================
; Open the file: Read only the user specified period
; ============================================================== 
  g= addfile ("yycompos.1DYsk7UH7u.nc", "r")

  wind = g->uwnd
  ;lat = g->lat
  ;lon = g->lon
  ;TIME = g->time
  printVarSummary(wind) ; variable overview
  sst1 =wind(:,:,:)

 printVarSummary(eof_ts) ; variable overview
  eof1 = eof_ts(0,:)
  print(eof_ts)

; =================================================================
; Correlations calculation
; =================================================================

   q = eof_ts(evn|0,time|:)

   y = wind(lat|:,lon|:,time|:)

    ; y&lat@units = "degrees_north"
    ; y&lon@units = "degrees_east"

  ccr = escorc(q, y)

 printVarSummary(ccr)

  ccr!0 = "lat" ; name dimensions
  ccr!1 = "lon"
  ccr&lat = y&lat ; assign coordinate values and
  ccr&lon = y&lon ; units attributes

 print(ccr) 
; =================================================================
; Extract the YYYYMM from the time coordinate 
; associated with eof_ts [same as SLP&time] 
; =================================================================

  yyyymm = cd_calendar(eof_ts&time,-2)/100  
;;yrfrac = yyyymm_to_yyyyfrac(yyyymm, 0.0); not used here

;============================================================
; PLOTS
;============================================================
  wks = gsn_open_wks("pdf","eof")
  gsn_define_colormap(wks,"GMT_drywet")
                                           ; choose colormap
  plot = new(neof,graphic)                ; create graphic array
                                          ; only needed if paneling
; EOF patterns

  res                      = True         
  res@gsnDraw              = False        ; don't draw yet
  res@gsnFrame             = False        ; don't advance frame yet

;---This resource not needed in V6.1.0
  res@gsnSpreadColors      = True         ; spread out color table

  res@gsnAddCyclic         = False        ; plotted dataa are not cyclic
 ; res@cnFillPalette        = "RedBlue"
  res@mpFillOn             = False        ; turn off map fill
  res@mpMinLatF            = latS         ; zoom in on map
  res@mpMaxLatF            = latN
  res@mpMinLonF            = lonW
  res@mpMaxLonF            = lonE

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

  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 labels

  yStrt                    = yyyymm(0)/100
  yLast                    = yyyymm(nyrs-1)/100
  resP@txString            = "PRECIPITATION: "+season+": "+yStrt+"-"+yLast

;*******************************************
; 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(wks,eof(n,:,:),res)
  end do
  gsn_panel(wks,plot,(/neof,1/),resP)     ; now draw as one plot

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

  rts           = True
  rts@gsnDraw   = False       ; don't draw yet
  rts@gsnFrame  = False       ; don't advance frame yet
  rts@gsnScale  = True        ; force text scaling               

; these four rtsources allow the user to stretch the plot size, and
; decide exactly where on the page to draw it.

  rts@vpHeightF = 0.40        ; Changes the aspect ratio
  rts@vpWidthF  = 0.85
  rts@vpXF      = 0.10        ; change start locations
  rts@vpYF      = 0.75        ; the plot


  rts@tiYAxisString = "AMPLITUDES"                    ; y-axis label      

  rts@gsnYRefLine           = 0.              ; reference line   
  rts@gsnXYBarChart         = True            ; create bar chart 
  rts@gsnAboveYRefLineColor = "blue"           ; above ref line fill red
  rts@gsnBelowYRefLineColor = "red"          ; below ref line fill blue

; panel plot only resources
  rtsP                      = True            ; modify the panel plot
  rtsP@gsnMaximize          = True            ; large format
  rtsP@txString             = "PRECIPITATION: "+season+": "+yStrt+"-"+yLast

  year = yyyymm/100

; create individual plots
  do n=0,neof-1
     rts@gsnLeftString  = "EOF "+(n+1)
     rts@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%"
     plot(n) = gsn_csm_xy (wks,year,eof_ts(n,:),rts)
  end do
  gsn_panel(wks,plot,(/neof,1/),rtsP)     ; now draw as one plot
;*******************************************
;correlation plot
;*******************************************
  wks = gsn_open_wks("X11","globalmap")
  gsn_define_colormap(wks,"gui_default") ; choose colormap

  rescn = True
  rescn@cnFillOn = True
  ;rescn@gsnDraw = False ; don't draw
  ;rescn@gsnFrame = False ; don't advance frame

  rescn@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
  rescn@cnMinLevelValF = -0.4 ; set min contour level
  rescn@cnMaxLevelValF = 0.4 ; set max contour level
  rescn@cnLevelSpacingF = 0.2 ; set contour spacing

  rescn@lbOrientation = "Vertical" ; vertical label bar

;---This resource defaults to True in NCL V6.1.0
  rescn@lbLabelAutoStride = True ; optimal label stride

  rescn@gsnSpreadColors = True ; use full range of colors
  ;rescn@gsnSpreadColorEnd = -3 ; don't use added gray

  rescn@mpCenterLonF = 180. ; center plot at 180

  rescn@gsnAddCyclic = True

     plot = gsn_csm_contour_map(wks,ccr,rescn)

end