diff --git a/binning_v_given_c.py b/binning_v_given_c.py
index 8d97ca2..7de70b8 100644
--- a/binning_v_given_c.py
+++ b/binning_v_given_c.py
@@ -7,8 +7,12 @@ import argparse
 import numpy as np
 
 
+program_description = '''\
+Compute conditional mean of v given condition c by data binning
+'''
+
 # Commandline argument parser
-parser = argparse.ArgumentParser()
+parser = argparse.ArgumentParser(description=program_description)
 parser.add_argument("-x", "--xindex", help="sampling x index", type=int, required=True)
 parser.add_argument("-v", "--variable", help="file containing variable to calculate conditional mean", required=True)
 parser.add_argument("-n", "--nbins", help="number of bins per condtion c interval", default=45, )
@@ -36,7 +40,7 @@ print ("Computing conditional mean at {}".format(xidx))
 
 def cmean_binning (c, v, nbins=100, boundary=0.1, quiet=False):
     '''
-    Compute conditional mean of v given condition c by binning
+    Compute conditional mean of v given condition c by data binning
     '''
     c_bins = np.hstack(
         (np.linspace(0, boundary, nbins),
diff --git a/concatenate_cmean.py b/concatenate_cmean.py
new file mode 100644
index 0000000..59a324e
--- /dev/null
+++ b/concatenate_cmean.py
@@ -0,0 +1,49 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+import sys
+import os
+import argparse
+import numpy as np
+import dnstool
+
+
+program_description = '''\
+concatenate all results (1~nx) from binning.
+'''
+
+# Commandline argument parser
+parser = argparse.ArgumentParser(description=program_description)
+parser.add_argument("-c", "--case", help="target case name", required=True)
+parser.add_argument("-v", "--variable", help="conditional mean variable name", required=True)
+args = parser.parse_args()
+params = vars(args)
+
+
+# Parameters 
+cases = dnstool.case_library()
+casename = params["case"]
+case = cases[casename]
+
+nx, ny, nz = case.shape
+
+vname = params["variable"]
+output_name = "cmean_{}_given_c".format(vname)
+
+
+batch_cmean = [np.load('./cmean_{}_given_c_{:03d}.npz'.format(vname, i)) for i in range(nx)]
+
+cstar = np.asarray([d['cstar'] for d in batch_cmean])
+cmean = np.asarray([d['cmean'] for d in batch_cmean])
+bin_edges = np.asarray([d["bins"]  for d in batch_cmean])
+bin_counts = np.asarray([d["count"] for d in batch_cmean])
+
+
+# Save result
+arr_dict = {}
+arr_dict["cstar"] = cstar
+arr_dict["cmean"] = cmean
+arr_dict["bins"] =  bin_edges
+arr_dict["count"] = bin_counts
+
+np.savez(output_name, **arr_dict)
diff --git a/extract_all.py b/extract_all.py
index 921e708..cbcc768 100644
--- a/extract_all.py
+++ b/extract_all.py
@@ -4,8 +4,14 @@ import numpy as np
 import dnstool
 from pycompact import CompactScheme
 
+
+program_description = '''\
+Read all DNS data files and extract u, c and ddx(c) to single array.
+New array dimension is (nx, time, ny, nz)
+'''
+
 # Commandline argument parser
-parser = argparse.ArgumentParser()
+parser = argparse.ArgumentParser(description=program_description)
 parser.add_argument("-c", "--case", help="target case name", required=True)
 args = parser.parse_args()
 params = vars(args)
diff --git a/fluctuation_product_u_ddxc.py b/fluctuation_product_u_ddxc.py
index 03344fe..4bfefdf 100644
--- a/fluctuation_product_u_ddxc.py
+++ b/fluctuation_product_u_ddxc.py
@@ -4,18 +4,20 @@ import numpy as np
 import dnstool
 from pycompact import CompactScheme
 
+
+program_description = '''\
+Calculate product of conditional fluctuations of u and ddx(c).
+'''
+
 # Commandline argument parser
-parser = argparse.ArgumentParser()
+parser = argparse.ArgumentParser(description=program_description)
 parser.add_argument("-c", "--case", help="target case name", required=True)
 args = parser.parse_args()
 params = vars(args)
 
-casename = params["case"]
-
 cases = dnstool.case_library()
-
+casename = params["case"]
 case = cases[casename]
-
 nx, ny, nz = case.shape
 
 cs = CompactScheme(nx, ny, nz, False, True, True, 4, 2, 2)
@@ -26,24 +28,32 @@ ufile = "u.dat"
 ddxcfile = "ddxc.dat"
 resfile = "uddxc.dat"
 
-c = np.memmap(cfile, mode="r", dtype=np.double).reshape((512,-1,256,256))
-u = np.memmap(ufile, mode="r", dtype=np.double).reshape((512,-1,256,256))
-ddxc = np.memmap(ufile, mode="r", dtype=np.double).reshape((512,-1,256,256))
-storage = np.memmap(resfile, mode='w+', dtype=np.double, shape=(nx, len(case.data_files), ny, nz))
+c = np.memmap(cfile, mode="r", dtype=np.double).reshape((nx,-1,ny,nz))
+u = np.memmap(ufile, mode="r", dtype=np.double).reshape((nx,-1,ny,nz))
+ddxc = np.memmap(ddxcfile, mode="r", dtype=np.double).reshape((nx,-1,ny,nz))
+storage = np.memmap(resfile, mode='w+', dtype=np.double, shape=c.shape)
 
+batch_cmean_ddxc_c = np.load('./cmean_ddxc_given_c.npz')
+print(*batch_cmean_ddxc_c.keys())
+
+cstar_ddxc_c = batch_cmean_ddxc_c['cstar']
+cmean_ddxc_c = batch_cmean_ddxc_c['cmean']
+
+batch_cmean_u_c = np.load('./cmean_u_given_c.npz')
+
+cstar_u_c = batch_cmean_u_c['cstar']
+cmean_u_c = batch_cmean_u_c['cmean']
 
 for xidx in range(nx):
 
-    u_cstar = 
-    u_mean = 
-    u_mean_at_x = np.interp (c[xidx], u_cstar, u_mean)
-    u_at_x = u[xidx] 
+    print (datetime.datetime.now(), xidx)
+
+    u_mean_at_x = np.interp (c[xidx], cstar_u_c[xidx], cmean_u_c[xidx])
+    u_at_x = u[xidx]
     u_flux_at_x = u_at_x - u_mean_at_x
 
-    ddxc_cstar = 
-    ddxc_mean = 
-    ddxc_mean_at_x = np.interp (c[xidx], ddxc_cstar, ddxc_mean)
-    ddxc_at_x = u[xidx] 
+    ddxc_mean_at_x = np.interp (c[xidx], cstar_ddxc_c[xidx], cmean_ddxc_c[xidx])
+    ddxc_at_x = u[xidx]
     ddxc_flux_at_x = ddxc_at_x - ddxc_mean_at_x
 
     storage[xidx,:,:,:] = u_flux_at_x * ddxc_flux_at_x
diff --git a/run_in_case_dir.sh b/run_in_case_dir.sh
index 25deb58..880803f 100755
--- a/run_in_case_dir.sh
+++ b/run_in_case_dir.sh
@@ -42,6 +42,9 @@ echo sbatch --array=0-511 -J $1-u-c $SCRIPT_ROOT/sbatch-job-binning-ddxc-array
 export BIN_TARGET=ddxc.dat
 echo sbatch --wait --array=0-511 -J $1-ddxc-c $SCRIPT_ROOT/sbatch-job-binning-ddxc-array
 
+echo time python $SCRIPT_ROOT/concatenate_cmean.py -c $1 -v u
+echo time python $SCRIPT_ROOT/concatenate_cmean.py -c $1 -v ddxc
+
 
 #=============================================================================#
 # Calculate product of conditional fluctuations of u, and ddx(c)
@@ -58,6 +61,8 @@ echo Batch Running - binning_v_given_c.py
 export BIN_TARGET=uddxc.dat
 echo sbatch --wait --array=0-511 -J $1-uddxc-c $SCRIPT_ROOT/sbatch-job-binning-ddxc-array
 
+echo time python $SCRIPT_ROOT/concatenate_cmean.py -c $1 -v uddxc
+
 
 #=============================================================================#
 echo Done