From f345c661ff2495ad34f5d379b40c9e3a9c7b328e Mon Sep 17 00:00:00 2001
From: weihe <amvdw95@gmail.com>
Date: Fri, 17 Mar 2017 08:59:16 +0000
Subject: [PATCH] revised additive scheme
---
da/stilt/stilt.co2.simu.2015.Sep.fpbc.r | 114 +++++++++++-------------
1 file changed, 54 insertions(+), 60 deletions(-)
diff --git a/da/stilt/stilt.co2.simu.2015.Sep.fpbc.r b/da/stilt/stilt.co2.simu.2015.Sep.fpbc.r
index e7324882..fe40eb0a 100755
--- a/da/stilt/stilt.co2.simu.2015.Sep.fpbc.r
+++ b/da/stilt/stilt.co2.simu.2015.Sep.fpbc.r
@@ -24,7 +24,7 @@ procs= as.integer(args[1])
#inputs and outputs
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="")
conf=as.matrix(read.table(fn,header=FALSE))
@@ -52,10 +52,10 @@ ztop=0 #upper vertical bound for influence projection, i
#if ztop set to zero, *surface* influence will be calculated
#set up an equivalent domain among footprints,fluxes, but not CO2 boundary
-ncol2=66 #number of pixels in x directions in grid
-nrow2=40 #number of pixels in y directions in grid
-LLLon2=(-129) #lower left corner of grid
-LLLat2=22 #lower left corner of grid
+ncol2=120 #66 #number of pixels in x directions in grid
+nrow2=70 #40 #number of pixels in y directions in grid
+LLLon2=-169.5 #(-129) #lower left corner of grid
+LLLat2=10.5 #22 #lower left corner of grid
ResLon2=1 #resolution in degrees longitude
ResLat2=1 #resolution in degrees latitude
@@ -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
{
- scalefacarr1[,,i+1] = 1
+ scalefacarr1[,,i+1] = 0
scalefacarr_bc1[,i+1] = 0
}
if (b1>=b0) #b1>b0
{
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)
- scalefacarr1[,,i+1] = scalefac
+ 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 * 1e6 # to umol/m2/s
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
print(scalefacarr_bc1[,i+1])
@@ -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=""))
- scalefac <- ncvar_get(ncf,"parametermap",start=c(52,113),count=c(ncol2,nrow2))
- scalefacarr2[,,i+1] = scalefac
+ scalefac <- ncvar_get(ncf,"parametermap",start=c(11,101),count=c(ncol2,nrow2))
+ scalefacarr2[,,i+1] = scalefac * 1e6 # to umol/m2/s
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
print(scalefacarr_bc2[,i+1])
@@ -279,7 +265,7 @@ for(mm in 1:length(pfbfns))
footp=load.ncdf(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))
@@ -355,8 +341,11 @@ for(mm in 1:length(pfbfns))
#endptsdate=ncvar_get(footnc,"endptsdate")
latarr=endpts[,3]
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
#dd=as.character(endptsdate)
#yrlist=substring(dd,1,4)
@@ -522,8 +511,8 @@ for(mm in 1:length(pfbfns))
{
inxd=ceiling(hh/24)
inxh=hh-(inxd-1)*24
- gppflux[,,ntimes-hh+1]=gpp[inxd,52:117,113:152,inxh]
- recflux[,,ntimes-hh+1]=rec[inxd,52:117,113:152,inxh]
+ gppflux[,,ntimes-hh+1]=gpp[inxd,11:130,101:170,inxh]
+ recflux[,,ntimes-hh+1]=rec[inxd,11:130,101:170,inxh]
}
inxd=ceiling(hh/24)
@@ -532,19 +521,19 @@ for(mm in 1:length(pfbfns))
#print(hh)
#print(inxh)
#print(dim(ocn))
- ocnflux[,,ntimes-hh+1]=ocn[inxd,52:117,113:152,inxh]
- fosflux[,,ntimes-hh+1]=fos[inxd,52:117,113:152,inxh]
- firflux[,,ntimes-hh+1]=fir[inxd,52:117,113:152,inxh]
+ ocnflux[,,ntimes-hh+1]=ocn[inxd,11:130,101:170,inxh]
+ fosflux[,,ntimes-hh+1]=fos[inxd,11:130,101:170,inxh]
+ firflux[,,ntimes-hh+1]=fir[inxd,11:130,101:170,inxh]
if(bioflux_flag == "SiBCASA")
{
- gppflux[,,ntimes-hh+1]=gpp[inxd,52:117,113:152,inxh]
- recflux[,,ntimes-hh+1]=rec[inxd,52:117,113:152,inxh]
+ gppflux[,,ntimes-hh+1]=gpp[inxd,11:130,101:170,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')
}
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
@@ -564,7 +553,7 @@ for(mm in 1:length(pfbfns))
neeflux[,,]=recflux[,,]-gppflux[,,] # nee*1e6 for SiBCASA?
if(bioflux_flag == "SiBCASA")
- neeflux[,,]=(recflux[,,]-gppflux[,,])*1e6
+ neeflux[,,]=(recflux[,,]-gppflux[,,])*1e6 # to umol/m2/s
if(bioflux_flag == "CT_OPT")
neeflux[,,]= bioflux[,,]*1e6 # for scaling CT flux
@@ -614,14 +603,14 @@ for(mm in 1:length(pfbfns))
n = 24-floor(as.numeric(hr) + as.numeric(mi)/60)
if (nchar(dd,type='width') ==10)
- {
+ {
n=24
#flog.info('No hr and mi, set n to 24',name='logger.b')
- }
+ }
#print("XXXn")
#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+n)]
# flog.info('n: %s, endtime: %s',n,endtime,name='logger.b')
@@ -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')
- neeflux1[,,hh] = neeflux[,,hh]*(scalefacarr2[,,i]) #delete "t" transform
+ neeflux1[,,hh] = neeflux[,,hh]+(scalefacarr2[,,i]) #delete "t" transform
dims=dim(scalefacarr_bc2)
#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))
{
#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])
+ neeflux1[,,hh] = neeflux[,,hh]+(scalefacarr1[,,i])
dims=dim(scalefacarr_bc1)
#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))
else #all scaling within second lag
{
for (hh in 1:ntimes)
- neeflux1[,,hh] = neeflux[,,hh]*(scalefacarr2[,,i])
+ neeflux1[,,hh] = neeflux[,,hh]+(scalefacarr2[,,i])
bc1[,hh] = scalefacarr_bc2[,i]
neefluxarr[,,,i] = neeflux1[,,]
@@ -894,31 +883,33 @@ for(mm in 1:length(pfbfns))
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
- if(latarr[i]<20 && lonarr[i]>= (-125) && lonarr[i] <= (-55))
+ counter_height = counter_height + 1
+ if(latarr[i]< 10.5 && lonarr[i]>= (-169.5) && lonarr[i] < (-50.5))
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
- 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
- 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
}
-
- if (typeflag == "surface") # surface footprints
+ # surface footprints
+ if (typeflag == "surface")
{
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
- 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
- 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
- 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
- }
+ }
#if (lonarr[p] < -125 || lonarr[p] > -55 || latarr[p] < 20 || latarr[p] > 75)
@@ -988,14 +979,14 @@ for(mm in 1:length(pfbfns))
}
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,]
len=length(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('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('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))
for(i in 1:nsam)
{
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==99)
+ if(i==1)
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])
--
GitLab