Index: /branches/eam_branches/ipp-20211108/ppStack/src/Makefile.am
===================================================================
--- /branches/eam_branches/ipp-20211108/ppStack/src/Makefile.am	(revision 41907)
+++ /branches/eam_branches/ipp-20211108/ppStack/src/Makefile.am	(revision 41908)
@@ -29,4 +29,6 @@
 	ppStackPrepare.c	\
 	ppStackConvolve.c	\
+	ppStackLoopByPercent.c  \
+	ppStackCombinePercent.c \
 	ppStackCombinePrepare.c	\
 	ppStackCombineInitial.c	\
Index: /branches/eam_branches/ipp-20211108/ppStack/src/ppStack.h
===================================================================
--- /branches/eam_branches/ipp-20211108/ppStack/src/ppStack.h	(revision 41907)
+++ /branches/eam_branches/ipp-20211108/ppStack/src/ppStack.h	(revision 41908)
@@ -115,4 +115,16 @@
 bool ppStackReadoutFinalThread(psThreadJob *job // Job to process
     );
+
+// Thread entry point for ppStackReadoutPercent
+bool ppStackReadoutPercentThread(psThreadJob *job);
+
+bool ppStackReadoutPercent(const pmConfig *config,
+			   pmReadout *outRO,
+			   const psArray *readouts,
+			   const psVector *mask,
+			   const psVector *weightings,
+			   const psVector *exposures,
+			   const psVector *addVariance);
+
 // Perform median stacking for background
 bool ppStackCombineBackground(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config);
Index: /branches/eam_branches/ipp-20211108/ppStack/src/ppStackCombineFinal.c
===================================================================
--- /branches/eam_branches/ipp-20211108/ppStack/src/ppStackCombineFinal.c	(revision 41907)
+++ /branches/eam_branches/ipp-20211108/ppStack/src/ppStackCombineFinal.c	(revision 41908)
@@ -2,5 +2,5 @@
 
 // This is the doomsday switch.
-// #define TESTING                         // Enable test output
+// # define TESTING                         // Enable test output
 
 //MEH -- adhoc addition to blank mask border of final stack since rejection different/none on order of KERNEL.SIZE with overlap in CombineInitial 
Index: /branches/eam_branches/ipp-20211108/ppStack/src/ppStackCombineInitial.c
===================================================================
--- /branches/eam_branches/ipp-20211108/ppStack/src/ppStackCombineInitial.c	(revision 41907)
+++ /branches/eam_branches/ipp-20211108/ppStack/src/ppStackCombineInitial.c	(revision 41908)
@@ -9,10 +9,4 @@
     psAssert(options, "Require options");
     psAssert(config, "Require configuration");
-
-/*     if (!options->convolve) { */
-/*         // No need to do initial combination when we haven't convolved */
-/*         // XXX either allocate inspect and rejected here, or do not require them downstream */
-/*         return true; */
-/*     } */
 
     psTimerStart("PPSTACK_INITIAL");
