diff --git a/README.md b/README.md
index 56550051..283d1a0c 100644
--- a/README.md
+++ b/README.md
@@ -51,5 +51,23 @@ Running Unit Tests
 
 	python run_unit_tests.py
 
+The modification by Xin Song
+------------------
+* Colour alarm for traces with very large traveltimes
+
+When the measured relative traveltime after MCCC exceeds the traveltime thresholds, the corresponding traces will show different colours. This feature is to avoid some high risky traces that might have very large traveltimes. It can have 3 traveltime thresholds maximum.
+
+This feature is set in the `ttconfig.py` in either `aimbat/src/pysmo/aimbat/` or either local directory
+
+For example in `ttconfig.py`:
+```
+thresholds = 2 2.5 3
+colorthresholds = green magenta red
+```
+It means the colour of trace will become green, magenta, or red if the traveltime of this trace is over 2s, 2.5s, or 3s after MCCC.
+
+* Perform the MCCC without changing the polarity of traces by setting `xcorr_func` in the `ttconfig.py`:
+
+`xcorr_func = xcorr_full_polarity`
 
 
diff --git a/src/pysmo/aimbat/algmccc.py b/src/pysmo/aimbat/algmccc.py
index 82f29db9..d5de3718 100755
--- a/src/pysmo/aimbat/algmccc.py
+++ b/src/pysmo/aimbat/algmccc.py
@@ -26,7 +26,7 @@
 from time import strftime, tzname
 from optparse import OptionParser
 from ttconfig import MCConfig 
-from sacpickle import loadData, saveData, SacDataHdrs, taper, taperWindow, windowIndex, windowData 
+from sacpickle import loadData, saveData, SacDataHdrs, taper, taperWindow, windowIndex, windowData
 from qualsort import initQual, seleSeis
 
 def getOptions():
