Commit f345c661 authored by weihe's avatar weihe
Browse files

revised additive scheme

parent 0e7df5f6
...@@ -24,7 +24,7 @@ procs= as.integer(args[1]) ...@@ -24,7 +24,7 @@ procs= as.integer(args[1])
#inputs and outputs #inputs and outputs
stiltfile = paste("stilt_",procs,".rc", sep="") stiltfile = paste("stilt_",procs,".rc", sep="")
curdir="/Storage/CO2/wei/test_4p_bc/ctdas-stilt-newBCfp-air+sur-flux+BC-revisedR-locsur/exec/da/stilt/" curdir="/Storage/CO2/wei/test_additive_fluxfac/ctdas-stilt-base-covL1000new-moreairsites-airlocR/exec/da/stilt/"
fn=paste(curdir, stiltfile,sep="") fn=paste(curdir, stiltfile,sep="")
conf=as.matrix(read.table(fn,header=FALSE)) conf=as.matrix(read.table(fn,header=FALSE))
...@@ -52,10 +52,10 @@ ztop=0 #upper vertical bound for influence projection, i ...@@ -52,10 +52,10 @@ ztop=0 #upper vertical bound for influence projection, i
#if ztop set to zero, *surface* influence will be calculated #if ztop set to zero, *surface* influence will be calculated
#set up an equivalent domain among footprints,fluxes, but not CO2 boundary #set up an equivalent domain among footprints,fluxes, but not CO2 boundary
ncol2=66 #number of pixels in x directions in grid ncol2=120 #66 #number of pixels in x directions in grid
nrow2=40 #number of pixels in y directions in grid nrow2=70 #40 #number of pixels in y directions in grid
LLLon2=(-129) #lower left corner of grid LLLon2=-169.5 #(-129) #lower left corner of grid
LLLat2=22 #lower left corner of grid LLLat2=10.5 #22 #lower left corner of grid
ResLon2=1 #resolution in degrees longitude ResLon2=1 #resolution in degrees longitude
ResLat2=1 #resolution in degrees latitude ResLat2=1 #resolution in degrees latitude
...@@ -119,23 +119,16 @@ for(i in 0:(nsam-1)) #parameters.000.2010010100_2010011100.nc ...@@ -119,23 +119,16 @@ for(i in 0:(nsam-1)) #parameters.000.2010010100_2010011100.nc
if (b1<b0) #b1==b0, revise on Aug 5,2015 if (b1<b0) #b1==b0, revise on Aug 5,2015
{ {
scalefacarr1[,,i+1] = 1 scalefacarr1[,,i+1] = 0
scalefacarr_bc1[,i+1] = 0 scalefacarr_bc1[,i+1] = 0
} }
if (b1>=b0) #b1>b0 if (b1>=b0) #b1>b0
{ {
ncf <- nc_open(paste(samdir,"parameters.",ii,".",tb1,"00","_",tb2,"00",".nc",sep="")) ncf <- nc_open(paste(samdir,"parameters.",ii,".",tb1,"00","_",tb2,"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) scalefac <- ncvar_get(ncf,"parametermap",start=c(11,101),count=c(ncol2,nrow2)) #real52:117,113:152,start=c(52,113),count=c(66,40)
scalefacarr1[,,i+1] = scalefac scalefacarr1[,,i+1] = scalefac * 1e6 # to umol/m2/s
scalefac <- ncvar_get(ncf,"parametervalues_bc") scalefac <- ncvar_get(ncf,"parametervalues_bc")
#dims=dim(scalefac)
#flog.info("dims of BC factor 1 array %s", dims[1], name='logger.b')
print("TTT")
print(scalefac)
print(scalefacarr_bc1[,i+1])
scalefacarr_bc1[,i+1] = scalefac scalefacarr_bc1[,i+1] = scalefac
print(scalefacarr_bc1[,i+1]) print(scalefacarr_bc1[,i+1])
...@@ -143,16 +136,9 @@ for(i in 0:(nsam-1)) #parameters.000.2010010100_2010011100.nc ...@@ -143,16 +136,9 @@ for(i in 0:(nsam-1)) #parameters.000.2010010100_2010011100.nc
} }
ncf <- nc_open(paste(samdir,"parameters.",ii,".",tb2,"00","_",tb3,"00",".nc",sep="")) 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)) scalefac <- ncvar_get(ncf,"parametermap",start=c(11,101),count=c(ncol2,nrow2))
scalefacarr2[,,i+1] = scalefac scalefacarr2[,,i+1] = scalefac * 1e6 # to umol/m2/s
scalefac <- ncvar_get(ncf,"parametervalues_bc") scalefac <- ncvar_get(ncf,"parametervalues_bc")
#dims=dim(scalefac)
#flog.info("dims of BC factor 2 array %s", dims[1], name='logger.b')
print("TTT")
print(scalefac)
print(scalefacarr_bc2[,i+1])
scalefacarr_bc2[,i+1] = scalefac scalefacarr_bc2[,i+1] = scalefac
print(scalefacarr_bc2[,i+1]) print(scalefacarr_bc2[,i+1])
...@@ -279,7 +265,7 @@ for(mm in 1:length(pfbfns)) ...@@ -279,7 +265,7 @@ for(mm in 1:length(pfbfns))
footp=load.ncdf(fn) footp=load.ncdf(fn)
footnc=nc_open(fn) footnc=nc_open(fn)
foot=ncvar_get(footnc,"foot1",start=c(41,12,1),count=c(ncol2,nrow2,-1)) foot=ncvar_get(footnc,"foot1",start=c(1,1,1),count=c(ncol2,nrow2,-1))
#print(mean(foot)) #print(mean(foot))
...@@ -355,8 +341,11 @@ for(mm in 1:length(pfbfns)) ...@@ -355,8 +341,11 @@ for(mm in 1:length(pfbfns))
#endptsdate=ncvar_get(footnc,"endptsdate") #endptsdate=ncvar_get(footnc,"endptsdate")
latarr=endpts[,3] latarr=endpts[,3]
lonarr=endpts[,4] lonarr=endpts[,4]
hgtarr=endpts[,5] + endpts[,6] #agl+grdht #hgtarr=endpts[,5] + endpts[,6] #agl+grdht
magl=endpts[,5] # meters above ground level
masl=endpts[,5] + endpts[,6] # meters above sea level
hgtarr=masl
# analyze the dates to recognize how many data to read # analyze the dates to recognize how many data to read
#dd=as.character(endptsdate) #dd=as.character(endptsdate)
#yrlist=substring(dd,1,4) #yrlist=substring(dd,1,4)
...@@ -522,8 +511,8 @@ for(mm in 1:length(pfbfns)) ...@@ -522,8 +511,8 @@ for(mm in 1:length(pfbfns))
{ {
inxd=ceiling(hh/24) inxd=ceiling(hh/24)
inxh=hh-(inxd-1)*24 inxh=hh-(inxd-1)*24
gppflux[,,ntimes-hh+1]=gpp[inxd,52:117,113:152,inxh] gppflux[,,ntimes-hh+1]=gpp[inxd,11:130,101:170,inxh]
recflux[,,ntimes-hh+1]=rec[inxd,52:117,113:152,inxh] recflux[,,ntimes-hh+1]=rec[inxd,11:130,101:170,inxh]
} }
inxd=ceiling(hh/24) inxd=ceiling(hh/24)
...@@ -532,19 +521,19 @@ for(mm in 1:length(pfbfns)) ...@@ -532,19 +521,19 @@ for(mm in 1:length(pfbfns))
#print(hh) #print(hh)
#print(inxh) #print(inxh)
#print(dim(ocn)) #print(dim(ocn))
ocnflux[,,ntimes-hh+1]=ocn[inxd,52:117,113:152,inxh] ocnflux[,,ntimes-hh+1]=ocn[inxd,11:130,101:170,inxh]
fosflux[,,ntimes-hh+1]=fos[inxd,52:117,113:152,inxh] fosflux[,,ntimes-hh+1]=fos[inxd,11:130,101:170,inxh]
firflux[,,ntimes-hh+1]=fir[inxd,52:117,113:152,inxh] firflux[,,ntimes-hh+1]=fir[inxd,11:130,101:170,inxh]
if(bioflux_flag == "SiBCASA") if(bioflux_flag == "SiBCASA")
{ {
gppflux[,,ntimes-hh+1]=gpp[inxd,52:117,113:152,inxh] gppflux[,,ntimes-hh+1]=gpp[inxd,11:130,101:170,inxh]
recflux[,,ntimes-hh+1]=rec[inxd,52:117,113:152,inxh] recflux[,,ntimes-hh+1]=rec[inxd,11:130,101:170,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') #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") if(bioflux_flag == "CT_OPT")
bioflux[,,ntimes-hh+1]=bio[inxd,52:117,113:152,inxh] bioflux[,,ntimes-hh+1]=bio[inxd,11:130,101:170,inxh]
} }
# replace NA values with 0, as NA was used as zero values for datasets # replace NA values with 0, as NA was used as zero values for datasets
...@@ -564,7 +553,7 @@ for(mm in 1:length(pfbfns)) ...@@ -564,7 +553,7 @@ for(mm in 1:length(pfbfns))
neeflux[,,]=recflux[,,]-gppflux[,,] # nee*1e6 for SiBCASA? neeflux[,,]=recflux[,,]-gppflux[,,] # nee*1e6 for SiBCASA?
if(bioflux_flag == "SiBCASA") if(bioflux_flag == "SiBCASA")
neeflux[,,]=(recflux[,,]-gppflux[,,])*1e6 neeflux[,,]=(recflux[,,]-gppflux[,,])*1e6 # to umol/m2/s
if(bioflux_flag == "CT_OPT") if(bioflux_flag == "CT_OPT")
neeflux[,,]= bioflux[,,]*1e6 # for scaling CT flux neeflux[,,]= bioflux[,,]*1e6 # for scaling CT flux
...@@ -614,14 +603,14 @@ for(mm in 1:length(pfbfns)) ...@@ -614,14 +603,14 @@ for(mm in 1:length(pfbfns))
n = 24-floor(as.numeric(hr) + as.numeric(mi)/60) n = 24-floor(as.numeric(hr) + as.numeric(mi)/60)
if (nchar(dd,type='width') ==10) if (nchar(dd,type='width') ==10)
{ {
n=24 n=24
#flog.info('No hr and mi, set n to 24',name='logger.b') #flog.info('No hr and mi, set n to 24',name='logger.b')
} }
#print("XXXn") #print("XXXn")
#print(n) #print(n)
datevalid=dateall[1]+0.5*3600-(0:(ntimes-1))*3600 ##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] #time backward
#ixflux=(1:ntimes)[(datevalid<=xx)][1:(endtime+n)] #ixflux=(1:ntimes)[(datevalid<=xx)][1:(endtime+n)]
# flog.info('n: %s, endtime: %s',n,endtime,name='logger.b') # flog.info('n: %s, endtime: %s',n,endtime,name='logger.b')
...@@ -688,7 +677,7 @@ for(mm in 1:length(pfbfns)) ...@@ -688,7 +677,7 @@ for(mm in 1:length(pfbfns))
# flog.info("hh: %s, neeflux[,,hh]=%f, scalefacarr2 = %f ",hh, neeflux[30,30,hh], scalefacarr2[30,30,i], name='logger.b') # 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 neeflux1[,,hh] = neeflux[,,hh]+(scalefacarr2[,,i]) #delete "t" transform
dims=dim(scalefacarr_bc2) dims=dim(scalefacarr_bc2)
#flog.info("dims of scalefacarr_bc2: %s, %s", dims[1], dims[2], name='logger.b') #flog.info("dims of scalefacarr_bc2: %s, %s", dims[1], dims[2], name='logger.b')
...@@ -701,7 +690,7 @@ for(mm in 1:length(pfbfns)) ...@@ -701,7 +690,7 @@ for(mm in 1:length(pfbfns))
{ {
#if(i==1) #if(i==1)
# flog.info("hh: %s, neeflux[,,hh]=%f, scalefacarr1 = %f",hh, neeflux[30,30,hh], scalefacarr1[30,30,i],name='logger.b') # 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]) neeflux1[,,hh] = neeflux[,,hh]+(scalefacarr1[,,i])
dims=dim(scalefacarr_bc1) dims=dim(scalefacarr_bc1)
#flog.info("dims of scalefacarr_bc1: %s, %s", dims[1], dims[2], name='logger.b') #flog.info("dims of scalefacarr_bc1: %s, %s", dims[1], dims[2], name='logger.b')
...@@ -718,7 +707,7 @@ for(mm in 1:length(pfbfns)) ...@@ -718,7 +707,7 @@ for(mm in 1:length(pfbfns))
else #all scaling within second lag else #all scaling within second lag
{ {
for (hh in 1:ntimes) for (hh in 1:ntimes)
neeflux1[,,hh] = neeflux[,,hh]*(scalefacarr2[,,i]) neeflux1[,,hh] = neeflux[,,hh]+(scalefacarr2[,,i])
bc1[,hh] = scalefacarr_bc2[,i] bc1[,hh] = scalefacarr_bc2[,i]
neefluxarr[,,,i] = neeflux1[,,] neefluxarr[,,,i] = neeflux1[,,]
...@@ -894,31 +883,33 @@ for(mm in 1:length(pfbfns)) ...@@ -894,31 +883,33 @@ for(mm in 1:length(pfbfns))
k=25 k=25
} }
if (typeflag == "aircraft" && k0 > 3000) # aircraft footprints # notice that in obspack.py,aircraft observations which are lower than 3 km have been excluded from event ids. Here take into account of particles higher than 3 km
# aircraft footprints
if (typeflag == "aircraft" && magl[p] > 3000)
{ {
counter_height = counter_height +1 counter_height = counter_height + 1
if(latarr[i]<20 && lonarr[i]>= (-125) && lonarr[i] <= (-55)) if(latarr[i]< 10.5 && lonarr[i]>= (-169.5) && lonarr[i] < (-50.5))
counter_bc1 = counter_bc1 + 1 counter_bc1 = counter_bc1 + 1
if(latarr[i]>75 && lonarr[i]>= (-125) && lonarr[i] <= (-55)) if(latarr[i]> 80.5 && lonarr[i]>= (-169.5) && lonarr[i] < (-50.5))
counter_bc2 = counter_bc2 + 1 counter_bc2 = counter_bc2 + 1
if(lonarr[i]< (-125) && latarr[i] >= 20 && latarr[i] <= 75) if(lonarr[i]< (-169.5) && latarr[i] >= 10.5 && latarr[i] < 80.5)
counter_bc3 = counter_bc3 + 1 counter_bc3 = counter_bc3 + 1
if(lonarr[i]> (-55) && latarr[i] >= 20 && latarr[i] <= 75) if(lonarr[i]> (-50.5) && latarr[i] >= 10.5 && latarr[i] < 80.5)
counter_bc4 = counter_bc4 + 1 counter_bc4 = counter_bc4 + 1
} }
# surface footprints
if (typeflag == "surface") # surface footprints if (typeflag == "surface")
{ {
counter_height = counter_height +1 counter_height = counter_height +1
if(latarr[i]<20 && lonarr[i]>= (-125) && lonarr[i] <= (-55)) if(latarr[i]< 10.5 && lonarr[i]>= (-169.5) && lonarr[i] < (-50.5))
counter_bc1 = counter_bc1 + 1 counter_bc1 = counter_bc1 + 1
if(latarr[i]>75 && lonarr[i]>= (-125) && lonarr[i] <= (-55)) if(latarr[i]> 80.5 && lonarr[i]>= (-169.5) && lonarr[i] < (-50.5))
counter_bc2 = counter_bc2 + 1 counter_bc2 = counter_bc2 + 1
if(lonarr[i]< (-125) && latarr[i] >= 20 && latarr[i] <= 75) if(lonarr[i]< (-169.5) && latarr[i] >= 10.5 && latarr[i] < 80.5)
counter_bc3 = counter_bc3 + 1 counter_bc3 = counter_bc3 + 1
if(lonarr[i]> (-55) && latarr[i] >= 20 && latarr[i] <= 75) if(lonarr[i]> (-50.5) && latarr[i] >= 10.5 && latarr[i] < 80.5)
counter_bc4 = counter_bc4 + 1 counter_bc4 = counter_bc4 + 1
} }
#if (lonarr[p] < -125 || lonarr[p] > -55 || latarr[p] < 20 || latarr[p] > 75) #if (lonarr[p] < -125 || lonarr[p] > -55 || latarr[p] < 20 || latarr[p] > 75)
...@@ -988,14 +979,14 @@ for(mm in 1:length(pfbfns)) ...@@ -988,14 +979,14 @@ for(mm in 1:length(pfbfns))
} }
dims=dim(bcoutarrtmp) dims=dim(bcoutarrtmp)
flog.info('dims of bcoutarrtmp : %s, %s',dims[1],dims[2],name='logger.b') #flog.info('dims of bcoutarrtmp : %s, %s',dims[1],dims[2],name='logger.b')
bcoutarr = FP_BC_c1*bcoutarrtmp[1,]+FP_BC_c2*bcoutarrtmp[2,]+FP_BC_c3*bcoutarrtmp[3,]+FP_BC_c4*bcoutarrtmp[4,] bcoutarr = FP_BC_c1*bcoutarrtmp[1,]+FP_BC_c2*bcoutarrtmp[2,]+FP_BC_c3*bcoutarrtmp[3,]+FP_BC_c4*bcoutarrtmp[4,]
len=length(bcoutarr) len=length(bcoutarr)
#print(bcoutarr) #print(bcoutarr)
flog.info('FP_BC_c1,2,3,4:%f,%f,%f,%f',FP_BC_c1,FP_BC_c2,FP_BC_c3,FP_BC_c4,name='logger.b') flog.info('FP_BC_c1,2,3,4:%f,%f,%f,%f',FP_BC_c1,FP_BC_c2,FP_BC_c3,FP_BC_c4,name='logger.b')
flog.info('value of bcoutarr[2],[98] : %s %s',bcoutarr[2],bcoutarr[98],name='logger.b') flog.info('value of bcoutarr[1] : %s',bcoutarr[1], name='logger.b')
flog.info('length of bcoutarr : %s',len,name='logger.b') flog.info('length of bcoutarr : %s',len,name='logger.b')
#flog.info('FP_BC_c based on lat, lon and altitude: %s',FP_BC_c,name='logger.b') #flog.info('FP_BC_c based on lat, lon and altitude: %s',FP_BC_c,name='logger.b')
...@@ -1035,10 +1026,13 @@ for(mm in 1:length(pfbfns)) ...@@ -1035,10 +1026,13 @@ for(mm in 1:length(pfbfns))
for(i in 1:nsam) for(i in 1:nsam)
{ {
deltaco2arr[i]=neeoutarr[i]+ocnout+fosout+firout deltaco2arr[i]=neeoutarr[i]+ocnout+fosout+firout
if(typeflag == "aircraft")
fsimuarr[i]=(fbou+bcoutarr[i])*1e-6
if(typeflag == "surface")
fsimuarr[i]=(fbou+bcoutarr[i]+deltaco2arr[i])*1e-6
fsimuarr[i]=(fbou+bcoutarr[i]+deltaco2arr[i])*1e-6 if(i==1)
if(i==99)
flog.info("fsimuarr[i] = %f, fbou = %f, bcoutarr[i] = %f, neeoutarr[i] = %f, deltaco2arr[i] = %f",fsimuarr[i]*1e6,fbou,bcoutarr[i],neeoutarr[i],deltaco2arr[i],name='logger.b') flog.info("fsimuarr[i] = %f, fbou = %f, bcoutarr[i] = %f, neeoutarr[i] = %f, deltaco2arr[i] = %f",fsimuarr[i]*1e6,fbou,bcoutarr[i],neeoutarr[i],deltaco2arr[i],name='logger.b')
outlist<-c(outlist, fsimuarr[i]) outlist<-c(outlist, fsimuarr[i])
......
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