Index: /branches/eam_branches/ipp-20211108/ppStack/src/ppStackCombinePercent.c
===================================================================
--- /branches/eam_branches/ipp-20211108/ppStack/src/ppStackCombinePercent.c	(revision 41908)
+++ /branches/eam_branches/ipp-20211108/ppStack/src/ppStackCombinePercent.c	(revision 41908)
@@ -0,0 +1,56 @@
+#include "ppStack.h"
+
+// #define TESTING
+
+bool ppStackCombinePercent(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config)
+{
+    psAssert(stack, "Require stack");
+    psAssert(options, "Require options");
+    psAssert(config, "Require configuration");
+
+    psTimerStart("PPSTACK_PERCENT");
+
+    bool status;                    // Status of read
+    int numChunk;                   // Number of chunks
+    for (numChunk = 0; true; numChunk++) {
+        ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, 0);
+        if (!status) {
+            // Something went wrong
+            psError(psErrorCodeLast(), false, "Unable to read chunk %d", numChunk);
+            return false;
+        }
+        if (!thread) {
+            // Nothing more to read
+            break;
+        }
+
+        psThreadJob *job = psThreadJobAlloc("PPSTACK_PERCENT_COMBINE"); // Job to start
+        psArrayAdd(job->args, 1, thread);
+        psArrayAdd(job->args, 1, options);
+        psArrayAdd(job->args, 1, config);
+        if (!psThreadJobAddPending(job)) {
+            return false;
+        }
+    }
+
+    if (!psThreadPoolWait(true, true)) {
+        psError(psErrorCodeLast(), false, "Unable to do initial combination.");
+        return false;
+    }
+
+    ppStackMemDump("initial");
+
+#ifdef TESTING
+    ppStackWriteImage("combined_image_initial.fits", NULL, options->outRO->image, config);
+    ppStackWriteImage("combined_mask_initial.fits", NULL, options->outRO->mask, config);
+    ppStackWriteImage("combined_variance_initial.fits", NULL, options->outRO->variance, config);
+
+    pmStackVisualPlotTestImage(options->outRO->image, "combined_image_initial.fits");
+#endif
+
+    if (options->stats) {
+      // psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_INITIAL", PS_META_REPLACE, "Time to make initial stack", psTimerMark("PPSTACK_INITIAL"));
+    }
+
+    return true;
+}
Index: /branches/eam_branches/ipp-20211108/ppStack/src/ppStackCombinePrepare.c
===================================================================
--- /branches/eam_branches/ipp-20211108/ppStack/src/ppStackCombinePrepare.c	(revision 41907)
+++ /branches/eam_branches/ipp-20211108/ppStack/src/ppStackCombinePrepare.c	(revision 41908)
@@ -31,7 +31,5 @@
       options->expRO = pmReadoutAlloc(expCell); //Output readout
     }
-/*     else { */
-/*       options->expRO = NULL; */
-/*     } */
+
     int bkg_r0 = 0,bkg_c0 = 0;
     int bkg_nC = 0,bkg_nR = 0;
@@ -53,6 +51,6 @@
       
       options->bkgRO = pmReadoutAlloc(bkgCell); // BKG readout
-      //      if (!pmHDUGenerateForFPA(options->bkgRO->parent->parent->parent)) {
       options->bkgRO->parent->parent->parent->hdu = pmHDUAlloc(NULL);
+
       if (!options->bkgRO->parent->parent->parent->hdu) {
 	fprintf(stderr,"failed to generate a HDU for this thing.\n");
@@ -68,6 +66,5 @@
 									ro->parent->parent->parent->concepts);
 
-    }
-    else {
+    } else {
       options->bkgRO = NULL;
     }
@@ -97,7 +94,4 @@
         return false;
       }
-      
-
-      
     }
 
Index: /branches/eam_branches/ipp-20211108/ppStack/src/ppStackLoop.c
===================================================================
--- /branches/eam_branches/ipp-20211108/ppStack/src/ppStackLoop.c	(revision 41907)
+++ /branches/eam_branches/ipp-20211108/ppStack/src/ppStackLoop.c	(revision 41908)
@@ -31,4 +31,10 @@
     ppStackMemDump("prepare");
     if (options->quality) return true; // Can't do anything else
+
+    // if we want to skip the rejection steps, and instead use a percentile clipping:
+    if (options->clipPercent) {
+      bool result = ppStackLoopByPercent (config, options);
+      return result;
+    }
 
     // Convolve inputs
@@ -88,5 +94,5 @@
         pmCellFreeData(options->cells->data[i]);
     }
-    //MEH -- must uncomment back out -- unclear is should be moved after pixel rejection
+    // MEH -- must uncomment back out -- unclear is should be moved after pixel rejection
     psFree(stack);
 
@@ -126,9 +132,11 @@
         return false;
       }
