Commit b6e203fc authored by weihe's avatar weihe
Browse files

bug fixed for BC parameter updating

parent 53f485fe
......@@ -25,7 +25,7 @@ procs= as.integer(args[1])
#inputs and outputs
stiltfile = paste("stilt_",procs,".rc", sep="")
curdir="/Storage/ivar/ctdas-stilt-moresites-aircraft-final/exec/da/stilt/"
curdir="/Storage/CO2/wei/ctdas-stilt-cofiltered-bioSiBCASA-bcCTNA-bcaircraft/exec/da/stilt/"
fn=paste(curdir, stiltfile,sep="")
conf=as.matrix(read.table(fn,header=FALSE))
......@@ -38,7 +38,7 @@ outdir = conf[6,3]
# data choice flag
bioflux_flag = "SiBCASA" #"SiB3", "SiBCASA", "CT_OPT"
boundary_flag= "CT" #"CT", "EMP"
boundary_flag= "CTNA" #"CTNA", "CTE", "EMP"
# optimization option (flux only, or flux + boundary)
opt_para_flag = "FLUX" # "FLUX_BOU","FLUX"
......@@ -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(0., dim=c(nsam))
scalefacarr_bc2=array(0., dim=c(nsam))
scalefacarr_bc1=array(NA, dim=c(nsam))
scalefacarr_bc2=array(NA, dim=c(nsam))
# make CTDAS to generate domain for North America, read data partly?
flog.info("Reading scaling factor files", name='logger.b')
......@@ -115,16 +115,8 @@ 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] = 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)
scalefacarr1[,,i+1] = 1
scalefacarr_bc1[i+1] = 0
}
if (b1>=b0) #b1>b0
{
......@@ -149,19 +141,25 @@ for(i in 0:(nsam-1)) #parameters.000.2010010100_2010011100.nc
# from corresponding period for optimzation
#---------------------------------------------------------------------------------
# according to b1, b2, decide which months
y1=substring(b2,1,4)
y2=substring(b3,1,4)
m1=substring(b2,6,7)
m2=substring(b3,6,7)
uyr=unique(c(y1,y2))
umon=unique(c(m1,m2))
flog.info("Determining ready-use footprint files", name='logger.b')
fns=NULL
for(mm in 1:length(umon))
for(yy in 1:length(uyr))
{
pfbpath=paste(footpdir,year,"/",umon[mm],"/",sep="") #"stilt2010x01x01x18x59x33.4057Nx081.8334Wx00285.nc"
tmp=list.files(pfbpath,pattern="nc", all=T)
fns=c(fns,tmp)
flog.info("Footprints in year %s",uyr[yy],name='logger.b')
for(mm in 1:length(umon))
{
pfbpath=paste(footpdir,uyr[yy],"/",umon[mm],"/",sep="") #"stilt2010x01x01x18x59x33.4057Nx081.8334Wx00285.nc"
tmp=list.files(pfbpath,pattern="nc", all=T)
fns=c(fns,tmp)
}
}
# 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")
......@@ -198,13 +196,25 @@ flog.info(log, name='logger.b')
#flog.info("Start convolving for all footprint files", name='logger.b')
newfile=T #output results into a new file
#
#read EMP boundary monthly into memory
#
if(boundary_flag == "EMP")
{
bouf=paste(bounddir,y1,"/",m1,"/","CO2.v201209_v1.init.stilt.",y1,m1,".txt",sep="")
print(bouf)
val<-as.matrix(read.table(bouf,header=TRUE))
}
for(mm in 1:length(pfbfns))
{
#--------------------------------------------
# STEP 1: Read footprints
#--------------------------------------------
yr=substring(pfbfns[mm],6,9)
mo=substring(pfbfns[mm],11,12)
fn=paste(footpdir,year,"/",mo,"/",pfbfns[mm],sep="")
fn=paste(footpdir,yr,"/",mo,"/",pfbfns[mm],sep="")
footp=load.ncdf(fn)
footnc=nc_open(fn)
......@@ -298,7 +308,11 @@ for(mm in 1:length(pfbfns))
# 2-1: Read fluxes and boundary data
#--------------------------------------------------------------
ndays=length(udaystring)
bouarr=array(0,dim=c(ndays,120,90,25,8)) #boundary : lon,lat,hgt,time
if(boundary_flag=="CTNA")
bouarr=array(0,dim=c(ndays,120,90,25,8)) #boundary : lon,lat,hgt,time
if(boundary_flag=="CTE")
bouarr=array(0,dim=c(ndays,360,180,25,8))
#biospheric fluxes
......@@ -339,9 +353,12 @@ for(mm in 1:length(pfbfns))
mn=substr(datestr,5,6)
dy=substr(datestr,7,8)
bou=load.ncdf(paste(bounddir,"CT2013B.molefrac_glb3x2_",yr,"-",mn,"-",dy,".nc",sep=""))
# co2(date, level, lat, lon) , ocn_flux_opt(date, lat, lon)
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)
if(boundary_flag=="CTE")
bou=load.ncdf(paste(bounddir,"3d_molefractions_1x1_",yr,mn,dy,".nc",sep=""))
bgf=load.ncdf(paste(bgfdir,"CT2013B.flux1x1.",yr,mn,dy,".nc",sep=""))
if(bioflux_flag == "SiB3")
......@@ -368,13 +385,22 @@ for(mm in 1:length(pfbfns))
if(bioflux_flag == "CT_OPT")
bio[d,,,]=bgf$bio.flux.opt
bouarr[d,,,,]=bou$co2
if(boundary_flag=="CTNA")
{
bouarr[d,,,,]=bou$co2
remove(list=c("bou"))
}
if(boundary_flag=="CTE")
{
bouarr[d,,,,]=bou$co2*1e6
remove(list=c("bou"))
}
ocn[d,,,]=bgf$ocn.flux.opt # here we can read less, so change the sizes of ocnarr
fos[d,,,]=bgf$fossil.flux.imp
fir[d,,,]=bgf$fire.flux.imp
remove(list=c("bou","bgf"))
remove(list=c("bgf"))
}
#--------------------------------------------------------------
......@@ -415,10 +441,10 @@ 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))
......@@ -500,6 +526,7 @@ for(mm in 1:length(pfbfns))
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')
#--------------------------------------------------------------
......@@ -583,10 +610,10 @@ for(mm in 1:length(pfbfns))
else #all scaling within second lag
{
for (hh in 1:tdim)
neeflux1[,,hh] = neeflux[,,hh]*(scalefacarr2[,,i])
neeflux1[,,hh] = neeflux[,,hh]*(scalefacarr2[,,i])
bc1[hh] = scalefacarr_bc2[i]
neefluxarr[,,,i] = neeflux1[,,]
neefluxarr[,,,i] = neeflux1[,,]
bcarr[,i] = bc1
}
}
......@@ -612,14 +639,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)
bcoutarr[i]=tail(bcarr[,i],n=1)
if(ident=="2010x01x22x19x27x40.0500Nx105.0040Wx00300")
flog.info("neeoutarr[i] = %f",neeoutarr[i], name='logger.b')
......@@ -656,33 +683,44 @@ for(mm in 1:length(pfbfns))
#--------------------------------------------------------------
# 2-4: Calculate boundary
#--------------------------------------------------------------
# boundary domain : lat 22.5~61.5, lon -128.5~-63.5
latmin=-90 #22.5
latmax=90 #61.5
lonmin=-180 #-128.5
lonmax=180 #-63.5
npts=dim(endpts)[1]
# calulate concentrations
pbou=0
nval=0
dmin=min(as.integer(substr(udaystring,7,8)))
for(p in 1:npts)
{
dy=0
goal=as.integer(daylist[p])
for(uu in 1:ndays)
{
dstr=udaystring[uu]
ddy=substr(dstr,7,8)
#case 1: data from CT2013B
if(boundary_flag == "CTNA" || boundary_flag == "CTE")
{
# boundary domain : lat 22.5~61.5, lon -128.5~-63.5
latmin=-90 #22.5
latmax=90 #61.5
lonmin=-180 #-128.5
lonmax=180 #-63.5
npts=dim(endpts)[1]
# calulate concentrations
pbou=0
nval=0
dmin=min(as.integer(substr(udaystring,7,8)))
for(p in 1:npts)
{
dy=0
goal=as.integer(daylist[p])
for(uu in 1:ndays)
{
dstr=udaystring[uu]
ddy=substr(dstr,7,8)
if(ddy<=goal)
dy=dy+1
if(ddy<=goal)
dy=dy+1
}
if(boundary_flag == "CTNA")
{
i=ceiling((lonarr[p]-lonmin)/3.0)
j=ceiling((latarr[p]-latmin)/2.0)
}
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
......@@ -705,24 +743,32 @@ for(mm in 1:length(pfbfns))
nval=nval+1
}
}
fbou=pbou/nval
}
fbou=pbou/nval
#case 2: data from Arlyn's emperical curtain
if(boundary_flag == "EMP")
{
#bouf=paste(bounddir,year,"/",mon,"/","CO2.v201209_v1.init.hysplit.",year,mon,".txt",sep="")
#print(bouf)
#val<-as.matrix(read.table(bouf,header=TRUE))
for(p in 1:dim(val)[1])
if(val[p,2]==eventid)
{
fbou=as.numeric(val[p,3])
break
}
}
#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])
#
# }
# if(xxleft<b2)
# fbouarr[i] = fbou + (scalefacarr_bc1[i])
# else
# fbouarr[i] = fbou + (scalefacarr_bc2[i])
#}
#################################################
# final results and output to files
#################################################
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment