diff --git a/da/stilt/stilt.co2.simu.2015.Sep.bc.r b/da/stilt/stilt.co2.simu.2015.Sep.bc.r
index 78fc7efcaa624ba71ebfd6778857af367dbd7c75..6de6fdb76b56df804f520caa8f98a2c1a163d02b 100755
--- a/da/stilt/stilt.co2.simu.2015.Sep.bc.r
+++ b/da/stilt/stilt.co2.simu.2015.Sep.bc.r
@@ -1,4 +1,4 @@
-#------------------------------------------------------------------------------------------------
+#i------------------------------------------------------------------------------------------------
 # CO2 concentration simulation based on WRF-STILT/HYSPLIT-NAM12 footprints and
 # SiB3/4 biosphere fluxes, CarbonTracker background fluxes and boundary conditions
 # 
@@ -70,7 +70,7 @@ library(futile.logger)
 
 logfile = paste("stilt_",procs,".log", sep="")
 flog.appender(appender.file(logfile), name='logger.b')
-flog.info("Launch WRF-STILT-based forward simulations %s", "file", name='logger.b')
+flog.info("======Launch WRF-STILT-based forward simulations======", name='logger.b')
 
 flog.info("The biosphere flux is from %s", bioflux_flag, name='logger.b')
 flog.info("The boundary is from %s", boundary_flag, name='logger.b')
@@ -204,7 +204,7 @@ newfile=T #output results into a new file
 if(boundary_flag == "EMP")
 {
     bouf=paste(bounddir,y1,"/",m1,"/","CO2.v201209_v1.init.stilt.",y1,m1,".txt",sep="")
-    print(bouf)
+    #print(bouf)
     val<-as.matrix(read.table(bouf,header=TRUE))
 }
 
@@ -221,6 +221,8 @@ for(mm in 1:length(pfbfns))
   footnc=nc_open(fn)
   foot=ncvar_get(footnc,"foot1",start=c(41,12,1),count=c(ncol2,nrow2,-1))
 
+  #print(mean(foot))
+  
   # get particles
   #part=footp$partfoot 
   #partna=footp$partfootnames
@@ -255,6 +257,8 @@ for(mm in 1:length(pfbfns))
     
     # get the time stamp for each foot step, going back time given by "endtime"
     xx=ISOdatetime(inityr,initmo,initdy,inithr,initmi,0,tz="UTC")
+    print("XXXxx?")
+    print(xx)
     yy=xx+(-foottimes*3600)
     cyy=as.character(yy)
     yrlist=substring(cyy,1,4)
@@ -268,7 +272,9 @@ for(mm in 1:length(pfbfns))
     udaystring=unique(daystring)
     yrmonstring=paste(yrlist,monlist, sep="")
     uyrmonstring=unique(yrmonstring)
-    
+   
+    print("XXXudaystring 1st from ident name")
+    print(udaystring)
     #current month
     sibyr=substring(uyrmonstring[1],1,4)
     sibmon=substring(uyrmonstring[1],5,6)	
@@ -292,24 +298,26 @@ for(mm in 1:length(pfbfns))
     hgtarr=endpts[,5] + endpts[,6]  #agl+grdht
     
     # analyze the dates to recognize how many data to read
-    dd=as.character(endptsdate)
-    yrlist=substring(dd,1,4)
-    monlist=substring(dd,6,7)
-    daylist=substring(dd,9,10)
-    hrlist=substring(dd,12,13)
-    milist=substring(dd,15,16)
+    #dd=as.character(endptsdate)
+    #yrlist=substring(dd,1,4)
+    #monlist=substring(dd,6,7)
+    #daylist=substring(dd,9,10)
+    #hrlist=substring(dd,12,13)
     
     # get unique months and days
-    daystring=paste(yrlist,monlist,daylist, sep="")
-    udaystring=unique(daystring)
-    yrmonstring=paste(yrlist,monlist, sep="")
-    uyrmonstring=unique(yrmonstring)
-    
+    #daystring=paste(yrlist,monlist,daylist, sep="")
+    #udaystring=unique(daystring)
+    #yrmonstring=paste(yrlist,monlist, sep="")
+    #uyrmonstring=unique(yrmonstring)
+
+    #print("daylist=") #it's wrong to use daylist independently, but to use daystring (yyyymmdd)
+    #print(daylist)
+
     #--------------------------------------------------------------
     # 2-1: Read fluxes and boundary data
     #--------------------------------------------------------------
     ndays=length(udaystring)
-    
+
     if(boundary_flag=="CTNA")
         bouarr=array(0,dim=c(ndays,120,90,25,8))  #boundary : lon,lat,hgt,time
     