-    }
-    else {
+    } else {
       // Since we haven't convolved, I believe we do need to normalize here.
-      //MEH -- see below for comment on ppStackCombineFinal and bscaleOffset
-      if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, true, true, true)) {
+      // MEH -- see below for comment on ppStackCombineFinal and bscaleOffset
+      // last 4 parameters are: safe, normalize, grow, bscaleoffset
+      // EAM : looks like we do NOT need to normalize here eithr because ppStackConvolve is still called and
+      // the inputs are normalized
+      if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, false, true, true)) {
         psError(psErrorCodeLast(), false, "Unable to perform final combination.");
         psFree(stack);
@@ -154,5 +162,6 @@
       return false;
     }
-    ppStackFileActivation(config,PPSTACK_FILES_BKG ,false);    
+    ppStackFileActivation(config, PPSTACK_FILES_BKG, false);    
+
     // Photometry
     psTrace("ppStack", 1, "Photometering stacked image....\n");
@@ -184,5 +193,5 @@
         return false;
     }
-    //MEH -- also must uncomment back out..
+    // MEH -- also must uncomment back out..
     psFree(stack);
     psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS"));
Index: /branches/eam_branches/ipp-20211108/ppStack/src/ppStackLoop.h
===================================================================
--- /branches/eam_branches/ipp-20211108/ppStack/src/ppStackLoop.h	(revision 41907)
+++ /branches/eam_branches/ipp-20211108/ppStack/src/ppStackLoop.h	(revision 41908)
@@ -7,4 +7,5 @@
     ppStackOptions *options             // Options for stacking
     );
+
 // Median only loop.
 bool ppStackMedianLoop(
@@ -12,4 +13,10 @@
     ppStackOptions *options             // Options for stacking
     );
+
+// stack by percentile range
+bool ppStackLoopByPercent(
+  pmConfig *config,
+  ppStackOptions *options
+  );
 
 // Setup
@@ -92,4 +99,6 @@
 
 
+bool ppStackCombinePercent(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config);
+
 bool ppStackUpdateHeader(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config);
 bool ppStackJPEGs(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config);
