Average over a 3-D array, clipping unusual deviants.
   Calculate the average value of an array, or calculate the average
   value over one dimension of an array as a function of all the other
   CCD data processing
   array = 3-D input array.  May be any type except string.  This cube of
              data is MODIFIED by this routine.  What's left on return may
              or may not be useful.
   BAD       - Bad flag cube that matches the input 3-D array.  0 means the
                  value is not known to be bad.  Anything else and the pixel
                  should not be used in the average.  If not provided, the
                  bad flag array is initially set to all good.  This cube
                  is modified to add any additional pixels that were flagged
                  in the process of creating the final image.
   CLEANMEAN - array (or name of fits file) that should be used to filter
                  outliers.  Normally, this image is computed from the median
                  of the image stack.  But, if you happen to have better
                  information provide it via this keyword.  The image size
                  must be consistent with the cropping region.  If the images
                  are to be scaled the mean of this image is not used in any
                  way so you are free to provide a normalized image if you
                  wish.  This feature is typically used when combining
                  flat field images.

   CLEANSIG - array (or name of fits file) that should be used as the sigma
                  image array.  By default, this is computed from a clipped
                  sqrt(cleanmean) which is normally the median of the stack.
                  If you provide this array you must also provide CLEANMEAN
                  or this keyword will be ignored.  Note that the units of
                  CLEANSIG and CLEANMEAN should match.  In other words, if
                  you do provide a normalized clean mean, then CLEANSIG
                  should bear the stamp of that normalization as well.  This
                  feature is typically used when combining flat field images.

   JUSTMEDIAN - Flag, if set will stop processing after the initial median
                  combination of the image cube and this will be the final
                  answer.  However, this is only honored if the NORMALIZE and
                  PRESCALE flags are both set.  This should only be used if
                  the intrinsic noise level of the result can never drop
                  below 1 DN.

   SCALE - 4 element vector which, if provide, defines the region of the
           array dimensions that are used to scale the mean
           of the arrays before combining.  If combined in this
           manner, the arrays are combined weighted by the means.
           If the region contains more than 20000 pixels, then only this
           many randomly chosen pixels will be used to derive the robust
           mean of the scaling region.

   NORMALIZE - Flag, if set and SCALE used, leaves the output average
                 normalized by the SCALE region.

   SILENT - Flat, if set will suppress all messages to screen.

   THRESH - Threshold, in units of a standard deviation, to flag and thus
               remove outliers.  The default is 3.0 sigma.

   NOISEMIN - Noise floor to use when combining images.  Typically this would
                be set to the readout noise so that data with unusually low
                noise will be treated more realistically. (default=1)

   NONEWBADFLAGS - Flag, if set will suppress searching for new bad pixels
                in the stack of images.  If you are providing a bad mask
                then it will be used as is.  If you are not providing a bad
                mask this will cause avgclip to return a straight average
                of the stack.

   NTOPCLIP - (default=0) How many of the brightest pixels to removed
                 before the robust statistics.  (see

   average - 2-D array that is the robust average of the stack.
   1992 Dec 30, Marc W. Buie, Lowell Observatory, cloned from AVG
      and added average sigma clipping.
   95/03/10, MWB, extensive re-write to optimize.
   97/06/19, MWB, added SILENT keyword
   2000/09/22, MWB, added THRESH keyword
   2006/03/18, MWB, fixed problem with bias frame averaging when using
                      overscan correction.
   2006/04/04, MWB, fixed problem with flat averaging caused by bias fix

   2006/06/29 Peter L. Collins, Lowell Observatory experiments on noise array
   2006/06/30 PLC, added NOISEMETHOD and READNOISE keywords.
   2006/07/14, MWB, final checkout and reincorporation into library.
   2006/07/25, MWB, slight change to the normalization scaling process.  Also
                  added CLEANMEAN and CLEANSIG keywords.
   2006/07/27, MWB, added JUSTMEDIAN keyword
   2007/08/08, MWB, added BAD keyword
   2007/09/04, MWB, fixed problem with CLEANMEAN and CLEANSIG
   2007/09/07, MWB, fixed more errors with handling of incoming bad arrays
                      also added the NONEWBADFLAGS keyword and pass BAD
                      through to medarr_mwb
   2011/12/16, MWB, added NTOPCLIP keyword
   2012/03/28, MWB, removed /nozero on allocation of output image, this
                      would cause problems for pixels that were all bad and
                      return completely unpredictible values.