@@ -454,9 +462,9 @@ for(mm in 1:length(pfbfns))
         inxd=ceiling(hh/24)
         inxh=ceiling( (hh-(inxd-1)*24)/3 )   #every three hours the same
        
-        print(hh)
-        print(inxh)
-        print(dim(ocn))
+        #print(hh)
+        #print(inxh)
+        #print(dim(ocn))
         ocnflux[,,ntimes-hh+1]=ocn[inxd,52:117,113:152,inxh]
         fosflux[,,ntimes-hh+1]=fos[inxd,52:117,113:152,inxh]
         firflux[,,ntimes-hh+1]=fir[inxd,52:117,113:152,inxh]
@@ -475,28 +483,7 @@ for(mm in 1:length(pfbfns))
     gppflux[,,]=replace(gppflux[,,],is.na(gppflux[,,]),0)
     recflux[,,]=replace(recflux[,,],is.na(recflux[,,]),0)
         
-
-    #for(hh in ntimes2:1)
-    #{
-    #  dateall2[ntimes2-hh+1]=time2[hh]
-    #  inxd=ceiling(hh/8)
-    #  inxh=hh-(inxd-1)*8
-    #  ocnflux[,,ntimes2-hh+1]=ocn[inxd,52:117,113:152,inxh]  #transform matrix, delete this 
-    #  fosflux[,,ntimes2-hh+1]=fos[inxd,52:117,113:152,inxh]  
-    #  firflux[,,ntimes2-hh+1]=fir[inxd,52:117,113:152,inxh]   
-     
-    #  if(bioflux_flag == "SiBCASA")
-    #  {
-    #      gppflux[,,ntimes2-hh+1]=gpp[inxd,52:117,113:152,inxh]
-    #      recflux[,,ntimes2-hh+1]=rec[inxd,52:117,113:152,inxh]
-    #  }
-
-    #  if(bioflux_flag == "CT_OPT")
-    #      bioflux[,,ntimes2-hh+1]=bio[inxd,52:117,113:152,inxh]
-
-    #}
-     
-    # replace NA values with 0, as NA was used as zero values for datasets
+    
     ocnflux[,,]=replace(ocnflux[,,],is.na(ocnflux[,,]),0)
     fosflux[,,]=replace(fosflux[,,],is.na(fosflux[,,]),0)
     firflux[,,]=replace(firflux[,,],is.na(firflux[,,]),0)
@@ -514,8 +501,6 @@ for(mm in 1:length(pfbfns))
      if(bioflux_flag == "CT_OPT")
          neeflux[,,]= bioflux[,,]*1e6  # for scaling CT flux
 
-    #flog.info("bioflux_flag = %s", bioflux_flag, name='logger.b')
-    #flog.info("boundary_flag = %s", boundary_flag, name='logger.b')
     #flog.info("nee flux, max=%f, min=%f, mean=%f", max(neeflux),min(neeflux),mean(neeflux), name='logger.b')
 
     #--------------------------------------------------------------
@@ -525,8 +510,7 @@ for(mm in 1:length(pfbfns))
     #**********************************
     # biosphere fluxes convolving
     #**********************************
-    # hourly flux
-    datevalid=dateall[1]+0.5*3600-(0:(ntimes-1))*3600
+    # hourly fxx-0.5*24*3600-ateall[1]+0.5*3600-(0:(ntimes-1))*3600
     ixflux=(1:ntimes)[(datevalid<=xx)][1:endtime]   #time backward
 
     if(bioflux_flag == "SiB3")
@@ -549,8 +533,11 @@ for(mm in 1:length(pfbfns))
     # involve scaling factors to neeflux[,,ixflux], partly, for example as follow
     # xx = "2010-01-24 18:10:00 UTC" || b1 = "2010-01-19 12:00:00 UTC", b3 = "2010-01-29 12:00:00 UTC"
     
-    xxleft=xx-ndays*24*3600
-   
+    #xxleft=xx-ndays*24*3600
+    xx=ISOdatetime(inityr,initmo,initdy,inithr,initmi,0,tz="UTC")
+    yy=xx+(-foottimes*3600)
+    xxleft=min(yy)
+
     #flog.info("xxleft = %s, b2=%s, b2-xxleft=%s", xxleft, b2, b2-xxleft,name='logger.b')
    
 
@@ -651,6 +638,11 @@ for(mm in 1:length(pfbfns))
     fosout=sum(fostmp,na.rm=TRUE) 
     firout=sum(firtmp,na.rm=TRUE)
 
+    print("XXX Line 654")
+    print(mean(foot))
+    print(mean(fosflux))
+    print(fosout)
+
     if(bioflux_flag == "CT_OPT")
     {
         biotmp=foot*bioflux[,,ixflux]*1e6
@@ -673,60 +665,101 @@ for(mm in 1:length(pfbfns))
       lonmax = 180   #-63.5
       npts=dim(endpts)[1]   
     
-      # calulate concentrations
+      # calculate concentrations
       pbou=0
       nval=0
-      dmin=min(as.integer(substr(udaystring,7,8)))
+      #dmin=min(as.integer(substr(aystring,7,8)))
+      #dmax=max(as.integer(substr(udaystring,7,8)))
+      dmin=min(endptsdate)
+      dmax=max(endptsdate)
+ 
+      xx=ISOdatetime(inityr,initmo,initdy,inithr,initmi,0,tz="UTC")
+      yy=xx+(-foottimes*3600)
+
+      print("XXXudaystring for selected flux period")
+      print(udaystring)
+      print("XXXdmin particles")
+      print(dmin)
+      print("XXXdmax particles")
+      print(dmax)
+      print("XXXdobs time")
+      print(xx)
+
+      dd1=as.character(endptsdate)
+      yrlist=substring(dd1,1,4)
+      monlist=substring(dd1,6,7)
+      daylist=substring(dd1,9,10)
+      hrlist=substring(dd1,12,13)
+      milist=substring(dd1,15,16)
+
       for(p in 1:npts) 
       { 
-        dy=0
-        goal=as.integer(daylist[p])
-        for(uu in 1:ndays) 
+        #dy=0
+        #goal=as.integer(daylist[p])
+        dy=ceiling( (as.integer(difftime(xx,endptsdate[p],units='hours')) )/24 )+1
+        #for(uu in 1:ndays) 
+        #{
+        #  dstr=udaystring[uu]
+        #  ddy=as.integer(substr(dstr,7,8))
+         
+         #print(ddy)
+          #print(goal)
+
+         # if(ddy<=goal)
+         #    dy=dy+1
+        #}
+
+        #print("XXXdy=")
+        #print(dy)
+
+        if(boundary_flag == "CTNA")
         {
-          dstr=udaystring[uu]
-          ddy=substr(dstr,7,8)
-        
-          if(ddy<=goal)
-            dy=dy+1
+           i=ceiling((lonarr[p]-lonmin)/3.0)
+           j=ceiling((latarr[p]-latmin)/2.0)
+        }
+        if(boundary_flag == "CTE")
+        {
+           i=ceiling((lonarr[p]-lonmin)/1.0)
+           j=ceiling((latarr[p]-latmin)/1.0)
         }
-     
-      if(boundary_flag == "CTNA")
-      {
-         i=ceiling((lonarr[p]-lonmin)/3.0)
-         j=ceiling((latarr[p]-latmin)/2.0)
-      }
-      if(boundary_flag == "CTE")
-      {
-         i=ceiling((lonarr[p]-lonmin)/1.0)
-         j=ceiling((latarr[p]-latmin)/1.0)
-      }
 
-      # height matching
-      k0=hgtarr[p]    #go upward stair stategy, one by one
-      k=1
-      for(l in 1:25) 
-      {
-        if(k0 > levelhgt[l] && k0 <= levelhgt[25])
-          k=l+1    
-        if(k0 > levelhgt[25])
-          k=25
-      }
+        # height matching
+        k0=hgtarr[p]
+        k=1
+        for(l in 1:25) 
+        {
+          if(k0 > levelhgt[l] && k0 <= levelhgt[25])
+             k=l+1    
+          if(k0 > levelhgt[25])
+             k=25
+        }
       
-      # 3-hourly matching /agt/alt
-      hpos=ceiling((as.integer(hrlist[p])+0.000001+as.integer(milist[p])/60.0 )/3.0)
+        # 3-hourly matching /agt/alt
+        hpos=ceiling((as.integer(hrlist[p])+1e-6+as.integer(milist[p])/60.0 )/3.0)
       
-      tt=bouarr[dy,i,j,k,hpos]
+        print("dy")
+        print(dy)
+        print("hpos")
+        print(hpos)
+
+        tt=bouarr[dy,i,j,k,hpos]
       
-      if(length(tt)==1)   #why sometimes we get a unnormal array?
-      {
-        pbou=pbou+tt  #sum  
-        nval=nval+1
-      }
+        if(length(tt)==1)   #why sometimes we get a unnormal array?
+        {
+          pbou=pbou+tt  #sum  
+          nval=nval+1
+        }
       
-     }
-     fbou=pbou/nval
+      }
+      print("pbou")
+      print(pbou)
+      print(nval)
+      fbou=pbou/nval
     }
 
+    print("XXXfbou")
+    print(fbou)
+
     #case 2: data from Arlyn's emperical curtain
     if(boundary_flag == "EMP")
     {