load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   ; High Level
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"    ; plot interfaces
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;**************************************************************** 

begin

a = addfile("/ptmp/notaro/run/ho.foam1-5lpj.ctl-pi-test.annual.sst.6351034.nc","r")   ; data file
neval = 5   ; # eigenvectors to compute
data = a->T(:,0,{-40:40},{210:350}) ;change variable name and specify the area;;
;;;;;;;;;;;;;;;;;;;;;do not change for most of cases;;;;;;;;;;;;;;;;;;;;;;;;;;
  data@_FillValue = 1e+35  ; missing value
  var = data(lat|:,lon|:,time|:)   ; reorder so time is last dim
  var@_FillValue = 1e+35 ; missing value

  var=dim_rmvmean(var)   ; remove mean
  dims=dimsizes(var)     ; get dim sizes for variable
  lat=var&lat            ; assign latitudes to variable
  do nl=0,dims(0)-1                      ; wgt by sqrt(cos(latitude))
    var(nl,:,:)  = var(nl,:,:)*sqrt(cos(lat(nl)*0.0174532))  
  end do

;;;;;;;;;;;;;

eof=eofcov_pcmsg_Wrap(var,neval,95.)  ; calculate eof
eof@_FillValue=1e+35   ; missing value
eof_ts = eofcov_ts_Wrap(var,eof) ; calculate time series of eof
eof_ts@_FillValue=1e+35 ; missing value

;;;;;;;;;;;;;;;;;

