PRO cleandata,  model=model

;; removes spikes from data, replaced with continuum

;; use:  cleandata
;;       cleandata,/model

; COMMON fileblock
COMMON fitsblock
COMMON procblock

; Define processed spectra
specf = spec
specr = spec

; Compute mean spectrum
specm = fltarr(nwave)
FOR iw = 0,nwave-1 DO specm(iw) = mean(specf(iw,*))

; Bale out here if this is a models and cleaning was not requried
if (keyword_set(model)) then return


if (n_elements(nwomit) EQ 0) THEN nwomit = 0

FOR iw = nwomit,nwave-1 DO BEGIN

; Replace bad pixels by continuum
  FOR it = 0,nspec-1 DO BEGIN
   IF specf(iw,it) GT 1.5 THEN specf(iw,it) = 1.0
   IF specf(iw,it) LT 0.3 THEN specf(iw,it) = 1.0
  ENDFOR

; Compute mean spectrum
  specm(iw) = total(specf(iw,*))/nspec

; Compute residual spectrum (clean / mean)
  FOR it = 0,nspec-1 DO BEGIN
   specr(iw,it) = specf(iw,it)/specm(iw)
  ENDFOR

ENDFOR


; Replace omitted spectrum by continuum
FOR iw = 0,nwomit-1 DO BEGIN
  FOR it = 0,nspec-1 DO BEGIN
    specr(iw,it) = 1.0
  ENDFOR
ENDFOR



END