function n802pp1,file,data=data,txt=txt,bin=bin,array=array,bi2=bi2,nored=nored,keepred=keepred,noisy=noisy,nofilename=nofilename,lun=lun,lbexp=lbexp

; Purpose: convert ERA (reduced) gaussian n80 data into pp-format
; Category: io data-conversion
;
; Note: distinguishes between daily fields (reads time info from file name)
;       and monthlies (will adjust time) according to file name ending in "avg"
; Note: expects input files to be named as for the output of "wmc_grib.pl" in
;       order to deduce time and variable info (probably wmc_grib1.pl by now,
;     	and don't forget EDITION=1)
;
; Options:
;	data=data - hand the routine the data, in which case it won't read it from file
;	/txt - input file is in text format
;	/bin - input file is in binary format (native 4-byte reals) (default)
;	/bi2 - input is in wmc_grib.pl "bin2" format
;	/array - return data interpolated, not wrapped up as pp-field
;        array=2 return original data (not interpolated onto rectangular grid)
;	/nored - input is 320*160, *not* reduced-grid
;       /keepred - input is reduced and we *dont* want it interpolated onto regular grid
;	/nofilename - the filename is not in the required format; don't attempt to deduce the
;	              field names/dates from it. Using this option probably means that
;		      you have forgotten to pass the files through wmc_grib1.pl and you are
;		      trying to read raw grib which will fail...
;       lun=lun - variable to keep the lun in
;	lbexp=lbexp - set lbexp to something other than ERA, eg 'eop' or number
@comm_error

if (keyword_set(noisy)) then print,'File: '+file

; Read in the numbers of points around each longitude circle (this is ERA-n80
; specific) and the latitudes of the rows (this too)
nlon=reform(readfromfile('$WMCDATA/n80.nlon.ERA'),160)
lats=reform(readfromfile('$WMCDATA/n80.lats.ERA'),160)

; Note the total number of points we expect
if (keyword_Set(nored)) then np=320l*160l else np=total(nlon)

; Read in the data itself (unless data supplied in call)
if (n_elements(data) eq 0) then begin
  ; May need to gunzip
  if (strpos(file,'.gz') ne -1) then begin
    file1=file
    spawn,'echo '+file+'|sed "s/.*\///"',file
    file='/tmp/'+file(0)
    spawn,'gunzip -c '+file1+'>'+file
  endif

  if (keyword_set(txt)) then $
    data=readfromfile(file) $
  else if (keyword_set(bi2)) then begin
    openr,lun,file,/get
    np1=1l & ref=0. & scale=1.
    readu,lun,np1,scale,ref
    if (keyword_set(noisy)) then print,'np, np1, scale, ref: ',np,np1,scale,ref
    if (np ne np1) then message,'whoops! np != np1: '+shtstr(np)+' '+shtstr(np1)
    data=bytarr(np*2) & readu,lun,data
    i=lindgen(np)
    data=(data(i*2)*256.+data(i*2+1))*scale+ref
    free_lun,lun
  endif else begin
    if (lun ne 0) then $
      message,/cont,'trying to keep reading from open lun '+shtstr(lun) $
    else $
      openr,lun,file,/get 
    junk=fstat(lun) 
    if (junk.size/4 ne np) then message,'File size/4 != total(nlon); maybe set /nored?: '+shtstr(junk.size/4)+' '+shtstr(np),/cont
    data=fltarr(np) & readu,lun,data
    if (not eof(lun)) then $
      message,/cont,'not at eof: will try to read another field' $
    else begin
      free_lun,lun
      lun=0
    endelse
  endelse

  if (n_elements(file1) ne 0) then begin
    spawn,'rm '+file,junk
    file=file1
 endif

endif

if (not keyword_set(keepred)) then begin
 
; J counts the points we've processed; data1 is the output array
  j=0
  data1=fltarr(320,160)

; Now loop over each row, putting data from "data" into "data1"
  if (keyword_Set(nored)) then $
    data1=reform(data,320,160) $
  else for i=0,n_elements(nlon)-1 do begin
    ind=indgen(nlon(i))+j
; Interpolate onto a 320-point grid, which is the resolution at the equator
; Note we wrap "data" around the longitude circle to allow interpolation
; past the "last" point back to the "first" (0E) point.
    d=interpolate([data(ind),data(ind(0))],indgen(320)*float(nlon(i))/320)
    data1(*,i)=d
    j=j+nlon(i)
  endfor

endif

; --------------------------------------------------------
; Past this point we become UKMO (pp-field) specific
; We may want to stop here and just return the data1 array
; --------------------------------------------------------
if (keyword_set(array)) then begin
  if (array eq 1) then return,data1			; Make sure /keepred not set!
  if (array eq 2) then return,data			; Would be sensible to set /keepred in this case
endif

; Wrap up "data1" as a pp-field. Note that the y-coords (lats) are supplied
; vi "y=lats" since they are irregular (gaussian)
if (keyword_set(keepred)) then begin
  y=fltarr(4,n_elements(lats))
  y(0,*)=lats
  y(1,*)=0
  y(2,*)=float(nlon)
  y(3,*)=float(nlon)
  pp=makeppfield(data,y=y,noisy=keyword_set(noisy))
  pp.lbhem=0
endif else begin
  dx=360./320.
  pp=makeppfield(data1,zx=-dx,dx=dx,y=lats,noisy=keyword_set(noisy))
  pp.lbhem=0
endelse

; Set up info from file name
if (not keyword_Set(nofilename)) then gofilename2header,pp,file

; Set lbexp as other than ERA, if required
if (keyword_set(lbexp)) then begin
  if (datatype(lbexp) ne 7) then pp.lbexp=lbexp $
  else pp.lbexp=pp_guess_lbexp(lbexp)
endif

; Return the field
return,pp
end

; -------------------------------------------

function n802pp,files1,nofilename=nofilename,_extra=e

@comm_error

if (not keyword_set(nofilename)) then begin
  if (n_elements(files1) eq 1) then $
    spawn,'ls '+files1,files $
  else $
    files=files1
  n=n_elements(files)
endif else begin
  files=''
  n=nofilename			; This is usually 1
endelse

lun=0
i=0
n_extra=0
while (i lt n+n_extra) do begin 
  if (i eq 0) then $
    pps=replicate(n802pp1(files(0),nofilename=nofilename,lun=lun,_extra=e),n) $
  else $
    pps(i)=n802pp1(files(i-n_extra),nofilename=nofilename,lun=lun,_extra=e)
  if (lun ne 0) then begin
    pps=[pps,pps(0)]
    n_extra=n_extra+1
  endif
  i=i+1
endwhile

return,pps
end