; red noise curve, 10% signif, 90% signif 
  spec1 = specx_anal(eof_ts(0,:),1,7,0.10)      ; AC spectrum of eof1 time series
  splt1 = specx_ci(spec1, 0.1, 0.9)             ; confid interval for eof 1 time series
  spec2 = specx_anal(eof_ts(1,:),1,7,0.10)      
  splt2 = specx_ci(spec2, 0.1, 0.9)             
  spec3 = specx_anal(eof_ts(2,:),1,7,0.10)
  splt3 = specx_ci(spec3, 0.1, 0.9)            
  spec4 = specx_anal(eof_ts(3,:),1,7,0.10)
  splt4 = specx_ci(spec4, 0.1, 0.9)           
  spec5 = specx_anal(eof_ts(4,:),1,7,0.10)
  splt5 = specx_ci(spec5, 0.1, 0.9)          

 wks   = gsn_open_wks("ps" ,"eof_sst")        ; open ps file
 plot  = new ( 12, "graphic")                ; create graphic array

 res                     = True            ; plot mods desired
 res@gsnDraw              = False          ; don't draw
 res@gsnFrame             = False           ; don't advance frame
 res@gsnAddCyclic                = False         ; regional data
 res@mpMinLatF    = -40.0                ; only plot 20N to 90N
 res@mpMaxLatF    = 40.0
 res@mpCenterLonF         = 180.               ; center plot at 180
 res@mpMinLonF = 210.
 res@mpMaxLonF = 350.
 res@cnLevelSelectionMode = "ExplicitLevels"  ; manual contour levels
 gsn_define_colormap(wks,"wxpEnIR")
 res@cnFillOn             = True            ; turn on color
 res@cnLevels = (/-.091,-.081,-.071,-.061,-.051,-.041,-.031,-.021,-.011,-.001,.001,.011,.021,.031,.041,.051,.061,.071,.081,.091/)
 res@cnLineDashPatterns=(/2,2,2,2,2,2,2,2,2,0,0,0,0,0,0,0,0,0/)
 res@cnFillColors         = (/-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,52,52,52,52,72,72,72,72,72/)
 res@tmXBOn               = False           ; eliminate bottom labels
 res@tmXTOn               = False           ; eliminate bottom labels
 res@tmYROn               = False           ; eliminate bottom labels
 res@cnLinesOn            = True           ; no contour lines
 res@cnMonoLineDashPattern = False
 res@cnLineLabelsOn = True
 res@gsnCenterString   = "Unrot EOF1 %Var=" + sprintf("%4.1f", eof@pcvar(0))
 res@gsnStringFontHeightF = .04
 res@gsnLeftString = ""
 res@gsnRightString = ""
 res@mpLandFillColor="black"
 res@mpOceanFillColor=-1 
 res@cnInfoLabelOn   = False   
 res@cnLineLabelPlacementMode = "computed"
 res@lbLabelBarOn = False
 res@pmLabelBarDisplayMode = False
 plot(0) = gsn_csm_contour_map_ce(wks,eof(0,:,:), res)

 delete(res@gsnCenterString)
 res@gsnCenterString   = "Unrot EOF2 %Var=" + sprintf("%4.1f", eof@pcvar(1))
 plot(3) = gsn_csm_contour_map_ce(wks,eof(1,:,:), res)

 delete(res@gsnCenterString)
 res@gsnCenterString   = "Unrot EOF3 %Var=" + sprintf("%4.1f", eof@pcvar(2))
 plot(6) = gsn_csm_contour_map_ce(wks,eof(2,:,:), res)

 delete(res@gsnCenterString)
 res@gsnCenterString   = "Unrot EOF4 %Var=" + sprintf("%4.1f", eof@pcvar(3))
 plot(9) = gsn_csm_contour_map_ce(wks,eof(3,:,:), res)

  res1=True
  res1@gsnDraw=False
  res1@gsnFrame=False
  res1@gsnStringFontHeightF = .034
  res1@gsnCenterString = "EOF1 Series"          ; yaxis
  res1@tiXAxisString = "Frequency (cycles/yr)"  ; xaxis
  res1@tiYAxisString = "Variance"                  ; yaxis

  res1@trYLog = False
  res1@trXLog = True                 ; log scaling
  res1@trXMinF  = 0.01                ; min value on y-axis
  res1@trXMaxF  =  0.3                ; max value on y-axis
  res1@vpHeightF = 0.35
  res1@trYMinF  = min(splt1(0,:))                ; min value on y-axis
  res1@trYMaxF  = max(splt1(0,:))                ; max value on y-axis
  res1@tmXBLabelsOn   = "True"
  res1@tmXBMode = "Explicit"
  res1@tmXBValues =  (/0.01, 0.02, 0.03, 0.05,0.1, 0.2, 0.3/)
  res1@tmXBLabels = (/"0.01","0.02","0.03","0.05","0.1","0.2","0.3"/)
  res1@xyLineColors        = (/"foreground","black","black","black"/)
  res1@xyDashPatterns      = (/0,0,1,1/)           ; Dash patterns
  res1@tmYROn="False"
  res1@tmXTOn="False"
  plot(1) = gsn_csm_xy(wks,spec1@frq,splt1,res1)

  res1@gsnCenterString = "EOF2 Series"          ; yaxis
  res1@trYMinF  = min(splt2(0,:))                ; min value on y-axis
  res1@trYMaxF  = max(splt2(0,:))                ; max value on y-axis
  plot(4) = gsn_csm_xy(wks,spec2@frq,splt2,res1)

  res1@gsnCenterString = "EOF3 Series"          ; yaxis
  res1@trYMinF  = min(splt3(0,:))                ; min value on y-axis
  res1@trYMaxF  = max(splt3(0,:))                ; max value on y-axis
  plot(7) = gsn_csm_xy(wks,spec3@frq,splt3,res1)

  res1@gsnCenterString = "EOF4 Series"          ; yaxis
  res1@trYMinF  = min(splt4(0,:))                ; min value on y-axis
  res1@trYMaxF  = max(splt4(0,:))                ; max value on y-axis
  plot(10) = gsn_csm_xy(wks,spec4@frq,splt4,res1)

  res1@gsnCenterString = "EOF1 Series"          ; yaxis
  res1@trYMinF  = min(splt1(0,:))                ; min value on y-axis
  res1@trYMaxF  = max(splt1(0,:))                ; max value on y-axis
  res1@trYLog = True
  res1@trXLog = False                 ; log scaling
  res1@tmXBValues =  (/0.01,0.05,0.1, 0.15,0.2,0.25,0.3/)
  res1@tmXBLabels = (/"0.01","0.05","0.1","0.15","0.2","0.25","0.3"/)
  plot(2) = gsn_csm_xy(wks,spec1@frq,splt1,res1)

  res1@gsnCenterString = "EOF2 Series"          ; yaxis
  res1@trYMinF  = min(splt2(0,:))                ; min value on y-axis
  res1@trYMaxF  = max(splt2(0,:))                ; max value on y-axis
  plot(5) = gsn_csm_xy(wks,spec2@frq,splt2,res1)

  res1@gsnCenterString = "EOF3 Series"          ; yaxis
  res1@trYMinF  = min(splt3(0,:))                ; min value on y-axis
  res1@trYMaxF  = max(splt3(0,:))                ; max value on y-axis
  plot(8) = gsn_csm_xy(wks,spec3@frq,splt3,res1)

  res1@gsnCenterString = "EOF4 Series"          ; yaxis
  res1@trYMinF  = min(splt4(0,:))                ; min value on y-axis
  res1@trYMaxF  = max(splt4(0,:))                ; max value on y-axis
  plot(11) = gsn_csm_xy(wks,spec4@frq,splt4,res1)

;******************************************
  pres             = True                 ; mod panel plot
  pres@gsnDraw=True
  pres@gsnFrame=True
  pres@gsnMaximize = True                        ; blow up plot
  pres@gsnPaperWidth=8.5
  pres@gsnPaperHeight=11.0
  pres@txString   = "SST Surrounding the Amazon (IntVeg 400Yrs)"
  pres@txFontHeightF = .012
  pres@gsnPaperOrientation="portrait"
  pres@gsnPaperMargin=0.1
  pres@gsnPanelLeft = 0.04
  pres@gsnPanelRight = 0.96
  pres@gsnPanelBottom = .05
  pres@gsnPanelTop = .95
  pres@gsnPanelYWhiteSpacePercent = 2
  pres@gsnPanelXWhiteSpacePercent = 5
  gsn_panel(wks,plot,(/4,3/),pres)        ; create panel plot

end