fn  = "zo_hist_1850-2005_ens_1-100.nc" ; define filename 

   in  = addfile(fn,"r")                                 ; open netcdf file
                           ; zo:[time|1872] x [ens|100] x [depth|1] x [lat|45] x [lon|90]
   x   = in->zo(:,:,0,:,:) 
   printVarSummary(x)      ; (time,ens,lat,lon) => (1872,100,45,90)
                           ;   0    1   2   3
   dimx = dimsizes(x)
   ntim = dimx(0)
;; nens = dimx(1)
;; nlat = dimx(2)
;; mlon = dimx(3)

;******************************************************************
;  Compute the ensemble statistics at each time step and grid point
;******************************************************************

   xStat = dim_stat4_n(x,1)   ; (4,ntim,nlat,mlon)
   printVarSummary(xStat)
   print("---")

   xAvg = xStat(0,:,:,:)
   xStd = xStat(1,:,:,:)
   xSkw = xStat(2,:,:,:)
   xKur = xStat(3,:,:,:)
                                 ; meta data
   xAvg@long_name = "Ens Wave Hgt: Avg"
   xAvg@units     = x@units
   xStd@long_name = "Ens Wave Hgt: Std Dev"
   xStd@units     = x@units
   xSkw@long_name = "Ens Wave Hgt: Skewness"
   xKur@long_name = "Ens Wave Hgt: Kurtosis"

   copy_VarCoords(x(:,0,:,:), xAvg) 
   copy_VarCoords(x(:,0,:,:), xStd) 
   copy_VarCoords(x(:,0,:,:), xSkw) 
   copy_VarCoords(x(:,0,:,:), xKur) 
                                 ; variable overview
   printVarSummary(xAvg)
   printMinMax(xAvg,0)
   print("---")
   printVarSummary(xStd)
   printMinMax(xStd,0)
   print("---")
   printVarSummary(xSkw)
   printMinMax(xSkw,0)
   print("---")
   printVarSummary(xKur)
   printMinMax(xKur,0)
   print("---")
     
;***************************************************************
;--- PLOTS
;***************************************************************

   wks  = gsn_open_wks ("png","EnsWaveStat")                ; send graphics to PNG file

   LAT = 45    ; arbitrary
   LON = 170

; ***********************************************
; create new 'time' array for use on the plot
; ***********************************************
  time    = xAvg&time
  Time    = cd_calendar(time,0)
  yyyy    = Time(:,0)
  mm      = Time(:,1)
  timePlt = yyyy + (mm-1)/12.
  ntStrt  = ind(timePlt.eq.1850.0)
  ntLast  = ntim-1

  res             = True                   ; plot mods desired
  res@gsnMaximize = False                  ; don't advance frame yet      

  res@vpHeightF   = 0.4                    ; change aspect ratio of plot
  res@vpWidthF    = 0.8                  
  res@vpXF        = 0.15                   ; start plot at x ndc coord 
  
  res@trXMinF  =  1840  ; yyyy(ntStrt)=1850
  res@trXMaxF  =  2010  ; yyyy(ntLast)=2005

  res@tiYAxisString    = ""        ; no  y-axis label [if desired]      
  res@tiMainString     = xAvg@long_name+": ["+LAT+","+LON+"]"
  plot = gsn_csm_xy (wks,timePlt(ntStrt:),xAvg(ntStrt:,{LAT},{LON}),res)

  res@tiMainString     = xStd@long_name+": ["+LAT+","+LON+"]"
  plot = gsn_csm_xy (wks,timePlt(ntStrt:),xStd(ntStrt:,{LAT},{LON}),res)

  res@gsnYRefLine      = 0.0               ; create a reference line 
  res@gsnYRefLineColor = "blue"

  res@tiMainString     = xSkw@long_name+": ["+LAT+","+LON+"]"
  plot = gsn_csm_xy (wks,timePlt(ntStrt:),xSkw(ntStrt:,{LAT},{LON}),res)

  res@tiMainString     = xKur@long_name+": ["+LAT+","+LON+"]"
  plot = gsn_csm_xy (wks,timePlt(ntStrt:),xKur(ntStrt:,{LAT},{LON}),res)