Index: /branches/eam_branches/ipp-20211108/ppStack/src/ppStackLoopByPercent.c
===================================================================
--- /branches/eam_branches/ipp-20211108/ppStack/src/ppStackLoopByPercent.c	(revision 41908)
+++ /branches/eam_branches/ipp-20211108/ppStack/src/ppStackLoopByPercent.c	(revision 41908)
@@ -0,0 +1,154 @@
+#include "ppStack.h"
+
+bool ppStackLoopByPercent(pmConfig *config, ppStackOptions *options)
+{
+    assert(config);
+
+    // Convolve inputs (this must be called even if the convolution is turned off because
+    // it performs the normalization of the inputs).  
+    psTrace("ppStack", 1, "Convolving inputs to target PSF....\n");
+    if (!ppStackConvolve(options, config)) {
+        psError(psErrorCodeLast(), false, "Unable to convolve images.");
+        return false;
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 2: Generate Convolutions and Save: %f sec", psTimerClear("PPSTACK_STEPS"));
+    ppStackMemDump("convolve");
+    if (options->quality) return true; // Can't do anything else
+
+    // Define threading elements
+    ppStackThreadData *stack = ppStackThreadDataSetup(options, config, true);
+    if (!stack) {
+        psError(psErrorCodeLast(), false, "Unable to initialise stack threads.");
+        return false;
+    }
+
+    // Prepare for combination (define outputs, generate output cells)
+    if (!ppStackCombinePrepare("PPSTACK.OUTPUT", "PPSTACK.OUTPUT.EXP", options->doBackground ? "PPSTACK.OUTPUT.BKGMODEL" : NULL, PPSTACK_FILES_STACK, stack, options, config)) {
+        psError(psErrorCodeLast(), false, "Unable to prepare for combination.");
+        psFree(stack);
+        return false;
+    }
+
+    // Initial combination
+    psTrace("ppStack", 1, "Initial stack of convolved images....\n");
+    if (!ppStackCombinePercent(stack, options, config)) {
+        psError(psErrorCodeLast(), false, "Unable to perform initial combination.");
+        psFree(stack);
+        return false;
+    }
+    ppStackMemDump("percent");
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 3: Make Initial Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
+
+    // Done with stack inputs for now
+    for (int i = 0; i < options->num; i++) {
+        pmCellFreeData(options->cells->data[i]);
+    }
+
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 5: Final Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
+    ppStackMemDump("final");
+
+    // Update Header
+    // Do this before performing photometry so that the cmf header gets all of the information.
+    if (!ppStackUpdateHeader(stack, options, config)) {
+        psError(psErrorCodeLast(), false, "Unable to update header.");
+        psFree(stack);
+        return false;
+    }
+
+    // Generate median background stack here.
+    if (options->doBackground && !ppStackCombineBackground(stack, options, config)) {
+      psError(psErrorCodeLast(), false, "Unable to generate median of background images.");
+      psFree(stack);
+      return false;
+    }
+    ppStackFileActivation(config, PPSTACK_FILES_BKG, false);    
+
+    // Photometry
+    psTrace("ppStack", 1, "Photometering stacked image....\n");
+    if (!ppStackPhotometry(options, config)) {
+        psError(psErrorCodeLast(), false, "Unable to perform photometry.");
+        return false;
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS"));
+    ppStackMemDump("photometry");
+
+    // Create JPEGS
+    if (!ppStackJPEGs(stack, options, config)) {
+        psError(psErrorCodeLast(), false, "Unable to make jpegs.");
+        psFree(stack);
+        return false;
+    }
+    // Assemble Stats
+    if (!ppStackStats(stack, options, config)) {
+        psError(psErrorCodeLast(), false, "Unable to assemble statistics.");
+        psFree(stack);
+        return false;
+    }
+
+   // Clean up
+    psTrace("ppStack", 2, "Cleaning up after combination....\n");
+    if (!ppStackCleanupFiles(stack, options, config, PPSTACK_FILES_STACK, PPSTACK_FILES_PHOT, true)) {
+        psError(psErrorCodeLast(), false, "Unable to clean up.");
+        psFree(stack);
+        return false;
+    }
+    psFree(stack);
+    psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS"));
+    ppStackMemDump("cleanup");
+
+    // Unconvolved stack --- it's cheap to calculate, compared to everything else!
+    // XXX unconvolved stack is currently using the convolved mask!  oops!
+    if (options->convolve) {
+        // Start threading
+        ppStackThreadData *stack = ppStackThreadDataSetup(options, config, false);
+        if (!stack) {
+            psError(psErrorCodeLast(), false, "Unable to initialise stack threads.");
+            return false;
+        }
+
+        // Prepare for combination
+        if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", NULL, PPSTACK_FILES_UNCONV,
+                                   stack, options, config)) {
+            psError(psErrorCodeLast(), false, "Unable to prepare for combination.");
+            psFree(stack);
+            return false;
+        }
+
+	// generate the unconvolved stack. NOTE: this one must be normalized since the inputs have not been
+        psTrace("ppStack", 2, "Stack of unconvolved images....\n");
+        if (!ppStackCombinePercent(stack, options, config)) {
+            psError(psErrorCodeLast(), false, "Unable to perform unconvolved combination.");
+            psFree(stack);
+            return false;
+        }
+        psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Unconvolved Stack: %f sec", psTimerClear("PPSTACK_STEPS"));
+        ppStackMemDump("unconv");
+
+	// Update Header
+	if (!ppStackUpdateHeader(stack, options, config)) {
+	    psError(psErrorCodeLast(), false, "Unable to update header.");
+	    psFree(stack);
+	    return false;
+	}
+	// Clean up unconvolved stack
+	psTrace("ppStack", 2, "Cleaning up after unconvolved stack....\n");
+	if (!ppStackCleanupFiles(stack, options, config, PPSTACK_FILES_UNCONV, PPSTACK_FILES_NONE, false)) {
+	    psError(psErrorCodeLast(), false, "Unable to clean up.");
+	    psFree(stack);
+	    return false;
+	}
+	psFree(stack);
+    }
+    psFree(options->cells); options->cells = NULL;
+
+    // Finish up
+    psTrace("ppStack", 1, "Finishing up....\n");
+    if (!ppStackFinish(options, config)) {
+        psError(psErrorCodeLast(), false, "Unable to finish up.");
+        return false;
+    }
+    ppStackMemDump("finish");
+
+    return true;
+}
+
Index: /branches/eam_branches/ipp-20211108/ppStack/src/ppStackMatch.c
===================================================================
--- /branches/eam_branches/ipp-20211108/ppStack/src/ppStackMatch.c	(revision 41907)
+++ /branches/eam_branches/ipp-20211108/ppStack/src/ppStackMatch.c	(revision 41908)
@@ -524,4 +524,5 @@
         psFree(conv);
     } else {
+	// this branch is done if options->convolve is FALSE
         // Match the normalisation
         float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation
Index: /branches/eam_branches/ipp-20211108/ppStack/src/ppStackOptions.c
===================================================================
--- /branches/eam_branches/ipp-20211108/ppStack/src/ppStackOptions.c	(revision 41907)
+++ /branches/eam_branches/ipp-20211108/ppStack/src/ppStackOptions.c	(revision 41908)
@@ -56,4 +56,7 @@
     options->num = 0;
     options->quality = 0;
+
+    options->clipPercent = false;
+
     options->psf = NULL;
     options->sumExposure = NAN;
Index: /branches/eam_branches/ipp-20211108/ppStack/src/ppStackOptions.h
===================================================================
--- /branches/eam_branches/ipp-20211108/ppStack/src/ppStackOptions.h	(revision 41907)
+++ /branches/eam_branches/ipp-20211108/ppStack/src/ppStackOptions.h	(revision 41908)
@@ -16,4 +16,7 @@
     int num;                            // Number of inputs
     int quality;                        // Bad data quality flag
+
+    bool clipPercent;                   // use percentile range to clip?
+
     // Prepare
     pmPSF *psf;                         // Target PSF
Index: /branches/eam_branches/ipp-20211108/ppStack/src/ppStackReadout.c
===================================================================
--- /branches/eam_branches/ipp-20211108/ppStack/src/ppStackReadout.c	(revision 41907)
+++ /branches/eam_branches/ipp-20211108/ppStack/src/ppStackReadout.c	(revision 41908)
@@ -19,4 +19,8 @@
                                          weightings, exposures, addVariance);
     thread->busy = false;
