function gridit,nx,ny,xlon,ylat,statlon,statlat,statval,statweight,$
  nstat=nstat,statflag=statflag,minv=minv,maxv=maxv
;
; **** THIS VERSION HAS A WEIGHTING FOR EACH STATION, THAT IS USED WHEN
; AVERAGING A NUMBER OF STATIONS INTO A SINGLE BOX ****
; ### AND GIVES MIN AND MAX IN EACH CELL ###
; Returns a gridded result from a list of station values
; where the grid is defined by xlon and ylat (centre of boxes)
; and where the stations are defined by statlon and statlat.
; If statval is 2D, its first dimension represents a sequence of fields to
; create.  The number of stations contributing to each grid-box can also
; be returned if required.  If more than one station contributes to a box,
; then the values are simply averaged.  If no stations in a box, NaN is
; returned.  Missing data must be stored as NaN to be accounted for.
; If statflag is given, then any station with a non-zero flag is ignored.
;
;-------------------------------------------------------------------
;
; Check dimensions, and ensure a 2D input array
;
n1=n_elements(statlon)
n2=n_elements(statlat)
statsize=size(statval)
ndim=statsize(0)
if (ndim lt 1) and (ndim gt 2) then message,'statval must be 1D or 2D'
if ndim eq 1 then begin
  print,'1D input'
  nt=1
  ns=statsize(1)
  stat2d=reform(statval,nt,ns)
  print,nt,ns
  help,stat2d
endif else begin
  print,'2D input'
  nt=statsize(1)
  ns=statsize(2)
  stat2d=statval
  print,nt,ns
  help,stat2d
endelse
if (ns ne n1) or (ns ne n2) then message,'Incompatible station arrays'
if n_elements(statweight) eq 0 then statweight=fltarr(ns)+1.
;print,statweight
;
; Define arrays
;
fdnew=dblarr(nt,nx,ny)
fdmin=fltarr(nt,nx,ny)  &  fdmin(*,*,*)=1.E10
fdmax=fltarr(nt,nx,ny)  &  fdmax(*,*,*)=-1.E10
nstat=fltarr(nt,nx,ny)
wstat=dblarr(nt,nx,ny)
;
; Compute grid-box edges
;
dx=xlon(1)-xlon(0)
xlonedge=[xlon(0)-0.5*dx,xlon+0.5*dx]
;print,'xlonedge'
;print,xlonedge
dy=ylat(1)-ylat(0)
ylatedge=[ylat(0)-0.5*dy,ylat+0.5*dy]
;print,'ylatedge'
;print,ylatedge
;
; Start filling the arrays
;
for istat = 0 , ns-1 do begin
;  print,'Station:',istat
  ;
  ; Ignore it if requested to do so
  ;
  iwant=1
  if n_elements(statflag) gt 0 then begin
    if statflag(istat) ne 0 then iwant=0
  endif
  if iwant eq 1 then begin
;    print,'Using it'
    ;
    ; Locate it's grid-box
    ;
    xlist=where( (xlonedge(0:nx-1) lt statlon(istat)) and $
                 (xlonedge(1:nx) ge statlon(istat)),nbox1)
    ylist=where( (ylatedge(0:ny-1) ge statlat(istat)) and $
                 (ylatedge(1:ny) lt statlat(istat)),nbox2)
    if (nbox1 gt 1) or (nbox2 gt 1) then message,'Oops!'
    if (nbox1 eq 1) and (nbox2 eq 1) then begin
;      print,'Location:',xlon(xlist),ylat(ylist)
      ;
      ; Check each time value separately
      ;
      for itime = 0 , nt-1 do begin
        if finite(stat2d(itime,istat)) then begin
          fdnew(itime,xlist,ylist)=fdnew(itime,xlist,ylist)+$
                                   stat2d(itime,istat)*statweight(istat)
          nstat(itime,xlist,ylist)=nstat(itime,xlist,ylist)+1.
          wstat(itime,xlist,ylist)=wstat(itime,xlist,ylist)+statweight(istat)
          fdmin(itime,xlist,ylist)=min([fdmin(itime,xlist,ylist),$
                                        stat2d(itime,istat)])
          fdmax(itime,xlist,ylist)=max([fdmax(itime,xlist,ylist),$
                                        stat2d(itime,istat)])
        endif
      endfor
      ;
    endif
    ;
  endif
  ;
endfor
;
; Average the totals
;
for i = 0 , nx-1 do begin
  for j = 0 , ny-1 do begin
    for itime = 0 , nt-1 do begin
      if nstat(itime,i,j) gt 0. then begin
        fdnew(itime,i,j)=fdnew(itime,i,j)/wstat(itime,i,j)
        if nstat(itime,i,j) eq 1. then begin
          fdmin(itime,i,j)=!values.f_nan
          fdmax(itime,i,j)=!values.f_nan
        endif
      endif else begin
        fdnew(itime,i,j)=!values.f_nan
        fdmin(itime,i,j)=!values.f_nan
        fdmax(itime,i,j)=!values.f_nan
      endelse
    endfor
  endfor
endfor
;
nstat=reform(nstat)
minv=reform(fdmin)
maxv=reform(fdmax)
return,float(reform(fdnew))
;
end