latest revisions from ivar. fixed scalefac alignment and added log statements

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

@@ -24,7 +24,7 @@ procs= as.integer(args[1])
 #inputs and outputs
 stiltfile = paste("stilt_",procs,".rc", sep="")
 
-curdir="/Storage/CO2/wei/ctdas-stilt-cofiltered-bioSiBCASA-bcCTNA-bcaircraft-fixedstilt/exec/da/stilt/"
+curdir="TYPE PATH HERE/exec/da/stilt/"
 fn=paste(curdir, stiltfile,sep="")
 conf=as.matrix(read.table(fn,header=FALSE))
 
@@ -166,10 +166,13 @@ for(yy in 1:length(uyr))
 # read needed event ids from sample_coordinates files
 ncf <- nc_open(paste(samdir,"sample_coordinates_",tb2,"00","_",tb3,"00","-",procs,".nc",sep="")) 
 eventidarr <- ncvar_get(ncf,"obs_num")
+obs_id <- ncvar_get(ncf,"obs_id")
 
 # filter before loop
 pfbfns=NULL
 fn=NULL
+eventids=NULL
+obsids=NULL
 for(mm in 1:length(fns))  
 {  
   
@@ -179,6 +182,8 @@ for(mm in 1:length(fns))
     if(eid==eventidarr[i])
     {
       pfbfns = c(pfbfns,fns[mm])
+      obsids = c(obsids,obs_id[i])
+      eventids = c(eventids,eventidarr[i])
       break
     }
   }
@@ -213,6 +218,8 @@ for(mm in 1:length(pfbfns))
   #--------------------------------------------
   # STEP 1: Read footprints
   #--------------------------------------------
+ flog.info('obsids: %s',obsids[mm], name='logger.b') 
+flog.info('eventids: %s',eventids[mm], name='logger.b')  #June 23 
   yr=substring(pfbfns[mm],6,9)
   mo=substring(pfbfns[mm],11,12)
   fn=paste(footpdir,yr,"/",mo,"/",pfbfns[mm],sep="")
@@ -389,6 +396,7 @@ for(mm in 1:length(pfbfns))
           gpp[d,,,]=biof$gpp
           rec[d,,,]=biof$resp
           remove(list=c("biof"))
+          flog.info('d: %s, yr: %s, mn: %s, dy: %s',d,yr,mn,dy,name='logger.b')
       }
 
       #flog.info("mean gpp=  %f", mean(gpp[d,,,]), name='logger.b')
@@ -478,6 +486,7 @@ for(mm in 1:length(pfbfns))
         {
            gppflux[,,ntimes-hh+1]=gpp[inxd,52:117,113:152,inxh]
            recflux[,,ntimes-hh+1]=rec[inxd,52:117,113:152,inxh]
+           flog.info('gppflux: %s, gpp: %s, %s, %s, %s',gppflux[30,30,ntimes-hh+1],gpp[inxd,82,143,inxh],ntimes-hh+1,inxd,inxh,name='logger.b')
         }
 
         if(bioflux_flag == "CT_OPT")
@@ -519,23 +528,50 @@ for(mm in 1:length(pfbfns))
     #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))