+
+    thread->status = job->results ? PPSTACK_THREAD_SUCCESS : PPSTACK_THREAD_FAILURE;
+
+    psAssert(job->results, "Stacking failed.");
 
     return job->results ? true : false;
@@ -48,4 +52,6 @@
     thread->busy = false;
 
+    thread->status = status ? PPSTACK_THREAD_SUCCESS : PPSTACK_THREAD_FAILURE;
+
     psAssert(status, "Stacking failed.");
 
@@ -53,4 +59,35 @@
 }
 
+bool ppStackReadoutPercentThread(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+
+    psArray *args = job->args;          // Arguments
+    ppStackThread *thread = args->data[0]; // Thread
+    ppStackOptions *options = args->data[1]; // Options
+    pmConfig *config = args->data[2];   // Configuration
+
+    pmReadout *outRO = options->outRO;  // Output readout
+    psVector *mask = options->inputMask; // Mask for inputs
+    psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image
+    psVector *exposures = options->exposures;   // Exposure times for each image
+    psVector *addVariance = options->matchChi2; // Additional variance when rejecting
+
+    // fprintf (stderr, "start thread %d : scan range %d - %d\n", thread->entry, thread->firstScan, thread->lastScan);
+
+    bool status = ppStackReadoutPercent(config, outRO, thread->readouts, mask, weightings, exposures, addVariance);
+
+    thread->busy = false;
+
+    thread->status = status ? PPSTACK_THREAD_SUCCESS : PPSTACK_THREAD_FAILURE;
+
+    // fprintf (stderr, "finished thread %d, status %d : scan range %d - %d\n", thread->entry, thread->status, thread->firstScan, thread->lastScan);
+
+    psAssert(status, "Stacking failed.");
+
+    return status;
+}
+
+//////////////////////////////////// 
 
 bool ppStackInspect(psThreadJob *job)
