; ==============================================================
; 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.
sign= 1.

; casename = "b40.1850.track1.1deg.006"
; datadir = "/ccsm/ocn/yeager/amoc_ts/"+casename+"/eofs/"
; netcdf = "MOC.700-1299.eof.atl.nc"

  casename = "taiesm"  ;
  datadir = "/work/yiwen89419/CCSM_ncl/"
;  netcdf = "./f09.B-PI.tn15.cmip6.j01.re/yearly/MOCtot.500-700.eof.atl.lz10+.nodt.nc"
;  netcdf = "MOCtot.500-700.eof.atl.nc"              ;detrend & smooth
;  netcdf = "MOCtot.500-700.eof.atl.lz10+.nodt.nc"  ;not detrend & smooth  [default]
;  netcdf = "MOCtot.500-700.eof.atl.lz10+.nc"       ;detrend & smooth(x)
;  netcdf = "MOCtot.500-700.eof.atl.nodt.nc"         ;not detrend & smooth(x) 
;  netcdf = "paper_data/sm5/MOCtot.957-1108.eof.atl.nodt.nc"  ; test change yr
   netcdf = "./G_JRA55_test1_vdc013l/5_6th_cycle/MOCtot.245-366.eof.atl.nodt.nc" ;

  fin= datadir+netcdf
  sfx = get_file_suffix(fin,0)
;  psout = "Fig2_AMOCeofs"
  psout = "G_JRA55_test1_vdc013l_EOF"

  
  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(0,:,:) = eof(0,:,:)*sign
  eofts(0,:) = eofts(0,:)*sign

  mocz = mocz/100.
  mocz = mocz/1000.
  mocz@units = "km"
  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!)
;jave = 1 ; number of periodograms to average (if < 3, no smoothing!) [default]
pct = 0.10  ; percent taper

;=================================================;
; Create plot
;=================================================;
  wks  = gsn_open_wks("ps",psout)             ; open a ps file
  gsn_define_colormap(wks,"gsdtol")
  blk    = NhlNewColor(wks,0.,0.,0.)          ; add gray to map
  grey    = NhlNewColor(wks,0.5,0.5,0.5)      ; add gray to map
  ltgrey    = NhlNewColor(wks,0.8,0.8,0.8)    ; add gray to map
  wt    = NhlNewColor(wks,1.,1.,1.)           ; 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 = blk
  res@cnMissingValFillPattern = 0
  res@lbLabelAutoStride    = True		  ; control labelbar labels
  res@cnLinesOn            = False                ; turn off contour lines
  res@cnLineLabelsOn 	   = False		; turn the line labels off
  res@gsnSpreadColors      = True                 ; use full colormap
  res@gsnSpreadColorStart  = 32                   ; start at color 10
  res@gsnSpreadColorEnd    = 8                   ; end at color 96
  res@gsnLeftStringFontHeightF = 0.025
  res@gsnRightStringFontHeightF = 0.025
  res@lbLabelBarOn        = False       ; Turn off labelbar
  res@cnInfoLabelOn       = False       ; Turn off informational label
  res@cnLevelSelectionMode = "ManualLevels"	; manually set the contour levels
;  res@cnMinLevelValF	   = 0.0
;  res@cnMaxLevelValF	   =  1.0
  res@cnMinLevelValF      = -2.0
  res@cnMaxLevelValF      =  2.0
  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.
  res@tmXTOn = False
  res@tmYROn = False
  res@tmXBLabelsOn = True
  res@tmXBMode = "Explicit"
  res@tmXBValues = (/-30.,0.,30.,60.,90./)
  res@tmXBLabels = (/"30~S~o~N~S","0~S~o~N~N","30~S~o~N~N","60~S~o~N~N","90~S~o~N~N"/)
 
  lineres = True
  lineres@gsnDraw              = False
  lineres@gsnFrame             = False
  lineres@cnFillOn             = False                 ; turn on color fill
  lineres@cnLinesOn            = True                ; turn off contour lines
  lineres@cnLineLabelFontHeightF = 0.018 
  lineres@cnLineLabelsOn       = True               ; turn the line labels off
  lineres@cnLevelSelectionMode = "ManualLevels"  ; manually set the contour levels
  lineres@cnMinLevelValF	   = -1.0
  lineres@cnMaxLevelValF	   =  1.0
  lineres@cnLevelSpacingF	   = 0.1
  lineres@cnMonoLineThickness = False
  lineres@cnMonoLineDashPattern = False
  lineres@cnLineThicknesses    = (/1,1,1,1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1/)
  lineres@cnLineDashPatterns   = (/1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0/)
  lineres@cnInfoLabelOn = False




;************************************************
; 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@gsnLeftStringFontHeightF = 0.025
  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.
  resxy@tmXTOn = False
  resxy@tmYROn = False

;************************************************
; 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@gsnLeftStringFontHeightF = 0.025
spec@gsnRightStringFontHeightF = 0.025
spec@xyLineThicknesses = (/2.,2.,1.,1./)
spec@xyDashPatterns = (/0,0,1,1/)
spec@xyLineColors = (/blk,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.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@gsnLeftString = "a. AMOC EOF1"
  res@gsnRightString = sprintf("%3.1f", eof@pcvar(0))+"%"
  tmp = eof(0,:,:)
  tmp = abs(eof(0,:,:))
  plot(0) = gsn_csm_contour(wks,tmp,res)
  ploto = gsn_csm_contour(wks,eof(0,:,:),lineres)
  overlay (plot(0), ploto)
  delete(tmp)

  resxy@gsnLeftString = "b. PC1 Time Series"
; resxy@gsnCenterString = ""
  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.)
  spec@gsnLeftString = "c. PC1 Power Spectrum"
  spec@gsnCenterString = ""
  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@gsnLeftString = "d. AMOC EOF2"
  res@gsnRightString = sprintf("%3.1f", eof@pcvar(1))+"%"
  res@tiXAxisString = "Latitude"
  tmp = eof(1,:,:)
  tmp = abs(eof(1,:,:))
  plot(3) = gsn_csm_contour(wks,tmp,res)
  ploto3 = gsn_csm_contour(wks,eof(1,:,:),lineres)
  delete(tmp)
  overlay (plot(3), ploto3)
  resxy@gsnLeftString = "e. PC2 Time Series"
  resxy@tiXAxisString = "Time (years)"
  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.)
  spec@gsnLeftString = "f. PC2 Power Spectrum"
  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@txFontHeightF = 0.01
  pres@gsnPanelLabelBar  = False
  pres@lbLabelAutoStride = True
  pres@gsnMaximize     = True                 ; fill up the page
  pres@gsnPanelRowSpec  = True         ; Specify plots per row
  pres@gsnPanelRight = 0.95
  gsn_panel(wks,plot,(/3,3/),pres)
  end