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("/Users/notaro/cru/Prc23_19001998.nc","r")   ; open data file

 lat = a->lat
 lon = a->lon
 times=ispan(1900,1998,1)    ; create time array

 prcann = a->PrcAnn ; read in precip as mm

 prcann!0="time"    ; assign attributes
 prcann!1="lat"
 prcann!2="lon"
 prcann&time=times
 prcann&lat=lat
 prcann&lon=lon
 prcann@_FillValue=-999.    ; missing value is -999.

pcp=new((/72,96,99/),float)
pcp=prcann(lat|:,lon|:,time|:)   ; reorder dimensions
pcp!0="lat"
pcp!1="lon"
pcp!2="time"
pcp&lat=lat
pcp&lon=lon
pcp&time=times
pcp@_FillValue=-999.

pcp=wgt_runave(pcp,(/.25,.5,.25/),-1)  ; here, I temporally smoothed or filtered the data twice
pcp=wgt_runave(pcp,(/.25,.5,.25/),-1)

;;;;;;;;

var = pcp(43:63,0:95,0:98)    ; select area of 20-70N only
var@_FillValue = -999.
var!0="lat"
var!1="lon"
var!2="time"
var&lat=lat(43:63)
var&lon=lon(0:95)
var&time=times

var=dim_standardize(var,0)
do nl=0,20                      ; wgt by sqrt(cos(latitude))
  var(nl,:,:)  = var(nl,:,:)*sqrt(cos(var&lat(nl)*0.0174532))  
end do

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

eof=eofcov_pcmsg_Wrap(var,5,80.)  ; calculate eof
eof@_FillValue=-999.
eof_ts = eofcov_ts_Wrap(var,eof) ; calculate time series of eof
eof_ts@_FillValue=-999.
print(max(eof))
print(min(eof))
print(max(eof_ts))
print(min(eof_ts))

do it=0,98    ; print out eof time series
  print(sprinti("%4i",it+1900)+" "+sprintf("%9.4f",eof_ts(0,it))+" "+sprintf("%9.4f",eof_ts(1,it))+" "+ \
        sprintf("%9.4f",eof_ts(2,it))+" "+sprintf("%9.4f",eof_ts(3,it))+" "+sprintf("%9.4f",eof_ts(4,it)))
end do

eof_ts=wgt_runave(eof_ts,(/1.,1.,1./),-1)  ; smooth time series before plotting
;;;;;;;;;;;;;;;;;

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

 res                     = True            ; plot mods desired
 res@gsnDraw              = False          ; don't draw
 res@gsnFrame             = False           ; don't advance frame
 res@mpShapeMode  = "FreeAspect"
 res@vpWidthF      = 0.9
 res@vpHeightF     = 0.2
 res@cnRasterModeOn = True
 res@gsnAddCyclic                = False         ; regional data
 res@mpMinLatF    = 20.0                ; only plot 20N to 90N
 res@mpMaxLatF    = 70.0
 res@mpOutlineBoundarySets = "AllBoundaries"
 res@cnLevelSelectionMode = "ExplicitLevels"  ; manual contour levels
 gsn_define_colormap(wks,"rainbow+white")
 res@cnFillOn             = True            ; turn on color
 res@cnLevels = (/-.05,-.001,.001,.05/)
 res@cnFillColors         = (/237,192,238,140,50/)
 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 = False
 res@gsnCenterString   = "Unrot EOF1 %Var=" + sprintf("%4.1f", eof@pcvar(0))
 res@gsnStringFontHeightF = .04
 res@gsnLeftString = ""
 res@gsnRightString = ""
 res@cnInfoLabelOn   = False   
 res@cnLineLabelPlacementMode = "computed"
 res@lbLabelBarOn = True
 res@cnFillDrawOrder      = "PreDraw"
 res@cnLineDrawOrder = "PreDraw"
 res@mpOceanFillColor     = "white"
 res@mpLandFillColor = -1
 res@mpInlandWaterFillColor = "white"
 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(1) = gsn_csm_contour_map_ce(wks,eof(1,:,:), res)

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

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

;******************************************
  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@gsnPaperOrientation="portrait"
  pres@gsnPaperMargin=0.02
  pres@gsnPanelLeft = 0.02
  pres@gsnPanelRight = 0.98
  pres@gsnPanelBottom = .02
  pres@gsnPanelTop = .98
  pres@gsnPanelYWhiteSpacePercent = 2
  pres@gsnPanelXWhiteSpacePercent = 5
  gsn_panel(wks,plot,(/4,1/),pres)        ; create panel plot


  resxy          = True                              ; xy plot mods desired
  resxy@gsnDraw  = False                             ; don't draw yet
  resxy@gsnFrame = False                             ; don't advance frame yet
  resxy@vpWidthF = 0.80                              ; plot width
  resxy@vpHeightF= 0.50                              ; plot height
  resxy@tiYAxisString=""
  do ne=0,3
    resxy@gsnLeftString    = "Covariance"           ; titles
    resxy@gsnCenterString  = "EOF "+(ne+1)
    resxy@gsnRightString   = "%Var=" + sprintf("%4.1f", eof@pcvar(ne))
    plot(ne) = gsn_csm_xy(wks,times,eof_ts(ne,:), resxy)
  end do
  gsn_panel(wks,plot(0:3),(/4,1/),False)          ; draw: 1-across, 3-down

end