@@ -145,5 +182,5 @@
 
     int num = readouts->n;              // Number of inputs
-    psArray *stack = psArrayAlloc(num); // Array for stacking
+    psArray *stackData = psArrayAlloc(num); // Array of data to be stacked
 
     for (int i = 0; i < num; i++) {
@@ -160,12 +197,19 @@
         }
 
-        stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i],
+	// stackData is an array of pmStackData structures
+        stackData->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i],
                                           addVariance->data.F32[i]);
     }
 
-    if (!pmStackCombine(outRO, NULL, stack, maskBad, maskSuspect, maskBlank, kernelSize, iter,
+    if (!pmStackCombine(outRO, NULL, stackData, maskBad, maskSuspect, maskBlank, kernelSize, iter,
                         combineRej, combineSys, combineDiscard, useVariance, safe, nminpix, false)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
-        psFree(stack);
+        psFree(stackData);
+	
+	// XXX EAM : if pmStackCombine fails above this return will cause the thread to block.
+	// The failure status results in a psThread.fault.  the psThreadLauncher function, which
+	// each background thread runs as a loop, will sleep until the fault is cleared externally
+	// by the handler.  psThreadPoolWait() catches, counts, and clears these faults, but
+	// ppStackThreadRead does not.  
         return false;
     }
@@ -175,5 +219,5 @@
     psArray *reject = psArrayAlloc(num);  // List of pixels rejected
     for (int i = 0; i < num; i++) {
-        pmStackData *data = stack->data[i]; // Data for this image
+        pmStackData *data = stackData->data[i]; // Data for this image
         if (!data) {
             continue;
@@ -186,5 +230,5 @@
         reject->data[i] = psMemIncrRefCounter(data->reject);
     }
-    psFree(stack);
+    psFree(stackData);
 
     //MEH change to trace
@@ -247,5 +291,5 @@
 
     int num = readouts->n;              // Number of inputs
-    psArray *stack = psArrayAlloc(num); // Array for stacking
+    psArray *stackData = psArrayAlloc(num); // Array of data to be stacked
 
     // We have rejection from a previous combination: combine without flagging pixels to inspect
@@ -272,5 +316,5 @@
                                              addVariance ? addVariance->data.F32[i] : NAN);
         data->reject = rejected ? psMemIncrRefCounter(rejected->data[i]) : NULL;
-        stack->data[i] = data;
+        stackData->data[i] = data;
 
         //MEH -- apply bscale offset before norm   
@@ -288,8 +332,8 @@
     }
 
-    if (!pmStackCombine(outRO, expRO, stack, maskBad, maskSuspect, maskBlank, 0, iter, combineRej,
+    if (!pmStackCombine(outRO, expRO, stackData, maskBad, maskSuspect, maskBlank, 0, iter, combineRej,
                         combineSys, combineDiscard, useVariance, safe, nminpix, rejected)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
-        psFree(stack);
+        psFree(stackData);
         return false;
     }
@@ -307,5 +351,5 @@
     expChip->data_exists = true;
 
-    psFree(stack);
+    psFree(stackData);
 
     //MEH change to trace
@@ -316,2 +360,87 @@
     return true;
 }
+
+bool ppStackReadoutPercent(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
+			       const psVector *mask, const psVector *weightings, const psVector *exposures,
+			       const psVector *addVariance)
+{
+    assert(config);
+    assert(outRO);
+    assert(readouts);
+    assert(mask && mask->n == readouts->n && mask->type.type == PS_TYPE_VECTOR_MASK);
+    assert(weightings && weightings->n == readouts->n && weightings->type.type == PS_TYPE_F32);
+    assert(addVariance && addVariance->n == readouts->n && addVariance->type.type == PS_TYPE_F32);
+    static int sectionNum = 0;          // Section number; for debugging outputs
+
+    // fprintf (stderr, "starting ReadoutPercent, %d\n", sectionNum);
+
+    // Get the recipe values
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
+    psAssert(recipe, "We've thrown an error on this before.");
+
+    bool mdok;                          // Status of MD lookup
+    psF64 minRange = psMetadataLookupF32(&mdok, recipe, "COMBINE.MIN.RANGE"); // Combination threshold
+
+    char defaultBlankStr[16] = "BLANK";
+
+    psString maskBadStr     = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask for bad
+    psString maskSuspectStr = psMetadataLookupStr(&mdok, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
+    psString maskBlankStr   = psMetadataLookupStr(&mdok, recipe, "MASK.BLANK"); // Name of bits to set for empty pixels
+    if (!maskBlankStr) {
+      maskBlankStr = psMetadataLookupStr(&mdok, recipe, "MASK.BAD"); // Old name for MASK.BLANK
+    }
+    if (!maskBlankStr) {
+      maskBlankStr = defaultBlankStr; // this is statically allocated above
+    }
+
+    psImageMaskType maskBad     = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+    psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits
+    psImageMaskType maskBlank   = pmConfigMaskGet(maskBlankStr, config); // Bits to mask for bad pixels
+
+    int num = readouts->n;              // Number of inputs
+    psArray *stackData = psArrayAlloc(num); // Array of data to be stacked
+
+    for (int i = 0; i < num; i++) {
+        pmReadout *ro = readouts->data[i];
+        if (!ro || mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
+            // Bad image
+            continue;
+        }
+
+        // Ensure there is a mask, or pmStackCombine will complain
+        if (!ro->mask) {
+            ro->mask = psImageAlloc(ro->image->numCols, ro->image->numRows, PS_TYPE_IMAGE_MASK);
+            psImageInit(ro->mask, 0);
+        }
+
+        stackData->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i],
+					      addVariance->data.F32[i]);
+    }
+
+    // XXX rename minRange to rejectFraction
+    // XXX add expmaps
+    if (!pmStackCombineByPercentile(outRO, NULL, stackData, minRange, maskBad, maskSuspect, maskBlank)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
+        psFree(stackData);
+        return false;
+    }
+
+    outRO->data_exists = true; // output readout
+    outRO->parent->data_exists = true; // output cell
+    outRO->parent->parent->data_exists = true; // output chip
+
+# if 0
+      pmCell *expCell = expRO->parent;    // Exposure cell
+      pmChip *expChip = expCell->parent;  // Exposure chip
+      expRO->data_exists = true;
+      expCell->data_exists = true;
+      expChip->data_exists = true;
+# endif
+
+    psFree(stackData);
+    sectionNum++;
+
+    return true;
+}
+
+
Index: /branches/eam_branches/ipp-20211108/ppStack/src/ppStackSetup.c
===================================================================
--- /branches/eam_branches/ipp-20211108/ppStack/src/ppStackSetup.c	(revision 41907)
+++ /branches/eam_branches/ipp-20211108/ppStack/src/ppStackSetup.c	(revision 41908)
@@ -18,4 +18,6 @@
     options->convolve = psMetadataLookupBool(NULL, recipe, "CONVOLVE"); // Convolve images?
     options->doBackground = psMetadataLookupBool(NULL, recipe, "BACKGROUND.MODEL");
+
+    options->clipPercent  = psMetadataLookupBool(NULL, recipe, "COMBINE.PERCENT"); // use percentile range to clip?
 
     if (!psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) {
Index: /branches/eam_branches/ipp-20211108/ppStack/src/ppStackThread.c
===================================================================
--- /branches/eam_branches/ipp-20211108/ppStack/src/ppStackThread.c	(revision 41907)
+++ /branches/eam_branches/ipp-20211108/ppStack/src/ppStackThread.c	(revision 41908)
@@ -19,4 +19,6 @@
     thread->firstScan = 0;
     thread->lastScan = 0;
+    thread->entry = -1; // this gets assigned after allocation
+    thread->status = PPSTACK_THREAD_NEW;
     return thread;
 }
@@ -124,5 +126,7 @@
             readouts->data[j] = ro;
         }
