; ================================================
; Lauren 24 July 2012
; test wgt_areaave with zeros or missing values
; ================================================

begin

latt = (/-90,-45,0,45,90/) ; dummy latitudes
rad  = 4.*atan(1.)/180.    ;
cost = cos(latt*rad)       ; cosine weights
;print(cost)

x       = new((/5,2/),float)  ; dummy fields
x2      = new((/5,2/),float)  ; dummy fields
x(:,0)  = (/0,1,1,1,0/)       ; with Zeros
x(:,1)  = (/0,2,3,0,0/)
x2(:,0) = (/-99,1,1,1,-99/)   ; with FillValue
x2(:,1) = (/-99,2,3,-99,-99/)

print(wgt_areaave(x,cost,1.0,0))
;(0)     1.414214

x2@_FillValue = -99
print(wgt_areaave(x2,cost,1.0,0))
;(0)     1.656854                 ; DIFFERS FROM PREVIOUS CALC ???

; Try with x(:,:) having zero as FillValue
x@_FillValue = 0
print(wgt_areaave(x,cost,1.0,0))
;(0)     1.656854                 ; Same as x2 but why should they differ at all ?
; Somehow handles missing values differently but why should it ?
; But zero might not always be a fill value ...

; Use calc as stated on http://www.ncl.ucar.edu/Document/Functions/Built-in/wgt_areaave.shtml
; x_ave = sum(x*wgty*wgtx)/sum(wgty*wgtx)

print(((sum(x(:,0)*cost)/sum(cost))+(sum(x(:,1)*cost)/sum(cost)))/2)
;(0)     1.414214
print(((sum(x2(:,0)*cost)/sum(cost))+(sum(x2(:,1)*cost)/sum(cost)))/2)
;(0)     1.414214

print(sum(x2(:,0)*cost)/sum(cost))
;(0)      1
print(sum(x2(:,1)*cost)/sum(cost))
;(0)     1.828427

print(sum(x(:,0)*cost)/sum(cost))
;(0)      1
print(sum(x(:,1)*cost)/sum(cost))
;(0)     1.828427

; This calc shows that they are the same ?

end