; ==============================================================
; 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 = "G_JRA55_test1_vdc013l"
;  datadir = "/work/yiwen89419/CCSM_ncl/G_JRA55_test1_vdc013l/"
;  netcdf = "MOCtot.306-366.eof.atl.lz10+.nodt.nc"

;  casename = "b_20ct_ctrl_1"  ;EOF1*(-1)
;  datadir = "/work/yiwen89419/CCSM_ncl/"
;  netcdf = "./Monthly/MOCtot.1921-2000.eof.atl.lz10+.nodt.nc"

  casename = "taiesm"  ;EOF1*(-1)
  datadir = "/work/yiwen89419/CCSM_ncl/"
  netcdf = "./f09.B-PI.tn15.cmip6.j01.re/yearly/MOCtot.500-700.eof.atl.lz10+.nodt.nc"
 
;  sign = -1.

  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./1000.
  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 = 1  ; 0 == remove mean 
jave = 5 ; 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
  gsn_define_colormap(wks,"BlueWhiteOrangeRed")         ; 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    = 255                   ; 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	   = -1.4
  res@cnMaxLevelValF	   =  1.4
  res@cnLevelSpacingF	   = 0.1
  res@gsnMaximize           = True         ; enlarge plot 
  res@trYReverse           = True       ; reverses y-axis
  res@gsnYAxisIrregular2Linear = True
  res@gsnXAxisIrregular2Linear = True
  res@tiYAxisString= "depth (km)"
  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           = True         ; enlarge plot 
  resxy@tiYAxisString= ""
  resxy@tiXAxisString= ""
  resxy@xyLineColors      = (/"black"/)          ; 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 = "Frequency x Power (Sv~S~2~N~)"
spec@tiXAxisString = ""
spec@gsnMaximize           = True         ; 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@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.5
 spec@trYMinF = 0.0
 spec@trYMaxF = 1.

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

lres = True
lres@gsLineDashPattern = 0
lres@gsLineColor = "grey"
nxlines = dimsizes(spec@tmXBValues) - 2
specxlines1 = new(nxlines,graphic)
specxlines2 = new(nxlines,graphic)


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

  res@gsnLeftString = "ci = "+flt2string(res@cnLevelSpacingF)+" Sv"
  res@gsnRightString = "EOF 1,   "+sprintf("%4.2f",eof@pcvar(0) ) + "%"
  plot(0) = gsn_csm_contour(wks,eof(0,:,:),res)

  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*conform(spectrum,dof@frq,(/1/)),spec)
  delete(tmp)
  delete(dof)
  delete(spectrum)
  getvalues plot(2)                        ; retrieve some of the plot resources
     "tmYLValues"  : tmYLValues          ; values used by NCL at major tick marks
     "tmYLMinorValues"  : tmYLMinorValues          ; values used by NCL at major tick marks
  end getvalues
  do i=0,nxlines-1
    x0 = spec@tmXBValues(i+1)
    specxlines1(i) = gsn_add_polyline(wks,plot(2),(/x0,x0/),(/spec@trYMinF,spec@trYMaxF/),lres)
  end do
  nyMajlines = dimsizes(tmYLValues)
  nyMinlines = dimsizes(tmYLMinorValues)
  specyMajlines1 = new(nyMajlines,graphic)
  specyMinlines1 = new(nyMinlines,graphic)
  do i=0,nyMajlines-1
    y0 = tmYLValues(i)
    specyMajlines1(i) = gsn_add_polyline(wks,plot(2),(/spec@trXMinF,spec@trXMaxF/),(/y0,y0/),lres)
  end do
  do i=0,nyMinlines-1
    y0 = tmYLMinorValues(i)
    specyMinlines1(i) = gsn_add_polyline(wks,plot(2),(/spec@trXMinF,spec@trXMaxF/),(/y0,y0/),lres)
  end do

  delete(tmYLValues)
  delete(tmYLMinorValues)
  delete(nyMajlines)
  delete(nyMinlines)

  res@gsnRightString = "EOF 2,   "+sprintf("%4.2f",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 = ""
  spec@tiXAxisString = "Period (years)"
  plot(5)=gsn_csm_xy(wks,dof@frq,spectrum*conform(spectrum,dof@frq,(/1/)),spec)
  delete(tmp)
  delete(dof)
  delete(spectrum)
  getvalues plot(5)                         ; retrieve some of the plot resources
     "tmYLValues"  : tmYLValues          ; values used by NCL at major tick marks
     "tmYLMinorValues"  : tmYLMinorValues          ; values used by NCL at major tick marks
  end getvalues
  do i=0,nxlines-1
    x0 = spec@tmXBValues(i+1)
    specxlines2(i) = gsn_add_polyline(wks,plot(5),(/x0,x0/),(/spec@trYMinF,spec@trYMaxF/),lres)
  end do
  nyMajlines = dimsizes(tmYLValues)
  nyMinlines = dimsizes(tmYLMinorValues)
  specyMajlines2 = new(nyMajlines,graphic)
  specyMinlines2 = new(nyMinlines,graphic)
  do i=0,nyMajlines-1
    y0 = tmYLValues(i)
    specyMajlines2(i) = gsn_add_polyline(wks,plot(5),(/spec@trXMinF,spec@trXMaxF/),(/y0,y0/),lres)
  end do
  do i=0,nyMinlines-1
    y0 = tmYLMinorValues(i)
    specyMinlines2(i) = gsn_add_polyline(wks,plot(5),(/spec@trXMinF,spec@trXMaxF/),(/y0,y0/),lres)
  end do


;************************************************
; 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