; ==============================================================
; Calculate AMOC EOFs 
; ==============================================================
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
; ==============================================================
; User defined parameters that specify region of globe to consider
; ==============================================================
yearadd = 0.
sign = 1.

  casename = "b40.1850.track1.1deg.006"
  datadir = "/work/yiwen89419/CCSM_ncl/paper_data/"
  netcdf = "./sm5/MOCtot.957-1108.eof.atl.lz10+.nodt.nc"
; netcdf = "./sm5/MOCtot.957-1108.eof.atl.nodt.nc"
 
  fin= datadir+netcdf
  sfx = get_file_suffix(fin,0)
  psout = sfx@fBase
  
  rEarth = 6.37122e8 		; Earth radius in m

; ==============================================================
  f      = addfile (fin, "r")
  eof    = f->eof
  eofts  = f->eofts
  mocz   = f->moc_z
  lat    = f->lat_aux_grid
  time   = f->time

  time = time + yearadd
  eof = eof*sign
  eofts = eofts*sign

  mocz = mocz/100.
  mocz@units = "m"
  eof&moc_z = mocz
  eof&moc_z@units = "m"
  dims   = dimsizes(time)
  nt   = dims(0)

; spectra = new((/4,nt/),float,eofts@_FillValue)
; pcts!0 = "ncurve"
; pcts!1 = "time"
; pcts&time = time

  eof@long_name = ""

;=================================================;
; spectra options
;=================================================;
opt = 0  ; 0 == remove mean 
jave = 1 ; number of periodograms to average (if < 3, no smoothing!)
pct = 0.10  ; percent taper

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

;************************************************
; resource list for EOF contours
;************************************************
  res                      = True                 ; plot mods desired
  res@gsnDraw		   = False		
  res@gsnFrame		   = False		
  res@cnFillOn             = True                 ; turn on color fill
  res@cnMissingValFillColor = "grey"
  res@cnMissingValFillPattern = 0
  res@lbLabelAutoStride    = True		  ; control labelbar labels
  res@cnLinesOn            = True                ; turn off contour lines
  res@cnLineLabelsOn 	   = True		; turn the line labels off
  res@gsnSpreadColors      = True                 ; use full colormap
  res@gsnSpreadColorStart  = 2                   ; start at color 10
  res@gsnSpreadColorEnd    = 17                   ; end at color 96
  res@lbLabelBarOn        = False       ; Turn off labelbar
  res@cnInfoLabelOn       = False       ; Turn off informational label
  res@cnLevelSelectionMode = "ManualLevels"	; manually set the contour levels
  res@cnMinLevelValF	   = -2.0
  res@cnMaxLevelValF	   =  2.0
  res@cnLevelSpacingF	   = 0.1
  res@gsnMaximize           = False         ; enlarge plot 
  res@trYReverse           = True       ; reverses y-axis
  res@gsnYAxisIrregular2Linear = True
  res@gsnXAxisIrregular2Linear = True
  res@tiYAxisString= "depth (m)"
  res@tiXAxisString= ""
  res@gsnCenterString = ""
  res@trXMinF = -30.
  res@trXMaxF = 90.

;************************************************
; resource list for EOF PC time series
;************************************************
  resxy = True
  resxy@gsnDraw            = False
  resxy@gsnFrame           = False
  resxy@gsnMaximize           = False         ; enlarge plot 
  resxy@tiYAxisString= ""
  resxy@tiXAxisString= ""
  resxy@xyLineColors      = (/"blue"/)          ; change line color
  resxy@xyMarkLineModes   = (/"lines"/)
  resxy@xyLineThicknesses = (/1/)
  resxy@xyDashPatterns   = (/0/)      ; make all lines solid
  resxy@trXMinF  = min(time)                   ; min value on y-axis
  resxy@trXMaxF  = max(time)                   ; max value on y-axis
  resxy@trYMinF  = -3.                   ; min value on y-axis
  resxy@trYMaxF  = 3.                   ; max value on y-axis
  resxy@tmYLMode = "Manual"
  resxy@tmYLTickSpacingF = 1.

;************************************************
; resource list for EOF PC time series spectra
;************************************************
spec = True
spec@gsnDraw = False
spec@gsnFrame = False
spec@tiYAxisString = "Variance/(unit freq.)"
spec@tiXAxisString = ""
spec@gsnMaximize           = False         ; enlarge plot 
spec@gsnLeftString = ""
spec@gsnRightString = ""
spec@gsnCenterString = ""
spec@xyLineThicknesses = (/2.,2.,1.,1./)
spec@xyDashPatterns = (/0,0,1,1/)
spec@xyLineColors = (/"black","red","red","red"/)
spec@tmXTOn              = False
spec@tmYROn              = False
spec@tmXBOn              = True
spec@tmXBMode            = "Explicit"
 spec@xyXStyle = "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.5
;spec@trYMaxF = 120.
;spec@trYMinF = 0.0
spec@trYMaxF = 20.

;************************************************
; Draw text on plot using plot coordinates.
;************************************************
  txres               = True                     ; text mods desired
  txres@txFontHeightF = 0.03                     ; font smaller. default big


;***********************************************
; Create array to hold plots.
  plot = new(6,graphic)

  res@gsnLeftString = "ci = "+flt2string(res@cnLevelSpacingF)+" Sv"
  res@gsnRightString = "EOF 1,   "+flt2string(eof@pcvar(0))+"%"
  plot(0) = gsn_csm_contour(wks,eof(0,:,:),res)
; label = "ci = "+flt2string(res@cnLevelSpacingF)
; text = gsn_add_text(wks,plot(0),label,70.,4500.,txres)
; label = "Sv"
; text = gsn_add_text(wks,plot(0),label,70.,5000.,txres)

  resxy@gsnCenterString = "PC 1 time series"
  plot(1) = gsn_csm_xy(wks,time,eofts(0,:),resxy)
  tmp = eofts(0,:)
  dof = specx_anal(tmp(ind(.not.ismissing(tmp))),opt,jave,pct)
  spectrum = specx_ci(dof,0.95,0.99)
  spec@gsnCenterString = "Power Spectrum"
  plot(2)=gsn_csm_xy(wks,dof@frq,spectrum,spec)
  delete(tmp)
  delete(dof)
  delete(spectrum)

  res@gsnRightString = "EOF 2,   "+flt2string(eof@pcvar(1))+"%"
  plot(3) = gsn_csm_contour(wks,eof(1,:,:),res)
  resxy@gsnCenterString = "PC 2 time series"
  plot(4) = gsn_csm_xy(wks,time,eofts(1,:),resxy)
  tmp = eofts(1,:)
  dof = specx_anal(tmp(ind(.not.ismissing(tmp))),opt,jave,pct)
  spectrum = specx_ci(dof,0.95,0.99)
  spec@gsnCenterString = ""
  plot(5)=gsn_csm_xy(wks,dof@frq,spectrum,spec)
  delete(tmp)
  delete(dof)
  delete(spectrum)

;************************************************
; multipanel stuff
;************************************************
  pres = True
; pres@txString   = casename
  pres@txFontHeightF = 0.01
  pres@gsnPanelLabelBar  = False
  pres@lbLabelAutoStride = True
  pres@gsnMaximize     = True                 ; fill up the page
  pres@gsnPanelRowSpec  = True         ; Specify plots per row
  gsn_panel(wks,plot,(/3,3/),pres)
  end