bug fix BC update parameters

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

@@ -37,7 +37,7 @@ samdir  = conf[5,3]
 outdir  = conf[6,3]
 
 # data choice flag
-bioflux_flag = "SiB3"  #"SiB3", "SiBCASA", "CT_OPT"
+bioflux_flag = "SiBCASA"  #"SiB3", "SiBCASA", "CT_OPT"
 boundary_flag= "CT"       #"CT", "EMP"
 
 # optimization option (flux only, or flux + boundary)
@@ -99,8 +99,8 @@ tb3 = paste(substring(b3,1,4),substring(b3,6,7),substring(b3,9,10),sep="")
 
 scalefacarr1=array(NA, dim=c(ncol2,nrow2,nsam))
 scalefacarr2=array(NA, dim=c(ncol2,nrow2,nsam))
-scalefacarr_bc1=array(NA, dim=c(nsam))
-scalefacarr_bc2=array(NA, dim=c(nsam))
+scalefacarr_bc1=array(0., dim=c(nsam))
+scalefacarr_bc2=array(0., dim=c(nsam))
 # make CTDAS to generate domain for North America, read data partly?
 
 flog.info("Reading scaling factor files", name='logger.b')
@@ -115,8 +115,16 @@ for(i in 0:(nsam-1))    #parameters.000.2010010100_2010011100.nc
   
   if (b1<b0)   #b1==b0, revise on Aug 5,2015
   {
-    scalefacarr1[,,i+1] = 1
-    scalefacarr_bc1[i+1] = 1
+    #scalefacarr1[,,i+1] = 1
+    #scalefacarr_bc1[i+1] = 0.0
+
+    ncf <- nc_open(paste(samdir,"parameters.",ii,".",tb2,"00","_",tb3,"00",".nc",sep=""))
+    scalefac <- ncvar_get(ncf,"parametermap",start=c(52,113),count=c(ncol2,nrow2))   #real52:117,113:152,start=c(52,113),count=c(66,40)
+    scalefacarr1[,,i+1] = scalefac
+    scalefac <- ncvar_get(ncf,"parametervalues_bc")
+    scalefacarr_bc1[i+1] = scalefac
+    nc_close(ncf)
+
   }
   if (b1>=b0)  #b1>b0
   {
@@ -407,8 +415,12 @@ for(mm in 1:length(pfbfns))
     neeflux=array(NA, dim=c(ncol2,nrow2,tdim))
     neeflux1=array(NA, dim=c(ncol2,nrow2,tdim))
     neefluxarr=array(NA, dim=c(ncol2,nrow2,tdim,nsam))
+    bc1=array(0., dim=c(tdim))
+    bcarr=array(0., dim=c(tdim,nsam))
+
 
     neeoutarr=array(NA, dim=c(nsam))
+    bcoutarr=array(NA, dim=c(nsam))
     fbouarr=array(NA, dim=c(nsam))
     fsimuarr=array(NA, dim=c(nsam))
    
@@ -555,6 +567,7 @@ for(mm in 1:length(pfbfns))
                  flog.info("neeflux[,,hh]=%f, scalefacarr1 = %f ",neeflux[,,hh], scalefacarr1[,,i], name='logger.b')
                 
                neeflux1[,,hh] = neeflux[,,hh]*(scalefacarr1[,,i])   #delete "t" transform
+               bc1[hh] = scalefacarr_bc1[i]
              }
         for (hh in (diff+1):tdim)
             if(hh<=tdim && hh>=diff+1)
@@ -562,15 +575,19 @@ for(mm in 1:length(pfbfns))
                 if(i==1 && ident=="2010x01x22x19x27x40.0500Nx105.0040Wx00300")
                    flog.info("neeflux[,,hh]=%f, scalefacarr2 = %f",neeflux[,,hh], scalefacarr2[,,i],name='logger.b')
                 neeflux1[,,hh] = neeflux[,,hh]*(scalefacarr2[,,i])
+                bc1[hh] = scalefacarr_bc2[i]
             }
         neefluxarr[,,,i] = neeflux1[,,]  
+        bcarr[,i]=bc1
       }      
       else  #all scaling within second lag
       {
         for (hh in 1:tdim) 
           neeflux1[,,hh] = neeflux[,,hh]*(scalefacarr2[,,i])    
+          bc1[hh] = scalefacarr_bc2[i]
         
         neefluxarr[,,,i] = neeflux1[,,]   
+        bcarr[,i] = bc1
       }
     }
    
@@ -595,12 +612,14 @@ for(mm in 1:length(pfbfns))
      # delta co2 on NEE for all ensemble members
      for(i in 1:nsam)
      {
-        neetmp=foot*neefluxarr[,,ixflux,i]   
+        neetmp=foot*neefluxarr[,,ixflux,i]  
         if(ident=="2010x01x22x19x27x40.0500Nx105.0040Wx00300")
             flog.info("neetmp = %f",neetmp, name='logger.b')
 
         neeout=sum(neetmp,na.rm=TRUE)
         neeoutarr[i]=neeout
+
+        bcoutarr[i]=tail(bcarr[,i],n=1)    
       
         if(ident=="2010x01x22x19x27x40.0500Nx105.0040Wx00300")
             flog.info("neeoutarr[i] = %f",neeoutarr[i], name='logger.b')
@@ -620,7 +639,7 @@ for(mm in 1:length(pfbfns))
     ocntmp=foot*ocnflux[,,ixflux]*1e6
     fostmp=foot*fosflux[,,ixflux]*1e6
     firtmp=foot*firflux[,,ixflux]*1e6
-    
+
     ocnout=sum(ocntmp,na.rm=TRUE)
     fosout=sum(fostmp,na.rm=TRUE) 
     firout=sum(firtmp,na.rm=TRUE)
@@ -690,19 +709,19 @@ for(mm in 1:length(pfbfns))
     fbou=pbou/nval
 
 
-    for (i in 1:nsam)
-    {
-      if(xxleft<b2)  #flux with scaling factors
-      {
-        fbouarr[i] = fbou + (scalefacarr_bc1[i])
-        
-      }
-      else
-      {
-        fbouarr[i] = fbou + (scalefacarr_bc2[i])
-        
-      }
-    }
+    #for (i in 1:nsam)
+    #{
+    #  if(xxleft<b2)  #flux with scaling factors
+    #  {
+    #    fbouarr[i] = fbou + (scalefacarr_bc1[i])
+    #    
+    #  }
+    #  else
+    #  {
+    #    fbouarr[i] = fbou + (scalefacarr_bc2[i])
+    #    
+    #  }
+    #}
     
     #################################################
     # final results and output to files
@@ -714,8 +733,8 @@ for(mm in 1:length(pfbfns))
     for(i in 1:nsam)
     { 
       deltaco2=neeoutarr[i]+ocnout+fosout+firout 
-      #fsimuarr[i]=(fbou+deltaco2)*1e-6        
-      fsimuarr[i]=(fbouarr[i]+deltaco2)*1e-6
+      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')