fixed more

da/stilt/stilt.co2.simu.2015.Sep.bc.r

-#i------------------------------------------------------------------------------------------------
+#------------------------------------------------------------------------------------------------
 # CO2 concentration simulation based on WRF-STILT/HYSPLIT-NAM12 footprints and
 # SiB3/4 biosphere fluxes, CarbonTracker background fluxes and boundary conditions
 # 
@@ -259,20 +259,21 @@ for(mm in 1:length(pfbfns))
     xx=ISOdatetime(inityr,initmo,initdy,inithr,initmi,0,tz="UTC")
     print("XXXxx?")
     print(xx)
-    yy=xx+(-foottimes*3600)
+    yy=xx+(-foottimes*3600)   #reversed order
     cyy=as.character(yy)
     yrlist=substring(cyy,1,4)
     monlist=substring(cyy,6,7)
     daylist=substring(cyy,9,10)
     hrlist=substring(cyy,12,13)
     milist=substring(cyy,15,16)
-    
+   
     # get unique months and days
     daystring=paste(yrlist,monlist,daylist, sep="")
-    udaystring=unique(daystring)
+    udaystring=unique(daystring) 
     yrmonstring=paste(yrlist,monlist, sep="")
     uyrmonstring=unique(yrmonstring)
    
+    udaystring=rev(udaystring)
     print("XXXudaystring 1st from ident name")
     print(udaystring)
     #current month
@@ -310,6 +311,7 @@ for(mm in 1:length(pfbfns))
     #yrmonstring=paste(yrlist,monlist, sep="")
     #uyrmonstring=unique(yrmonstring)
 
+
     #print("daylist=") #it's wrong to use daylist independently, but to use daystring (yyyymmdd)
     #print(daylist)
 
@@ -359,7 +361,10 @@ for(mm in 1:length(pfbfns))
       yr=substr(datestr,1,4)
       mn=substr(datestr,5,6)
       dy=substr(datestr,7,8)
-    
+
+      print("#####################")
+      print(datestr)
+
       if(boundary_flag=="CTNA")
           bou=load.ncdf(paste(bounddir,"CT2013B.molefrac_glb3x2_",yr,"-",mn,"-",dy,".nc",sep=""))
           #co2(date, level, lat, lon) , ocn_flux_opt(date, lat, lon)
@@ -510,10 +515,26 @@ for(mm in 1:length(pfbfns))
     #**********************************
     # biosphere fluxes convolving
     #**********************************
-    # 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")
+    # hourly flux
+    #print("XXXdateall")
+    #print(dateall)  #reversed order
+    
+    # note: ixflux should not be from 1:ntimes, should be more accurate in hours
+    dd=as.character(min(yy))
+    #yrlist=substring(dd,1,4)
+    #monlist=substring(dd,6,7)
+    #daylist=substring(dd,9,10)
+    hr=substring(dd,12,13)
+    mi=substring(dd,15,16)
+
+    n = ceiling(as.numeric(hr) + as.numeric(mi)/60)
+    print("XXXn") 
+    print(n)
+    datevalid=dateall[1]+0.5*3600-(0:(ntimes-1))*3600
+    #ixflux=(1:ntimes)[(datevalid<=xx)][1:endtime]   #time backward
+    ixflux=(1:ntimes)[(datevalid<=xx)][n:(endtime+n-1)]
+   
+    if(bioflux_flag == "SiB3")   #footprint is backward
     {
         gpptmp=foot*gppflux[,,ixflux]
         recotmp=foot*recflux[,,ixflux]
@@ -638,10 +659,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)
+    #print("XXX Line 640")
+    #print(dim(foot))
+    #print(dim(fosflux))  #tailed by line 520
+    #print(ixflux)
+    #print(fosout)
 
     if(bioflux_flag == "CT_OPT")
     {
@@ -691,12 +713,13 @@ for(mm in 1:length(pfbfns))
       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])
-        dy=ceiling( (as.integer(difftime(xx,endptsdate[p],units='hours')) )/24 )+1
+        #dy=ceiling( (as.integer(difftime(xx,endptsdate[p],units='hours')) )/24 )+1 #if boundary reversed
+        dy=ceiling( (as.integer(difftime(endptsdate[p],min(yy),units='hours')) )/24 )+1
         #for(uu in 1:ndays) 
         #{
         #  dstr=udaystring[uu]
@@ -742,7 +765,7 @@ for(mm in 1:length(pfbfns))
         print("hpos")
         print(hpos)
 
-        tt=bouarr[dy,i,j,k,hpos]
+        tt=bouarr[dy,i,j,k,hpos]  #notice bouarr not be reversed
       
         if(length(tt)==1)   #why sometimes we get a unnormal array?
         {