undef("taylor_diagram_stats")
function taylor_diagram_stats(x[*][*]:numeric, y[*][*]:numeric, weight:numeric \
                             , caseLabel[1]:string, varLabel[1]:string, opt[1]:logical)
;
; Compute Taylor statistics
; Arguments:
;    x        - reference array
;    y        - comparison array
;    weight   - array with weights: 1 means no area averaging; 
;                                [*] and [*][*]; most commonly cosine(latitude)
;    caseLable- identifier for (possible) later use in plotting
;    varLable - identifier for (possible) later use in plotting
;    opt      - optional: currently not used
;
local dimx, dimy, dimw, rankw, wgt2d, sumw, x_avg, y_avg, x_var, y_var, ratio, cc, BIAS, bias, tay_stats
begin
  dimx  = dimsizes(x)
  dimy  = dimsizes(y)
  dimw  = dimsizes(weight) 
  rankw = dimsizes(dimw)

; error checking

  if (all(dimx.ne.dimy)) then
      print("taylor_diagram_stats: Fatal:x and y dimensions are not the same size")
      print("                      x="+dimx+"   y="+dimy)
      exit
  end if

  if (rankw.gt.2) then
      print("taylor_diagram_stats: Fatal: weight can be a scalar [1], [*] or [*][*]")
      print("                      rankw: "+rankw)
      exit
  end if

  if (rankw.eq.2 .and. .not.all(dimx.eq.dimw)) then
      print("taylor_diagram_stats: Fatal: weight[*][*] must have the same dimensions as x[*][*]")
      print("                      dimx: "+dimx+"    dimw="+dimw)
      exit
  end if

  if (isatt(x,"units") .and. isatt(y,"units") .and. x@units.ne.y@units) then
      print("taylor_diagram_stats: Warning: units mismatch")
      print("                      x@units: "+x@units )
      print("                      y@units: "+y@units )
  end

; make generic weight array (convenience)

  if (rankw.eq.1) then 
      if (dimw.eq.1) then                     ; must be scalar
          wgt2d = conform(x, weight, -1)      ; expand scalar
      else      
          wgt2d = conform(x, weight, 0)       ; [*] ==>[*][*]
      end if
  else
      wgt2d = weight
  end if

  sumw  = sum(wgt2d)                           ; sum of weights (scalar)

; compute stats

  x_avg = wgt_areaave2(x, wgt2d,  0)           ; weighted areal means (scalar) 
  y_avg = wgt_areaave2(y, wgt2d,  0)
                                             
  x_var  = sum(wgt2d*(x - x_avg)^2)/sumw       ; wgted areal variances
  y_var  = sum(wgt2d*(y - y_avg)^2)/sumw

  ratio  = sqrt( y_var / x_var )               ; variance ratio (sqrt)
  cc     = pattern_cor(x, y, wgt2d, 0)         ; centered, weighted pattern cor
  BIAS   = y_avg - x_avg                       ; raw bias 

  tay_stats = new( 3, "double", getVarFillValue(x)) 

  tay_stats(0) = ratio 
  tay_stats(1) = cc 
  if (x_ref.ne.0) then
      tay_stats(2) = (BIAS/x_avg)*100 
  end if

  tay_stat!0 = "taylor"
  tay_stats@long_name = "Taylor Diagram Stats: ratio (0), cc (1), % bias (2)"
  if (isatt(x,"long_name")) then
      tay_stats@x_long_name = x@long_name
  end if
  if (isatt(y,"long_name")) then
      tay_stats@y_long_name = y@long_name
  end if
  tay_stats@caseLabel = caseLabel
  tay_stats@varLabel  = varLabel
  return(tay_stats)

end