;===============================================;
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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
; ================================================;
begin
  ;; data files

  casename = "b40.1850.track1.1deg.006"
  datadir = "/work/yiwen89419/CCSM_ncl/paper_data/"
;  f1 = "./sm5/MOCtot.957-1108.eof.atl.lz10+.nodt.nc"
  f1 = "./sm5/MOCtot.957-1108.eof.atl.nodt.nc"

  year0 = 957.
  year1 = 1109.

  psout = "spectrum.AMOCsm5pc1."+casename+"."+flt2string(year0)+"-"+flt2string(year1-1)

  nc2     = addfile(f1,"r")
  eofts   = nc2->eofts
  time    = nc2->time

  tmpr = eofts(0,:)
  dims = dimsizes(tmpr)
  nt = dims(0)

  wgt  = (/1.0, 1.0, 1.0, 1.0, 1.0/)  ; 5-year boxcar
  wgt  = wgt/sum(wgt)             ; normalize

  curves = new((/1,nt/),float,eofts@_FillValue)
  curves!0 = "ncurve"
  curves!1 = "time"
  curves&time = time
  curves(0,:) = (/ tmpr /)

;=================================================;
; Create plot
;=================================================;
  wks  = gsn_open_wks("ps",psout)             ; open a ps file
; gsn_define_colormap(wks,"BkBlAqGrYeOrReViWh200")
  gsn_define_colormap(wks,"wgne15")
; i    = NhlNewColor(wks,0.8,0.8,0.8)      ; add gray to map

;************************************************
; resource list for time series spectra
;************************************************
opt = 1  ; 0 == remove mean, 1 = detrend
jave = 3 ; number of periodograms to average (if < 3, no smoothing!)
pct = 0.10  ; percent taper
spec = True
spec@gsnDraw = False
spec@gsnFrame = False
spec@tiYAxisString = "Variance/(unit freq.)"
spec@tiXAxisString = "Period (years)"
spec@gsnMaximize           = True         ; enlarge plot
spec@gsnLeftString = ""
spec@gsnRightString = ""
spec@gsnCenterString = ""
spec@xyLineThicknesses = (/1,1,1,1/)
spec@xyDashPatterns = (/0,0,1,1/)
spec@xyLineColors = (/"black","grey","grey","grey"/)
spec@tmXTOn              = False
spec@tmYROn              = False
spec@tmXBOn              = True
spec@tmXBMode            = "Explicit"
 spec@xyXStyle = "Log"
 spec@xyYStyle = "Log"
 spec@tmXBValues = (/(1./nt),0.005,0.01,0.02,0.05,0.1,0.2,0.5/)
 low = ""
 low = nt
 spec@tmXBLabels = (/low,"200","100","50","20","10","5","2"/)
 spec@trXMinF = 1./nt
;  spec@trXMaxF = 0.2
 spec@trXMaxF = 0.5
 spec@trYMinF = 1.
 spec@trYMaxF = 100.
;spec@trYMinF = 0.
;spec@trYMaxF = 120.

  tmp = curves(0,:)
  dof = specx_anal(tmp(ind(.not.ismissing(tmp))),opt,jave,pct)
  spectrum = specx_ci(dof,0.95,0.99)
  spec@gsnCenterString = casename+": AMOC PC1 Timeseries "+flt2string(year0)+"-"+flt2string(year1-1)
  plot1=gsn_csm_xy(wks,dof@frq,spectrum,spec)
  delete(spectrum)

; tmp = curves(0,:)
; jave = 3 ; number of periodograms to average (if < 3, no smoothing!)
; dof = specx_anal(tmp(ind(.not.ismissing(tmp))),opt,jave,pct)
; spectrum = specx_ci(dof,0.95,0.99)
; spec@gsnCenterString = casename+": AMOC PC1 Timeseries "+flt2string(year0)+"-"+flt2string(year1-1)
; spec@xyLineThicknesses = (/3,1,1,1/)
; spec@xyLineColors = (/"blue","grey","grey","grey"/)
; plot2=gsn_csm_xy(wks,dof@frq,spectrum,spec)
; delete(spectrum)
; overlay(plot1,plot2)


  draw(plot1)
  frame(wks)
end