-        threads->data[i] = ppStackThreadAlloc(readouts);
+        ppStackThread *thread = ppStackThreadAlloc(readouts);
+	thread->entry = i;
+        threads->data[i] = thread;
         psFree(readouts);               // Drop reference
     }
@@ -157,6 +161,8 @@
     while (true) {
         // check for any groups which can read data
+	// NOTE: threads->n is set (in ppStackThreadDataSetup) to be 1 more than the argument to -threads
         for (int j = 0; j < threads->n; j++) {
             ppStackThread *thread = threads->data[j];
+	    // fprintf (stderr, "grab thread %d, scan: %d - %d : busy: %d read: %d status: %d\n", thread->entry, thread->firstScan, thread->lastScan, thread->busy, thread->read, thread->status);
             if (thread->read) {
                 continue;
@@ -166,4 +172,5 @@
             thread->lastScan = stack->lastScan + rows;
             stack->lastScan = thread->lastScan;
+	    thread->status = PPSTACK_THREAD_READ;
 
             psArray *readouts = thread->readouts;
@@ -172,4 +179,5 @@
 
             psTrace("ppStack", 2, "Reading data for chunk %d into group %d....\n", numChunk, j);
+            // fprintf (stderr, "Reading data for chunk %d into group %d (thread %d, scan: %d - %d)....\n", numChunk, j, thread->entry, thread->firstScan, thread->lastScan);
             for (int i = 0; i < readouts->n; i++) {
                 pmReadout *ro = readouts->data[i]; // Input readout
@@ -177,4 +185,6 @@
                     continue;
                 }
+
+		// fprintf (stderr, "read data for readout %d\n", i);
 
                 // override the recorded last scan
@@ -232,4 +242,5 @@
 
             thread->read = thread->busy = true;
+	    thread->status = PPSTACK_THREAD_RUN;
             return thread;
         }
@@ -239,7 +250,10 @@
         for (int j = 0; j < threads->n; j++) {
             ppStackThread *thread = threads->data[j];
+	    // fprintf (stderr, "test thread %d, scan: %d - %d : busy: %d, read: %d, status: %d\n", thread->entry, thread->firstScan, thread->lastScan, thread->busy, thread->read, thread->status);
+
             if (thread->busy) {
                 continue;
             }
+	    thread->status = PPSTACK_THREAD_NEW;
             thread->read = false;
             wait = false;
@@ -283,4 +297,10 @@
     }
 
+    {
+        psThreadTask *task = psThreadTaskAlloc("PPSTACK_PERCENT_COMBINE", 3);
+        task->function = &ppStackReadoutPercentThread;
+        psThreadTaskAdd(task);
+        psFree(task);
+    }
     return;
 }
Index: /branches/eam_branches/ipp-20211108/ppStack/src/ppStackThread.h
===================================================================
--- /branches/eam_branches/ipp-20211108/ppStack/src/ppStackThread.h	(revision 41907)
+++ /branches/eam_branches/ipp-20211108/ppStack/src/ppStackThread.h	(revision 41908)
@@ -11,4 +11,6 @@
     int firstScan;                      // First row of the chunk to be read for this group
     int lastScan;                       // Last row of the chunk to be read for this group
+    int entry;				// thread number for tracking progress
+    int status;				// thread number for tracking progress
 } ppStackThread;
 
@@ -47,3 +49,6 @@
 
 
+enum {PPSTACK_THREAD_NEW, PPSTACK_THREAD_READ, PPSTACK_THREAD_RUN, PPSTACK_THREAD_SUCCESS, PPSTACK_THREAD_FAILURE};
+
+
 #endif