-    dd=as.character(xx)
-    #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 = 24-floor(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)]
-    print(ixflux)
-    
+#    # 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 = floor(as.numeric(hr) + as.numeric(mi)/60)
+#    n = 24-n
+#    print("XXXn") 
+#    print(n)
+#    datevalid=dateall[1]+0.5*3600-(0:(ntimes-1))*3600
+#    #ixflux=(1:ntimes)[(datevalid<=xx)][1:endtime]   #time backward
+#    #ixflux=(ntimes:1)[(datevalid<=xx)][(endtime+n-1):n]
+#    #ixflux=(ntimes:1)[(datevalid<=xx)][(endtime+n-1):(n)]
+#    ixflux<-c(n:endtime+n-1)
+##    ixflux=(1:ntimes)[(datevalid<=xx)][n:(endtime+n-1)]
+
+
+     #dd=as.character(min(yy))
+     dd=as.character(xx)
+     #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 = 24-floor(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)][1:(endtime+n)]
+     ixflux <-array(n:(endtime+n-1),dim=c(240))
+
+    flog.info('ixflux endtime: %s,%s,%s',endtime,n,as.numeric(hr)+(as.numeric(hr)/60),name='logger.b')
+
+    for (ttt in 1:240)  
+      flog.info('ixflux: %s,%s',ttt,ixflux[ttt],name='logger.b')
+
+
+
+ 
     if(bioflux_flag == "SiB3")   #footprint is backward
     {
         gpptmp=foot*gppflux[,,ixflux]
@@ -571,30 +607,33 @@ for(mm in 1:length(pfbfns))
     {
       if(xxleft<b2)  #flux with scaling factors located in first lag
       {    
-        diff=abs(as.numeric(difftime(b2,xxleft,units='hours')))
+        diff=abs(as.numeric(difftime(b2,xxleft,units='hours')))+24
 
         #if(bioflux_flag == "SiBCASA" || bioflux_flag == "CT_OPT")
         #    diff = diff/3.0
         diff=ceiling(diff)
-        
+        align=24*ceiling((ntimes-diff)/24)#-(ntimes-diff) 
+        if (i==1)
+           flog.info('diff: %s, align: %s',diff,align,name='logger.b')
         #if(ident=="2010x01x22x19x27x40.0500Nx105.0040Wx00300")
         #   flog.info("after diff = %f", diff, name='logger.b')
 
-        for (hh in 1:(ntimes-diff)) 
-            if(hh<=ntimes-diff && hh>=1)
+        for (hh in 1:(align)) 
+            if(hh<=align && hh>=1)
             {
-              #if(i==1 && ident=="2010x01x22x19x27x40.0500Nx105.0040Wx00300")              
-              #   flog.info("neeflux[,,hh]=%f, scalefacarr1 = %f ",neeflux[,,hh], scalefacarr1[,,i], name='logger.b')
+              if(i==1)              
+                 flog.info("hh: %s, neeflux[,,hh]=%f, scalefacarr2 = %f ",hh, neeflux[30,30,hh], scalefacarr2[30,30,i], name='logger.b')
+                
                 
                neeflux1[,,hh] = neeflux[,,hh]*(scalefacarr2[,,i])   #delete "t" transform
                bc1[hh] = scalefacarr_bc2[i]
              }
 
-        for (hh in ( (ntimes-diff+1):ntimes ) )
-            if(hh<=ntimes && hh>=ntimes-diff+1)
+        for (hh in ( (align+1):ntimes ) )
+            if(hh<=ntimes && hh>=align+1)
             {
-                #if(i==1 && ident=="2010x01x22x19x27x40.0500Nx105.0040Wx00300")
-                #   flog.info("neeflux[,,hh]=%f, scalefacarr2 = %f",neeflux[,,hh], scalefacarr2[,,i],name='logger.b')
+                if(i==1)
+                   flog.info("hh: %s, neeflux[,,hh]=%f, scalefacarr1 = %f",hh, neeflux[30,30,hh], scalefacarr1[30,30,i],name='logger.b')
                 neeflux1[,,hh] = neeflux[,,hh]*(scalefacarr1[,,i])
                 bc1[hh] = scalefacarr_bc1[i]
             }
@@ -630,8 +669,10 @@ for(mm in 1:length(pfbfns))
      for(i in 1:nsam)
      {
         neetmp=foot*neefluxarr[,,ixflux,i]   
-        #if(ident=="2010x01x22x19x27x40.0500Nx105.0040Wx00300")
-        #    flog.info("neetmp = %f",neetmp, name='logger.b')
+        neetmp2=foot*neefluxarr[,,6:245,i]
+
+        if(i==1)
+            flog.info("neetmp = %f, %f",sum(neetmp,na.rm=TRUE), sum(neetmp2,na.rm=TRUE),name='logger.b')
 
         neeout=sum(neetmp,na.rm=TRUE)
         neeoutarr[i]=neeout
@@ -673,6 +714,7 @@ for(mm in 1:length(pfbfns))
         bioout=sum(biotmp,na.rm=TRUE)
     }
     
+    foot2=foot 
     remove(list=c("ocnflux","fosflux","firflux","ocntmp","fostmp","firtmp","foot","dateall"))
     gc()
     
@@ -767,6 +809,9 @@ for(mm in 1:length(pfbfns))
         print("hpos")
         print(hpos)
 
+        flog.info('p3:, %s, pbou: %s, dy: %s, i: %s, j: %s, k: %s, hpos: %s',p,pbou,dy,i,j,k,hpos,name='logger.b')
+
+
         tt=bouarr[dy,i,j,k,hpos]  #notice bouarr not be reversed
       
         if(length(tt)==1)   #why sometimes we get a unnormal array?
@@ -820,12 +865,46 @@ for(mm in 1:length(pfbfns))
       fsimuarr[i]=(fbou+bcoutarr[i]+deltaco2)*1e-6        
        #fsimuarr[i]=(fbouarr[i]+deltaco2)*1e-6
 
-      #if(ident=="2010x01x22x19x27x40.0500Nx105.0040Wx00300")
-      #    flog.info("fsimuarr[i] = %f",fsimuarr[i], name='logger.b')
+      if(i==1)
+          flog.info("fsimuarr[i] = %f, fbou = %f, neeoutarr[i] = %f, bcoutarr[i] = %f",fsimuarr[i]*1e6,fbou, neeoutarr[i],bcoutarr[i], name='logger.b')
 
       outlist<-c(outlist, fsimuarr[i])
     }
 
+
+
+#   if (eventids[mm] ==289997 || eventids[mm] ==300326)
+#    {
+#    output_fname=paste("/Storage/ivar/run_test_allsites_bcopt_stiltout_wei/exec","/","lefdata_",eventids[mm],".nc",sep="")
+#
+#   if (file.exists(output_fname)== FALSE)
+#    
+#    {
+#
+#    xdim <- ncdim_def( 'Lon', 'degree', 1:ncol2 )
+#    ydim <- ncdim_def( 'Lat', 'degree', 1:nrow2 )
+#    timedim <- ncdim_def( 'time', 'degree', 1:240 )
+#    ensdim <- ncdim_def( 'ensemble', 'degree', 1:100 )
+#
+#
+#
+#     mv <- 0 # missing value
+#     var_nee <- ncvar_def( name="nee", units="-", dim=list(xdim,ydim,timedim,ensdim), longname="nee",missval=mv )
+#     var_neeco2 <- ncvar_def( name="neeco2", units="-", dim=list(ensdim), longname="neeco2",missval=mv )
+#     var_foot <- ncvar_def( name="foot", units="-", dim=list(xdim,ydim,timedim), longname="foot",missval=mv )
+#
+#
+#    ncid_new <- nc_create( output_fname, list(var_nee,var_neeco2,var_foot))
+#    ncvar_put( ncid_new,var_nee, neefluxarr, start=c(1,1,1,1), count=c(ncol2,nrow2,240,100))
+#    ncvar_put( ncid_new,var_neeco2, neeoutarr, start=c(1), count=c(100))
+#    ncvar_put( ncid_new,var_foot, foot2, start=c(1,1,1), count=c(ncol2,nrow2,240))
+#
+#    nc_close( ncid_new )
+#
+#    }
+#    }
+
+
     #add component varibles:gpp,reco,ocn,fos,fir,boundary
     #outlist<-c(outlist,gppout,recoout,ocnout,fosout,firout,fbou,bioout)
     if(bioflux_flag == "SiB3" || bioflux_flag == "SiBCASA")