@@ -279,7 +279,7 @@ def WriteFileWithDelay(mcpara, solist, solution, outvar, outcc, t0_times, delay_
 	line1 = 'station, mccc delay,    std,    cc coeff,  cc std,   pol   , t0_times  , delay_times\n'
 	ofile.write( line0 )
 	ofile.write( line1 )
-	fmt = ' {0:<9s} {1:9.4f} {2:9.4f} {3:>9.4f} {4:>9.4f} {5:4d}  {6:<s}  {7:9.4f}  {8:9.4f}\n'
+	fmt = ' {0:<9s} {1:11.8f} {2:9.4f} {3:>9.4f} {4:>9.4f} {5:4d}  {6:<s}  {7:9.4f}  {8:9.4f}\n'
 
 	selist_LonLat = zeros(shape=(len(solist),2))
 	nsta = len(solist)
@@ -308,7 +308,7 @@ def WriteFileWithDelay(mcpara, solist, solution, outvar, outcc, t0_times, delay_
 
 	return selist_LonLat, delay_times
 
-def WriteFileOriginal(mcpara, solist, solution, outvar, outcc, itmean):
+def WriteFileOriginal(mcpara, solist, solution, outvar, outcc, itmean, invdata):
 	ofilename = "original"+mcpara.mcname
 	kevnm = mcpara.kevnm
 	delta = mcpara.delta
@@ -326,7 +326,7 @@ def WriteFileOriginal(mcpara, solist, solution, outvar, outcc, itmean):
 	line1 = 'station, mccc delay,    std,    cc coeff,  cc std,   pol\n'
 	ofile.write( line0 )
 	ofile.write( line1 )
-	fmt = ' {0:<9s} {1:9.4f} {2:9.4f} {3:>9.4f} {4:>9.4f} {5:4d}  {6:<s}\n'
+	fmt = ' {0:<9s} {1:11.8f} {2:9.4f} {3:>9.4f} {4:>9.4f} {5:4d}  {6:<s}\n'
 
 	selist_LonLat = zeros(shape=(len(solist),2))
 	nsta = len(solist)
@@ -351,9 +351,11 @@ def WriteFileOriginal(mcpara, solist, solution, outvar, outcc, itmean):
 	# write phase and event
 	ofile.write( 'Phase: {0:8s} \n'.format(mcpara.phase) )
 	ofile.write( mcpara.evline + '\n' )
+        for i in range(len(invdata)):
+                ofile.write( str(invdata[i]) + '\n')
 	ofile.close()	
 
-def corrite(solist, mcpara, reftimes, solution, outvar, outcc):
+def corrite(solist, mcpara, reftimes, solution, outvar, outcc, invdata):
 	""" Write output file, set output time picks.
 	"""
 	ofilename = mcpara.mcname
@@ -370,13 +372,14 @@ def corrite(solist, mcpara, reftimes, solution, outvar, outcc):
 	itmean = mean(reftimes)
 	for i in range(nsta):
 		wt = itmean + solution[i,0]
+		solist[i].sethdr('user5', solution[i,0])
 		solist[i].sethdr(wpick, wt)
 	t0_times = [sacdh.gethdr('t0') for sacdh in solist]
 	delay_times = [(solution[i,0]-(t0_times[i]-itmean)) for i in xrange(nsta)]
 
 	selist_LonLat, delay_times = WriteFileWithDelay(mcpara, solist, solution, outvar, outcc, t0_times, delay_times, itmean)
-	WriteFileOriginal(mcpara, solist, solution, outvar, outcc, itmean)
-	return selist_LonLat, delay_times
+	WriteFileOriginal(mcpara, solist, solution, outvar, outcc, itmean, invdata)
+	return selist_LonLat, delay_times, itmean
 
 def mccc(gsac, mcpara):
 	""" Run MCCC. 
@@ -422,7 +425,7 @@ def mccc(gsac, mcpara):
 	outvar = sqrt(sum(resmatrix**2)/2/(nsta*(nsta-1)/2))
 	outcc = mean(ccmean)
 
-	selist_LonLat, delay_times = corrite(solist, mcpara, reftimes, solution, outvar, outcc)
+	selist_LonLat, delay_times, itmean = corrite(solist, mcpara, reftimes, solution, outvar, outcc, invdata)
 
 	# set wpick as ipick for deleted ones and array stack
 	for sacdh in delist:
@@ -430,7 +433,7 @@ def mccc(gsac, mcpara):
 	stkdh = gsac.stkdh
 	stkdh.sethdr(wpick, stkdh.gethdr(ipick))
 
-	return solution, selist_LonLat, delay_times
+	return solution, selist_LonLat, delay_times, itmean
 
 
 def eventListName(evlist='event.list', phase='S', isol='PDE'):
diff --git a/src/pysmo/aimbat/pickphase.py b/src/pysmo/aimbat/pickphase.py
index 2029f306..03cbcee4 100755
--- a/src/pysmo/aimbat/pickphase.py
+++ b/src/pysmo/aimbat/pickphase.py
@@ -204,7 +204,8 @@ def labelStation(self):
 			cc = sacdh.gethdr(hdrcc)
 			sn = sacdh.gethdr(hdrsn)
 			co = sacdh.gethdr(hdrco)
-			slab += 'qual={0:4.2f}/{1:.1f}/{2:4.2f}'.format(cc, sn, co)
+                        mc = sacdh.gethdr('user5')
+			slab += 'qual={0:5.3f}/{1:4.2f}/{2:.1f}/{3:4.2f}'.format(mc, cc, sn, co)
 		trans = transforms.blended_transform_factory(axpp.transAxes, axpp.transData)
 		font = FontProperties()
 		font.set_family('monospace')
@@ -235,7 +236,10 @@ def changeColor(self):
 		""" Change color of a seismogram based on selection status. 
 		"""
 		if self.sacdh.selected:
-			col = self.opts.pppara.colorwave
+                        if self.color == 'gray':
+                                col = self.opts.pppara.colorwave
+                        else:
+                                col = self.color
 		else:
 			col = self.opts.pppara.colorwavedel
 		setp(self.stalabel, color=col)
@@ -424,6 +428,10 @@ def initIndex(self):
 		delist = self.gsac.delist
 		nsel = len(selist)
 		ndel = len(delist)
+                if hasattr(self.gsac, 'itmean'):
+                        self.itmean = self.gsac.itmean
+                else:
+                        self.itmean = 0.
 		maxsel, maxdel = opts.maxnum
 		pagesize = maxsel + maxdel
 		aipages, ayindex, aybases, ayticks = indexBaseTick(nsel, ndel, pagesize, maxsel)
@@ -451,13 +459,25 @@ def plotWave(self):
 		# get colors from sacdh.selected
 		colsel = opts.pppara.colorwave
 		coldel = opts.pppara.colorwavedel
+                threshd = opts.pppara.thresholds
+                colthd = opts.pppara.colorthresholds
 		colors = [[None,] * npsel , [None,] * npdel]
 		for j in range(2):
 			for k in range(nsede[j]):
+                                #plists[j][k].users[3] = plists[j][k].thdrs[3] - self.itmean
+                                dt = plists[j][k].users[5]
 				if plists[j][k].selected:
-					colors[j][k] = colsel
+                                        if abs(dt) > 10000 or abs(dt) < threshd[0]:
+                                                colors[j][k] = colsel
+                                        elif threshd[0] <= abs(dt) < threshd[1]:
+                                                colors[j][k] = colthd[0]
+                                        elif threshd[1] <= abs(dt) < threshd[2]:
+                                                colors[j][k] = colthd[1]
+                                        if threshd[2] <= abs(dt) < 10000:
+                                                colors[j][k] = colthd[2]
 				else:
 					colors[j][k] = coldel
+
 		# plot
 		pps = []
 		for j in range(2):
diff --git a/src/pysmo/aimbat/qualctrl.py b/src/pysmo/aimbat/qualctrl.py
index fe01a4a1..4e946142 100755
--- a/src/pysmo/aimbat/qualctrl.py
+++ b/src/pysmo/aimbat/qualctrl.py
@@ -167,7 +167,7 @@ def setLabels(self):
 		trans = transforms.blended_transform_factory(axpp.transAxes, axpp.transData)
 		font = FontProperties()
 		font.set_family('monospace')
-		axpp.text(1.025, 0, ' '*8+'qual= CCC/SNR/COH', transform=trans, va='center', 
+		axpp.text(1.025, 0, ' '*8+'qual= MCCC/CCC/SNR/COH', transform=trans, va='center', 
 			color='k', fontproperties=font)
 
 	def plotStack(self):
@@ -1109,7 +1109,7 @@ def mccc(self, event):
 		mcpara.evline = evline
 		mcpara.mcname = mcname
 		mcpara.kevnm = gsac.kevnm
-		solution, solist_LonLat, delay_times = mccc(gsac, mcpara)
+		solution, solist_LonLat, delay_times, itmean = mccc(gsac, mcpara)
 		self.gsac.solist_LonLat = solist_LonLat
 		self.gsac.delay_times = delay_times
 
diff --git a/src/pysmo/aimbat/ttconfig.py b/src/pysmo/aimbat/ttconfig.py
index 2f3801aa..9b2a06c2 100644
--- a/src/pysmo/aimbat/ttconfig.py
+++ b/src/pysmo/aimbat/ttconfig.py
@@ -80,6 +80,8 @@ def __init__(self):
 		self.srate = config.getfloat('sacplot', 'srate')
 		self.tapertype = config.get('signal', 'tapertype')
 		self.taperwidth = config.getfloat('signal', 'taperwidth')
+                self.thresholds = [ float(val) for val in config.get('sacplot', 'thresholds').split() ]
+                self.colorthresholds = config.get('sacplot', 'colorthresholds').split()
 
 
 
@@ -109,6 +111,9 @@ def __init__(self):
 		self.pickstyles = config.get('sacplot', 'pickstyles').split()
 		self.minspan = config.getint('sacplot', 'minspan')
 		self.fstack = config.get('iccs', 'fstack')
+                self.thresholds = [ float(val) for val in config.get('sacplot', 'thresholds').split() ]
+                self.colorthresholds = config.get('sacplot', 'colorthresholds').split()
+
 
 
 class CCConfig:
diff --git a/src/pysmo/aimbat/ttdefaults.conf b/src/pysmo/aimbat/ttdefaults.conf
index bfe36eeb..3e476835 100644
--- a/src/pysmo/aimbat/ttdefaults.conf
+++ b/src/pysmo/aimbat/ttdefaults.conf
@@ -12,8 +12,8 @@ alphatwfill = 0.2
 alphatwsele = 0.6
 # plot picks t0-t5
 npick = 6
-pickcolors = kmrcgyb
-pickstyles = --  :
+pickcolors = rkmcgyb
+pickstyles = -  :
 
 ### Fig size and rect for seismogram [x0, y0, dx, dy]
 figsize = 8 10
@@ -24,6 +24,10 @@ minspan = 5
 ### sample rate. read first file if srate<0
 srate = -1
 
+### thresholds for relative traveltime after mccc
+thresholds = 3.5 4.0 5.0
+colorthresholds =  green magenta red
+
 ##############################################
 [sachdrs]
 ### time window
@@ -47,14 +51,14 @@ qweights = 0.3333 0.3333 0.3333
 srate = -1 
 ### cross-correlation 
 xcorr_modu = xcorrf90
-xcorr_func = xcorr_fast
+xcorr_func = xcorr_full_polarity
 shift = 10
-maxiter = 10
+maxiter = 20
 ### convergence criterion
 ### type: coef for correlation coefficient, resi for residual
 ### epsi: small number
 convtype = coef
-convepsi = 0.001
+convepsi = 0.0001
 ### stack weight: coef for correlation coefficient, otherwise even weight
 stackwgt = coef
 
@@ -67,7 +71,7 @@ srate = -1
 ### output file: "$evdate".mc if mc
 ofilename = mc
 xcorr_modu = xcorrf90
-xcorr_func = xcorr_faster
+xcorr_func = xcorr_full_polarity
 shift = 10
 extraweight = 1000.
 ### lsqr solution without weighting;
diff --git a/src/pysmo/aimbat/xcorr.f90 b/src/pysmo/aimbat/xcorr.f90
index fc377e7b..326687d1 100644
--- a/src/pysmo/aimbat/xcorr.f90
+++ b/src/pysmo/aimbat/xcorr.f90
@@ -84,6 +84,49 @@ end subroutine xcorr_full
 !!!---------------------------------------------------------------------!!!
 
 
+!!!---------------------------------------------------------------------!!!
+subroutine xcorr_full_polarity(x, y, n, shift, delay, ccmax, ccpol)
+! Cross-correlation of array x and y.
+! It does not correct the polarity, meaning it does only calculate the maximum value of cross correlation, not including minimum, which is different from xcorr_full
+! Return time shift, correlation coefficient, and polarity at maximum correlation.
+implicit none
+real(8) :: crosscorr
+integer :: n, shift, delay, ccpol
+real(8), dimension(0:n-1) :: x, y
+real(8) :: cc, ccmax, ccmin
+integer :: k, kmin
+!f2py intent(in)  :: x, y, n, shift
+!f2py intent(out) :: delay, ccmax, ccpol
+!f2py intent(hide):: c, k, kmin, ccmin
+
+shift = 1
+ccmax = 0
+ccmin = 0
+kmin = 0
+do k = -n+1,n-1,shift
+    cc = crosscorr(x,y,n,k)   
+    if (cc.gt.ccmax) then
+        ccmax = cc
+        delay = k
+    endif
+    if (cc.lt.ccmin) then
+        ccmin = cc
+        kmin = k
+    endif
+enddo
+! if (ccmax.gt.-ccmin) then
+!     ccpol = 1
+! else
+!     ccmax = -ccmin
+!     delay = kmin
+!     ccpol = -1
+ccpol = 1
+ccmax = ccmax/sqrt( dot_product(x,x) * dot_product(y,y) )
+
+end subroutine xcorr_full_polarity
+!!!---------------------------------------------------------------------!!!
+
+
 !!!---------------------------------------------------------------------!!!
 subroutine xcorr_fast(x, y, n, shift, delay, ccmax, ccpol)
 ! Fast cross-correlation using 1 level of coarse shift.
diff --git a/src/pysmo/aimbat/xcorr.py b/src/pysmo/aimbat/xcorr.py
old mode 100644
new mode 100755
index f0a9656b..7289e01a
--- a/src/pysmo/aimbat/xcorr.py
+++ b/src/pysmo/aimbat/xcorr.py
@@ -50,6 +50,28 @@ def _xcorr(x, y, cmode='full'):
 		delay -= (len(y)-1)/2
 	return delay, ccmax, ccpol
 
+def _xcorr_polarity(x, y, cmode='full'):
+	""" 
+	Cross-correlation of two 1-D arrays using 'full' or 'same' mode.
+		c[k] = sum_i x[i]*y[i+k]
+	Full mode: indices of array c --> j=0:nx+ny-2 <-- k=ny-1:-nx+1:-1 (j=ny-1-k)
+        It does not correct the polarity, meaning it does only calculate the maximum
+        value of cross correlaton, not including minimum, which is different from
+        _xcorr
+	Return time shift of y relative to x at maximum correlation and polarity.
+	"""
+	cc = correlate(x, y, cmode)
+	imax = argmax(cc)
+	imin = argmin(cc)
+	ccmax = cc[imax]
+	ccmin = cc[imin]
+        ccpol = 1
+        delay = len(y)-1-imax
+        ccmax =  ccmax / sqrt(dot(x,x)*dot(y,y))
+	if cmode == 'same':
+		delay -= (len(y)-1)/2
+	return delay, ccmax, ccpol
+
 def xcorr_full(x, y, shift=1):
 	""" 
 	Cross-correlation of two 1-D arrays using 'full' mode.
@@ -57,6 +79,14 @@ def xcorr_full(x, y, shift=1):
 	"""
 	return _xcorr(x, y, 'full')
 
+def xcorr_full_polarity(x, y, shift=1):
+        """ 
+	Cross-correlation of two 1-D arrays using 'full' mode.
+	Argument shift=1 is here only in order to make the same number of arguments for all xcorr functions.
+        Not correct the polarity
+	"""
+	return _xcorr_polarity(x, y, 'full')
+
 def xcorr_same(x, y):
 	""" 
 	Cross-correlation of two 1-D arrays using 'same' mode.